Las epidemias de enfermedades infecciosas son fenómenos sociales aparte de ser problemas importantes de salud pública y medicina. Una de sus características más determinantes es que su dinámica se ve fuertemente afectada por la heterogeneidad de las poblaciones en muchos niveles de integración, que comprenden elementos genéticos, fenotípicos, conductuales y económicos que inducen, a su vez, variabilidad e incertidumbre en el proceso de transmisión de la enfermedad. El modelo clásico de la epidemiología que ha cobrado en estos tiempos de pandemia cierta popularidad, es el modelo de Kermack y McKendrick publicado en el primer cuarto del siglo XX. En la Figura 1 se ilustra una típica solución de este modelo que será discutido más adelante, en la sección 4. Por el momento nos limitaremos a identificar algunas características básicas de esta solución. La gráfica roja es ejemplo de una curva epidémica, en este caso de prevalencia o porcentaje de casos activos diarios. La curva epidémica usada comúnmente en salud pública es la de incidencia o número de nuevos casos diarios. La curva negra es el número de personas recuperadas de la enfermedad y la azul el número de personas susceptibles, es decir, que todavía no han contraído la enfermedad y pueden enfermarse. La curva epidémica tiene un único máximo y tanto su crecimiento como su decrecimiento es exponencial. Fue por primera vez usada para describir la epidemia de peste bubónica de principios del siglo XX en la India por Kermack y McKendrick [KM27].
Anderson G. McKendrick fue un médico escocés nacido a finales del siglo XIX. Trabajó en el Sudán y en la India en donde estudió, de manera autodidáctica, matemáticas [Bac11]. Resultado de su experiencia y estudios fue el artículo [McK25]. En él, McKendrick presenta una primera elaboración de un modelo sobre epidemias que desarrolló más ampliamente con W.O. Kermack en una serie de artículos iniciada en 1927. William O. Kermack también fue escoces y nació hacia finales del siglo XIX. Una explosión en el laboratorio donde trabajaba lo dejó completamente ciego en 1924 [Bac11] lo que no le impidió iniciar con McKendrick, en 1927, una serie de artículos llamados Contributions to the Mathematical Theory of Epidemics que constituyen el fundamento de los modelos epidémicos determinísticos clásicos [KM27; KM32; KM33; KM37; KM38]. Cabe resaltar que, aparte de su trabajo en epidemias, McKendrick [McK25] presenta por primera vez uno de los modelos más básicos y generales de la demografía. Brevemente, pues en una sección posterior expandiremos el tema, el modelo de McKendrick postula que si $N(a,t) da$ es el número de individuos entre las edades $a$ y $a+da$ al tiempo $t$, y $\mu(a)$ es la probabilidad de muerte a la edad $a$, entonces la tasa de cambio de la población en función de su edad y el tiempo transcurrido, está dada por la ecuación
\[\frac{\partial}{\partial t} N(a,t)+\frac{\partial }{\partial a}N(a,t) =-\mu(a) N(a,t).\]Esta ecuación es importante, entre otras cosas, porque incorpora la heterogeneidad poblacional inherente a los procesos demográficos y a la transmisión de enfermedades infecciosas: la edad. En este caso, la mortalidad es específica de la edad y no un parámetro constante general. A lo largo de las secciones de este texto, usaremos la formulación de McKendrick para calcular indicadores importantes de transmisión.
Las enfermedades infecciosas se transmiten mediante algún proceso de contacto entre individuos infecciosos y aquellos que son susceptibles a la enfermedad. Los modelos más básicos, como los que discutiremos aquí, presuponen que todos los individuos son iguales en todos los aspectos relevantes para la transmisión y que lo único que los diferencia es el estado con respecto a la enfermedad en que se hallan. Estos modelos son compartamentales. La población en riesgo se divide en clases disjuntas de individuos típicamente clasificados como susceptibles, expuestos, infecciosos, inmunes pero estas clases se pueden ampliar con otras pertinentes. En la presente pandemia se han publicado muchos artículos donde se proponen compartimientos para personas hospitalizadas, asintomáticas, en terapia intensiva, aisladas, etc. La proliferación de compartimientos conlleva un incremento significativo en el número de parámetros y, por consiguiente, a la necesidad de asignar valores adecuados a éstos. Por lo común, el número de parámetros en estas generalizaciones es muy elevado, lo que lleva al dilema de sobreajuste (overfitting). Una solución socorrida consiste en tomar valores de la literatura que se suponen «correctos» (usualmente el valor promedio encontrado en laboratorio o campo) y asignarlos directamente al parámetro en cuestión con el resultado de que el ajuste a los datos se centra en solo 3 o 4 parámetros de un total dos o tres veces mayor. Esta estrategia no resuelve el problema del overfitting, simplemente lo oculta. En particular, cuando el objetivo del modelo es predecir el número de individuos en compartimientos de interés, esta forma de reducción paramétrica se vuelve crítica pues las predicciones dependen altamente de la sensibilidad de las soluciones a los parámetros que arbitrariamente se han mantenido fijos. Afortunadamente siempre ha habido alternativas. A lo largo de la pandemia, para efectos de predicción, la metodología que ha mostrado mayor robustez y poder predictivo es la construcción de modelos «ensamble» o ensamblados, que combinan diferentes tipos de estrategias de modelación (modelos estadísticos, compartamentales, estocásticos, de agentes, inteligencia artificial) que informan las proyecciones pero considerando diferentes niveles de variabilidad e incertidumbre. El otro tipo de modelos que históricamente han sido muy informativos en epidemiología, y que en esta pandemia han mostrarlo serlo también, son los modelos más generales y simples desde el punto de vista estructural, basados en ecuaciones diferenciales. El objetivo de estos modelos no es la predicción de tendencias sino, más bien, la generación de escenarios e hipótesis explicativas de los fenómenos de relevancia. En esta contribución exploraremos, muy someramente, los resultados básicos existentes para cuatro de estos modelos: Kermack-McKendrick, McKendrick-von Foester, Sharpe-Lotka y el número reproductivo.
Alfred J. Lotka fue un demógrafo norteamericano nacido a finales del siglo XIX que contribuyó notablemente al desarrollo de lo que ahora llamamos dinámica de poblaciones. Encontró, simultánea e independientemente de Vito Volterra, las ecuaciones de dinámica poblacional para la interacción depredador-presa (de ahí que ambos compartan el crédito en el nombre de estas ecuaciones). En colaboración con F.R. Sharpe (de quien, desafortunadamente, no he podido encontrar mayores datos biográficos) desarrolló la idea de la distribución estable de edades asociada al crecimiento (o decrecimiento) exponencial de una población. La ecuación de Sharpe-Lotka relaciona la tasa de fertilidad y probabilidad de supervivencia, con la capacidad de crecimiento de una población. En esta sección revisaremos someramente los principios esenciales de este modelo siguiendo libremente a [Web84; Mur02; IP14]
Dada una especie biológica arbitraria de reproducción sexual, definamos a $\pi(a)$ como la probabilidad de sobrevivir a la edad $a$, $0\leq a\leq L$; $\pi$ es una función monótonamente decreciente, suave en $[0,L)$ que satisface
\[\pi \mayorque 0,\qquad \pi(0)=1,\qquad \lim_{a\to L}\pi(a) = 0.\]$L$ se denomina la esperanza de vida máxima en la población. Generalmente, para facilitar los cálculos, se supone que $L=\infty.$ La tabulación de $\pi(a)$ se conoce como la tabla de vida de una población. Podemos definir también
\[\pi(a,\tau)=\frac{\pi(\tau)}{\pi(a)},\qquad 0\leq a \leq \tau\leq L\]como la probabilidad de sobrevivir a la edad $\tau$ dado que se está vivo en la edad $a$. En consecuencia
\[1-\pi(a,\tau)\]es la probabilidad de que un individuo de edad $a$ muera en el intervalo $(a,\tau)$. Obsérvese que podemos definir la siguiente aproximación
\[\mu(a)h\approx 1-\pi(a, a+h)=\frac{ \pi(a)-\pi(a+h)}{\pi(a)}\]cuando $h\to 0$. De aquí resulta que, para un individuo de edad $a$, con un pequeño error de aproximación, $\mu(a)h$ nos da la probabilidad de que muera en el intervalo $(a, a+h)$. A $\mu(a)$ se le denomina la función de mortalidad; $\mu(a) da$ es la probabilidad de morir en el intervalo $(a, a+da)$. Usando la función de mortalidad, podemos reescribir $\pi(a)$ y $\pi(a,\tau)$ como
\begin{eqnarray*} \pi(a)&=&\exp(-\int^a_0 \mu(s)ds),\\ \pi(a,\tau)&=&\exp(-\int^{\tau}_a \mu(s)ds) \end{eqnarray*}La integral
\[\int^{a_2}_{a_1}\rho(a,t)da\]es el número de individuos con edades entre $a_1$ y $a_2$ al tiempo $t$. Similarmente
\[P(t)=\int_0^L\rho(a,t)da\]es la población total. A la función
\[\rho:[0,L) \times R^+\to R^+\]se le denomina la distribución de edades de la población; finalmente
\[B(t)=\rho(0,t)\]es la tasa de nacimientos al tiempo $t$. La relación entre la tasa de nacimientos $B(t)$ y la distribución de edades $\rho$ se da a través de la integral:
\[B(t) = \int_0^L\beta(a)\rho(a,t)da\]donde $\beta(a) da$ es el número esperado de nacimientos de un individuo de edad $a$ al tiempo $t$. A la función $\beta: [0,L)\to R^+$ le llamaremos la función de fertilidad. La dinámica de una población (el cambio del número de individuos en el tiempo) estructurada por edades, está completamente definida para toda $t \mayorque 0$ por $\beta$, por $\pi$ y por la llamada distribución de edades inicial
\[\rho(a,0)=\phi(a).\]Para entender esta relación con mayor intuición, usaremos ahora un ejemplo ilustrativo. Supongamos que iniciamos un estudio demográfico amplio en una población humana al tiempo $t=0$, que continuamos por muchos años. Los sujetos de $70$ años de edad en $t=10$, es decir, 10 años después de iniciado el estudio, tenían $60$ años cuando el estudio se inició en $t=0$. Por lo tanto
\[\rho(70,10)=\rho(60,0)\pi(60,70)=\phi(60)\pi(60,70).\]Siempre que $t\leq a$, es decir, siempre que el tiempo transcurrido desde el inicio del estudio sea menor o igual a la edad biológica del sujeto observado, tenemos
\[\rho(a,t)=\phi(a-t)\pi(a-t,a), \qquad t\leq a.\]Ahora tomemos a los sujetos de $10$ años de edad al tiempo $t=70$. Estos sujetos nacieron cuando el estudio llevaba $t=60$ años y sobrevivieron hasta llegar a los 10 años de edad. Por lo tanto
\[\rho(10,70)=B(60)\pi(10).\]Siempre que $0\leq a\leq t$, es decir, que el tiempo de observación es mayor o igual a la edad biológica del sujeto observado, tenemos
\[\rho(a,t)=B(t-a)\pi(a), \qquad 0\leq a\leq t.\]Sustituyendo estas relaciones en la tasa de nacimientos
\[B(t) = \int_0^L\beta(a)\rho(a,t)da,\]obtenemos que si $t \menorque L$
\[B(t)=\int^t_0 \beta(a)\pi(a)B(t-a)da+\int_t^L\beta(a)\pi(a-t,a)\phi(a-t)da;\]además, si $L\leq t$
\[B(t) =\int_0^L \beta (a) \pi(a) B(t-a) da. \]Si definimos como cero las funciones de natalidad, supervivencia y distribución inicial cuando $a \mayorque L$ y definimos
\[\Phi(t)=\int_t^L\beta(a)\pi(a-t,a)\phi(a-t)da, \qquad k(a)=\beta(a)\pi(a),\]obtenemos la llamada Ecuación de Sharpe-Lotka
\[B(t)=\int^t_0k(a)B(t-a) da + \Phi(t),\]que es una ecuación integral lineal tipo Volterra para $B$. Conociendo $B$ se puede encontrar $\rho$, lo que la relaciona directamente con la ecuación discutida en la siguiente sección. La ecuación de Sharpe-Lotka tiene una solución única $B$ en $R^+$ si y sólo si $r$ es la solución única real de
\begin{equation}\label{eq:char0} \int^L_0\exp(-ra)\beta(a)\pi(a) da =1. \end{equation}A $r$ se le llama la tasa intrínseca de crecimiento poblacional.
Heinz von Foester fue un investigador austriaco, nacionalizado estadounidense, nacido a principios del siglo XX. Contribuyó al estudio del ciclo celular, para el que planteó una ecuación similar a la de McKendrick, aparentemente de manera independiente. El modelo que describiremos en esta sección es denominado, consecuentemente, como el modelo de McKendrick-von Foester. Este es un modelo de naturaleza demográfica y es importante porque muchas de las ideas y conceptos que tenemos en el estudio de enfermedades infecciosas provienen de él al igual que del modelo de Sharpe-Lotka revisado en la sección anterior. En este caso estaremos trabajando con cohortes de individuos. Una cohorte es un grupo de individuos con la misma edad. La «densidad poblacional» de estos individuos está dada por la función $\rho(a,t)$; $t+h$ unidades de tiempo después, su densidad será $\rho(a+h, t+h)$. Definimos la derivada de cohorte como
\[D\rho(a,t) =\lim_{h\to o}\frac{\rho(a+h,t+h)-\rho(a,t)}{h}\]que proporciona la tasa a la que la densidad de cada cohorte cambia en el tiempo. La mortalidad en la cohorte la calculamos como
\[\int^h_0\rho(a+h,t+h)\mu(a+h)dh.\]En consecuencia, el cambio en la densidad poblacional de una cohorte es el balance entre su envejecimiento y su mortalidad
\[\rho(a,t)-\rho(a+h, t+h) =\int^h_0\rho(a+h,t+h)\mu(a+h)dh.\]Dividiendo por $h$ y dejando $h\to 0$ obtenemos,
\[ D\rho +\mu\rho=0.\]Si $\rho$ es una función suave, entonces $D\rho= \rho_a+\rho_t$ y el modelo puede escribirse como
\begin{equation}\label{eq:KF} \frac{\partial}{\partial a} \rho(a,t)+\frac{\partial}{\partial t}\rho(a,t)-\mu(a,t) \rho(a,t)=0. \end{equation}El problema matemático de la dinámica poblacional estructurada por edades es el de encontrar una solución $\rho:[0,L) \times R^+ \to R^+$ de la ecuación de McKendrick-von Foester en $(0,L)\times (0,\infty)$ que satisfaga apropiadas condiciones iniciales y de frontera dadas por la tasa de nacimientos
\[ \rho(0,t) = \int_0^L\beta(a)\rho(a,t) da\]para $t \mayorque 0$ y la distribución inicial de edades
\[\rho(a,0) = \phi(a)\]para $0\leq a \menorque L$. La teoría y resultados sobre la existencia, unicidad y propiedades de las soluciones de esta ecuación son muy amplias. Por ahora nos concentraremos en un tipo especial de soluciones que permiten ligar la demografía con la epidemiología. Estamos interesados en encontrar soluciones de la ecuación de McKendrick-von Foester de la forma
\[\rho(a,t) = A(a)T(t).\]Este tipo de soluciones se denominan soluciones etarias estables. Sustituyendo esta expresión en la ecuación ( 2) obtenemos que
\[\frac{A'(a)}{A(a)}+\mu(a)=-\frac{\dot{T}(t)}{T(t)}=-p,\]donde $p$ es una constante. Resolviendo estas ecuaciones diferenciales ordinarias obtenemos:
\begin{eqnarray*} A(a) &= &C_1 \exp(-\int_0^a (\mu(s)+p)ds),\\ T(t)&=&C_2 e^{pt}, \end{eqnarray*}Sustituyendo la solución hallada en
\[ \rho(0,t) = \int_0^L\beta(a)\rho(a,t) da\]llegamos, después de un poco de álgebra, a la expresión
\[\int_0^L\beta(a)\pi(a) e^{-pa}da =1.\]La solución real única, que denotaremos por la letra $r$ de esta ecuación, es la tasa intrínseca de crecimiento poblacional que encontramos ya en nuestro estudio de la ecuación de Sharpe-Lotka. Podemos concluir entonces que las soluciones etarias estables tienen forma
\[\rho(a,t) =B_0 \pi(a) e^{-ra}e^{rt},\]donde $B_0$ es la tasa de nacimientos inicial. Finalmente el número
\[R = \int_0^L \beta(a)\pi(a) da,\]se denomina el número reproductivo y representa el número de nacimientos viables per cápita en la población. Este número fue generalizado años después en epidemiología dando lugar al número reproductivo básico.
El llamado modelo general de Kermack-McKendrick es desarrollado en [KM27] y es un complejo modelo integrodiferencial que toma en cuenta la edad de la infección en cada individuo. Más adelante abundaremos sobre este concepto también. Un caso especial del modelo, presentado como la ecuación (29) en la referencia antes citada, es el clásico modelo de transmisión epidemiológica. Definamos $S(t)$, $I(t)$ y $R(t)$ a los individuos en una población que son susceptibles, infecciosos o recuperados de alguna enfermedad, respectivamente. Entonces, la evolución de la epidemia queda descrita con las ecuaciones
\begin{eqnarray}\label{eq:KK} S'(t)&=&-\kappa S(t) I(t),\\\nonumber I'(t)&=& \kappa S(t) I(t)-\gamma I(t),\\ R'(t)&=&\gamma I(t),\nonumber \end{eqnarray}donde $\gamma$ es la tasa de recuperación. A la expresión $\kappa=\phi c$ se denomina tasa efectiva de contacto, con componentes $c$, la tasa de contactos per cápita por unidad de tiempo, y $\phi$, la probabilidad de que uno de los contactos resulte en infección en el hospedero susceptible. En este modelo el número de contactos per cápita por unidad de tiempo es constante, es decir, supone implícitamente que todos los individuos son iguales y que sus hábitos y conductas individuales no cambian debido a la edad, preferencias, sexo, lo que es una generalización bastante fuerte. Para efectos de comparación con una epidemia real y cercana, en la Figura 2 se muestra la curva epidémica de la Ciudad de México desde su inicio hasta noviembre del 2020. Su forma evidencia que los supuestos antes mencionados del modelo Kermack-McKendrick no son satisfechos y la pregunta necesaria ahora es saber por qué. Por el momento, nos centraremos en caracterizar algunas de las propiedades del modelo para después, en las conclusiones, comentar al respecto de lo recién dicho.
Observemos que, en ausencia de nuevos contactos infecciosos, con $I(0)=I_0$ hospederos infectados iniciales, tenemos $I'(t)=-\gamma I(t)$ y, en consecuencia, podemos escribir a la fracción de los individuos infecciosos iniciales que se han recuperado de la enfermedad al tiempo $t$, como
\[1-\frac{I(t)}{I_0}=1-\exp(-\gamma t)\]y, por lo tanto,
\[\frac{1}{\gamma}=\int_0^\infty t\gamma\exp (-\gamma t)dt\]es el tiempo medio de residencia en el estado infeccioso pues podemos tomar a $F(t)=1-\exp(-\gamma t)$ como la distribución (exponencial) de los tiempos de residencia en el compartimiento $I$. El término $\kappa S(t)I(t)$ es la incidencia y proporciona el número de casos nuevos por unidad de tiempo; por su lado $\kappa I$ es la fuerza de infección que describe el número de infecciones potenciales generadas en promedio por un individuo infeccioso típico. Las ecuaciones ( 3) han aparecido desde hace años en muchos libros de texto (eg. [Mar15; Het05; BVW08]) sobre biología matemática, modelación y epidemiología cuantitativa junto con sus muchas variantes. Siendo este modelo un modelo compartamental que, además, tiene la nolinealidad más común en el término de incidencia, es fácil añadirle compartimientos que, en el SARS-CoV-2, han sido aquellos para individuos hospitalizados, reportados, en terapia intensiva, por ejemplo. Añadir es fácil, parametrizar esos compartimientos e interpretarlos correctamente es otra cosa. La simplicidad estructural de los modelos Kermack-McKendrick no debe ser tomada a la ligera.
En esta sección revisaremos algunos resultados directos y clásicos del modelo Kermack-McKendrick ( 3). Para ello nos basaremos en las referencias [Bac11; Mar15; Het05; BVW08; MW20]. Para empezar definamos el parámetro que surge de la segunda ecuación en ( 3) cuando $S=S_0$ constante al inicio de la epidemia
\begin{equation}\label{eq:R0} {\cal R}_0=\kappa S_0\frac{1}{\gamma}. \end{equation}Este parámetro es el número de infecciones que un individuo infeccioso puede generar ($\kappa$) en una población completamente susceptible (i.e., $S_0\approx N$) durante el tiempo en que es infeccioso ($1/\gamma$). Este es el famoso, no hallo otra palabra mejor para describirlo, número reproductivo básico. La población total $N=S(t)+I(t)+R(t)$ es constante para todo tiempo y, en lo que resta de este texto, la normalizaremos de tal forma que $N=1$. En particular, si $S(0)=S_0$, $I(0)=I_0$ y $R(0)=R_0$ son cantidades no negativas, podemos ver que $S(t)$ es una función no creciente del tiempo y que $I(t)$ crece si ${\cal R}_0 \mayorque 1$ y decrece si ${\cal R}_0 \menorque 1$. Sumando la primera y segunda ecuación en ( 3) e integrando obtenemos
\[S(t)-S(0)+I(t)-I(0)=-\gamma\int^t_0I(s)ds,\]lo que significa que $I(t)$ debe estar acotada y, más aún, que $\int_0^{\infty}I(t)dt \menorque \infty.$ De la primera ecuación en ( 3) y usando las observaciones anteriores podemos obtener un resultado importante e informativo sobre las epidemias simplemente integrándola:
\[S(t)=S(0)\exp\Bigl(-\kappa \int_0^tI(s)ds\Bigr),\]y, por consiguiente,
\[\lim_{t \to {\infty}} S(t)=S_{\infty} \mayorque 0,\qquad \lim_{t\to \infty}I(t)=0,\]es decir, en una epidemia, al final, la proporción de infecciosos será cero (la epidemia terminará) pero siempre queda un conjunto de gente que no se enferma. Este fenómeno no tiene que ver con inmunidad sino con el carácter gregario de la población humana y la dinámica poblacional propia del patógeno. Tomando entonces el límite cuando $t\to\infty$ y observando que
\[\log \Bigl(\frac{S(t)}{S_0}\Bigr)=-\kappa \int_0^{t}I(s)ds,\]tenemos que
\[S_{\infty}=S(0)+I(0)+\frac{\gamma}{\kappa} \log\Bigl(\frac{S_{\infty}}{S_0}\Bigr).\]Usando las relaciones anteriores podemos ahora obtener una relación entre el total de casos $S_{\infty}-S_0$ y los casos acumulados al tiempo $t$ que denotaremos por $T_C$. Tenemos
\[T_C(t)=-(\log S(t)-\log S_0)=\frac{\gamma}{\kappa}\int^t_0I(s)ds,\]y
\[T_C(\infty)=\lim_{t\to\infty}T_C(t)=\log\frac{S_0}{S_{\infty}}=\frac{\gamma}{\kappa}\int^{\infty}_0I(s)ds\]de donde
\[S_{\infty}=S_0\exp\Bigl(-\frac{\kappa}{\gamma}T_C(\infty)\Bigr).\]Usando esta expresión, junto con las tres igualdades anteriores, obtenemos una versión de la llamada ecuación del tamaño final
\[S(0)\exp\Bigl(-\frac{\kappa}{\gamma}T_C(\infty)\Bigr)=S(0)+I(0)-\frac{\gamma}{\kappa}T_C(\infty),\]que nos permite estimar el número de casos totales $T_C(\infty)$ que se esperan conociendo el número reproductivo y la población inicial de susceptibles.
El número reproductivo básico es uno de los parámetros epidemiológicos más importantes pero también uno de los más malentendidos en estas épocas de la epidemia. Su equivocada interpretación ha generado, en diversos lugares alrededor del mundo, sesgos muy notorios en los modelos Kermack-McKendrick usados para hacer predicciones sobre la incidencia, muertes y hospitalizaciones. En esta sección nos basaremos en los resultados de [WL07; Par+19] para desarrollar la presentación. La definición del número reproductivo básico es muy simple y ya ha sido dada en la sección anterior justificándola con una fórmula derivada del modelo SIR. Es importante entender que la definición del número reproductivo es independiente de los modelos Kermack-McKendrick. Para efectos de claridad se presenta nuevamente:
Número reproductivo básico. El número reproductivo básico, denotado por ${\cal R}_0$, es el número de infecciones secundarias generadas en una población completamente susceptible por un individuo infeccioso típico o promedio, durante el tiempo en que puede transmitir la enfermedad.
Es importante resaltar que los modelos Kermack-McKendrick son, como todo modelo, aproximaciones de la realidad y están basados en suposiciones específicas que deben ser tomadas en cuenta cuando los modelos se aplican a problemas reales. Por ejemplo, el modelo SIR de Kermack-McKendrick supone que el tiempo durante el cual un individuo infeccioso puede transmitir la enfermedad es constante e igual al periodo de recuperación $1/\gamma$. Esta aproximación es correcta para algunas enfermedades pero, para otras, sobre todo aquellas cuyo tiempo de recuperación es largo, no lo es. Antes de ampliar este comentario, empecemos primero caracterizando este importante parámetro. Si tomamos la ecuación para $I$, suponemos que la población es completamente susceptible y estamos al inicio de la epidemia, con un número inicial de infecciosos $I_0$, para $t$ corto tenemos $I'(t)=(S_0\kappa-\gamma)I(t)$ cuya solución está dada por $I(t)=I_0\exp\Bigl(\gamma ({\cal R}_0-1)t\Bigr);$ es decir, el crecimiento inicial de la epidemia es exponencial y depende de la magnitud de ${\cal R}_0$: si ${\cal R}_0 \mayorque 1$ la epidemia crece, si ${\cal R}_0 \menorque 1$ decrece. Obsérvese también que si $r=\gamma ({\cal R}_0-1)$ es la tasa de crecimiento de la epidemia (ya discutida un par de veces en las secciones anteriores pero ahora referente a la prevalencia diaria o casos activos diarios), encontramos que
\begin{equation}\label{eq:r} {\cal R}_0=1+\frac{r}{\gamma}. \end{equation}Cuando ${\cal R}_0 \mayorque 1$, el crecimiento exponencial de la epidemia es muy sensible a la incertidumbre en el valor de ${\cal R}_0$ y en el número inicial de infectados.
Para lo que sigue vamos a regresar por un momento a las ecuación de McKendrick y Sharpe-Lotka que presentamos en secciones anteriores pero reformulada en términos de infecciones. Cuando nos enfermamos pasamos por varias etapas. Inicialmente el patógeno entra a nuestro organismo en muy pequeñas cantidades. Si el sistema inmune es exitoso defendiéndonos, no nos enfermaremos pero si no es así, la concentración del patógeno en la sangre crecerá a lo largo de horas o días hasta rebasar un umbral a partir del cual comenzamos a mostrar síntomas y en algunos casos, como en el SARS-CoV-2, podemos ya ser infecciosos. El incremento del patógeno en la sangre continúa hasta alcanzar un punto máximo después del cual la concentración de patógenos disminuye hasta extinguirse. En algún momento, antes de esta extinción, un enfermo deja de ser infeccioso. Es entonces posible caracterizar a un individuo enfermo en términos de la edad de su infección. Esta caracterización es importante porque la capacidad de infectar a otra persona dependerá de la edad de la infección. Sea entonces $i(a,t)da$ el número de individuos infecciosos existente al tiempo $t$ con edad de infección entre $a$ y $a+da$. La ecuación de McKendrick correspondiente a esta formulación es la siguiente
\[\frac{\partial}{\partial t} i(a,t)+\frac{\partial }{\partial a}i(a,t) =-\gamma(a) i(a,t).\]La tasa de recuperación, $\gamma=\gamma (a)$, depende ahora de la edad de infección. Las condiciones iniciales y de frontera de esta ecuación son
\[i(a,0)=\phi (a),\qquad i(0,t)=B(t)=\int_0^{\infty}\kappa (a) i(a,t)da,\]en donde $\kappa(a)$ representa la distribución de la transmisibilidad de individuos infectados con edad de infección $a$ o, equivalentemente, la tasa efectiva de contactos dependiente de la edad de infección y $\phi(a)$ se denomina la distribución inicial de edades de infección. La distribución $\kappa(a)$ generalmente tiene un máximo y tiende a cero cuando la edad de la infección se hace grande, es decir, la enfermedad termina o la persona se muere. Estamos interesados en encontrar soluciones de la ecuación diferencial parcial separables, como hicimos con la ecuación de McKendrick-von Foester, de la forma $i(a,t)=T(t)A(a)$. Sustituyendo en la ecuación, puede determinarse que estas soluciones son posibles siempre y cuando exista una única $r$ real tal que satisfaga la llamada ecuación característica
\begin{equation}\label{eq:char} \int_0^{\infty}e^{-ra}\kappa (a)e^{-\int_0^a\gamma(s)ds}da=1. \end{equation}Esta expresión es, como puede verificarse, la misma que obtuvimos del modelo de Sharpe-Lotka ( 1) y también en la del McKendrick-von Foester. El parámetro $r$ ahora se interpreta como la tasa de crecimiento intrínseca de la epidemia y está, valga la redundancia, intrínsecamente ligada al tiempo de doblamiento del número de casos definido como $T_d=\frac{\log 2}{r}$. La expresión
\[K(a)=\kappa (a)e^{-\int_0^a\gamma(s)ds}\]es el producto de la tasa de transmisibilidad $\kappa(a)$ y la probabilidad de ser todavía infeccioso a la edad de infección $a$ dada por $e^{-\int_0^a\gamma(s)ds}$. $K(a)$ la podemos interpretar entonces como la tasa de generación de infecciones secundarias de un individuo infeccioso de edad de infección $a$. El número reproductivo básico es simplemente el número total de casos secundarios generado por un individuo infeccioso durante todo el tiempo en el que puede infectar, es decir, está dado por la siguiente integral
\[{\cal R}_0=\int_0^{\infty}K(a)da.\]Recordemos que en el modelo SIR, la probabilidad de ser infeccioso a la edad $t$ es exponencial, $I(t)/I_0=\exp{(-\gamma t)}$, y la transmisibilidad es independiente de la edad de infección, $\kappa S_0$. Entonces
\[{\cal R}_0=\int_0^{\infty}K(a)da=\kappa S_0\int_0^{\infty}e^{-\gamma t }dt=\frac{\kappa S_0}{\gamma}\]que recupera la expresión ( 4) hallada antes para el modelo SIR. A $K(a)$ también se le llama la distribución del periodo generacional. Este es uno de los conceptos más importantes en el cálculo de ${\cal R}_0$ pero que, desafortunadamente, ha sido largamente desatendido en la comunidad de modelación matemática debido a la casi exclusiva atención a métodos de cálculo del número reproductivo sustentados en los modelos Kermack-McKendrick. Si definimos ahora
\[G(a)=\frac{K(a)}{{\cal R}_0},\]obtenemos que $0\leq G(a) \leq 1$ para toda $a$ y tenemos la función de densidad del periodo o intervalo generacional. Éste se define como el tiempo que pasa entre que una persona se contagia de la infección y que contagia a alguien. En términos de edad de infección, el intervalo generacional es la edad de infección a la que un individuo infectado infecta a alguien más. En consecuencia, dividiendo ( 6) por ${\cal R}_0$ obtenemos que el recíproco del número reproductivo es el resultado de la transformada de Laplace, o función generadora de momentos, de la distribución de intervalos generacionales evaluada en la tasa intrínseca de crecimiento $r$:
\begin{equation}\label{eq:laplace} M(r)=\int_0^{\infty}e^{-ra}G(a)da=\frac{1}{{\cal R}_0}. \end{equation}El intervalo generacional de un brote epidémico se calcula directamente de la incidencia, usando la fecha de inicio de síntomas como aproximación a la fecha de contagio, y haciendo uso extensivo de métodos de tamizaje y seguimiento de contactos. A esta aproximación se le llama intervalo serial. Importantemente, el promedio de la distribución del intervalo seriales o generacionales, no necesariamente coincide con el período de recuperación (o período de infecciosidad como también se le llama) que es unicamente una aproximación específica del modelo SIR Kermack-McKendrick ( 3). Podemos ahora fácilmente calcular el número reproductivo básico para diferentes distribuciones de intervalos seriales. Definamos a $T_g$ como el intervalo serial promedio. El lector puede confirmar, después de hacer unos cálculos directos, que para la distribución exponencial
\[G(a)=\gamma\exp(-\gamma t),\qquad \frac{1}{M(r)}={\cal R}_0=1+r/\gamma\]en donde $T_g=1/\gamma$ que coincide con la fórmula presentada en ( 5). Igualmente, usando una distribución $\delta$ centrada en el periodo de recuperación $1/\gamma$ obtenemos
\[G(a)=\delta_{1/\gamma},\qquad \frac{1}{M(r)}={\cal R}_0=e^{r/\gamma},\]y, finalmente, para una distribución normal
\[G(a)=\frac{1}{\sqrt{2\pi}\sigma}e^{(a-1/\gamma)^2/(2\sigma)^2},\qquad \frac{1}{M(r)}={\cal R}_0=e^{r/\gamma-\frac{1}{2}\sigma^2r^2}.\]Obsérvese, en este último caso, que una mayor varianza de la distribución, disminuye el valor del número reproductivo.
A manera de conclusión permítame el lector desarrollar un poco más las ideas en torno a los modelos Kermack-McKendrick. Arriba revisamos solamente el modelo SIR que es uno de los dos más simples (el otro es el modelo SIS pero la historia de éste está más ligada a otro famoso epidemiólogo, Ronald Ross en el estudio de enfermedades transmitidas por insectos vectores, particularmente paludismo que ya he tratado en otro lugar [Vel99]). El modelo Kermack-McKendrick que le sigue en complejidad al SIR es el SEIR que tiene la forma
\begin{eqnarray*} S'(t)&=&-\kappa S(t) I(t),\\\nonumber E'(t)&=&\kappa S(t) I(t)-\eta E(t),\\ I'(t)&=& \eta E(t)-\gamma I(t),\\ R'(t)&=&\gamma I(t),\nonumber \end{eqnarray*}donde el nuevo compartimiento $E$ contiene a los individuos que han sido infectados pero todavía no son infecciosos; el parámetro $1/\eta$ es el período de incubación de la enfermedad. En este modelo también se puede calcular la tasa de crecimiento de la epidemia linealizando alrededor del llamado punto libre de la enfermedad, que tiene coordenadas $(S^*,E^*,I^*,R^*)=(1,0,0,0)$ (recordando que la población total la hemos normalizado para que $N=1$), y calculando el eigenvalor dominante $r$. También en este modelo existe una relación entre $R_0$ y $r$ que es, naturalmente, un poco más compleja que la que hallamos en el modelo SIR. Para el modelo SEIR esta relación es
\[{\cal R}_0= \Bigl(1+\frac{r}{\eta}\Bigr)\Bigl( 1+\frac{r}{\gamma}\Bigr).\]La aproximación del tiempo generacional para el modelo SIR es $T_g=1/\gamma$; para el modelo SEIR la aproximación resulta ser $T_g=1/\gamma+1/\eta$. Cabe señalar que el modelo SEIR proporciona una aproximación conceptual más precisa del intervalo serial pues incorpora tanto al período de incubación como al período de recuperación en su aproximación. Para ilustrar la naturaleza aproximada de los modelos de Kermack-Mckendrick para estimar ${\cal R}_0$, usemos un ejemplo práctico. Se sabe [NLA20] que el intervalo serial para el SARS-CoV-2, calculado directamente de los datos (usando tamizajes y seguimiento de contactos) tiene una mediana de 5.7 días; sabemos también que el período de incubación, $1/\eta$, de este virus es de aproximadamente 5.1 días y que el período de recuperación, $1/\gamma$, puede variar entre 14 y 25 días, dependiendo de la gravedad de la enfermedad. Supongamos ahora que el tiempo de doblamiento de la epidemia es de una semana lo que nos permite calcular de manera independiente del modelo la tasa de crecimiento, $r$, que resulta aproximadamente igual a $0.14$. Con estos datos, usando $1/\gamma =14$ días, obtenemos que con el intervalo serial supuesto en el modelo SIR ${\cal R}_0=2.96$ mientras que para el SEIR ${\cal R}_0=5.7$. Usando el intervalo serial verdadero obtenemos, para el SIR ${\cal R}_0=1.84$. La sobreestimación del intervalo serial conlleva la sobreestimación de ${\cal R}_0$ que a su vez sobreestima $r$, adelanta erróneamente la fecha del máximo de incidencia y sobreestima el tamaño total de la epidemia (ver Figura 3).
En el caso de la influenza usaré [Cow+09] para ilustrar el ejemplo. Para influenza usaremos también la información de la página web de los CDC que indica que el período de recuperación para esta otra enfermedad es $1/\gamma=6$ días y el de incubación es $1/\eta=2$ días en promedio. El intervalo serial medio es $T_g=3.9$ [Cow+09]. En consecuencia usando el modelo SIR ${\cal R}_0=1.84$; usando el SEIR ${\cal R}_0=2.3$ y usando $T_g$, ${\cal R}_0=1.54$. El rango de variación del número reproductivo en este caso es bastante más estrecho que en el caso del SARS-CoV-2 y, por tanto, la aproximación proporcionada por los modelos Kermack-McKendrick es mucho mejor que para el caso del COVID-19.
Un modelo es una aproximación de la realidad para algún propósito bien definido, y esto nunca debe olvidarse. La presente epidemia de SARS-CoV-2 nos ha permitido apreciar el poder de los modelos matemáticos como herramientas estratégicas para entender la ciencia básica de su dinámica poblacional y para diseñar políticas públicas de mitigación y control. En particular, los modelos han ayudado a entender el papel de los medidas de mitigación en el crecimiento no exponencial de muchas curvas epidémicas en el mundo, del papel de los procesos de superdispersión en el disparo de brotes epidémicos, la efectividad del uso de cubrebocas, de la naturaleza dinámica del intervalo generacional, del papel de la destrucción de hábitat en el derrame (spillover) de zoonosis sobre la población humana, de la mejor estimación del período de incubación usando mediciones virales de pacientes, de las presiones de selección que pudieran generar mutantes de escape de las nuevas vacunas. Esta lista es una pequeña (y sesgada) muestra de importantes problemas teóricos y prácticos en los que los modelos matemáticos han tenido una contribución significativa. Sin embargo, la historia reciente del manejo de la epidemia también ha confirmado un hecho fundamental: este uso estratégico de la modelación matemática en salud pública y medicina es una actividad multidisciplinaria que requiere de un juicio crítico en la interpretación de los modelos tanto en el significado de sus parámetros como en la naturaleza y alcance de sus proyecciones. Así como en el mundo ha habido historias notables de éxito en la aplicación de modelos para predecir las tendencias poblacionales de la COVID-19, así también ha habido notables fallos, muchos ciertamente provocados por el uso político de estas herramientas, pero otros teniendo realmente sus raíces en un conocimiento incompleto de la teoría básica de la biología de una epidemia, en una serie de generalizaciones e interpretaciones que las hipótesis de los modelos tipo Kermack-McKendrick no justifican, y, finalmente, en el abuso de la estadística priorizando la estimación de parámetros para el ajuste de datos, usando una amplia variedad de herramientas, pero olvidando que una epidemia es, ante todo, un fenómeno biológico que es expresión de la dinámica poblacional del patógeno y de la conducta social humana, factores ambos que le imprimen una alta variabilidad, una inherente heterogeneidad y una significativa incertidumbre en su expresión fenomenológica tanto en tiempo como en espacio.
Agradezco al Comité Editorial de Motivos Matemáticos la invitación hecha para escribir esta pieza y los comentarios de Fernando Saldaña sobre la primera versión de este trabajo. Finalmente, agradezco el apoyo financiero de los proyectos DGAPA-PAPIIT UNAM IV100220 (especial COVID-19) y DGAPA-PAPIIT IN115720.
