TRABAJO DE FIN DE GRADO
CAMBIO DE LA CRITICIDAD EN MODELOS EPIDÉMICOS DEBIDO AL CONFINAMIENTO
Jaime García Lozano
Grado de Física Facultad de Ciencias
Año Académico 2020-21
CAMBIO DE LA CRITICIDAD EN MODELOS EPIDÉMICOS DEBIDO AL CONFINAMIENTO
Jaime García Lozano
Trabajo de Fin de Grado Facultad de Ciencias
Universidad de las Illes Balears
Año Académico 2020-21
Palabras clave del trabajo:
Modelos epidémicos, confinamiento, criticidad, modelo determinista, modelos estocásticos, SIR, SEIR
Nombre Tutor/Tutora del Trabajo Jose Javier Ramasco Sukia
Se autoriza la Universidad a incluir este trabajo en el Repositorio Institucional para su consulta en acceso abierto y difusión en línea, con fines exclusivamente académicos y de investigación
Autor Tutor Sí No Sí No
☒ ☐ ☒ ☐
Resumen
Este trabajo, motivado por la pandemia del COVID-19, pretende ser una introducción al modelo epidemiológico SIR y a su variante SEIR. Se desarrollarán ambos modelos tanto desde una visión determinista como estocástica y con ellos se analizarán los posibles regímenes en los que puede estar un brote epidémico, haciendo especial hincapié en el llamado régimen crítico, que puede aparecer al aplicar un confinamiento parcial en la población.
Índice
1. Introducción 3
2. Modelo determinista 4
2.1. SIR . . . 4
2.1.1. Ecuaciones . . . 5
2.1.2. Resolución numérica . . . 6
2.1.3. Posibles estados . . . 7
2.1.4. Confinamiento . . . 8
2.2. SEIR . . . 9
3. Modelos estocásticos 11 3.1. Gillespie . . . 11
3.1.1. SIR . . . 12
3.1.2. SEIR . . . 12
3.1.3. Confinamiento . . . 13
3.2. Binomiales . . . 13
3.2.1. SIR . . . 14
3.2.2. SEIR . . . 14
3.2.3. Confinamiento . . . 15
3.2.4. Relación entre probabilidad y tasa . . . 16
4. Resultados 17 4.1. Caso general . . . 17
4.2. Caso con confinamiento . . . 19
4.3. Diagrama de fase . . . 20
5. Limitaciones de los modelos 21 6. Conclusión 22 7. Apéndices 23 7.1. Cálculo de R0 mediante la Next Generation Matrix . . . 23
7.2. Método de ”binning” aplicado a los datos del Gillespie . . . 24
1. Introducción
El estudio y la comprensión de la evolución de una epidemia es un tema de gran interés, pues es un fenómeno que, en caso de descontrolarse, puede tener consecuencias sociales y económicas devastadoras. Se ha podido ver a lo largo de la historia con casos como la Peste Negra, la Gripe Española, el Ébola o más recientemente el COVID-19.
A fin de encontrar patrones en las epidemias y plantear estrategias de contención para mitigar sus efectos en la población, la elaboración de modelos epidemiológicos resulta de gran utilidad.
El tipo de epidemia que se tratará en este trabajo será la causada por una enfermedad infecciosa, microparasitaria y de transmisión directa ([1], capítulo 2). Ésta última sucede cuando la enfermedad se propaga por contacto cercano entre un infectado y una persona sana. En general, las enfermedades microparasitarias son de transmisión directa y sus patógenos no pueden sobrevivir durante mucho tiempo fuera de un huésped.
Cuando una epidemia de estas características está descontrolada, el número de infec- tados crece exponencialmente hasta que se alcanza la inmunidad de grupo, lo cual sucede cuando el patógeno no encuentra más individuos a los que transmitirse. A partir de ese punto, las infecciones empiezan a decaer hasta desaparecer. Estas dos etapas correspon- den al llamado régimen supercrítico. Si, por el contrario, desde un inicio el patógeno no consigue propagarse (debido a una baja infectividad, por ejemplo) la epidemia estará en un régimen subcrítico y dejará de haber contagios al cabo de poco tiempo. El número re- productivo básicoR0 es un parámetro que puede indicar el régimen en el que se encuentra el sistema, ya que representa el número de infecciones secundarias que surgen a partir de un caso primario. Si su valor es mayor que uno, cabe esperar que la enfermedad se propa- gue [2]. Es por ello que las medidas de contención buscan disminuir su valor efectivo para forzar al sistema a entrar en el régimen subcrítico reduciendo así la duración del brote y las posibles víctimas. Sin embargo, puede aparecer un tercer estado que se caracteriza por una evolución prácticamente constante de la curva de infectados durante un largo periodo de tiempo. Se trata del régimen crítico.
Figura 1: Evolución de las infecciones diarias en España durante la pandemia del COVID- 19 desde el 15 de marzo del 2020 hasta el 23 de mayo de 2021 [3]. Nótese que el confi- namiento más restrictivo se aplicó desde el 15 de marzo de 2020 hasta el 21 de junio de 2020 .
En la Figura 1 podemos ver los sucesivos brotes que han tenido lugar en España durante la pandemia del COVID-19 con la forma característica del régimen supercrítico:
un ascenso exponencial hasta que las medidas de contención hacen efecto seguido de un pronunciado descenso hasta el relajamiento de las restricciones donde vuelven a subir los contagios.
Durante el primer brote se adoptaron las medidas más drásticas al someterse a gran parte de la población a un confinamiento prácticamente total, lo que dio lugar a un apla- namiento de la curva de infecciones manteniéndose prácticamente constante hasta unos meses después del relajamiento de las restricciones (véase Figura 1). Estudios sugieren que este comportamiento corresponde al régimen crítico [4] y para analizarlo utilizan el modelo epidemiológico SIR (Susceptibles-Infectados-Recuperados) descrito por Kermack y McKendrick en 1927 [5].
En resumen, este trabajo tiene varios objetivos:
Por una parte pretende servir como introducción al modelo SIR y su variante SEIR.
Se deducirán las ecuaciones del modelo y se mostrará detalladamente su resolución numérica. Después de este análisis hecho desde una visión puramente determinista, se resolverá el mismo modelo pero desde un planteamiento estocástico. Se utilizará el algoritmo de Gillespie y su aproximación, el de las Binomiales (que reduce enor- memente la carga computacional), a fin de poder ver el efecto de la aleatoriedad en la evolución del brote.
Finalmente, se hará un estudio del régimen crítico donde se analizarán los posibles mecanismos que lo hacen aparecer. Para ello, se observará el efecto de un con- finamiento parcial aplicándolo tanto en el modelo determinista como en los dos estocásticos previamente mencionados.
2. Modelo determinista
2.1. SIR
Se considerará una única población sin demografía (no hay nacimientos ni muertes), una aproximación válida para infecciones cuyos tiempos de evolución sean mucho más pequeños que los de cambio poblacional. Además, se asumirá que el patógeno provoca síntomas durante un periodo corto (días o semanas) seguido de la inmunidad. De esta manera, la población se divide en tres grupos: los susceptibles serán aquellos que aún no han pasado la enfermedad, los infectados los portadores del patógeno y los recuperados los que hayan superado la enfermedad y queden inmunizados. Un individuo podrá pasar de un grupo a otro según el siguiente esquema:
Figura 2: Esquema de los compartimentos del modelo SIR.β yµson las tasas de infección y de recuperación respectivamente.
Un infectado se trasladará al grupo de recuperados tras haber superado el tiempo pro- medio en el que el patógeno se mantiene activo, mientras que la transición de susceptible a infectado dependerá de tres factores: de la cantidad de infectados que haya en ese mo- mento, de la interacción que puedan tener los individuos entre ellos y de la probabilidad de infección tras el contacto con un portador del patógeno.
2.1.1. Ecuaciones
Asumimos que todos los individuos interactúan igual que el resto en un escenario de población mezclada (”mean field”). Así, definimos k como la cantidad de contactos que puede tener una persona con otros individuos por unidad de tiempo,i como el número de infectados y N como el de individuos.
En un intervalo de tiempo pequeño entre t y t+δt el número de contactos con in- fectados que tendrá un individuo susceptible será kNi δt, siendo Ni la probabilidad de que un contacto esté infectado. Si definimos ccomo la probabilidad de que haya transmisión tras un contacto, entonces 1−c, será la probabilidad de que no haya infección. Así, la probabilidad de que un susceptible no se infecte tras haber tenido kNi δt contactos con infectados vendrá dada por:
1−q= (1−c)kNiδt (2.1)
De donde q es la probabilidad de que sí se infecte.
Definimosβ =−kln(1−c). Sustituyendo en (2.1) y despejando q nos queda:
q= 1−e−βNiδt (2.2)
Si ahora se desarrolla la exponencial y se coge el límite δt→0:
q ≈βdt i
N (2.3)
La expresión (2.3) representa la probabilidad per cápita de ser infectado. Por tanto, si hay s individuos susceptibles, su variación vendrá dada por:
ds ≈ −s·q→ ds
dt =−sβ i
N (2.4)
De donde β es la tasa de infección y el signo (-) nos indica que el número de susceptibles tiene que decrecer con el tiempo.
El sistema de ecuaciones diferenciales que determinará la evolución temporal de cada uno de los tres grupos, en función de sus respectivas proporciones (S = s/N, I = i/N, R =r/N), vendrá dado por:
dS
dt =−βSI (2.5)
dI
dt =βSI−µI (2.6)
dR
dt =µI (2.7)
Dondeµes la tasa de recuperación (su inverso la duración media de la infección). Nótese que S+I +R = 1.
Se trata de un sistema de ecuaciones diferenciales no lineales que resolveremos numé- ricamente.
2.1.2. Resolución numérica
El método de Runge-Kutta-Fehlberg nos permite resolver numéricamente ecuaciones diferenciales ordinarias de primer orden. Si reescribimos las ecuaciones del SIR (Ecs (2.5)- (2.7)):
f(S, I) =−βSI (2.8)
g(S, I) = βSI−µI (2.9)
p(I) =µI (2.10)
Aplicando el método iterativo de quinto orden, la proporción de susceptibles, infectados y recuperados para cada paso de tiempo vendrá dada por:
Sn+1 =Sn+h
16
135f0+ 6656
12825f2+ 28561
56430f3− 9
50f4+ 2 55f5
(2.11) In+1 =In+h
16
135g0+ 6656
12825g2+ 28561
56430g3− 9
50g4+ 2 55g5
(2.12) Rn+1 =Rn+h
16
135p0+ 6656
12825p2+28561
56430p3− 9
50p4+ 2 55p5
(2.13) De donde la h es el paso de integración (h = 0,1 en nuestro caso) y las f’s, g’s y p’s son las funciones intermedias:
f0 =f(S0, I0) g0 =g(S0, I0) p0 =p(I0)
(2.14)
rk1(a, b) =a+ h 4b
f1 =f(rk1(S0, f0), rk1(I0, g0)) g1 =g(rk1(S0, f0), rk1(I0, g0)) p1 =p(rk1(I0, g0))
(2.15)
rk2(a, b, c) =a+ 3h
32b+9h 32c
f2 =f(rk2(S0, f0, f1), rk2(I0, g0, g1)) g2 =g(rk2(S0, f0, f1), rk2(I0, g0, g1)) p2 =p(rk2(I0, g0, g1))
(2.16)
rk3(a, b, c, d) =a+1932h
2197 b−7200h
2197 c+7269h 2197 d f3 =f(rk3(S0, f0, f1, f2), rk3(I0, g0, g1, g2)) g3 =g(rk3(S0, f0, f1, f2), rk3(I0, g0, g1, g2)) p3 =p(rk3(I0, g0, g1, g2))
(2.17)
rk4(a, b, c, d, e) = a+439h
216 b−8hc+ 3680h
513 d− 845h 4104e f4 =f(rk4(S0, f0, f1, f2, f3), rk4(I0, g0, g1, g2, g3)) g4 =g(rk4(S0, f0, f1, f2, f3), rk4(I0, g0, g1, g2, g3)) p4 =p(rk4(I0, g0, g1, g2, g3))
(2.18)
rk5(a, b, c, d, e, j) =a−8h
27b+ 2hc−3544h
2656 d+ 1859h
4104 e− 11h 40 j f5 =f(rk5(S0, f0, f1, f2, f3, f4), rk5(I0, g0, g1, g2, g3, g4))
g5 =g(rk5(S0, f0, f1, f2, f3, f4), rk5(I0, g0, g1, g2, g3, g4)) p5 =p(rk5(I0, g0, g1, g2, g3, g4))
(2.19)
Donde las ecuaciones se han obtenido a partir de [6], capítulo 5.
Mediante las condiciones iniciales I0, R0 y S0 = 1−I0 −R0 se podrá conseguir la evolución temporal de S, I y R.
Es importante destacar que no hay ningún punto en el queI = 0, debido que la solución es asintótica (el modelo está pensado para poblaciones infinitas). Por tanto, detendremos la iteración en el momento en el que la proporción de infectados sea considerablemente pequeña (en I = 10−5, por ejemplo) (Véase Fig. 3.2).
2.1.3. Posibles estados
Reescribimos la ecuación (2.6) para ver el efecto de las condiciones iniciales sobre la evolución de la curva de infectados:
dI
dt =Iβ(S− µ
β) (2.20)
De esta manera, el sistema puede comenzar en tres regímenes distintos:
SiS0 > µ/β, entoncesdI/dt >0, que corresponde a la fase de expansión del régimen supercrítico, donde I(t) crece exponencialmente hasta que se alcanza la inmunidad de grupo. Se llega a este punto cuando la proporción de susceptibles ha disminuido lo suficiente como para que S(t) < µ/β y, en consecuencia, dI/dt < 0, pasando así a la fase de decrecimiento donde la infección acaba desapareciendo. Por tanto, la inmunidad de grupo se alcanzará cuando S =µ/β.
El estado subcrítico puede aparecer desde un inicio si S0 < µ/β. Como dI/dt <0, no se alcanzará ningún máximo y el patógeno desaparecerá al cabo de poco tiempo.
Finalmente, siS0 =µ/β,I(t)=constante, se tiene el régimen crítico, donde la curva de infectados se mantiene horizontal.
En resumen, la proporción de susceptibles iniciales debe superar el umbral µ/β para que el patógeno pueda propagarse. Este valor representa la tasa de eliminación relativa, mientras que su inverso se conoce como el número reproductivo básico R0 y representa el número promedio de casos secundarios que surgen a partir de un caso primario. Como inicialmente S0 ≈1, el valor de R0 también nos indica el régimen en el que se encuentra la epidemia: R0 > 1 (dI/dt > 0) supercrítico, R0 = 1 (dI/dt = 0) crítico y R0 < 1 (dI/dt <0) subcrítico (Véase Ec. (2.20)).
3.1: Inicio de la curva de infectados para distin- tos valores deR0y un valor fijo deµ= 0,16.
3.2: Evolución completa de los tres grupos del SIR. Se ha utilizado unaR0= 2,5 yµ= 0,16.
Figura 3: Ejemplos de las solución determinista en diferentes regímenes, para una pobla- ción de N = 1000 individuos con i0 = 10 infectados iniciales.
En la Figura 3.1 se enseña cómo aumenta o disminuye la proporción de infectados dependiendo del valor de R0. En los casos supercríticos, cuanto mayor es R0 más pro- nunciado es el crecimiento y antes se llega al pico de infectados, mientras que en los subcríticos, cuanto más pequeño es R0 más rápido decae y muere la infección.
Por otra parte, la Figura 3.2 muestra la evolución completa de un caso supercrítico.
Nótese que el máximo de infectados sucede cuando S = 0,4 que, en efecto, se trata del punto en el que se alcanza la inmunidad de grupo, donde S =µ/β = 1/R0 = 0,4.
2.1.4. Confinamiento
Las medidas de contención , como la vacunación o el confinamiento, actúan variando el número de individuos en cada uno de los grupos que hemos definido en nuestro modelo. La primera traslada a individuos susceptibles directamente al grupo de recuperados, mientras que la segunda disminuye el número de personas que participan en el sistema. Este trabajo se va a centrar en el segundo caso. Modificaremos el modelo SIR con el fin de intentar aproximarnos al efecto de un confinamiento parcial sobre la evolución de la epidemia.
Así, definimos Xs ∈ (0,1) como la fracción de la población que ha sido confinada y Pth (Prevalence threshold [4]) como el umbral de la curva de infectados a partir del cual se aplica esta medida.
El confinamiento dificulta al patógeno a propagarse debido a que parte de los indivi- duos susceptibles están aislados. Si planteamos las ecuaciones del SIR teniendo en cuenta esta situación, quedan:
dS
dt =−βS(1−Xs)I (2.21) dI
dt =βS(1−Xs)I−µI (2.22)
dR
dt =µI (2.23)
Reescribimos la ecuación (2.22):
dI
dt =Iβ(1−Xs)
"
S− µ
β(1−Xs)
#
(2.24)
El valor del número reproductivo básico efectivo será entonces:
R0ef = β(1−Xs)
µ =R0(1−Xs) (2.25)
Por tanto, el confinamiento aplicado en un régimen supercrítico puede tener varios efectos: si es lo suficientemente fuerte como para que R0ef < 1 , entonces dI/dt < 0 (Ec.
(2.24)), por lo que se conseguirá pasar a un régimen subcrítico y que la curva de infectados entre en fase descendente. Si, por el contrario,Ref0 >1, la inmunidad de grupo se alcanzará antes que en el caso en el que no se hubieran impuesto medidas de contención, ya que Ref0 < R0 y, en consecuencia 1/Ref0 > 1/R0. Por último, se puede dar la situación en la que Ref0 = 1, que daría lugar al régimen crítico. A partir de la expresión (2.25) se puede estimar el valor de Xs que haría que el sistema entrara en este régimen:
Xscrit= 1− µ
β = 1− 1
R0 (2.26)
En la práctica, se integrarán las ecuaciones imponiendo un caso supercrítico dejando evolucionar la solución hasta que la fracción de infectados alcance el umbral Pth, a partir del cual se impondrá el confinamiento y se detendrá la iteración de la misma manera que en el caso sin confinamiento (cuando I <10−5) (véase Figura 4).
Figura 4: Evolución de la curva aplicando diferentes confinamientos a partir de un mismo umbral Pth= 0,1. Escenario con N = 1000,i0 = 10, µ= 0,16 yβ = 2µ.
2.2. SEIR
Algunos patógenos, como la gripe o el COVID-19, suelen requerir un periodo de incu- bación en el que el huésped no presenta síntomas ni puede infectar al resto de individuos.
El modelo SEIR incluye un nuevo grupo formado por los que se han infectado y se en- cuentran en este periodo latente: los expuestos.
Figura 5: Esquema de los compartimentos del modelo SEIR. β, γ y µ son las tasas de infección, exposición y recuperación respectivamente.
γ es la tasa de exposición y su inverso el periodo medio de exposición. La transición de ”susceptible” a ”expuesto” tendrá la misma dinámica que en el SIR la de ”susceptible”
a ”infectado”, mientras que la de ”expuesto” a ”infectado” sucederá después de que el individuo en cuestión supere el periodo de incubación promedio 1/γ.
De esta manera, el sistema de ecuaciones diferenciales, con las mismas aproximaciones mencionadas en el SIR y teniendo en cuenta el confinamiento, viene dado por:
dS
dt =−βS(1−Xs)I (2.27) dE
dt =βS(1−Xs)I−γE (2.28)
dI
dt =γE−µI (2.29)
dR
dt =µI (2.30)
6.1: Modelo SIR. 6.2: Modelo SEIR.
Figura 6: Comparación SIR-SEIR, Escenario con N = 106, i0 = 10, e0 = 0, µ=γ = 0,16 y β = 2µ.
En la Figura 6 podemos ver cómo ambos modelos presentan una dinámica bastante distinta. El SIR muestra un máximo de infectados más alto y se alcanza antes que en el SEIR, lo que implica que el brote muera más rápidamente.
Podríamos resumir el efecto de añadir el compartimento de los expuestos en que la incidencia máxima disminuye, la duración del brote aumenta y el número de individuos que han pasado la enfermedad tras acabar el brote queda un poco por debajo que en el SIR.
Por otra parte, cuanto más grande esγ, más pequeño es el periodo medio de exposición 1/γ, por lo que en el límite γ → ∞, recuperamos el SIR.
3. Modelos estocásticos
Hasta ahora se ha planteado el modelo desde una visión totalmente determinista:
cada conjunto de parámetros iniciales siempre van a dar lugar a la misma solución. En la realidad, sin embargo, existe una fuerte componente de aleatoriedad. Si pudiéramos repetir una epidemia bajo las mismas condiciones iniciales, no tendrían por qué infectarse los mismos individuos en el mismo orden que la vez anterior. Para parametrizar dicha aleatoriedad, se podría añadir una función de ruido en las ecuaciones deterministas, pero este trabajo se limitará a plantear dos modelos estocásticos diferentes con individuos discretizados.
En los modelos estocásticos las fluctuaciones tienen un papel fundamental. Se vio cómo la solución determinista es asintótica y en ningún momento los infectados llegan a ser cero. En el caso estocástico, sin embargo, detenemos la iteración en el momento en el que se cumple esta condición, por lo que el hecho de que I = 0 en un punto determinado es únicamente consecuencia de las fluctuaciones.
Para cada simulación, el modelo estocástico da una solución diferente (véase Figura 7). Si repitiéramos el proceso infinitas veces con una población suficientemente grande y lo promediáramos, obtendríamos la solución determinista, pues esta es el promedio de todos los resultados estocásticos en el límite de una población infinita.
Figura 7: Comparación entre la solución determinista y tres posibles resultados del modelo de las Binomiales. Parámetros utilizados: i0 = 10; r0 = 0; µ= 0,16; β = 2µ; N = 1000.
3.1. Gillespie
El Método Directo de Gillespie [7] genera una posible solución para un sistema de ecuaciones estocástico en el que las tasas de transición entre estados son conocidas. Se basa en estimar el tiempo hasta el próximo evento (3.1) a partir de la tasa total acumulada ωT (la suma de las tasas de todos los posibles eventos) y de un número aleatorio ˆu entre
0 y 1:
∆ˆt=− 1
ωTln (ˆu) (3.1)
Donde el valor de ∆ˆt está extraído de una distribución exponencial de tiempos, lo que implica que el tiempo entre eventos está descrito por la dinámica de Poisson.
Por otra parte, la probabilidad de que un evento con tasaωj ocurra será:
P = ωj
ωT (3.2)
3.1.1. SIR
En el SIR las infecciones y recuperaciones ocurren con una tasa de iβNs eiµ respecti- vamente.
Escribimos el modelo en forma de algoritmo:
1. Introducimos las condiciones iniciales: i0 y r0 (s0 =N −i0−r0).
2. Calculamos la tasa total:
ωT =iµ+iβ s
N (3.3)
3. Si i≤0 ó ωT = 0, detenemos el bucle.
4. A partir de un número aleatorio ˆu∈[0,1] calculamos con (3.1) el tiempo que pasará hasta el siguiente evento.
5. Obtenemos la probabilidad de recuperación p con (3.2) y generamos otro número aleatorio ˆu∈[0,1].
6. Si ˆu≤phabrá una recuperación:
it+ ∆ˆt=i(t)−1 rt+ ∆ˆt=r(t) + 1 7. En caso contrario habrá una infección:
st+ ∆ˆt=s(t)−1 it+ ∆ˆt=i(t) + 1
8. Por último, hacemos avanzar el tiempo t→t+ ∆ˆt y volvemos al paso 2.
3.1.2. SEIR
Ahora se ha de añadir el evento de exposición al caso anterior, con tasa γe:
1. Introducimos las condiciones iniciales: i0,e0 y r0 (s0 =N −i0−r0−e0).
2. Calculamos la tasa total:
ωT =iµ+iβ s
N +γe (3.4)
3. Si i≤0 ó ωT = 0, detenemos el bucle.
4. A partir de un número aleatorio ˆu∈[0,1] calculamos con (3.1) el tiempo que pasará hasta el siguiente evento.
5. Calculamos la probabilidad de infección p1 con (3.2)
6. La probabilidad de pasar la etapa de exposiciónp2 será entonces p2 = ωγe
T +p1 7. Obtenemos otro número aleatorio ˆu∈[0,1].
8. Si ˆu≤p1 habrá una infección:
s(t+ ∆ˆt) =s(t)−1 et+ ∆ˆt=e(t) + 1 9. Si p1 ≤uˆ≤p2 un expuesto pasará a ser infeccioso:
e(t+ ∆ˆt) =e(t)−1 it+ ∆ˆt=i(t) + 1 10. En otro caso habrá una recuperación:
it+ ∆ˆt=i(t)−1 r(t+ ∆ˆt) =r(t) + 1
11. Por último, hacemos avanzar el tiempo t→t+ ∆ˆt y volvemos al paso 2.
3.1.3. Confinamiento
Tanto en el SIR como en el SEIR el efecto del confinamiento será una reducción de la tasa de infección. Así, si confinamos una fracción Xs de la población, ésta quedará como βNis(1−Xs).
Es sencillo modificar los algoritmos previamente explicados para representar esta si- tuación. Simplemente se ha de añadir la condición de que si la fracción de infectados supera o iguala un umbral Pth, entonces la tasa de infección se modifica multiplicándola por (1-Xs).
3.2. Binomiales
El modelo de Gillespie requiere un trabajo computacional considerable. A medida que se aumenta la población, el tiempo de ejecución del programa crece linealmente ([1], capí- tulo 6). Cuantos más susceptibles haya, más sucesos de infección o recuperación ocurrirán, provocando una disminución del ∆ˆt, lo que se traduce en el cálculo de más puntos y en que el tiempo avance más lentamente.
El modelo que se explicará a continuación [8] es una aproximación del Gillespie que consigue reducir la carga computacional. Ahora el tiempo será una variable discreta y en cada uno de sus puntos se calculará el número de nuevos susceptibles, infectados y recuperados.
3.2.1. SIR
Se considera que hay dos sucesos en un instante de tiempo concreto: que un número
∆i de susceptibles se infectan y que un número ∆r de infectados se recuperan. A cada paso de tiempo la evolución del número de individuos en cada uno de los tres grupos vendrá descrita por el siguiente esquema:
s→s-∆i i→i+∆i-∆r
r→r+∆r
Definimos Pβ y Pµ como la probabilidad de infección y recuperación por posible even- to respectivamente. Si tenemos s susceptibles e i infectados, la probabilidad de que un individuo se infecte tras haber tenidokcontactos seráP = NiPβk. En otras palabras, tene- mos una distribución binomial con n =S ensayos independientes y con una probabilidad p= Ni Pβk de que la infección sea exitosa. Por tanto, la cantidad de nuevos infectados ∆i será un número aleatorio extraído de dicha distribución. De esta manera, la probabilidad de infección irá variando en función de la cantidad de infectados que haya. Por otra parte, los nuevos recuperados se extraerán de una binomial con n =i y p=Pµ.
Escrito en forma de algoritmo:
1. Introducimos las condiciones iniciales: i0 y r0.(s0 =N −i0−r0).
2. Calculamos el número de nuevos infectados y recuperados:
∆i=Binomial(n =s, p= Ni Pβk)
∆r=Binomial(n=i, p=Pµ) 3. Redefinimos las variables:
s(t+δt) =s(t)−∆i i(t+δt) = i(t) + ∆i−∆r
r(t+δt)=r(t) + ∆r 4. Si i≤0, detenemos el bucle; si no, se vuelve al paso 2.
Se recuperará el Gillespie si los intervalos temporales son muy pequeños, donde las binomiales pasarían a ser poissonianas y los tiempos entre eventos exponenciales.
3.2.2. SEIR
En el caso de considerar que existe un tiempo de incubación en el que el huésped no infecta, se deberá modificar parcialmente el esquema anterior añadiendo el grupo de expuestos:
s→s-∆e e→e+∆e−∆i
i→i+∆i-∆r r→r+∆r
Definimos Pγ como la probabilidad de que un individuo expuesto pase al grupo de infectados. Los nuevos expuestos se calcularán de la misma manera que los infectados en el SIR, mientras que el incremento de estos últimos vendrán dados por una distribución binomial con n=e y p=Pγ.
El algoritmo será:
1. Introducimos las condiciones iniciales: e0, i0 y r0 (s0 =N −e0−i0−r0).
2. Calculamos el número de nuevos expuestos, infectados y recuperados:
∆e=Binomial(n =s, p= Ni Pβk)
∆i=Binomial(n=e, p=Pγ)
∆r=Binomial(n=i, p=Pµ) 3. Redefinimos las variables:
s(t+δt) =s(t)−∆e e(t+δt) =e(t) + ∆e−∆i
i(t+δt) = i(t) + ∆i−∆r r(t+δt)=r(t) + ∆r 4. Si i≤0, detenemos el bucle; si no, se vuelve al paso 2.
3.2.3. Confinamiento
En este modelo el confinamiento se aplica exactamente igual que en el Gillespie, con la diferencia de que actúa reduciendo el número de susceptibles.
Cuando la fracción de infectados llegue al umbral Pth y se confine a una fracción Xs de la población, ∆i en el caso del SIR y ∆e en el caso del SEIR pasarán a ser:
∆i=Binomial(n =s(1−Xs), p= NiPβk)
∆e=Binomial(n =s(1−Xs), p= NiPβk)
3.2.4. Relación entre probabilidad y tasa
Para las tasas β, µ y γ nos interesa saber su relación con las probabilidades Pβ, Pµ y Pγ a fin de poder simular la misma situación en los tres modelos.
Como µdt es la probabilidad de recuperarse en un tiempo dt, la probabilidad de no recuperarse en un intervalo Dt dividido en pequeñosdt será:
p= (1−µdt)Dtdt (3.5)
En el límite dt →0:
p≈e−µDt (3.6)
Por tanto, la probabilidad de recuperación vendrá dada por:
Pµ= 1−e−µDt (3.7)
Donde se ha impuesto que Dt= 1 día.
De la misma manera se obtienenPβ y Pγ.
En este modelo, la transición entre estados viene descrita por las probabilidades, mien- tras que en el determinista y en el Gillespie por las tasas, lo que provoca que el R0 del modelo de las Binomiales sea parcialmente distinto (véase Sección 7.1 para más detalles):
R0 = Pβ
Pµ (3.8)
4. Resultados
El escenario que analizaremos será el de una población con N = 106 individuos, en el que las tasas de exposición y recuperación serán γ = µ = 0,16 y la de infección β = 2µ (R0 = 2); con condiciones iniciales: i0 = 10, e0 = 0 y r0 = 0. En del modelo de las Binomiales obtenemos la probabilidad de recuperación Pµ a partir de (3.7). La probabilidad de infección será Pβ = 2Pµ (R0 = 2) y en el SEIR la de exposición Pγ =Pµ (Se ha considerado un número de contactos k = 1).
Se ejecutarán los tres modelos con este conjunto de parámetros. En el caso de los estocásticos se hará la mediana de 103 de sus posibles soluciones para el de las Binomiales y de 102 para el Gillespie (102 en todos los diagramas de fase) y se representará junto a sus percentiles del 5 % y del 95 %. (Nota: Los datos del Gillespie son continuos en el tiempo y para poder hacer la mediana de n simulaciones es necesario disponer de datos discretos. Para ello, se ha utilizado el método de ”binning”, véase Sección 7.2 para más detalles.)
4.1. Caso general
Determinista Gillespie Binomiales
8.1 8.2 8.3
8.4 8.5 8.6
Figura 8: Fila de arriba SIR. Fila de abajo SEIR. En los modelos estocásticos el área en color representa el conjunto de soluciones entre los percentiles 5 y 95.
En esta primera situación el modelo SIR determinista (Fig. 8.1) predice que un 80 % de la población pasará la infección (véase la proporción de recuperados final) y que el pico
de infectados sucederá ent ≈75 días (conI(t)≈0,18). En comparación, el SEIR (Fig. 8.4) retrasa el punto de incidencia máxima a lost ≈160 días y disminuye su valor aI(t)≈0,1.
Además, aumenta la duración de la epidemia, pero la proporción de recuperados final sigue siendo en torno al 80 %.
Vemos cómo efectivamente la mediana de las soluciones del modelo estocástico del Gillespie (Fig. 8.2 y 8.5) da lugar a la determinista. Cabe destacar el aumento de las fluctuaciones durante el crecimiento y decrecimiento de la curva de infectados.
Por último, el modelo de las Binomiales retrasa varios días el pico de infectados tanto en el SIR (Fig. 8.3) como en el SEIR (Fig. 8.6), lo que implica una mayor duración del brote. Sin embargo, predice correctamente la proporción de recuperados final y el valor del máximo de incidencia.
4.2. Caso con confinamiento
Repetimos el escenario anterior aplicando un confinamiento del 50 % de la población (Xs = 0,5) cuando la proporción de infectados alcance el umbral Pth = 0,05. Éste valor de Xs no se ha escogido aleatoriamente: como hemos impuesto que β = 2µ (Pβ = 2Pµ en el de las Binomiales), la Ec. (2.26) nos dice que el régimen crítico aparecerá siXs = 0,5.
Determinista Gillespie Binomiales
9.1 9.2 9.3
9.4 9.5 9.6
Figura 9: Fila de arriba SIR. Fila de abajo SEIR. En los modelos estocásticos el área en color representa el conjunto de soluciones entre los percentiles 5 y 95.
Si comparamos la Fig. 9 con la Fig. 8 se ve claramente que el efecto principal del con- finamiento es una reducción drástica de la proporción final de recuperados. Sin embargo, no afecta igual al SIR que al SEIR. En el primero menos del 30 % de los individuos pasa en algún momento la enfermedad, mientras que en el segundo lo hacen más del 40 % (véanse Fig. 9.1 y 9.4).
Otra consecuencia apreciable es la aparición del régimen crítico: la duración del brote es más larga que en el caso sin confinamiento -aunque este aumento es mucho más pro- nunciado en el SIR que en el SEIR- y, tras aplicar el confinamiento, la curva de infectados se mantiene prácticamente constante durante un largo periodo de tiempo.
Los resultados del Gillespie (Fig. 9.2 y 9.5) coinciden con el determinista, mientras que el de las Binomiales difiere en lo mismo que en el apartado anterior: el retraso del
pico implica que el confinamiento se aplica más tarde, lo que provoca una mayor duración que en el caso del Gillespie.
4.3. Diagrama de fase
Recordemos que el régimen crítico se caracteriza por una evolución constante de la curva de infectados durante un largo periodo de tiempo. Por ello, el dato que utilizaremos para estudiar la aparición de dicho régimen será el de la duración del brote tmáx.
Con el fin de analizar cuándo aparece la criticidad, consideraremos una lista de posibles valores de Xs y otra dePth. Para cada pareja Xs−Pth se determinará tmáx.
En el caso de los modelos estocásticos se hará la mediana de 100 simulaciones.
Determinista Gillespie Binomiales
10.1 10.2 10.3
10.4 10.5 10.6
Figura 10: Fila de arriba SIR. Fila de abajo SEIR. Comparación entre los diagramas de fase de los tres modelos. Nótese que la escala de los colores en los tres casos es muy diferente.
En el modelo determinista (Fig. 10.1 y 10.4) vemos cómo efectivamente el régimen crítico aparece en Xs = 0,5 y también, que cuanto antes se aplique el confinamiento más larga es la duración del brote (el máximo aparece en el Pth más pequeño). Es importante tener en cuenta que en este caso no nos interesa la duración del brote en sí, ya que
la solución determinista está pensada para poblaciones infinitas, lo que implica que sea asintótica (en ninguno de los casos el brote llega a desaparecer). Por tanto, al imponer un mismo límite (arbitrario) a partir del cual el brote se detiene en todas las simulaciones, nos permite comparar cuánto se tarda en llegar a dicho límite en cada caso.
En los modelos estocásticos también el régimen crítico aparece enXs= 0,5 (Fig. 10.2 y 10.5 Gillespie , 10.2 y 10.5 Binomiales) con la diferencia de que el tmáx más grande no corresponde a los valores más pequeños de Pth, lo que se puede deber a que muchas de las simulaciones en las que el confinamiento se aplica nada más empezar, el brote muere casi instantáneamente, provocando una disminución considerable de la mediana. Por otro lado, SIR y SEIR presentan diferencias en cuanto a la forma del diagrama: el máximo no aparece en los mismos Pth. Finalmente, a pesar de que el modelo de las Binomiales predice una duración más larga que el Gillespie, coincide en el punto crítico lo que lo vuelve adecuado para el estudio de esta situación.
5. Limitaciones de los modelos
Los modelos pueden ser realmente útiles a la hora hacer predicciones, pero hay que ser consciente de sus limitaciones a fin de saber en qué situaciones se pueden utilizar.
La primera pregunta que nos deberíamos hacer en nuestro caso es qué modelo entre el SIR y el SEIR es el más adecuado. Se elegirá uno u otro en función de si el periodo de incubación 1/γ es lo suficientemente corto como para poder despreciarlo (SIR) o si, por el contrario, su larga duración afecta a la dinámica de la epidemia (SEIR).
El siguiente punto a considerar es si se resuelve el caso mediante el modelo determi- nista o uno de los estocásticos. Para ello, es necesario plantearse si la dinámica de cada individuo puede tener una importante influencia en la evolución del brote (casi siempre es así) o si se puede despreciar, lo que depende del tamaño de la población N. Para poblaciones pequeñas, la aleatoriedad gana importancia y, por tanto, el modelo determi- nista (poblaciones infinitas) deja de ser válido. En la Fig. 11 se muestra cómo aumenta la incertidumbre en el caso de una población con pocos individuos.
11.1:N = 103 11.2:N = 106
Figura 11: Comparación entre dos simulaciones del SIR hechas con el modelo de las Binomiales con el fin de ver el aumento de las fluctuaciones en una población con una N pequeña. Parámetros utilizados: µ= 0,16; β = 2µ.
Por último, si se decide utilizar el modelo estocástico a fin de tener en cuenta las posibles fluctuaciones en la solución, se ha de optar entre el Gillespie o el de las Binomiales.
El segundo es una aproximación del primero y se ha mostrado en la Sección 4 las posibles diferencias con respecto a la solución teórica a las que puede dar lugar. Sin embargo, la enorme carga computacional que supone utilizar el Gillespie con una población grande (los diagramas de fase requerían del orden de 8 horas de ejecución, mientras que los del modelo de las Binomiales menos de 5 min) vuelve realmente complicado llevar a cabo determinados análisis. Por tanto, puede valer la pena sacrificar cierta precisión a fin de ganar eficacia computacional. No obstante, el modelo de las Binomiales deja de ser válido si las probabilidades Pβ, Pµ o Pγ son demasiado grandes, ya que podría darse el caso de que un individuo se infectara y se recuperara dentro del intervalo de tiempo discreto (el modelo sólo considera un evento por individuo en cada intervalo de tiempo).
6. Conclusión
En este trabajo, motivado por la pandemia del COVID-19, se ha buscado explicar de manera detallada los modelos epidemiológicos SIR y SEIR sin demografía, así como también su resolución tanto determinista como estocástica. Se ha planteado la situación de un confinamiento parcial modificando los modelos originales y se han visto los efectos que tiene sobre la evolución de la epidemia. El lector podrá reproducir cualquier apartado si así lo desea.
Se han comparado los tres casos (determinista, Gillespie y Binomiales) destacando las posibles diferencias que podían aparecer entre cada uno de ellos. Hemos visto que la mediana de las soluciones del Gillespie da lugar a la solución determinista, mientras que su aproximación (el de las Binomiales) retrasa el pico de infectados y aumenta la duración del brote, aunque predice correctamente el valor del máximo de incidencia y la cantidad de recuperados final.
Se ha analizado la aparición del régimen crítico en los tres modelos dependiendo de los parámetros del confinamiento Xs y Pth que eligiéramos en un caso supercrítico con R0 = 2 para una población de un millón de habitantes, con el resultado de que tanto en el SIR como en el SEIR la criticidad aparece en Xs = 0,5 para los tres casos, comprobando así la validez de la Ec. (2.26). Por otro lado, la gran similitud entre los diagramas de fase de los dos modelos estocásticos nos indica que el de las Binomiales puede ser una buena alternativa al Gillespie (teniendo en cuenta sus limitaciones), debido a que el trabajo computacional que requiere éste último puede dificultar estudios en los que se requieran un gran número de simulaciones con un tamaño de población grande.
Destacar que es importante el estudio de este régimen, ya que al caracterizarse por una evolución constante de la curva de infectados durante un largo periodo, significa que el relajamiento de las medidas de contención antes de tiempo darían lugar a otro brote.
En la situación ideal planteada, se podría pensar que un confinamiento parcial del 50 % de la población sería una buena opción para acabar con la epidemia: no se confina a todos los individuos consiguiendo no paralizar del todo la economía y, al reducir drásticamente el número de susceptibles, se logra erradicar al patógeno. Sin embargo, como hemos visto en la Fig. 10, eso sería un error, pues daría lugar al régimen crítico alargando la duración de la epidemia y provocando un rebrote en cuanto las autoridades decidieran desconfinar a la población.
El presente trabajo ha tratado una situación idealizada a la que se le puede añadir complejidad con el fin de ser más realista. Por tanto, puede servir como punto de par- tida para futuros estudios de la criticidad en los que, por ejemplo, se tenga en cuenta la demografía (nacimientos y muertes) o se considere un sistema de varias poblaciones
interconectadas. Por otro lado, también se puede introducir el efecto de la vacunación o utilizar modelos con otros compartimentos (SIS, SEIS, MSEIR...).
7. Apéndices
7.1. Cálculo de R
0mediante la Next Generation Matrix
El número reproductivo básico de un modelo epidemiológico se obtiene utilizando la llamada Next Generation Matrix [9]. En este apartado explicaremos su significado aplicándola en el modelo SEIR (fuentes que se han consultado [10],[11] y [12]).
Sean las ecuaciones de expuestos (2.28) e infectados (2.29) linealizadas en torno al punto de absorción (S ≈1,I << 1,E <<1):
dE
dt ≈βI−γE (7.1)
dI
dt =γE−µI (7.2)
Si definimos el vector X~ = (E, I)0, el sistema de ecuaciones anterior se puede escribir de la siguiente manera:
d ~X
dt = −γ β
γ −µ
!
X~ =
"
0 β 0 0
!
− γ 0
−γ µ
!#
X~ = (T −S)X~ (7.3) DondeT es la parte de transmisión y S la de transición.
Si resolvemos la Ec. (7.3), nos queda:
X~ =X~0exphT S−1−1tSi (7.4) T S−1es laNext Generation Matrix yR0se define como su autovalor dominante, ya que éste será el que determinará la dinámica de las curvas de los distintos compartimentos de la clase ”portadores del patógeno” (expuestos e infectados en nuestro caso). Así, el autovalor más grande de la matriz T S−1 será:
R0 = β
µ (7.5)
Quedando demostrado que el SIR y el SEIR tienen el mismo número reproductivo básico.
Para el SIR el cálculo deR0 es trivial, véase [11].
En el caso del modelo de las Binomiales, estamos discretizando el sistema de ecuaciones en intervalos en los que las transiciones ya no vienen determinadas por las tasas, sino por las probabilidades:
dE
dt ≈ E(t+Dt)−E(t)
Dt ≈PβI−PγE (7.6)
dI
dt ≈ I(t+Dt)−I(t)
Dt ≈PγE−PβI (7.7)
Por ello, elR0 en este modelo viene dado por:
R0 = Pβ
Pµ (7.8)
7.2. Método de ”binning” aplicado a los datos del Gillespie
El método de ”binning” permite agrupar datos continuos en ”bins” donde cada ”bin”
corresponde a un tiempo discreto.
A modo de ejemplo se mostrará cómo se ha aplicado en los datos del SIR:
1. Primeramente, modificaremos el algoritmo de la sección 3.1.1 añadiendo el cálculo de sbin, ibin y rbin ( inicialmente sbin =s0, ibin = i0 y rbin = r0). En cada iteración se han de redefinir e ir guardando los nuevos valores:
sbin →sbin+s ibin →ibin+i rbin →rbin+r
2. Al acabar la simulación, guardamos la duración máxima del brote tmáx, la lista de tiempo continuo y las de los valores de sbin, ibin y rbin .
3. Definimos los ”bins” como una lista de tiempos discretos: [0,1,2, ..., tmáx].
4. El dato del ”bin” n para una variable concreta ( sbin, ibin o rbin) será la suma de todos los valores de la variable en cuestión, cuyos tiempos continuos estén entre el
”bin”ny eln+1, dividida entre el número de estos valores. Por ejemplo, imaginemos que tenemos unos valores de sbin = [s1, s2..., sf] cuyos tiempos continuos son t = [0,5; 0,6;...; 0,9]. Como estos están entre 0 y 1, corresponden al ”bin” (0). Entonces los susceptibles en t= 0 serán:
s(0) =
Pf j=1sj
f (7.9)
De esta manera se consiguen datos de s, i y r discretizados con los que es posible calcular la mediana entre distintas simulaciones del Gillespie.
Referencias
[1] M. J. Keeling y P. Rohani, Modeling Infectious Diseases in Humans and Animals (Princeton University Press, 2008).
[2] Número reproductivo básico, https://www.ncbi.nlm.nih.gov/pmc/articles/
PMC6291769/.
[3] COVID-19 en España, https://cnecovid.isciii.es/covid19/#ccaa.
[4] F. Radicchi y G. Bianconi, «Epidemic plateau in critical susceptible-infected-removed dynamics with nontrivial initial conditions», Physical Review E102 (2020).
[5] W. O. Kermack y A. G. McKendrick, «A Contribution to the Mathematical Theory of Epidemics», Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character115, 700-721 (1927).
[6] P. L. DeVries, A First Course in Computational Physics, 1st (John Wiley & Sons, Inc., USA, 1993).
[7] D. T. Gillespie, «Exact Stochastic Simulation of Coupled Chemical Reactions»,The Journal of Physical Chemistry81, 2340-2361 (1977).
[8] D. Balcan, V. Colizza, B. Gonçalves, H. Hu, J. J. Ramasco y A. Vespignani, «Multis- cale mobility networks and the spatial spreading of infectious diseases»,Proceedings of the National Academy of Sciences106, 21484-21489 (2009).
[9] O. Diekmann y J. Heesterbeek, Mathematical epidemiology of infectious diseases:
model building, analysis and interpretation, English, Wiley series in mathematical and computational biology (John Wiley y Sons, United States, 2000).
[10] Notes on R0, http : / / web . stanford . edu / ~jhj1 / teachingdocs / Jones - on - R0.pdf.
[11] Compute R0 using Next Generation Operator, https://qrlssp.asu.edu/sites/
default/files/brn_mtbi_2016.pdf.
[12] The construction of next-generation matrices for compartmental epidemic models, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2871801/#RSIF20090386C2.
