Modelos compartimentales.

Los modelos compartimentales (SIS, SIR, SEIR y sus variantes) son formas simplificadas de representar la evolución de una epidemia en una población grande. Estos modelos separan la población en distintos estados o “compartimentos” y tratan de explicar la forma en la que los individuos de esa población evolucionan de unos a otros. En el caso del modelo SIR (también llamado “modelo de Kermack-McKendrick”) se asume que la población (de tamaño \(N\), constante) se divide en “susceptibles”, \(S(t)\), “infectados”, \(I(t)\) y “recuperados”, \(R(t)\). El modelo asume que la tasa de disminución de la proporción de individuos susceptibles es proporcional al producto de las proporciones de susceptibles e infectados. La constante de proporcionalidad es la “tasa de contagio”. También se asume que la tasa de paso de infectado a recuperado es constante y a tal constante se la conoce como “tasa de recuperación”. Los individuos recuperados adquieren “inmunidad”, es decir, no contagian ni pueden ser contagiados. Formalmente todo esto significa que las funciones \(S,I\) y \(R\) evolucionan de aucerdo con el sistema de EDO’s \[ \mbox{(SIR)}\quad\left\{ \begin{array}{ll} \dot{S}=&-\beta \frac{SI}N,\\ \dot{I}=&\beta \frac{SI}N-\gamma I,\\ \dot{R}=&\gamma I. \end{array} \right. \] El cociente entre la tasa de contagio y la tasa de recuperación es el parámetro conocido como \(R_0\) o “número de reproducción”. Este parámetro determina el comportamiento a largo plazo de la epidemia. Esencialmente, valores de \(R_0\) menores que uno conducen a un equilibrio libre de enfermedad, mientras que si ese valor es más que uno se alcanzará, a largo plazo, un equilibrio “endémico”, es decir, con una proporción no nula de individuos infectados. Una presentación completa de estas y otras propiedades del modelo SIR, así como de otros modelos compartimentales puede consultarse en [Li, M. (2018) “An Introduction to Mathematical Modeling of Infectious Diseases”, Springer].

Cuando se trata de modelar la evolución de una epidemia que, previsiblemente, puede tener un grave impacto sobre las condiciones de vida de una sociedad organizada, en principio, el modelo SIR parece demasiado simplista. Posiblemente se pueda aceptar que la tasa de recuperación se mantenga constante en el tiempo, pero en general las sociedades actúan para reducir las tasas de contagio, introduciendo medidas de aislamiento, o a través de vacunaciones, para conseguir erradicar o, al menos, limitar el impacto de la epidemia. Hay variantes del modelo SIR que tratan de adaptarse a esta realidad. Una posibilidad es introducir una función de “mitigación”, es decir, una función que modele la forma en que la tasa de contagio va disminuyendo en el tiempo, de manera que la evolución del sistema está dado por el sistema de EDO’s \[ \mbox{(SIRM)}\quad \left\{ \begin{array}{ll} \dot{S}=&-\beta(t) \frac{SI}N,\\ \dot{I}=&\beta(t) \frac{SI}N-\gamma I,\\ \dot{R}=&\gamma I, \end{array} \right. \] donde \(\beta(t)>0\) es una función que permita modelar de forma apropiada el descenso en la tasa de contagio. Hay muchas otras posibilidades, que permiten incluir parámetros demográficos (tasas de nacimiento, muerte, inmigración o emigración) o el efecto de vacunaciones. Con vistas a un estudio rápido de la evolución de la actual epidemia de COVID19 entendemos que, en principio, las variaciones demográficas a corto plazo tendrán un efecto despreciable. Por razones que trataremos de explicar mejor más adelante esto nos ha llevado a ignorar esas posibles generalizaciones. Usar modelos SIR con vacunación puede parecer inapropiado, puesto que no hay por el momento vacunas disponibles contra el COVID-19, pero, en realidad, el paso de “susceptible” a “recuperado” (que es lo que realmente se introduce en el modelo con vacunación) puede darse de forma inadvertida en los recuentos si, como creemos, una proporción importante de casos de COVID-19 no son detectados, por manifestarse únicamente en síntomas leves, y, quizá, puede producir la inmunización (el paso a “recuperados”) de una parte de la población. Podemos modelar esta transición directa mediante el sistema \[ \mbox{(SIRV)}\quad \left\{ \begin{array}{ll} \dot{S}=&-\beta \frac{SI}N-vS,\\ \dot{I}=&\beta \frac{SI}N-\gamma I,\\ \dot{R}=&\gamma I+vS, \end{array} \right. \] La \(v\) en la etiqueta anterior hace referencia a que modelamos esta forma natural de adquirir inmunidad, aunque advertimos que el modelo SIRV no es el modelo con vacunación presente en la literatura que hemos consultado sobre el tema. Todavía hay otras variantes. En el modelo SEIR se trata de incorporar la posibilidad de que los individuos que entran en contacto con infectados pasen un periodo de latencia antes de ser trasmisores de la enfermedad. A los individuos en este estado intermedio se le denomina “expuestos”. La dinámica del modelo SEIR se describe por las siguiente EDO’s: \[ \mbox{(SEIR)}\quad \left\{ \begin{array}{ll} \dot{S}=&-\beta \frac{SI}N,\\ \dot{E}=&\beta \frac{SI}N-\alpha E,\\ \dot{I}=&\alpha E-\gamma I,\\ \dot{R}=&\gamma I, \end{array} \right. \] En el modelo SEIR el parámetro \(\alpha\) o, más bien, \(\frac 1 \alpha\) representa el tiempo medio de recuperación de un infectado.

Los modelos anteriores pueden modificarse de forma que se tengan en cuenta varios de los efectos que queremos introducir. Por ejemplo, podemos considerar un modelo SIR con mitigación e inmunización natural al mismo tiempo. A este modelo lo denominaremos SIRMV. Igualmente podemos manejar un modelo SEIRM, con una tasa \(\beta(t)\) variable en el tiempo, o un modelo SEIRV en el que sea posible la transición del estado “expuesto” a “recuperado” sin pasar por el estado “infectado”.

Análisis de regresión no lineal con modelos compartimentales.

Con todo el atractivo que puedan tener los modelos compartimentales en términos de interpretabilidad, el ajuste a series de datos reales de cualquiera de ellos requiere la búsqueda de valores apropiados de los parámetros implicados. Una posible estrategia es la calibración de los parámetros a partir de estudios independientes. Por ejemplo, podemos estimar la tasa de contagio \(\beta\) a partir de los valores disponibles en estudios de epidemias similares o de la misma epidemia en otros lugares en los que su desarrollo esté más avanzado. Esta aproximación, sin embargo, es bastante problemática. En ocasiones los parámetros incluidos en un modelo actúan de formas aproximadamente proporcionales y se produce un fenómeno parecido al de la multicolinealidad en regresión lineal. Por otro lado los valores observados de los números de casos activos son enteros y parece más lógico interpretar que se trata de versiones perturbadas de la curva \(I(t)\). Por estas razones nos proponemos asjutar modelos de regresión no lineal sobre las series de casos activos, de forma que si \(\hat{I}_t\) representa el número de casos activos en el día \(t\) entonces \[\hat{I}_t=I_\theta(t)+\varepsilon_t,\quad t=1,\ldots,T,\] donde \(I_\theta(t)\) es la función dentro del modelo correspondiente con parámetro \(\theta\) y condiciones iniciales iguales a las observadas, mientras que \(\varepsilon_1,\ldots,\varepsilon_T\) son errores centrados e independientes. Por ejemplo, en el ajuste del modelo SIR tendremos \(\theta=(\beta,\gamma)\) y la función \(I_\theta(t)\) es la (parte correspondiente a la) solución del sistema (SIR) con parámetros \(\beta\) y \(\gamma\) y condición inicial \(S_{\theta}(0)=\hat{S}_0\), \(I_{\theta}(0)=\hat{I}_0\) y \(R_{\theta}(0)=\hat{R}_0\). Asumimos conocidos estos valores iniciales, igual que la serie de casos activos \(\hat{I}_t\).

La estimación de los parámetros puede hacerse con el método de mínimos cuadrados, tal como se describe en [Bates, D.M. y Watts, D. G. (1988) Nonlinear Regression Analysis and its Applications, Wiley]. El estimador de mínimos cuadrados es el minimizador de la suma de desviaciones cuadráticas entre los valores observados \(\hat{I}_t\) y los valores ajustados del modelo, es decir, \[\hat{\theta}=\underset{\theta}{\mbox{argmin}}V(\theta),\] donde \(V(\theta)=\sum_{t=1}^T (\hat{I}_t-I_\theta(t))^2\). El valor mínimo de \(V\) es la suma de cuadrados residual, \[\mbox{RSS}=V(\hat{\theta}).\] Bajo la hipótesis de varianza constante, \(\sigma^2\), de los errores el estimador \[\hat{\sigma}^2=\frac{\mbox{RSS}}{T-p}\] (donde \(p\) es el número de parámetros del modelo) es un estimador consistente de \(\sigma^2\). A partir de \(\hat{\theta}\) y \(\hat{\sigma}^2\) se pueden calcular regiones de confianza para \(\theta\) o regiones de predicción para la respuesta.

Estudio de la evolución de la epidemia COVID-19 en algunas provincias chinas a a partir de modelos compartimentales. Comparación con el caso de Castilla y León.

Estudiamos a continuación la evolución de la epidemia COVID19 en lugares en los que se encuentra en un estadio más avanzado desde el punto de vista de los modelos compartimentales. En particular nos fijamos en los datos de algunas provincias chinas. Las siguientes gráficas muestran la evolución de los grupos de casos acumulados (azul), recuperados (verde) y muertos (rojo). En algunas de ellas se observa un corte en laa series de casos totales (al menos en Gansu y Hubei). Otras corresponden a provincias con poca población (el caso de Macau) o con pocos casos y evolución incompleta. Por ello nos centramos en las provincias con mayor número de casos con la excepción de Hubei, es decir, Anhui, Chongquing, Guangdong, Henan, Hunan, Jiangxi y Zhejiang.

En los modelos aplicados el tamaño de a población se mantiene constante en el tiempo. La expansión del COVID-19 está siendo muy rápida y por ello no parece que la hipótesis de tamaño de población constante suponga una grave desviación de la realidad. No incluimos por ello tasas de nacimiento o muerte. Además, en nuestro estudio contamos todas las salidas del estado “infectado como recuperados”, aunque una parte de ellas se correspondan con fallecimientos, pero, en todo caso, conducen a la situación de no trasmisor de la enfermedad.

Para las provincias elegidas hemos probado los modelos compartimentales descritos en la sección anterior. En los modelos que incluyen función de mitigación, ésta es tipo logístico, es decir, de la forma \[\beta(t)=\frac{\beta}{1+e^{\mu(x-x_0)}},\quad t\geq 0.\] Los parámetros \(\beta\), \(\mu\) y \(x_0\) se estiman a partir de los datos. El parámetro \(\beta\) representa la tasa de trasmisión antes de aplicar medidas de aislamiento, mientras que \(\mu\) tiene que ver con la velocidad de descenso en la tasa de transmisión y \(x_0\) es el punto central en el proceso de decrecimiento (si el decrecimiento es rápido \(x_0\) marca el “punto de cambio”). En el ajuste del modelo SIRM el valor estimado de \(x_0\) es 0 en todas las provincias en estudio, por lo que hemos incluido también en nuestro estudio el modelo SIRM0 en el que este modelo está fijado en 0, al igual que en el resto de modelos con mitigación considerados. La siguiente tabla resume la suma de cuadrados residual correspondiente al ajuste de cada modelo en cada provincia.

## SIR SIRV SIRM SIRM0 SIRMV SEIR SEIRV SEIRM
## Anhui 1534063 62965.0 86125.3 86125.3 62965.0 8299048 1540102 42028.4
## Chongqing 479192 18143.6 25108.1 25108.1 18143.6 2436145 516055 11897.8
## Guangdong 1812950 57175.1 60904.8 60904.8 57175.1 11706428 1908532 53050.0
## Henan 1632947 89088.7 117814.4 117814.4 89088.7 11233035 1672684 54556.6
## Hunan 1361477 20269.4 34238.3 34238.3 20269.4 6563321 1404415 14175.8
## Jiangxi 1133859 46182.6 63404.1 63404.1 46182.5 7391570 1191529 28227.0
## Zhejiang 2271713 71555.8 97113.2 97113.2 71555.8 10953662 2349289 73915.6

Aunque es difícil valorar la diferencia cuantitativa entre valores para los distintos modelos, es claro que el ajuste de los modelos SIR, SEIR y SEIRV es mucho pero que el de los demás modelos. Entre los modelos de tres compartimentos el modelo SIRV es el ganador en todas las provincias. El modelo SIRM0 no pierde nada frente al modelo SIRM (en el que se incluye \(x_0\) como parámetro libre). SIRM0 funciona algo peor que SIRV, pero el ajuste es razonable (luego volvemos sobre esta afirmación). El modelo SIRMV incluye como submodelos a SIRV y a SIRM. Lo que nos dice la tabla es que si el efecto de inmunización natural está incluido en el modelo entonces la función de mitigación no ayuda a mejorar el ajuste a los datos. En resumen, a la vista de la tabla los modelos de dos compartimentos que mejor describen la evolución de la epidemia son SIRV y SIRM0, por ese orden.

¿Qué significa buen ajuste o mal ajuste con estos datos? No hemos usado ningún procedimiento formal de contraste de ajuste de los modelos a estos conjuntos de datos y no creemos que sea posible extraer ninguna conclusión especialmente fuerte por esta vía con series de esta longitud. Pero una comparación gráfica de los datos de casos activos (\(\hat{I}_t\)) frente a los modelos ajustados (\(I_{\hat{\theta}} (t)\)) puede servir para tener una idea aproximada. Mostramos algunos de estos gráficos a continuación.