Quédate en cajas: nociones de la modelación matemática por compartimentos


                        Luis Zavala Sansón


* Nota de divulgación para la publicación electrónica TODoS@CICESE.
** Documento PDF disponible en la edición del 20 de mayo de 2020.
*** Versión actual: 26 de mayo 2020

Introducción

En los últimos meses se ha hablado mucho sobre la modelación matemática de la pandemia de Covid-19 que está cambiando al mundo. Es quizá uno de los residuos positivos de la enfermedad, ya que ha despertado la curiosidad de muchos estudiantes, profesionales y público en general por entender, o al menos conocer un poco, la forma de hacer diagnósticos o predicciones sobre un fenómeno que nos tiene resguardados en casa. Afortunadamente, los componentes básicos de los modelos que mayormente se utilizan para representar las epidemias se pueden comprender con algunos conocimientos de matemáticas superiores a nivel licenciatura. Concretamente, de ecuaciones diferenciales ordinarias de primer orden. De hecho, los modelos a los que nos referiremos son aplicables en diferentes disciplinas de las ciencias naturales y sociales. Así que lo que se expone a continuación puede ser fácilmente digerible para los estudiantes de carreras científicas (física, química, biología, medicina), ingenierías o ciencias económicas (actuaría, matemáticas financieras, entre otras [1]).

El concepto abstracto fundamental es el de considerar que el objeto de estudio consiste en uno o varios compartimentos o cajas separadas. Cada caja tiene un "contenido", el cual cambia en el tiempo según su relación con las demás cajas; es decir, los compartimentos están interconectados unos con otros. Los modelos consisten en establecer la regla que determina el cambio en el tiempo del contenido de cada caja (lo cual incluye las reglas de interconexión). Para aclarar esto, vamos a discutir en forma inductiva modelos de una, dos y tres cajas, y al final hablaremos de la generalización de la idea.

Una caja

Supongamos que queremos saber cuántos habitantes habrá en un tiempo futuro en una región o país. La caja es el lugar en cuestión y el contenido de la caja es la población, que llamaremos P. El cambio de la población en el tiempo t está dado por su derivada dP/dt, con unidades de habitantes/tiempo. Un modelo básico para este problema consiste en la siguiente regla: el cambio de la población en cada unidad de tiempo (dP/dt) es proporcional a la población P. Se trata del modelo poblacional formulado en 1798 por el poco piadoso clérigo británico Thomas Malthus. Bajo este supuesto la población aumenta cada vez más rápido porque los habitantes que van engrosando la población aportan en conjunto más descendientes. La expresión matemática (que llamaremos ecuación de evolución) y su solución exponencial son las siguientes:



La constante a es la tasa de aumento neto de la población, con unidades de 1/tiempo (lo cual se deduce fácilmente al recordar que las unidades en cada término de cualquier ecuación deben ser las mismas). Este número es fundamental, ya que indica la rapidez con la que se añaden nuevos individuos a la población. La constante P0 es el número de habitantes al tiempo t=0.
Figura 1. Crecimiento de la población P en el modelo exponencial de Malthus (curva azul) para un valor inicial de P0=100 personas y una tasa de crecimiento a=0.04 en unidades arbitrarias de 1/tiempo. La curva roja corresponde al crecimiento lineal, con b=10 personas/tiempo.

La Figura 1 muestra la curva (azul) de crecimiento de la población. Vemos que P aumenta progresivamente más rápido, ya que se trata de un crecimiento exponencial, también llamado geométrico. Por comparación, hemos añadido en la misma gráfica un crecimiento lineal o aritmético partiendo de la misma población inicial. El modelo lineal y su solución son:



donde b es una constante que indica el número de habitantes por unidad de tiempo que se añaden a la población (con unidades de habitantes/tiempo). Dependiendo del valor de b, podemos ver que para tiempos pequeños la curva lineal puede crecer más rápido que la exponencial, pero posteriormente la curva exponencial la supera ampliamente. A principios del siglo XIX, el crecimiento Malthusiano generó la perspectiva de una catástrofe ante el supuesto aumento lineal del suministro de alimentos. Pero al poco tiempo comenzaron a aparecer modelos poblacionales que daban cuenta de factores limitantes, como veremos en el siguiente modelo de dos cajas.

Dos cajas

Un problema con dos compartimentos es el de la población de dos especies que conviven en un ecosistema dado: una es la presa (P) y la otra es la depredadora (D). Por lo tanto, ahora necesitamos las ecuaciones de evolución para dP/dt y dD/dt. Si inicialmente se tiene una población P bien establecida y aparece una especie depredadora D, ¿cómo cambian las poblaciones de ambas especies en el tiempo? Este problema se puede modelar con ecuaciones de tipo Lotka-Volterra, conocidas desde principios del siglo XX y generalizadas por el célebre matemático soviético A. Kolmogorov, pilar de las teorías sobre turbulencia en Mecánica de Fluidos. La versión más sencilla del modelo es suponer que la especie presa sucumbe ante la depredadora a una tasa c (con unidades de 1/tiempo), y que el aumento de la especie depredadora es esa misma tasa. Las ecuaciones de evolución se pueden escribir como:


donde N=P+D es la población total del ecosistema. El lado derecho de ambas ecuaciones tiene dos aspectos a destacar. Primero, vemos que es el mismo término pero con signo contrario, lo que significa que lo que pierde la especie P lo gana la especie D exactamente. De hecho, al sumar ambas ecuaciones vemos que dN/dt=0, es decir, el número total de individuos permanece constante. Segundo, se trata de un término no lineal (porque es el producto de ambas poblaciones PD). Esta no-linealidad  expresa la interacción entre las dos especies.

El hecho de que N sea constante implica que las expresiones (3) y (4) se pueden reescribir como la llamada ecuación diferencial logística, que aquí mostramos junto con sus respectivas soluciones:



donde P0 y D0 son las poblaciones iniciales de ambas especies. Ambas soluciones se llaman, por supuesto, función logística, un nombre bastante ilógico que propuso el matemático belga (de la región flamenca) Pierre Verhulst en la primera mitad del siglo XIX. El trabajo de Verhulst estaba enfocado en refinar el modelo de crecimiento Malthusiano.
Figura 2. Modelo logístico de presa-depredador (5)-(6). Los valores iniciales de las respectivas poblaciones son P0=999 y D0=1 y la tasa de interacción es c=0.15 en unidades arbitrarias de 1/tiempo.

Para apreciar la evolución de la presa y el depredador vamos a graficar las curvas logísticas en la Figura 2. Al tiempo t=0 la especie P tiene un valor inicial dado y la D un valor pequeño. Si vemos con detenimiento las ecuaciones de evolución, notaremos que para tiempos iniciales D crece exponencialmente a una tasa cP/N, mientras que P decrece a una tasa más lenta cD/N (este decrecimiento es menor porque la población D apenas acaba de llegar). Pero los papeles se invierten cuando D ha crecido y P disminuido: para tiempos largos D crece más lento porque hay menos individuos P.

Un modelo de presa-depredador más general es el siguiente:


Nótese que ahora la constante que mide el consumo de presas c tiene un significado diferente, pues tiene unidades de la tasa temporal/individuo. Además, es diferente al coeficiente d de la especie D. Por otro lado, hemos añadido una tasa de natalidad n de las presas y una mortalidad m de los depredadores (ambas con unidades de 1/tiempo). Estos términos lineales también se pueden interpretar como entradas o salidas de individuos que salen de la caja respectiva (por ejemplo, por migración).
Figura 3. Modelo de presa-depredador (7)-(8). Las poblaciones iniciales son P0=500 y D0=500. Los parámetros de interacción son c=0.001 y d=0.003 en unidades arbitrarias de 1/(tiempo X individuos). La tasa de natalidad de las presas es n=0.4 y la mortalidad de los depredadores m=0.3, ambas en unidades de 1/tiempo.

Existen soluciones y propiedades del sistema (7)-(8) que podemos discutir en otra ocasión. Por ahora vamos a graficar en la Figura 3 las curvas para las dos especies considerando algunos valores particulares de los parámetros. Lo que se observa es una oscilación de las poblaciones en el tiempo. Ambas comienzan con el mismo número de individuos. Al prinicipio los depredadores capturan presas, lo que hace crecer D (curva roja) y disminuir P (curva azul). Pero con el tiempo escasean las presas y D comienza a disminuir después de alcanzar un máximo. Al haber menos depredadores ahora crece P, y así sucesivamente. Ambas poblaciones están en un equilibrio oscilante. La gráfica incluye la suma de las especies N=P+D (curva negra), la cual ya no es una constante sino que también oscila en el tiempo. Lo relevante de este ejemplo es que el comportamiento de las soluciones depende críticamente de los parámetros utilizados (c,d,m,n). Si utilizamos valores diferentes podríamos obtener curvas que se comporten parecido pero con amplitudes y periodos distintos. O podríamos tener la extinción de una de las especies o incluso de las dos. El comportamiento de este modelo y sus muchas variantes ha sido y sigue siendo objeto de investigación desde que fue formulado hace más de 100 años.

Tres cajas

Ahora ya podemos entender uno de los modelos epidemiológicos más básicos. Se trata de un modelo de N personas contenidas en tres cajas: en la primera hay personas susceptibles (S) de enfermarse , es decir sanas; en la segunda hay personas infectadas (I) y en la tercera hay personas que dejan de estar infectadas (R), ya sea porque se recuperaron o porque murieron (vamos a suponer que se recuperaron). En este modelo los recuperados no regresan al grupo S, sino que se quedan en R. Evidentemente, la denominación del modelo es SIR. Las ecuaciones de evolución de cada compartimento son




Al sumar las tres ecuaciones vemos que dN/dt=0, por lo que el número de personas N=S+I+R permanece constante todo el tiempo.

Los términos que relacionan a susceptibles e infectados son no lineales (producto SI), similares a los que vimos en el modelo de presa-depredador. En efecto, en este caso los infectados "consumen" a los susceptibles. El parámetro β contiene la tasa de infección (con unidades de 1/tiempo). Epidemiológicamente se trata del número de contactos por individuo por unidad de tiempo multiplicado por la probabilidad de transmitir la enfermedad. Se puede advertir que es un número difícil de determinar, porque implica conocer la movilidad promedio de los individuos de la población para saber cuántos contactos tiene en promedio por unidad de tiempo. Además, se requiere conocer la probabilidad de contagio, lo cual puede ser difícil de saber en enfermedades de reciente aparición. Por otro lado, los infectados se reducen conforme las personas pasan al compartimento R a una tasa γ. Dicha tasa es el inverso del tiempo de recuperación y depende del tipo de enfermedad. El parámetro clave de los modelos de epidemias es el número básico de reproducción R0. En el modelo SIR, este número es R0=β/γ. Se trata de la comparación de la tasa de infección con respecto a la tasa de recuperación. Si R0>1 (β>γ) la epidemia se desarrolla, porque la razón de infección es más rápida que la de recuperación.


Figura 4. (a) Modelo SIR (9)-(10) con valores iniciales S0=300,000, I0=10 y R0=0 individuos. Los parámetros son β=0.5 y γ=0.17 (d-1). El tiempo final es 105 días (15 semanas). (b) Comparación de las curvas de infectados I para dos valores de β.

Vamos a ver un ejemplo tipo Covid utilizando las unidades de tiempo en días. Supongamos que β=1/2 (d-1). El tiempo de recuperación ronda los 6 días, por lo que γ=1/6 (d-1). De ambos parámetros resulta R0=3, considerado como elevado y que revela una enfermedad muy perniciosa. La Figura 4a muestra las curvas SIR para una población de 300,000 habitantes en la que inicialmente hay sólo 10 personas infectadas. El número de susceptibles S decae en el tiempo porque pasan a ser infectados I. Al poco tiempo comienza a haber recuperados R. La curva más relevante es la de los infectados (curva roja): aumenta gradualmente, alcanza un pico y despues decae casi completamente: la epidemia ha desaparecido. A su vez, los recuperados aumentan poco a poco. Al tiempo final hay un pequeño número de susceptibles que nunca se enfermó, casi nadie queda infectado y la mayoría son individuos recuperados.

Para mitigar una epidemia se puede tratar de reducir el número de contactos entre personas, es decir, β. Esta es la razón básica de las medidas de distanciamiento social que se han aplicado en casi todo el mundo para reducir los efectos del Covid-19. Si en nuestro ejemplo anterior reducimos β a 1/3 (d-1), la curva de infectados I queda como lo indica la Figura 4b (línea discontinua). El pico de infectados se reduce, aunque la curva se ensancha notoriamente. Sin embargo, ese es el efecto buscado, porque el objetivo es impedir que I sea muy alto en un tiempo corto, y así evitar la saturación de los servicios de salud. Con la curva "aplanada" sigue habiendo muchos infectados, pero ahora distribuidos en el tiempo.

Comentarios finales

Con lo expuesto anteriormente ya podemos apreciar la lógica de los modelos de cajas: se trata de representar al sistema bajo estudio en compartimentos interconectados, y de plantear la ecuación de evolución para cada caja. Las conexiones, las entradas y salidas de las cajas se modelan con parámetros que suponemos representan dichas interacciones. Hay que recordar que los compartimentos pueden contener grupos de individuos, organismos, sustancias, objetos o lo que el problema bajo estudio demande.

Figura 5
Figura 5. Modelo NPZ de plancton oceánico (ver texto).

Por ejemplo, un problema de tres compartimentos en el océano aborda la cantidad de nutrientes (N) disueltos en la superficie del mar, los cuales son consumidos por el fitoplancton (P), el cual es a su vez consumido por el zooplancton (Z). El esquema de la Figura 5 ilustra este tipo de modelo, comúnmente llamado NPZ. Las flechas horizontales indican la relación entre cada caja y las verticales se refieren a entradas y salidas de cada compartimento debido a los factores que se quieran modelar. En el ejemplo hay entrada y salida de nutrientes por procesos físicos (afloramientos o hundimientos de agua [2]). También hay pérdida de fitoplancton y zooplancton por mortalidad. Cada proceso se modela por al menos un parámetro (indicado con las letras s,f,b,g,m,n en el esquema).  Los modelos NPZ se pueden ampliar para considerar diferentes nutrientes o especies de plancton. También se pueden adaptar a modelos dinámicos que resuelvan el movimiento del océano en una región determinada [3]. Dichos movimientos introducen nuevos términos entre compartimentos (nuevas flechas en el esquema). Y pueden llegar a ser tan complicados como se quiera.

En los modelos epidemiológicos existe una amplísima gama que puede ser de varias cajas que representen, por ejemplo, infectados latentes (que no contagian durante cierto tiempo), hospitalizados, visitantes o migrantes, defunciones, y un largo etcétera. Como en los modelos poblacionales, las posibilidades son enormes, por lo que siguen siendo objeto de investigación en la epidemiología matemática [4]. Otra vertiente de investigación es la forma de determinar los parámetros del modelo. Por ejemplo, el valor del coeficiente β en varios modelos tipo SIR se puede tratar de ajustar mediante el uso de datos que se van presentando conforme la epidemia avanza. De esta forma, las curvas (y las predicciones del modelo) se van ajustando día con día [5].

Es importante recordar que los modelos de compartimentos son altamente idealizados, porque tratan de resumir mediante parámetros las interacciones que puede haber entre los elementos de las diferentes cajas, así como sus fuentes y sumideros. Debido a estas sobre-simplificaciones, es importante que los modelos conlleven un cálculo de incertidumbre o error que permita estimar su confiabilidad. Otra acción útil es la de reajustar los parámetros en los casos en los que vayan apareciendo datos del fenómeno. Es lo que ha sucedido con la epidemia actual bajo el punto de vista del modelo SIR: el conocimiento de la enfermedad permite estimar γ, mientras que los registros de personas infectadas son importantes (aunque siempre insuficientes) para estimar β. En cualquier caso, se debe ser cauteloso con los resultados en este y otros problemas. No se debe olvidar que en la modelación matemática es útil dividir al mundo para entenderlo, pero si no juntamos las piezas posteriormente, lo habremos destruido.

Referencias

[1] Existen varios sitios que exponen el problema en forma muy didáctica, incluyendo la página oficial del gobierno de México:
https://coronavirus.gob.mx/
https://institucional.us.es/blogimus/2020/03/covid-19-analisis-por-medio-de-un-modelo-seir/
https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology
https://www.idmod.org/docs/hiv/model-seir.html

[2] https://www.myroms.org/wiki/images/0/04/Ecosystem.pdf

[3] Zavala Sansón, L. & Provenzale, A., 2009. The effects of abrupt topography on plankton dynamics. Theoretical population biology, 76(4), pp.258-267.

[4] Zlojutro, A., Rey, D. \& Gardner, L., 2019. A decision-support framework to optimize border control for global outbreak mitigation. Scientific reports, 9(1), pp.1-14.

[5] https://covid19-projections.com/