El objetivo de este escrito es presentar los resultados de un algoritmo de ajuste de un modelo epidemiológico SIR a la evolución de la pandemia del coronavirus COVID-19 con número reproductivo variable (R0t). Este ajuste se ha aplicado a la actual difusión de la epidemia en diversos Estados y comunidades autónomas: Italia, España, Francia, Gran Bretaña, Estados Unidos, Alemania, Corea del Sur, Madrid, Cataluña, Euskadi y Navarra. El algoritmo es fácilmente aplicable a otras zonas.
El modelo SIR utilizado es el modelo clásico, sin dinámicas vitales, suponiendo que la población total (N0=S+I+R+F) se mantiene constante, sin nacimientos ni muertes más allá de las derivadas de la epidemia (F). Los parámetros de propagación del modelo considerados han sido los siguientes, apoyados en la bibliografía (más detalles sobre el modelo y algoritmo utilizados al final del texto):
- Tiempo de infección (Tinf), desde la infección hasta la recuperación, de 20 días (30%, 14 días; 55%, 21 días; 10%, 25 días; 5%, >25 días); el paciente es vector de infección durante todo este tiempo.
- Tiempo de defunción (Tm), desde la infección hasta el fallecimiento, de 15 días.
- Tasa de letalidad (FR) del 1% (se ha hecho un análisis de sensibilidad considerando FR = 0.1%).
Existen modelos S(E)IR on-line, como el Epidemic Calculator, que permiten modelizar la evolución de la epidemia fijados los parámetros (Tinf, Tm, FR) y el número reproductivo básico (R0), permitiendo a veces fijar un segundo número reproductivo a partir de una supuesta intervención social orientada a limitar la propagación de la epidemia (confinamiento, distanciamiento social, uso de mascarillas, etc.).
El número reproductivo varía a lo largo del tiempo (Rt). Este número reproductivo es también proporcional a la proporción de infectados sobre la población total, y puede escribirse:
Rt = R0 · S/N
R0 es el número reproductivo al inicio de la epidemia, dado que en el instante inicial S = N.
El objetivo de las intervenciones es situar Rt por debajo de 1. Si Rt < 1, entonces la epidemia alcanza el pico de infectados, y éstos comienzan a reducirse, hasta desaparecer al cabo de un tiempo. Si R0 no varía, el pico de infectados en la epidemia se puede calcular aplicando la condición:
Rt = R0 · S/N = 1
fracción de susceptibles, S/N = 1/R0
fracción de afectados, (1-S/N) = 1 – 1/R0
Este umbral es el que proporciona lo que se conoce como inmunidad de grupo, que sin embargo no debe confundirse con la proporción de afectados al final de la epidemia, ya que conforme cae la curva de infectados, la epidemia continúa sumando afectados (infectados que acabarán recuperándose, o falleciendo). Puede comprobarse que la proporción de afectados totales una vez apagada la epidemia, cuando R0 no varía y Rt solo queda limitada por la inmunidad de grupo, puede aproximarse por:
afectados finales, (1-Sf/N) = 1 – 4 · exp(–1.5·R0)
Lograr la inmunidad de grupo no es una estrategia aceptable para hacer frente a una epidemia, a no ser que R0 sea muy cercano a 1 (si R0 = 1.4, el 51% de la población acabará infectándose; si la tasa de letalidad es del 1%, el 0.5% de la población acabará falleciendo. Con R0>2.4, los afectados superarán el 80% y los fallecidos el 0.8% de la población total: 8 mil personas por millón de habitantes). Por ello, es necesario realizar intervenciones que reduzcan el valor de R0 a lo largo del tiempo:
Rt = R0t · S/N
donde R0t es el número reproductivo que varía a lo largo del tiempo. Este número reproductivo variable es, en principio, desconocido. En esta trabajo hemos desarrollado un algoritmo de ajuste que, dados los parámetros Tinf, Tm y FR, y la evolución de los fallecidos en un área o región concreta, ajusta la evolución del número reproductivo básico, R0t, a lo largo del tiempo. Para ello se utiliza un algoritmo de optimización que minimiza el error entre los fallecidos proporcionados por el modelo y los fallecidos registrados en las estadísticas oficiales, p. e. aquí y aquí.
A modo de ejemplo, la siguiente figura muestra el resultado de la modelización para Italia (N = 60.48 millones de personas).
Último registro de fallecidos: 2020-04-18
Puede observarse en la gráfica de la izquierda que R0t es cercano a 9 antes de las intervenciones sucesivas en Italia, que se inician en algunas zonas del norte del país a finales de febrero, y no se extienden a toda Italia hasta el 9 de marzo. En las últimas semanas de marzo las medidas de contención parecen limitar la propagación de la epidemia, tal y como se muestra en la gráfica de la derecha. A continuación se muestran los resultados más significativos de la simulación realizada para Italia, alimentada con datos disponibles hasta el 18 de abril:
Día del primer infectado: 26 de enero (01-26)
Pico de infectados: 1.24 millones de personas, el día 31 de marzo (03-31).
Situación el último día registrado, 18 de abril (04-18): 0.88 millones de infectados (1.5% de la población) y 1.65 millones de personas recuperadas (2.7% de la población).
Desde el pico de infectados, estos se han reducido el 28.6% en 18 días, o a un ritmo del 1.9% diario.
Según la modelización, en Italia el pico de infectados se habría alcanzado (31 de marzo) tres semanas después de extender el confinamiento a todo el país, y un mes después de los primeros confinamientos localizados.
Las siguientes figuras muestran los resultados del ajuste y modelización para España, Francia, Reino Unido, Estados Unidos, Alemania, Corea del Sur, y las comunidades autónomas de Madrid, Cataluña, Euskadi y Navarra.
Los resultados de los ajustes están accesibles aquí.
La siguiente tabla muestra los principales resultados para las simulaciones realizadas, suponiendo una tasa de letalidad del 1%:
FR = 1%
Del análisis de estos datos pueden derivarse las siguientes conclusiones:
-
La epidemia presenta un pico de infectados al menos dos semanas después de la puesta en marcha de las medidas más estrictas de contención (Italia, España, Francia, CC.AA.).
-
USA todavía no habría alcanzado el pico de infección, o estaría en él, y UK lo habría superado hace pocos días (04-11).
-
Las medidas de contención podrían necesitar de una semana para reducir los infectados por debajo de los recuperados (04-05 en Italia, 04-10 en España).
-
Una vez alcanzado el pico, los infectados se van reduciendo, pero lentamente (la mayoría de las regiones con ritmos medios diarios de entre el 2% y el 4%). En Italia los infectados se han reducido el 29% en 18 días; en España el 42% en 16 días; en Madrid más de la mitad en 20 días.
-
Tras un mes de contención, los infectados descienden, pero se mantienen por encima del 1% de la población (más del 1.5% en España). En las zonas más afectadas los infectados pueden superar en estos momentos el 2% (Euskadi, 2.2%), e incluso el 3% (Madrid, 3.2%) y el 4% (Navarra, 4.5%).
-
En estos momentos la epidemia habría afectado ya (infectados + recuperados) a entre el 4% (Italia, España, Francia) y casi el 11% de la población (Madrid, 10.9%). En una posición intermedia se encuentran Cataluña (4.9%), Euskadi (5.4%) y Navarra (7.8%).
-
Alemania, pese a presentar un número de fallecidos más reducido que otros países, podría tener en estos momentos más de 350 mil personas infectadas, y casi 300 mil recuperadas.
-
En Corea del Sur el pico de infección solo habría afectado a 10 mil personas, el 0.02% de la población, y podría haber sufrido repuntes el último mes. En estos momentos podrían quedar todavía varios miles de infectados en el país.
En el caso de Euskadi y Navarra, si la tasa de letalidad se situara en el entorno del 1%, con las actuales medidas de confinamiento la cantidad de infectados continuaría situándose en el orden de decenas de miles de personas, y con un número reproductivo cercano a la unidad, lo que exigiría una muy cuidadosa y vigilada planificación del proceso de desescalado de la contención, para evitar futuros repuntes de la epidemia, que presentan un riesgo evidente. En Navarra se podría haber dado un repunte de la epidemia las últimas semanas; se necesitarían más datos de la serie para confirmarlo.
Análisis de sensibilidad. FR=0.1%
A continuación se comparan los resultados para Italia, España y Euskadi suponiendo tasas de letalidad del 0.1% y del 1%.
Si la tasa de letalidad fuera la décima parte de lo supuesto inicialmente, la penetración de la epidemia se multiplicaría por diez. Pero con una tasa de letalidad del 0.1% todavía existiría margen de avance para la epidemia y para posibles repuntes. La población susceptible se habría reducido casi a la mitad con respecto a la situación si la tasa fuera del 1%, pero la población infectada se habría multiplicado por 10, en una proporción mucho mayor. Con tanta población infectada las posibilidades de repunte seguramente aumentarían, al ser mucho más difícil el control de los infectados (entre el 15 y el 20% de la población se mantendría contagiosa el 18 de abril). Si R0t aumentara por encima de 1.4 (el valor que fija la inmunidad de grupo en el 50%), la epidemia repuntaría. Si R0t se situara por encima de 2.4, la epidemia superaría el 80% de la población, y la letalidad superaría las 800 personas por millón de habitantes.
Anexo. Modelo SIR(F) con R0 variable y ajustable
En nuestro modelo epidemiológico se fijan los siguientes parámetros:
N : población
FR : tasa de letalidad
Tinf : tiempo de infección
Tm : tiempo desde la infección hasta el fallecimiento
f0 : lista de la evolución de fallecimientos acumulados, a lo largo de los días
La epidemia tiene un comportamiento exponencial en el arranque, y la evolución de los fallecidos se puede ajustar con una regresión exponencial:
b,a = polyfit([i for i in range(n)], log(f0[0:n]), 1)
donde n son primeros valores de la serie (n = 10 para todos los países salvo Corea, con n = 30).
Se puede comprobar que, suponiendo comportamiento exponencial, el número de días que el primer infectado se adelanta al primer fallecido se calcula como:
dPrimerInfectado = 1/b·log(Tm/FR·b·exp(a))
Las ecuaciones en diferencias del modelo SIR(F) son las siguientes:
iS = -R0t[i]/Tinf·I[i]·S[i]/N[i] # incremento de susceptibles
iI = R0t[i]/Tinf·I[i]·S[i]/N[i] – I[i]/Tinf # incremento de infectados
iR = I[i]/Tinf # incremento de recuperados
iF = FR·I[i]/Tm # incremento de fallecidos
S[i+1] = S[i] + iS # susceptibles
I[i+1] = I[i] + iI # infectados
R[i+1] = R[i] + iR # recuperados
F[i+1] = F[i] + iF # fallecidos
N[i+1] = N[i] – iF # población
El algoritmo de minimización del error (F-f0) hace uso de scipy.optimize.minimize.
El código está accesible aquí.