Modelos matemáticos en epidemiología

Traballo Fin de Grao Modelos matemáticos en epidemiología Alejandra Comesaña García 2020/2021 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA GRAO DE MATEMÁTICAS Traballo Fin de Grao Modelos matemáticos en epidemiología Alejandra Comesaña García Julio, 2021 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA Trabajo propuesto Área de Coñecemento: Análise Matemática Título: Modelos matemáticos en epidemioloxía Breve descrición do contido: A modelización matemática é unha ferramenta fundamental no ámbi- to da epidemoloxía, xa que permite analizar cuestións esenciais dunha epidemia como a incidencia, o contaxio ou a propagación. Isto asegura un mellor coñecemento que facilita a toma de decisións e establecemento de medidas oportunas para minimizar o seu impac- to. Por exemplo, resulta de gran utilidade para combatir e controlar o densenvolvemento das epidemias, permitindo predecir as posibles consecuencias de introducir medidas específicas. O obxectivo deste traballo é analizar distintos modelos epidemiolóxicos baseados principalmente en ecuacións diferenciais. Recomendacións Outras observacións iii Índice general Resumen viii Introducción xi 1. Ecuaciones diferenciales ordinarias 1 1.1. Sistemasautónomos ............................... 1 1.2. Estabilidad .................................... 3 2. Ecuaciones diferenciales con retardo 7 2.1. Conceptosbásicos................................. 8 2.2. Existenciayunicidad............................... 9 2.3. Estudiocualitativo................................ 13 2.3.1. Sistemas lineales autónomos . . . . . . . . . . . . . . . . . . . . . . . 14 2.3.2. La ecuación característica . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3.3. La ecuación escalar . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.3.4. Principio de estabilidad linealizada . . . . . . . . . . . . . . . . . . . 18 2.3.5. Estabilidad absoluta . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3. Modelos epidemiológicos 23 3.1. Introducción.................................... 23 3.2. Del modelo SIR al SEIR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.3. Modelos epidemiológicos con retardo . . . . . . . . . . . . . . . . . . . . . . 37 4. La pandemia del coronavirus en Italia 43 4.1. Simulaciónnumérica ............................... 44 4.1.1. Modelos sin medidas de contención . . . . . . . . . . . . . . . . . . . 45 4.1.2. Medidas de contención . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.2. Conclusiones ................................... 52 v vi ÍNDICE GENERAL A. Resultados de análisis complejo 55 B. Código MATLAB utilizado 57 Bibliografía 63 xiv INTRODUCCIÓN Capítulo 1 Ecuaciones diferenciales ordinarias En este capítulo se repasarán algunos conceptos y resultados acerca de las ecuaciones diferenciales ordinarias que utilizaremos a lo largo del trabajo. Todos ellos han sido estudiados en las asignaturas de Introducción a las Ecuaciones Diferenciales Ordinarias y Ecuaciones Diferenciales Ordinarias. Nos centraremos en los sistemas autónomos y, más concretamente, en el estudio cualitativo de los mismos, puesto que serán resultados de interés a la hora de estudiar los modelos que propondremos posteriormente. Una posible referencia para profundizar en los temas a tratar es [4]. 1.1. Sistemas autónomos La imposibilidad, en muchos casos, de obtener expresiones explícitas de las soluciones de los sistemas de ecuaciones no lineales hace que el estudio cualitativo de los mismos tenga especial interés. Es decir, el objetivo es el de obtener la máxima cantidad de información posible de las soluciones y, más concretamente, de su comportamiento y sus propiedades en base, simplemente, al propio sistema. Veremos entonces diferentes resultados que nos faciliten este estudio. Consideremos la ecuación x0=F(x),(1.1) donde x=x(t)yF:A⊂Rn−→ Rn, con Aun abierto de Rn. Definición 1.1. La ecuación (1.1) recibe el nombre de ecuación diferencial autónoma, puesto que la función Fno depende explícitamente de t. Observación 1.2.Supondremos que la función Fes localmente lipschitziana en Rn. De esta forma, tendremos garantizada, dada una condición inicial (t0, x0)∈R×A, existencia y unicidad de solución de (1.1) pasando por (t0, x0). 1 2CAPÍTULO 1. ECUACIONES DIFERENCIALES ORDINARIAS Observación 1.3.Además, si x(t)es una solución de (1.1) en Imax = (a, b), para todo c∈Rse tiene que x(t−c)es solución de la misma ecuación en (a+c, b +c). Es decir, las soluciones de una ecuación autónoma son, en cierto modo, invariantes por traslaciones temporales. Como consecuencia, para conocer las soluciones de la ecuación (1.1) podremos considerar t0= 0. Llamaremos Ix= (α, β)al intervalo maximal de definición de la solución que pasa por (0, x). Sea ϕx:Ix−→ Rnla solución maximal de (1.1) que pasa por (0, x), con x∈A. Definición 1.4. El conjunto γx={ϕx(t) : t∈Ix} recibe el nombre de órbita de x. De la definición anterior obtenemos que, para cada instante inicial considerado, la órbita es la proyección de la solución en Rn. Además, debido a que ϕxes solución de la ecuación, se tiene que el campo Fes tangente a la órbita en cada punto. Tenemos entonces definida una orientación de la curva acorde con el recorrido temporal. Observación 1.5.Para cada punto x0∈Aexiste una única órbita. Por lo tanto, dos órbitas de (1.1) o son la misma o no se intersecan. Teniendo en cuenta esta última observación, el conjunto de órbitas será disjunto. Definición 1.6. Llamaremos diagrama de fases oretrato de fases al conjunto definido como la unión disjunta de todas las órbitas (orientadas en el sentido del tiempo) que pasan por todos los puntos x∈A. Obviamente, debido a que la curva solución se definiría como (t, x(t)) ∈R×Rn, la ventaja de estudiar el conjunto de las órbitas es que hemos pasado de dimensión n+ 1 an. Por ello, el estudio cualitativo que se hace en los sistemas autónomos tiene mucho que ver con el estudio del diagrama de fases, y el objetivo será dar una descripción lo más precisa posible del mismo. Veamos ahora cuáles son los puntos del dominio de definición de Fque resultarán interesantes para su estudio en el diagrama de fases. Definición 1.7. Dado el sistema autónomo (1.1) y dado x0∈A, diremos que x0es una singularidad opunto de equilibrio si F(x0)=0. x0es regular si F(x0)6= 0. Observación 1.8.Si x0es una singularidad, entonces γx0={x0}. 1.2. ESTABILIDAD 3 1.2. Estabilidad Una vez definidos los conceptos básicos de los sistemas autónomos, pasemos a una parte fundamental del estudio cualitativo: la estabilidad de los puntos singulares. De aquí en adelante consideraremos x0=F(x),(1.2) con F:A⊂Rn−→ Rnuna función de clase C1yAun abierto. Denotaremos por ϕx:Ix⊂ R−→ Rnla solución maximal del sistema con condición inicial ϕx(0) = x. Suponemos que las definiciones de estabilidad, estabilidad asintótica e inestabilidad son conocidas, por lo que nos centraremos directamente en los criterios que faciliten el estudio. Una posible referencia para profundizar en los resultados de esta sección es [5, Lecture 23]. Sea, de aquí en adelante, x0una singularidad de F. Definición 1.9. Diremos que x0es un atractor si existe un entorno Wde x0tal que si x∈W, se tiene que ϕx(t)está definida para todo t≥0yl´ım t→∞ ||ϕx(t)−x0|| = 0. Definición 1.10. Se dice que x0es una singularidad hiperbólica si todos los autovalores de DF (x0)tienen parte real no nula, donde DF (x0)representa la matriz jacobiana de F en x0. Método de la primera aproximación. La ventaja de este método es que las hipótesis que se piden son muy fáciles de verificar. Sin embargo, dejaremos atrás el caso de que una singularidad sea estable pero no asintóticamente estable. Teorema 1.11. Sea x0una singularidad de F. Entonces si todos los autovalores de DF (x0)tienen parte real negativa, x0es asintóticamente estable. si algún autovalor de DF (x0)tiene parte real positiva, x0es inestable. Bajo la hipótesis de que la singularidad sea hiperbólica, se obtienen condiciones necesarias y suficientes. Corolario 1.12. Si x0es una singularidad hiperbólica, se verifican x0es asintóticamente estable si, y solamente si, todos los autovalores de DF(x0) tienen parte real negativa. 4CAPÍTULO 1. ECUACIONES DIFERENCIALES ORDINARIAS x0es inestable si, y solamente si, algún autovalor de DF (x0)tiene parte real positiva. La formalización del paso del Teorema al Corolario se basa en el Teorema de Hartman- Grobman, que se puede encontrar en [6]. Método directo de Liapunov. Como hemos dicho, con el método anterior no podemos saber si una singularidad es estable sin ser asintóticamente estable. El método de Liapunov que veremos a continuación no tendrá esa limitación. A mayores, podremos dar aproximaciones de las regiones de atracción. Sea F∈C1(A, Rn)yx0una singularidad de F. Definición 1.13. Sea V:D⊂A−→ R, con Dun abierto y x0∈D. Diremos que Ves una función de Liapunov asociada a F si se verifican las siguientes condiciones: 1. V∈C1(D),V(x0) = 0. 2. V(x)>0∀x∈D,x6=x0. 3. V∗(x) := h∇V(x), F(x)i ≤ 0∀x∈D. Además, se dice que una función de Liapunov es estricta si V∗(x)<0para todo x∈D, x6=x0. Observación 1.14.Se tiene que V∗(x) := h∇V(x), F(x)i=Pn i=1 ∂V ∂xi(x)Fi(x). A partir de la función de Liapunov, se pueden obtener los siguientes resultados. Teorema 1.15. Sea x0una singularidad para F, con Fen las hipótesis mencionadas. Se verifican las afirmaciones que siguen. Si existe una función de Liapunov relativa a x0, entonces x0es estable. Si existe una función de Liapunov estricta relativa a x0, entonces x0es asintóticamente estable. Teorema 1.16. Si existe una función V:D⊂A−→ R, con x0∈Dtal que 1. V∈C1(D),V(x0)=0, 2. V(x)>0∀x∈D,x6=x0, y 3. V∗(x)>0∀x∈D{x0}, entonces x0es inestable. 1.2. ESTABILIDAD 5 Como habíamos anticipado, una ventaja del método de Liapunov con respecto al de la primera aproximación es la posibilidad de dar estimaciones de las regiones de atracción. Formalicemos este concepto. Definición 1.17. Sea x0una singularidad estable para el sistema x0=F(x). Definimos la región de atracción de x0como R={y∈Atal que ϕy(t)está definida para todo t≥0yl´ım t→∞ ϕy(t) = x0}.(1.3) Definición 1.18. Dado el sistema x0=F(x), sea M⊂A. Decimos que Mes positivamente invariante para el sistema si ∀x∈M,ϕx(t)está definida para todo t≥0yγ+ x⊂M. Decimos que Mes negativamente invariante para el sistema si ∀x∈M,ϕx(t)está definida para todo t≤0yγ− x⊂M. Decimos que Mes invariante si ∀x∈M,ϕx(t)está definida para todo t∈Ry γx⊂M. Teorema 1.19. La región de atracción de una singularidad asintóticamente estable x0es positivamente invariante. Para finalizar este capítulo, enunciemos el siguiente resultado, que nos facilitará la tarea de encontrar una región de atracción para una singularidad x0. Teorema 1.20. Sea V:D−→ Rrelativa a x0para el sistema x0=F(x). Supongamos que el conjunto S={x∈D:V∗(x) = 0} no contiene ninguna solución del sistema excepto x0. Entonces, se verifican los enunciados que siguen. x0es asintóticamente estable Si K⊂Des un conjunto compacto y positivamente invariante, entonces Kpertenece a la región de atracción de x0. Si Ck, componente conexa del conjunto Sk={x∈D:V(x)≤k}, es compacta, entonces forma parte de la región de atracción de x0. 6CAPÍTULO 1. ECUACIONES DIFERENCIALES ORDINARIAS Capítulo 2 Ecuaciones diferenciales con retardo La teoría de las ecuaciones diferenciales con retardo (EDR) requiere un conocimiento más amplio sobre otras ramas del análisis que las ecuaciones diferenciales ordinarias. Sin ir más lejos, en el ejemplo que se propone a continuación como introducción, son necesarias diversas nociones de análisis complejo. Al final del trabajo se encuentra un apéndice con los resultados y definiciones necesarios de esta rama que aparecerán a lo largo del capítulo. Ejemplo 2.1. Para el caso de las ecuaciones diferenciales lineales con coeficientes constantes, una idea a la que se recurre es la de proponer soluciones del tipo x(t) = eλt y obtener los valores de λque anulan el polinomio característico asociado a la EDO. Por ejemplo, si consideramos la ecuación x0(t) = x(t), sabemos que su polinomio característico es P(λ) = λ−1. La única raíz del polinomio es λ= 1 y, por ello, una solución a la EDO es x(t) = et. Realmente, cualquier función de la forma x(t) = Cet, con C∈Res solución de la EDO, y hemos obtenido así todas las posibles soluciones reales de la ecuación de partida. Veamos ahora como cambia la situación cuando introducimos un retardo τ > 0y pasamos a trabajar en el plano complejo. La EDR resulta x0(t) = x(t−τ)y, reemplazando x(t) = eλt, obtenemos λx(t) = x(t−τ) = eλ(t−τ)=e−λτ x(t). Podemos observar que la ecuación obtenida no es un polinomio, si no una función P:C−→ C, definida como P(λ) = λ−e−λτ , de la cual debemos obtener sus raíces. Sabemos que λ= 0 no es una raíz, por lo que podemos tomar z= 1/λ y hacer que la ecuación a resolver sea ahora ze−τ/z = 1. Es aquí donde entra en juego el análisis complejo: la función f(z) = ze−τ/z tiene una singularidad esencial en z= 0 y, además, no se anula en C\{0}por lo que, en virtud del Teorema grande de Picard se deduce que f−1(1) tiene 7 8CAPÍTULO 2. ECUACIONES DIFERENCIALES CON RETARDO un número infinito de elementos. ¿Qué significa esto en nuestro problema? Que el número de soluciones de la ecuación ze−τ/z = 1 es infinito. Es decir, hemos obtenido un conjunto {λk}k∈N, con |λk| → ∞ de raíces características. Por tanto, la ecuación diferencial con retardo de partida tiene infinitas soluciones de la forma x(t) = eλkty, a mayores, cualquier combinación lineal en Cde estas funciones también lo será. Hemos obtenido entonces un espacio de soluciones de dimensión infinita. El ejemplo que acabamos de ver pone de manifiesto el hecho de que introducir un retardo supone un cambio sustancial en la resolución de las ecuaciones. Comencemos con la teoría de las ecuaciones diferenciales con retardo. Para las demostraciones de los resultados pueden consultarse [7, 8]. 2.1. Conceptos básicos Las ecuaciones diferenciales con retardo nos permitirán abordar el estudio de sistemas de forma más realista, puesto que en estas ecuaciones se involucran también los estados anteriores de estos sistemas. Si bien, como hemos visto, la resolución de las EDRs requiere técnicas más avanzadas que en el caso de las EDOs, el avance de la computación y la posibilidad de simular las soluciones de dichos sistemas han hecho que el estudio de los mismos resulte un problema muy interesante. Definición 2.2. Sea x(t)una función dependiente de la variable temporal t. Dada x(t−τ), decimos que τ∈Res un retardo discreto. Consecuentemente, dados los retardos τi,i= 1, . . . , m, definimos las ecuaciones diferenciales con retardo discreto como aquellas de la forma x0(t) = f(t, x(t), x(t−τ1), x(t−τ2), . . . , x(t−τm)),(2.1) con t∈R,x∈Rn,f:D⊂R×Rn(m+1) −→ Rn. Nos restringiremos al caso en que τsea un parámetro positivo fijo, pero podríamos también tomarlo negativo, puesto que basta pensar en un cambio de variable y volver al caso inicial. Definición 2.3. En adición a la ecuación (2.1), para determinar un problema de valor inicial se requiere que la solución x(t)satisfaga la siguiente condición en un intervalo temporal: x(t) = ϕ(t), t0−max{τ1, τ2, . . . , τm} ≤ t≤t0,(2.2) donde t0∈Ryϕ∈Rn. La función ϕrecibe el nombre de función inicial. 2.2. EXISTENCIA Y UNICIDAD 9 Observación 2.4.Si τi= 0,i= 1, . . . , m, entonces la ecuación (2.1) se reduce a una EDO y, además, la condición inicial (2.2) se convierte en x(t0) = ϕ(t0), es decir, hemos generalizado la noción de Problema de Cauchy que ya conocíamos. Observación 2.5.Además de los retardos discretos, existen también retardos del tipo Rt t−τk(t−s)x(s)ds =Rτ 0k(z)x(t−z)dz, con 0≤τ≤ ∞, que reciben el nombre de retardos distribuidos, puesto que reflejan la media ponderada de los retardos x(t−z), como se observa en el segundo término de la igualdad. Las ecuaciones que incluyen un retardo distribuido tienden a ser más realistas, pero también requieren del uso de herramientas más avanzadas y que no trataremos en este trabajo. Es por ello que asumiremos siempre que cuando se hable de EDR, se está considerando el caso de retardo discreto. 2.2. Existencia y unicidad El estudio de la existencia y unicidad de soluciones en las ecuaciones diferenciales con retardo discreto desciende del método de los pasos para la resolución de las mismas. Este método se lleva a cabo de forma iterativa. Escogido un intervalo temporal específico, se realiza un cambio de variable en tde forma que, en dicho intervalo, la EDR se convierte en una EDO que sabemos resolver. Una vez obtenida su solución, se considera el siguiente intervalo temporal, se realiza el nuevo cambio de variable y se toma como condición inicial la solución obtenida previamente. Pongamos un ejemplo. Ejemplo 2.6. Supongamos la EDR unidimensional x0(t) = αx(t−τ) x(t) = a, −τ≤t≤0,(2.3) con τ > 0el retardo discreto y ayαparámetros reales positivos. Trabajemos, en primer lugar, con ten el intervalo [0, τ]. Tomando r=t−τ, se tiene que −τ≤r≤0. Por lo tanto, x(t−τ) = x(r) = a. Sustituyendo en (2.3), x0(t) = αx(t−τ) = αx(r) = αa, 0≤t≤τ, de donde, resolviendo la EDO por variables separadas, x(t) = x(0) + Zt 0 aαds =a+aαt. (2.4) Determinemos ahora la solución en el intervalo [τ, 2τ]. Con el mismo cambio de variable, 16 CAPÍTULO 2. ECUACIONES DIFERENCIALES CON RETARDO o, equivalentemente, λv =L(expλv). Escribiendo ahora v=Pn j=1 vjej, con {ej}jla base canónica de Cn, obtenemos que L(expλv) = Pn j=1 vjL(expλej). Definimos entonces Lcomo la matriz n×n Lλ= (L(expλe1)|L(expλe2)|. . . |L(expλen)) = (Li(expλej)), donde Li(ϕ)es la componente i-ésima de L(ϕ). Entonces, se tiene que L(expλv) = Lλv, y podemos ver que x(t) = eλtves una solución no nula de (2.14) si λes solución de la ecuación característica det(λI −Lλ)=0.(2.17) Habitualmente nos referiremos a λ∈Ccomo una raíz característica. Entra en juego ahora el análisis complejo. Esto no hace más que reforzar una idea que ya hemos mencionado: las ecuaciones diferenciales con retardo necesitan de un conocimiento de otras ramas de las matemáticas mucho más profundo que el de las ecuaciones diferenciales ordinarias. Si echamos la vista atrás, esto ya fue ilustrado con el ejemplo con el que empezamos el capítulo. Lema 2.14. h(λ) = det(λI −Lλ)es una función entera. Lema 2.15. Dado σ∈R, existe un número finito de raíces características tal que Re(λ)> σ. Si hay un número infinito de raíces características distintas {λn}n, entonces Re(λn)→ −∞, n → ∞. Este último resultado implica la existencia de un σ∈Ry de un conjunto finito de raíces características cuya parte real será igual a σ, mientras que todas las demás tendrán parte real estrictamente menor que σ. Otra apreciación es que las raíces características complejas vendrán en pares conjugados. Formalmente: Proposición 2.16. Sea Luna función tal que L(C([−τ, 0],Rn)) ⊂Rn. Entonces, λes una raíz característica si, y solamente si, λlo es. Demos ahora el resultado central de este apartado, la estabilidad de la solución x= 0 para (2.14). Teorema 2.17. Supongamos que Re(λ)< µ para toda raíz característica λ. Entonces, existe K > 0tal que |x(t, ϕ)| ≤ Keµt||ϕ||, t ≥0, ϕ ∈C, (2.18) 2.3. ESTUDIO CUALITATIVO 17 donde x(t, ϕ)es la solución de (2.14) con condición inicial x0=ϕ. En particular, x= 0 es asintóticamente estable para (2.14) si Re(λ)<0para toda raíz característica, y es inestable si existe algún λtal que Re(λ)>0. Por lo tanto, hemos obtenido condiciones para la estabilidad asintótica y para la inestabilidad. Estas tienen que recordarnos a las que ya enunciamos en la sección anterior. Es natural entonces hacerse la siguiente pregunta: ¿hasta qué punto es relevante el retardo en el estudio de las raíces características? Demos una respuesta considerando el sistema con retardo z0(t) = Az(t) + Bz(t−τ) y el sistema sin retardo asociado z0(t)=(A+B)z(t). Queremos estudiar el tipo de relación que hay entre las raíces características de h(λ, τ) = det[λI −A−e−λτ B] = 0 y los autovalores de A+B h(λ, 0) = det[λI −A−B] = 0. La relación entre las raíces de las dos ecuaciones y, por tanto, la respuesta a la pregunta, se encuentra en el siguiente resultado. Teorema 2.18. Sean z1, z2,...zklos distintos autovalores de A+B,δ > 0ys∈Rtal que s < m´ıni=1,...,k Re(zi). Entonces, existe τ0>0tal que si 0< τ < τ0yh(z, τ) = 0 para algún z, entonces o bien Re(z)< s o|z−zi|< δ para algún i. Es decir, si el retardo es suficientemente pequeño, las raíces características del sistema con retardo, o bien están muy cerca de los autovalores de A+B, o bien tienen una parte real más pequeña que cualquiera de los autovalores de A+B. En otras palabras, podemos considerar que los retardos pequeños no afectan considerablemente en el sentido que explicaremos a continuación. Esto es, en caso de tener estabilidad asintótica para τ= 0, entonces podemos seguir teniéndola para retardos pequeños, simplemente tomando δlo suficientemente pequeño como para que la bola de radio δcentrada en cada uno de los autovalores de A+Bsiga estando en el semiplano negativo y escogiendo snegativo. Por otro lado, se puede ver que si tenemos inestabilidad para τ= 0 debida a una raíz simple positiva o una pareja de raíces complejas conjugadas con parte real positiva, la inestabilidad 18 CAPÍTULO 2. ECUACIONES DIFERENCIALES CON RETARDO permanece siempre que consideremos retardos pequeños. 2.3.3. La ecuación escalar Particularicemos ahora el estudio de las raíces características y, consecuentemente, de la estabilidad. Esto lo haremos porque, al igual que en los sistemas de ecuaciones ordinarias, un método clave para estudiar la estabilidad de las singularidades en un sistema no lineal es el de considerar el linealizado del sistema alrededor de la singularidad. Consideremos entonces la ecuación escalar x0(t) = Ax(t) + Bx(t−τ),(2.19) con A,B∈R. El caso matricial se razona de modo análogo. Tenemos que la ecuación característica es λ=A+Be−λτ .(2.20) Si z=λτ,α=Aτ,β=Bτ, obtenemos que la ecuación anterior se puede escribir como z=α+βe−z.(2.21) Nuestro objetivo ahora sería estudiar las raíces de esta ecuación. El desarrollo del mismo está en [8, pág.50]. Volviendo al problema original, la estabilidad del estado estacionario x= 0 para la ecuación escalar (2.19) depende de las raíces de la ecuación característica (2.20). Asumamos que A+B6= 0, puesto que, en caso contrario, λ= 0 sería una raíz. Teorema 2.19. Las siguientes afirmaciones se verifican para el sistema (2.19). 1. Si A+B > 0, entonces x= 0 es inestable. 2. Si A+B < 0yB≥A, entonces x= 0 es asintóticamente estable. 3. Si A+B < 0yB < A, entonces existe τ∗>0tal que x= 0 es asintóticamente estable para 0< τ < τ∗e inestable para τ > τ∗. En el caso (3), existe un par de raíces imaginarias puras en τ=τ∗= (B2−A2)−1/2cos−1(−A/B). 2.3.4. Principio de estabilidad linealizada Pasemos ahora a la determinación de la estabilidad de las singularidades de sistemas no lineales mediante el método del linealizado. 2.3. ESTUDIO CUALITATIVO 19 Consideremos el siguiente sistema no lineal x0(t) = f(xt).(2.22) Entonces x(t) = x0∈Rn,t∈Res una solución estacionaria para (2.22) si, y solamente si f(ˆx0) = 0, donde ˆx0∈Ces la función constante igual a x0. Si x(t)es una solución de (2.22) y x(t) = x0+y(t), entonces y(t)satisface y0(t) = f(ˆx0+yt).(2.23) Queremos entender el comportamiento de las soluciones de (2.22) que empiezan cerca de ˆx0y, para ello, es suficiente con entender el comportamiento de las soluciones de (2.23) para soluciones que empiezan cerca de y= 0. Asumiremos que f(ˆx0+ϕ) = L(ϕ) + g(ϕ), ϕ ∈C, (2.24) donde, como hemos considerado hasta ahora, L:C−→ Rnes una función lineal acotada yg:C−→ Rnverifica l´ım ϕ→0 |g(ϕ)| ||ϕ|| = 0. Esta última condición significa que, para todo ε > 0, existe δ > 0tal que ||ϕ|| ≤ δ=⇒ |g(ϕ)| ≤ ε||ϕ||. El sistema lineal z0(t) = L(zt)(2.25) recibe el nombre de ecuación linealizada alrededor del equilibrio ˆx0. Debemos verla en el espacio C=C([−τ, 0],Cn). Demos ahora el resultado central de esta sección, que relaciona los autovalores del linealizado con la estabilidad del punto singular. Teorema 2.20. Sea ∆(λ) = 0 la ecuación característica correspondiente a (2.25) y supongamos que −σ:= m´ax ∆(λ)=0 Re(λ)<0. 20 CAPÍTULO 2. ECUACIONES DIFERENCIALES CON RETARDO Entonces, ˆx0es un estado estacionario localmente asintóticamente estable de (2.22). Además, existe b > 0verificando ||ϕ−ˆx0|| < b =⇒ ||xt(ϕ)−ˆx0|| ≤ K||ϕ−ˆx0||e−σt/2, t ≥0.(2.26) Si Re(λ)>0para alguna raíz característica, entonces ˆx0es inestable. En concreto, podemos suponer el caso particular de la ecuación (2.22) x0(t) = F(x(t), x(t−τ)),(2.27) donde asumimos que F:D×D−→ Rnes de clase C1yD⊂Rnes un abierto. Si F(x0, x0)=0para algún x0∈D, entonces x(t) = x0, con t∈R, es un equilibrio. Entonces, f(ϕ) = F(ϕ(0), ϕ(−τ)), por lo que (2.24) se convierte en f(ˆx0+ϕ) = Aϕ(0) + Bϕ(−τ) + G(ϕ(0), ϕ(−τ)), donde A=fx(x0, x0)yB=fy(x0, x0). De aquí sigue el que el sistema linealizado para x=x0correspondiente a (2.27) es x0(t) = Ax(t) + Bx(t−τ),(2.28) caso que ya hemos estudiado. 2.3.5. Estabilidad absoluta Anticipándonos a lo que nos sucederá en los modelos epidemiológicos con retardo, trataremos el caso en el que la ecuación característica es de la forma p(λ) + q(λ)e−τλ = 0,(2.29) donde pyqson polinomios con coeficientes en Ryτ > 0es el retardo. Este caso se da, por ejemplo, en sistemas bidimensionales con det(B)=0, algo muy común en caso de que el sistema tenga un solo argumento retardado. Además, habitualmente, pes de grado mayor que q. Enunciemos el siguiente resultado. Proposición 2.21. Sean p,qpolinomios con coeficientes reales. Supongamos que p(λ)6= 0,Re(λ)≥0, |q(iy)|<|p(iy)|,0≤y < ∞, 2.3. ESTUDIO CUALITATIVO 21 l´ım |λ|→∞ |q(λ)/p(λ)|= 0. Entonces, Re(λ)<0para toda raíz λy para todo τ≥0. La conclusión de la proposición recibe el nombre de estabilidad absoluta, puesto que la estabilidad se mantiene, bajo esas hipótesis, para cualquier retardo. Por último, simplificaremos las hipótesis del resultado anterior para los casos en que q(λ)es constante. Esto será útil en los casos el los que, al linearizar el sistema de interés, la dependencia de λdesaparezca en los términos que multiplican a e−λτ . Corolario 2.22. Sea pun polinomio con coeficientes reales y coeficiente principal 1. Sea q=cuna constante. Si todas las raíces de pson reales y negativas y |p(0)|>|c|, o p(λ) = λ2+aλ +b,a, b > 0, y •b > |c|ya2≥2b, o •ap(4b−a2)>2|c|ya2<2b, entonces Re(λ)<0para toda raíz λy para todo τ≥0. Criterios de Routh-Hurwitz Para terminar con el marco teórico, y ahora que hemos visto que la estabilidad de los sistemas de ecuaciones diferenciales con retardo también mantiene una estrecha relación con los autovalores del jacobiano, nos interesa establecer condiciones para decidir cuándo la parte real de los autovalores de la matriz es negativa. En este caso, recordemos, tenemos estabilidad asintótica. Debido a que los sistemas sobre los que trabajaremos serán de orden 2 o 3, incluimos a continuación los criterios que nos facilitarán el estudio de los mismos. Proposición 2.23. (Routh-Hurwitz en dimensión 2) Sea A∈ M2×2. Sus autovalores vienen dados por la ecuación λ2−tr(A)λ+ det(A)=0, donde tr(A)indica la traza de la matriz y det(A)su determinante. Se tiene que Re(λ)<0 para todo λsi, y solamente si, tr(A)<0ydet(A)>0. Proposición 2.24. (Routh-Hurwitz en dimensión 3) Sea la ecuación λ3+a1λ2+a2λ+a3= 0. Se tiene que Re(λ)<0para todo λsi, y solamente si, a1>0,a3>0ya1a2−a3>0. 22 CAPÍTULO 2. ECUACIONES DIFERENCIALES CON RETARDO Estos resultados se pueden generalizar al caso de una ecuación de grado n. Un estudio detallado se encuentra en [9, sec. B.1]) Capítulo 3 Modelos epidemiológicos 3.1. Introducción El objetivo de este capítulo será el de aplicar los conceptos que hemos visto hasta ahora en un área que tiene una relación muy estrecha con la situación actual: la epidemiología. Esta disciplina médica se encarga del estudio de la distribución, frecuencia y factores determinantes en el desarrollo de las enfermedades que afectan a distintas poblaciones humanas. Por lo tanto, se tratarán diversos modelos más o menos adecuados para la enfermedad que abordaremos en el próximo capítulo: la causada por el SARS-CoV-2, comúnmente llamada COVID-19. Debemos tener en cuenta que no solamente la enfermedad a estudiar es determinante, si no que también se debe especificar la población sobre la que se realiza el estudio, puesto que las diferencias entre los diversos continentes, países o incluso zonas dentro de los mismos, tanto a nivel de desarrollo como de distribución de la población influyen de modo considerable en la manera de evolucionar que tiene la enfermedad. Sin ir más lejos, y teniendo en cuenta que en el siguiente capítulo trabajaremos siempre con los datos de Italia, es conocido que Lombardía fue la región con más casos de COVID-19, sobre todo al inicio de la pandemia. También es la más densamente poblada del país, y sus tres aeropuertos están entre los cuatro más importantes de Italia. Teniendo en cuenta las características de la enfermedad considerada, es obvio que las grandes afluencias de gente contribuyen a la expansión de la misma. Concluimos, por tanto, que debemos siempre especificar la población en la que se realiza el estudio. Cuando se plantea la modelización matemática de un problema real, debemos decidir en primer lugar cuál es el enfoque que queremos darle: uno probabilístico o uno determinístico. En el primero de ellos se asume que existe una componente de azar, mientras que, en el segundo, se supone que los datos de partida se conocen con certeza. En nuestro caso, escogeremos un enfoque determinístico. Dentro del mismo, asumiremos que la variable 23 24 CAPÍTULO 3. MODELOS EPIDEMIOLÓGICOS temporal es continua. Cabe notar que no es la única elección que podemos hacer, basta pensar en que los fenómenos que hacen variar las poblaciones tengan lugar en tiempos discretos y obtener así sistemas de ecuaciones en diferencias. Por último, también podemos plantearnos la continuidad de las funciones que rigen las poblaciones. Es decir, plantearnos de nuevo si las variables que representan el número de individuos enfermos, sanos, etc. son discretas o continuas. Asumiremos su continuidad (es decir, no estableceremos restricciones sobre las soluciones). Por tanto, en resumen, los modelos epidemiológicos que estudiaremos en esta sección se basarán en dos hipótesis fundamentales: la hipótesis determinista del modelo y la continuidad de las variables. Dentro de la modelización epidemiológica, en este trabajo trataremos siempre modelos compartimentales, es decir, dividiremos al conjunto total de la población en distintos conjuntos y estudiaremos las relaciones que se establezcan entre ellos. Los compartimentos serán característicos del modelo a estudiar y serán más o menos específicos dependiendo del mismo. Dos posibles referencias para el estudio de modelos epidemiológicos son [10, 11]. Demos paso ahora al estudio de modelos sin retardo, y veremos una forma natural de introducir el mismo posteriormente. 3.2. Del modelo SIR al SEIR El modelo SIR epidémico El modelo más conocido en epidemiología es el SIR. A pesar de que hay otros más simples (SI, SIS), el SIR, desarrollado a principios del siglo XX por William Kermack y A. G. McKendrick en A contribution to the mathematical theory of epidemics, [12], es el modelo más famoso dentro del estudio epidemiológico. Este divide a la población en 3 compartimentos, S,I,R, donde S=S(t)son los individuos susceptibles de padecer la enfermedad, I=I(t)los infectados por la misma y, por último, R=R(t)son los recuperados. Debemos observar que, una vez se supera la enfermedad, un individuo permanece de ahí en adelante en R. Es decir, de forma esquemática tenemos: Sβ −→ Iα −→ R En base a esto, podemos establecer el punto de inicio de nuestro estudio en el siguiente sistema de ecuaciones diferenciales, siendo Nla población total:        S0=−β NSI, I0=β NSI −αI, R0=αI. (3.1) 3.2. DEL MODELO SIR AL SEIR 25 No debemos olvidarnos de que, siguiendo el modelo de McKendrick y Kermack, estamos suponiendo que Nse mantiene constante, por lo que no tenemos en cuenta ninguna dinámica de poblaciones. ¿Qué significa el sistema (3.1)? Por un lado, αyβson coeficientes del modelo, que se determinan empíricamente en cada enfermedad. Sobre el primero de ellos, 1 αse interpreta como el tiempo medio de duración de la enfermedad, es decir, el tiempo medio que un individuo permanece en el compartimento I. Sobre β, este se interpreta como el número medio de contactos "efectivos" entre una persona susceptible y otra infectada por unidad de tiempo. Por lo tanto, β Nes la tasa de contagio de la enfermedad. En lo que respecta al significado de las ecuaciones, se tiene que la primera de ellas, S0=−β NSI expresa la variación de individuos susceptibles que pasan a ser infectados. Esto se debe a que, para que haya un contagio, tienen que suceder dos cosas: la primera de ellas es que haya un encuentro entre un individuo en Sy otro en I, hecho representado en el producto SI (Ley de acción de masas), mientras que la segunda de ellas es que ese contacto de lugar a un contagio. Esto se determina a partir de β Nde acuerdo con el significado que le acabamos de dar. Obviamente, con respecto al término −αI que aparece en la segunda ecuación y, en vista de lo que αrepresenta, es obvio que expresa el número de individuos que pasan al compartimento Runa vez superado el tiempo medio de la enfermedad. Realizaremos la siguiente simplificación del modelo: si consideramos s(t) = S(t)/N,i(t) = I(t)/N yr(t) = R(t)/N, el sistema que se obtiene es el siguiente.        s0=−βsi, i0=βsi −αi, r0=αi. (3.2) En lugar de considerar el número de individuos en cada clase, ahora estamos trabajando con proporciones dentro de la población total. Es por ello que s+i+r= 1, ya que Nes constante. Renombrando las variables s−→ S,i−→ Iyr−→ R, obtenemos el sistema sobre el que trabajaremos:        S0=−βSI, I0=βSI −αI, R0=αI. (3.3) Esto lo haremos siempre de forma implícita en el resto de modelos. Pasemos ahora al estudio del sistema. Tenemos que S+I+R= 1, y, por tanto, el número de individuos en Rlo conocemos si tenemos determinado el valor de Sy de I. 32 CAPÍTULO 3. MODELOS EPIDEMIOLÓGICOS S I R D Suponiendo de nuevo que Nes constante, se tendrá S+I+R+D=N. Veamos qué significa el esquema anterior. Por un lado, la población en Sse moverá a Imediante una tasa βque mantiene el mismo significado que hasta ahora. La diferencia viene después. Se introduce en este modelo una tasa de mortalidad de la enfermedad p. En este caso, no toda la población infectada se recupera, si no que solo lo hace una parte de ellos. Coherentemente, el sistema de ecuaciones asociado a este modelo y ya normalizado resulta:            S0=−βSI, I0=βSI −αI, D0=pαI, R0= (1 −p)αI. (3.9) Obviamente, las condiciones iniciales no difieren de las consideradas hasta ahora, S(0) = S0≥0,I(0) = I0≥0,R(0) = R0≥0yD(0) = D0≥0. Debemos notar que el añadir un compartimento hace que el sistema reducido (que podemos definir, puesto que ahora S+I+R+D= 1) no tendrá dos ecuaciones, si no tres. Se tiene el siguiente sistema.        S0=−βSI, I0=βSI −αI, D0=pαI. (3.10) Podemos considerar ahora la región Td={(S, I, D)∈R3:S≥0, I ≥0, D ≥ 0, S +I+D≤1}.Tdes positivamente invariante para el sistema, y siguiendo el razonamiento que hemos hecho con los modelos anteriores, la solución de (3.10) y, consecuentemente, de (3.9), existe, es única y está definida para todo t≥0. Estudiemos entonces los puntos de equilibrio de (3.10) y su comportamiento en Td. Igualando a cero las ecuaciones, se tiene que los equilibrios se obtienen cuando I= 0. Por lo tanto podemos ya sacar la primera conclusión del modelo. Incluso introduciendo muertes por enfermedad, independientemente del valor de los parámetros αyβy de la tasa de mortalidad p, se tiene que, cuando t→ ∞, la enfermedad tiende a desaparecer. Recordemos que esto es lo que sucedía en el caso SIR epidémico. El objetivo ahora será el de estudiar qué sucede en Tdcon SyDantes de que la enfermedad desaparezca. De nuevo, 3.2. DEL MODELO SIR AL SEIR 33 tenemos que Ses decreciente para todo t≥0, y Des creciente para todo t≥0. En lo que respecta a I, de forma totalmente análoga al modelo SIR endémico, crece hasta que S=α βy después decrece. Por lo tanto, R0=β αy se tendrá que si R0>1tendremos una situación de epidemia, mientras que en el caso R0<1, no la tenemos. Hemos obtenido el mismo número básico de reproducción que en modelo de partida de este capítulo. El objetivo ahora será el de estudiar la estabilidad del punto de equilibrio (S∞,0, D∞), con S∞= l´ım t→∞ S(t)yD∞= l´ım t→∞ D(t), en Td. Las conclusiones que saquemos del mismo se podrán extender al punto de equilibrio del sistema completo (S∞,0, R∞, D∞)libre de enfermedad, donde R∞= l´ım t→∞ R(t). Estudiemos el linealizado del sistema reducido, obteniendo el siguiente jacobiano: A=    −βI −βS 0 βI βS −α0 0pα 0     .(3.11) Debido a que en el caso SIR epidémico habíamos estudiado el linealizado del sistema completo y a la estructura de la matriz A, se tiene que los autovalores son los mismos (el término pα no aparece en det(A−λI)). Por tanto, sustituyendo previamente en (S∞,0, D∞), o bien λ= 0, o bien λ=βS∞−α. En caso de que βS∞−α < 0, las trayectorias tienden al punto de equilibrio, mientras que si βS∞−α > 0, las trayectorias primero se alejan de los puntos de equilibrio, alcanzando el número máximo de infectados hasta que βS∞−α < 0. Hemos obtenido las mismas conclusiones que en el caso del SIR epidémico. ¿Qué otro enfoque podemos darle a un modelo que tenga en cuenta las muertes por enfermedad? Pensemos en el caso ya estudiado en que consideramos una tasa de nacimientos y muertes. Teniendo en cuenta que nuestro objetivo es siempre el de terminar estudiando cómo se comportan los modelos en la enfermedad COVID-19 y considerando que esta, viendo las consecuencias de la vacunación, tiene un tiempo de vida bastante limitado, podemos asumir que el número total de la población no varía. Es decir, ni las muertes provocadas por la enfermedad ni los nacimientos y muertes naturales tienen un impacto relevante en N. Es por ello por lo que estamos asumiendo hipótesis que realmente no podríamos pasar por alto en otro tipo de enfermedades que llevan presentes en la población muchos años. Por lo tanto, el otro punto de vista que podemos aportar para introducir las muertes por enfermedad es ampliar el modelo estudiado con tasa de nacimientos y muertes, y equilibrarlos de tal forma que la tasa de nacimientos sea igual a la suma de la tasa de muertes de forma natural y la tasa de mortalidad de la enfermedad a estudiar. De esta manera, obtenemos que, nuevamente, la población se mantiene constante. El esquema del modelo es el siguiente: 34 CAPÍTULO 3. MODELOS EPIDEMIOLÓGICOS ↓µ+pαI S−→ I−→ R ↓µS ↓pαI ↓µI ↓µR El sistema de ecuaciones asociado al modelo resulta:        S0=−βSI +µ(1 −S) + pαI, I0=βSI −pαI −µI −(1 −p)αI =βSI −µI −αI, R0= (1 −p)αI −µR. (3.12) Los parámetros se mantienen de acuerdo con todo lo visto hasta el momento, y es obvio que es el sistema más sofisticado que hemos visto hasta ahora. Se añaden las condiciones iniciales S(0) = S0≥0,I(0) = I0≥0,R(0) = R0≥0para tener el sistema bien definido y se considera la región positivamente invariante Tenf ={(S, I)∈R2:S≥0, I ≥ 0, S +I≤1}para el sistema reducido (S0=−βSI +µ(1 −S) + pαI, I0=βSI −µI −αI. (3.13) La región Tenf es la misma que Tend. Veamos que también las conclusiones que se obtienen son las mismas, simplemente cambiando el punto de equilibrio endémico. Calculemos en primer lugar los equilibrios del sistema. Por un lado, existe el equilibrio libre de enfermedad (1,0). Por otro lado, tenemos un equilibrio endémico del sistema reducido, dado por Se=1 σ, Ie=µ(σ−1) β−σpα ,con σ=β α+µ.(3.14) Siguiendo con el procedimiento realizado hasta ahora, podemos asegurar que R0=β α+µ. Simplemente basta fijarse en que, de nuevo, S(t)es decreciente para todo t(debido a que S≤1) y que el crecimiento o decrecimiento de I(t)depende de R0. Por lo tanto, concluimos que el hecho de añadir muertes por enfermedad no altera el valor de R0. Estudiemos la estabilidad de los puntos de equilibrio en función a los autovalores del jacobiano del linealizado. Se obtiene: A= −βI −µ−βS +pα βI βS −α−µ!.(3.15) Por un lado, sustituyendo en (1,0),det(A−λI) = (−µ−λ)(β−µ−α−λ). Por tanto, se tiene que, bajo la hipótesis R0<1, los autovalores de Ason ambos negativos y, por lo tanto, (1,0) es asintóticamente estable para el sistema reducido. En el caso (Se, Ie), se obtiene que tr(A)<0ydet(A) = −βµ(σ−1) β−σpα (−(α+µ) + pα)>0si 3.2. DEL MODELO SIR AL SEIR 35 R0>1. Entonces, por el Criterio de Routh-Hurwitz (Proposición 2.23), la parte real de todos los autovalores es negativa y el equilibrio es asintóticamente estable. Con todos los modelos desarrollados hasta ahora hemos dado solución, al menos parcialmente, a dos de las carencias del SIR propuesto por Kermack y McKendrick: la población no tiene por qué ser siempre la misma, si no que se pueden introducir individuos nuevos y se pueden eliminar otros (a pesar de que hemos trabajado siempre bajo la hipótesis de que el número de individuos de la población se mantiene constante), y también se ha tenido en cuenta la mortalidad del virus. Pasemos ahora a una cuestión que nos permitirá tener en cuenta el tiempo de incubación del virus. El modelo SEIR Centrémonos en el modelo SEIR. La motivación de elegirlo se debe a que, en primer lugar, de la enfermedad a tratar en el capítulo siguiente conocemos que existe un tiempo de latencia en el que el individuo no es contagioso, es decir, el individuo permanece en el compartimento Ede expuestos antes de pasar a I(la novedad de este modelo con respecto al SIR). En segundo lugar, este modelo se presta a introducir de modo natural un retardo temporal, puesto que la clase Epuede ser sustituida con un retardo temporal en I. Una discusión interesante versará sobre las diferencias entre introducir un nuevo compartimento Eo un retardo en I. Vayamos entonces con el estudio del método SEIR. Partamos de la situación ilustrada en el esquema: ↓µ S−→ E−→ I−→ R ↓µS ↓µE ↓µI ↓µR Podemos entonces escribir el sistema:            S0=µ−βSI −µS, E0=βSI −(µ+κ)E, I0=κE −(µ+α)I, R0=αI −µR, (3.16) donde βestá definido como hasta ahora, µes la tasa de natalidad y mortalidad, 1/α es el tiempo medio de periodo infeccioso y 1/κ es el tiempo medio de incubación. Con la elección de que la tasa de natalidad y mortalidad coincidan, estamos suponiendo, de nuevo, que N es constante. Además, añadiremos al sistema (3.16) las condiciones iniciales S(0) = S0≥0, E(0) = E0≥0,I(0) = I0≥0yR(0) = R0≥0. De acuerdo con lo que ya hemos dicho, S(t) + E(t) + I(t) + R(t) = 1. 36 CAPÍTULO 3. MODELOS EPIDEMIOLÓGICOS Pasemos al estudio cualitativo. Por un lado, en el momento en que µ > 0, se tiene que el conjunto [0,∞)4es invariante para el sistema y, puesto que S+E+I+R= 1, podemos considerar el sistema reducido        S0=µ−βSI −µS, E0=βSI −(µ+κ)E, I0=κE −(µ+α)I. (3.17) El conjunto [0,∞)3resulta ahora invariante para este sistema. Además, podemos ser más precisos y considerar el compacto T=(S, E, I)∈R3:S, E, I ≥0, S +E+I≤1, que es positivamente invariante para el sistema (3.17). De aquí se deduce la existencia global y unicidad de la solución para el sistema reducido y, consecuentemente, para el sistema completo. El objetivo ahora es estudiar la estabilidad de los puntos de equilibrio. Para ello, igualamos las ecuaciones a 0. Por un lado, se obtiene un estado libre de enfermedad, (1,0,0) del sistema reducido, que se extiende a (1,0,0,0) en el sistema completo. Por otro lado, se alcanza un equilibrio endémico dado por Se=1 σ, Ee=µ β(σ−1), Ie=(α+µ)µ βκ (σ−1),con σ=κβ (κ+µ)(α+µ).(3.18) Se tiene, como hasta ahora, que R0=σ. Basta de nuevo con fijarse en que Ses decreciente y que los intervalos de crecimiento de Ey de Iestán determinados en base a σ. Obviamente, para que el estado de equilibrio endémico esté dentro de la región de interés, necesitamos que σo, equivalentemente, R0sea mayor que 1. Una vez determinados los puntos de equilibrio del sistema (todavía trabajando con el reducido), estudiemos la estabilidad. Siempre mediante el método de la primera aproximación, partimos con el jacobiano de la linealización del sistema. Se tiene que A=    −βI −µ0−βS βI −(µ+κ)βS 0κ−(µ+α)     .(3.19) Para el punto (1,0,0) se tiene que, por el Criterio de Routh-Hurwitz (Proposición 2.24), si R0<1, la parte real de todos los autovalores es negativa y, por tanto, el equilibrio es asintóticamente estable. Por el mismo resultado, para el equilibrio endémico (Se, Ee, Ie)se tiene que, en caso de que R0>1, todos los autovalores tienen parte real negativa y, de nuevo, el equilibrio es asintóticamente estable. 3.3. MODELOS EPIDEMIOLÓGICOS CON RETARDO 37 Veamos ahora que, bajo la hipótesis R0<1,(1,0,0) atrae a todo el compacto T. Para ello, consideramos la siguiente función de Liapunov: V(S, E, I) = κE + (κ+µ)I. (3.20) Se tiene que V(1,0,0) = 0,V(S, E, I)>0para todo punto distinto del equilibrio y V∗(S, E, I) = βIE S−1 σ≤0, puesto que σ < 1en este caso. Por tanto, Ves una función de Liapunov asociada al sistema. Además, viendo que la única solución de V∗(S, E, I)=0 en la región de interés es el punto de equilibrio (1,0,0), y que Tes compacto y positivamente invariante, se tiene que el punto de equilibrio libre de enfermedad atrae todo T. No solo el equilibrio (1,0,0) verifica esta propiedad, si no que, en el caso R0>1, se tiene también que el equilibrio endémico atrae a toda la región T. La demostración de este resultado y su generalización están en [15]. Para finalizar, fijémonos en el número básico de reproducción R0=κβ (κ+µ)(α+µ). Si µ= 0, se tiene que R0=β α, es decir, volvemos al caso del SIR endémico. Sin embargo, excluyendo el caso µ= 0, se tiene que si el tiempo de latencia es muy grande, entonces κ se hará más pequeño. En este caso, R0decrecerá y, por lo tanto, se reducirá también la posibilidad de que haya una epidemia. Por tanto, las enfermedades con tiempos de latencia muy grandes son menos propensas a convertirse en epidemias. 3.3. Modelos epidemiológicos con retardo Hemos terminado la sección anterior solucionando una de las carencias del modelo SIR introduciendo el compartimento de expuestos. Nos formulamos la siguiente pregunta: ¿podemos evitar la creación de un nuevo compartimento y, consecuentemente, evitar añadir una nueva ecuación al modelo, pero solucionando igualmente el problema? Teniendo en cuenta el marco teórico que hemos proporcionado en el segundo capítulo del trabajo, todo hace pensar que sí. El compartimento Ese convertirá ahora en un retardo temporal fijo y discreto con respecto a la variable I. Este retardo corresponderá con el tiempo de latencia de la enfermedad. Viendo que en el modelo SEIR consideramos Nconstante, pero también una tasa de nacimientos y muertes (que suponemos igual), heredaremos esta elección para el modelo que definiremos a continuación. Posibles referencias para la estabilidad de los equilibrios en este modelo son [16, 17, 18]. Por un lado, eliminando el compartimento de expuestos, volvemos al modelo SIR. Por lo tanto, el que propondremos a continuación podría recibir el nombre de modelo SIR con retardo discreto. A nivel esquemático, y sin tener en cuenta el retardo, se tiene la siguiente 38 CAPÍTULO 3. MODELOS EPIDEMIOLÓGICOS dinámica compartimental: ↓µ S−→ I−→ R ↓µS ↓µI ↓µR En efecto, esta es la misma que la del modelo SIR con tasa de nacimientos y muertes que fue estudiada en la sección anterior. El sistema asociado al modelo resulta        S0(t) = −βS(t)I(t−τ) + µ(1 −S(t)), I0(t) = βS(t)I(t−τ)−αI(t)−µI(t), R0(t) = αI(t)−µR(t). (3.21) El significado de este modelo es obvio si hemos entendido los anteriores. Los coeficientes mantienen el mismo significado que hasta ahora, seguimos teniendo S+I+R= 1 (la población total se mantiene constante) y la única diferencia recae sobre el parámetro τ que, como hemos explicado, es el tiempo de latencia. Un individuo contrae la enfermedad tras un contacto, pero permanece un tiempo τsin contagiar a otros. Definamos las condiciones iniciales como S(θ) = ϕ1(θ),I(θ) = ϕ2(θ)yR(θ) = ϕ3(θ), de forma que ϕi(θ)sea continua en [−τ, 0] yϕi(θ)≥0,i∈ {1,2,3}. Una vez con el problema definido, y como S+I+R= 1, sabemos que la solución de Res conocida si tenemos determinadas las soluciones de Sy de I. Por lo tanto, pasaremos ya al estudio del sistema reducido siguiente: (S0(t) = −βS(t)I(t−τ) + µ(1 −S(t)), I0(t) = βS(t)I(t−τ)−αI(t)−µI(t).(3.22) Comprobemos, en primer lugar, que este sistema tiene solución, es única y, además, dado que las condiciones iniciales son no negativas, la solución también lo es. Escribamos el sistema como X0=F(X, Xτ), donde X= (S, I)TyXτ= (Sτ, Iτ)T. Por un lado, se tiene que F= (F1, F2)T, con F1(X, Xτ) = −βSIτ+µ(1−S)yF2(X, Xτ) = βSIτ−αI −µI. Es obvio que cada componente de Fes continua y, por tanto, Flo es. Además, derivando se tiene: ∂F1 ∂S =−βIτ−µ, ∂F1 ∂I = 0,∂F2 ∂S =βIτ,∂F2 ∂I =−α−µ. Todas ellas son continuas y, por lo tanto, tenemos garantizada, en virtud del Teorema 2.7, existencia y unicidad de solución del problema reducido (y, consecuentemente, del inicial). A mayores, podemos ver que si S= 0, entonces F1(X, Xτ) = µ≥0, y si I= 0, entonces F2(X, Xτ) = βSIτ≥0. Se tiene, por el Teorema 2.10, que las soluciones son no negativas en todo el dominio de definición. 3.3. MODELOS EPIDEMIOLÓGICOS CON RETARDO 39 Pasemos ahora al estudio cualitativo de este sistema. Para ello, en primer lugar calculemos los puntos de equilibrio. Por definición, necesitamos que estos sean estacionarios, es decir, que no varíen con el tiempo. Por tanto, podemos igualar a 0 las ecuaciones ignorando el retardo. Obtenemos así dos equilibrios: uno libre de enfermedad (1,0) y otro endémico dado por Se=1 σ, Ie=µ(σ−1) β,con σ=β α+µ.(3.23) Considerando la región Tτ=(S, I)∈R2:S, I ≥0, S +I≤1, que es la que nos interesa estudiar, tendremos que el equilibrio endémico pertenece al interior de Tτsi, y solamente si, σ > 1. Esta elección de σno es casual, y, como hasta ahora, coincidirá con el número básico de reproducción del modelo. Pasemos ahora a estudiar la estabilidad de los equilibrios. Para ello, hagamos el linealizado del sistema. Por un lado, y en lo que respecta a las variables sin retardo, tenemos: J= −βIτ−µ0 βIτ−α−µ!.(3.24) Por otro lado, derivando con respecto a las variables con retardo: Jτ= 0−βS 0βS !.(3.25) Comencemos por el equilibrio libre de enfermedad. En él, se tiene que el sistema linealizado es el siguiente: X0(t) = JX(t) + JτX(t−τ),(3.26) donde J= −µ0 0−α−µ! y Jτ= 0−β 0β!. La ecuación característica es h(λ, τ) = det[λI −J−e−λτ Jτ]. Si todas las raíces de h, es decir, todas las soluciones de h(λ, τ)=0verifican que Re(λ)<0, se tiene que el punto de equilibrio es asintóticamente estable. Escribamos la matriz A=λI −J−e−λτ Jτ: A= λ+µ βe−λτ 0λ+α+µ−βe−λτ !. 40 CAPÍTULO 3. MODELOS EPIDEMIOLÓGICOS Obtenemos que h(λ, τ) = det(A)=(λ+µ)(λ+α+µ−βe−λτ ).(3.27) Por tanto, λ1=−µ < 0yλ2es la solución de λ=−α−µ+βe−λτ . Como µ > 0, para la estabilidad asintótica nos podremos centrar en λ2. Usando el Teorema 2.19, se tiene que, bajo la condición β α+µ<1, el equilibrio es asintóticamente estable. Por lo tanto, juntando las dos conclusiones sobre los autovalores, y teniendo en cuenta esa hipótesis, la asintótica estabilidad está garantizada. Ahora bien, la restricción obtenida es equivalente a decir σ < 1. En resumen, tenemos que si σ < 1, existe un equilibrio asintóticamente estable que es libre de enfermedad. Tiene sentido ya lo que anticipábamos: σcoincide con el número básico de reproducción R0. En caso de que se cumpla R0>1, tendremos otro equilibrio, el endémico. Estudiemos su estabilidad asintótica siguiendo el método que hemos usado para el punto (1,0). En (Se, Ie)el sistema linealizado sigue siendo X0(t) = JX(t) + JτX(t−τ),(3.28) donde ahora J= −µσ 0 µ(σ−1) −α−µ! y Jτ= 0−β σ 0β σ!. La ecuación característica asociada al linealizado alrededor del equilibrio endémico es g(λ, τ) = det[λI −J−e−λτ Jτ]. Estudiemos las raíces de g. Ahora B=λI −J−e−λτ Jτ resulta: B= λ+µσ β σe−λτ −µ(σ−1) λ+α+µ−βe−λτ !. Entonces, g(λ, τ) = det(B)=(λ+µσ)(λ+α+µ−β σe−λτ ) + µ(σ−1)β σe−λτ .(3.29) Entra aquí en juego la estabilidad absoluta que tratamos en el capítulo previo. La función g(λ, τ)puede escribirse también como g(λ, τ)=(λ+µσ)(λ+α+µ) + β σ(−µ−λ)e−λτ .(3.30) 3.3. MODELOS EPIDEMIOLÓGICOS CON RETARDO 41 Poniendo p(λ) = (λ+µσ)(λ+α+µ)yq(λ) = β σ(−µ−λ), comprobemos que se verifican las hipótesis del Teorema 2.21. p(λ)=0no tiene soluciones con Re(λ)≥0. Esto podemos verlo directamente, puesto que p(λ) = (λ+µσ)(λ+α+µ), y todos los parámetros son positivos. Tenemos que |p(iy)|=p(y2+ (µσ)2)(y2+ (µ+α)2)y|q(iy)|=qβ2 σ2(µ2+y2). Teniendo en cuenta que α+µ=β σy que estamos trabajando bajo la hipótesis de que σ > 1(solo así este punto de equilibrio está en la región de interés), se verifica que |q(iy)|<|p(iy)|para todo y∈[0,∞). l´ım|λ|→∞,Re(λ)≥0|q(λ)/p(λ)|= 0, puesto que el grado del polinomio p(λ)es mayor que el de g(λ). Hemos comprobado que se satisfacen las hipótesis del Teorema 2.21. Por ello, Re(λ)<0 para toda raíz λde la ecuación característica (3.29). En resumen, si σ > 1o, equivalentemente, R0>1, el equilibrio endémico (Se, Ie)es asintóticamente estable. Para finalizar, recordar que al principio de la sección consideramos el sistema reducido. Estas conclusiones las podemos extender al sistema completo y decir que el equilibrio (1,0,0) es asintóticamente estable para el sistema si R0<1(la enfermedad desaparece) y el equilibrio endémico (Se, Ie, Re), donde Re= 1 −Se−Ie, es asintóticamente estable para el sistema si R0>1. ¿Qué hemos ganado con respecto al modelo SEIR? La primera diferencia notoria es que, siempre bajo la hipótesis de mantener un número de población constante, el no introducir un nuevo compartimento nos permite reducir el sistema y enfocar nuestro estudio en un espacio bidimensional. Por otro lado, y ya centrándonos en el ámbito de la epidemiología, el número básico de reproducción, que hace de umbral entre la situación de epidemia o no, no depende del tiempo de latencia si se introduce un retardo, mientras que sí lo hacía cuando introducíamos el compartimento de expuestos. Esto se debe a que R0lo estamos siempre calculando en base a los equilibrios del sistema. Cuando consideramos modelos con retardo como el que hemos estudiado, buscamos soluciones estacionarias, es decir, constantes con respecto al tiempo. Por ello, el retardo no influye y el número básico de reproducción no lo tiene en cuenta. 48 CAPÍTULO 4. LA PANDEMIA DEL CORONAVIRUS EN ITALIA ello, podemos realizar lo siguiente: teniendo en cuenta que el tiempo medio que un individuo permanece en el compartimento de expuestos antes de pasar al de infectados es de 5 días, tenemos que los datos de nuevos infectados de los 5 días posteriores al inicial habrán estado en Eel 24 de febrero, nuestro día de partida. Entonces, sumando y dividiendo por el total de la población, podemos considerar E0= 1.489 ×10−5. Como hasta ahora, I0= 3.66 ×10−6. Por lo tanto, S0= 1 −I0−E0. Se obtiene: Figura 4.5: Modelo SEIR. Y se tiene la Figura 4.6, donde se muestran las gráficas que representan la población infectada con respecto al tiempo. (a) Solución numérica. (b) Datos reales. Figura 4.6: Comparación de la población infectada con el modelo SEIR. Con respecto a la Figura 4.5, podemos ver que el modelo tiende al equilibrio endémico, 4.1. SIMULACIÓN NUMÉRICA 49 hecho coherente con lo que hemos visto teóricamente. Viendo la Figura 4.6, sorprende que la gráfica se haya desplazado tanto hacia la derecha. A pesar de que era algo esperable por el significado del parámetro κ(para ser un individuo en I, es necesario pasar un tiempo 1 κen E), intentaremos dar más adelante una respuesta a este desplazamiento. Siguiendo con la misma figura, podemos observar que la distribución de los nuevos infectados en el tiempo que nos proporciona el modelo SEIR es más parecida a la que obtenemos mediante los datos reales. En este caso, la concavidad de la curva es menos pronunciada que en los casos que habíamos visto. Esto es consecuencia de añadir el compartimento E, ya que los contagios no se producen de forma inmediata, como sí lo hacían en el modelo SIR. SIR con retardo Como último modelo antes de pasar a tener en cuenta las medidas de control, realizamos una simulación del modelo SIR con retardo. En cuanto a los parámetros del modelo, podemos reutilizar los calculados para el SIR con muerte por enfermedad (pno aparece en el R0de ese modelo). Por tanto, α= 0.05,µ= 0.01 yβ= 0.3. Obviamente, τ= 5. A diferencia de los modelos anteriores, la determinación de la condición inicial no es inmediata. Necesitamos conocer el comportamiento del sistema en los tiempos entre [−τ, 0] = [−5,0]. A diferencia del modelo anterior, en el que podíamos saber qué pasaba en el tiempo inicial conociendo los datos de los 5 días siguientes, en el caso del retardo discreto se realiza un salto a estados pasados. Estos datos no están disponibles, así que simplemente supondremos que los datos iniciales ϕSyϕIson constantes en [−5,0] e iguales a las condiciones iniciales S0,I0que hemos considerado hasta ahora. Como solución alternativa, en caso de conocer los nuevos infectados por día de los últimos 5 antes del inicial, podríamos interpolar esos puntos mediante un polinomio y establecerlo como condición inicial. Figura 4.7: Modelo SIR con retardo. 50 CAPÍTULO 4. LA PANDEMIA DEL CORONAVIRUS EN ITALIA (a) Solución numérica. (b) Datos reales. Figura 4.8: Comparación de la población infectada introduciendo un retardo. En primer lugar vemos que el sistema tiende al equilibrio endémico, a pesar de que no lo haga de una forma tan veloz como en otros modelos sin retardo. De nuevo, observando la Figura 4.8, existe un traslado temporal de la función a la derecha que sabríamos que se tendría, pero que es bastante grande. Como hemos dicho antes, intentaremos dar una respuesta posteriormente. 4.1.2. Medidas de contención Pasemos ahora a introducir al modelo SIR con muertes por enfermedad el efecto de las restricciones. Recordando la ecuación (4.2), tenemos que determinar los valores de h(t)y k. Por un lado, daremos un ksuficientemente grande de modo que la población cumple las reglas. Pongamos, por ejemplo, k= 50. Por otro lado, y según las pequeñas nociones que hemos dado de las restricciones en Italia durante la primera ola, se tiene que el 9 de marzo (día 13) se establecen unas medidas que se endurecen el 23 del mismo mes (día 27). Podemos entonces definir la función h(t)a trozos como sigue: h(t) =        0t≤13, 0.2 13 < t ≤27, 0.4 27 < t ≤150. (4.3) Es decir, hasta el día 13 no hay medidas, del día 13 al 27 las medidas son pequeñas y, a partir del día 27, son más severas. Cabe destacar que estamos ignorando que a partir del 4 de mayo las medidas volvieron a relajarse. Vamos a introducir esta modificación en el modelo SIR con muertes por enfermedad. Usemos las condiciones iniciales I0= 3.66×10−6, 4.1. SIMULACIÓN NUMÉRICA 51 S0= 1 −I0y los parámetros α= 0.05,p= 0.08 yβ0= 0.3. Tenemos, entonces β(t)=0.3(1 −h(t)) (1 −0.04I(t))50 ,(4.4) con h(t)definida en (4.3). Los resultados obtenidos se muestran en las Figuras 4.9 y 4.10 Figura 4.9: Modelo SIR completo modificando β. (a) Solución numérica. (b) Datos reales. Figura 4.10: Comparación de la población infectada modificando β. En la Figura 4.9 se aprecia ya una notable diferencia con el modelo que habíamos considerado previamente sin modificar β. Mientras que en la Figura 4.3 el porcentaje de la población llegaba a alcanzar en 50 % en su punto más alto, en esta apenas alcanza el 30 %. Por lo tanto uno de nuestros objetivos está conseguido. En la comparación que se realiza en la Figura 4.10, vemos que esta se adecúa mucho más a la curva de los datos reales. De hecho, aparece la concavidad de la que se hablaba en los modelos en los que se introducían 52 CAPÍTULO 4. LA PANDEMIA DEL CORONAVIRUS EN ITALIA tiempos de latencia. Ciñéndonos a los datos, hemos reducido el valor máximo de infectados, aproximadamente, a la mitad. Por último, observamos que el crecimiento de la función de infectados es mucho más lento al principio que en el caso tratado previamente. A pesar de que no haya medidas hasta el tiempo t= 27 (por ello hemos puesto h(t) = 0 cuando t varía entre 0 y 17), la modificación de βresulta β(t)=0.3 (1 −0.04I(t))50. Como I(t)>0 para todo t, nuestro parámetro β0= 0.3comienza a disminuir incluso aunque no haya restricciones impuestas, debido al miedo y al autoaislamiento de la población por los casos críticos conocidos. 4.2. Conclusiones Tratemos en primer lugar las simulaciones que hemos hecho sin tener en cuenta las medidas de control existentes. Todas ellas tienen un punto en común: la diferencia en el número máximo de personas infectadas en comparación con los datos reales. Esto refuerza la idea que dimos al principio y que motivó la modificación sobre β. Por otro lado, debemos observar que en el momento en que introducimos los tiempos de latencia, tanto mediante el compartimento de expuestos como mediante el retardo, el tiempo en el que se alcanzaba el número máximo de infectados proporcionado por el modelo se distanció del tiempo real en que esto ocurría. Este desajuste puede deberse a que estamos considerando el inicio de la enfermedad. El desconocimiento de la existencia de personas asintomáticas y la falta de material sanitario llevaron a que solo se les hiciese un test a aquellas personas con síntomas. Esto implica, por un lado, que no sabemos cuántos infectados había realmente (lo que hace que nos planteemos si la diferencia entre el máximo número de infectados proporcionado por los modelos y el que obtenemos mediante los datos oficiales es tan grande como parece) y, además, que no tenemos determinado qué pasaba en el tiempo anterior al dato inicial del 24 de febrero. Por lo tanto, la condición inicial del sistema con retardo (que, recordemos, hemos supuesto constante), podría, como ya hemos dicho, perfeccionarse. Quizás, sobre todo si se tratan situaciones relativas al inicio de la pandemia, un enfoque más interesante sería el de considerar los hospitalizados por la enfermedad y no los individuos infectados. En lo relativo a las medidas de control, hemos visto que estas sirven para reducir el número de personas infectadas y también para hacer que los procesos infecciosos de la enfermedad estén más repartidos temporalmente, evitando así picos de incidencia y, consecuentemente, colapsos hospitalarios como hemos visto en el caso de esta pandemia. Debemos fijarnos en que el número real de infectados comienza a crecer antes de lo que lo hace la solución numérica del sistema modificado. Esto tiene una explicación, y es que, a 4.2. CONCLUSIONES 53 pesar de que el 9 de marzo se endurecen las medidas oficialmente en el conjunto del país, un error informático provocó que esta noticia estuviese disponible un día antes, provocando así un enorme movimiento de personas hacia el sur de Italia. Esto contribuyó a que la enfermedad se extendiera por todo el país y también que el número de contactos en ese día fuese muy alto. Es decir, en contra de lo que hemos supuesto, β0no siempre disminuyó, si no que alcanzó un pico el día previo a las primeras restricciones. Este razonamiento puede dar respuesta al retraso temporal del modelo con respecto a los datos reales. De todos modos, la modificación que se ha realizado depende de parámetros bastante subjetivos, como pueden ser la dureza de las medidas o la responsabilidad de la población. Como comentario final, podríamos decir que para realizar un estudio más riguroso de la situación sanitaria deberíamos añadir compartimentos en los que se tuvieran en cuenta, mediante los distintos nuevos parámetros, las restricciones vigentes en cada momento. 54 CAPÍTULO 4. LA PANDEMIA DEL CORONAVIRUS EN ITALIA Apéndice A Resultados de análisis complejo En este apéndice se incluyen definiciones y resultados de análisis complejo necesarios para la correcta comprensión del trabajo. Para un estudio completo, puede consultarse [24]. Definición A.1. Sea f: Ω ⊂C−→ C, con Ωun abierto. Decimos que fes holomorfa en a∈Ωsi existe y es finito el siguiente límite: l´ım h→0 f(a+h)−f(a) h. Decimos también que fes holomorfa en Ωsi lo es en todo punto de Ω. Definición A.2. Una función entera es una función holomorfa con dominio de definición C. Las propiedades de las funciones enteras no triviales son las siguientes: 1. Cada raíz característica tiene orden finito. 2. Existe, como máximo, una cantidad numerable de raíces características. 3. El conjunto de raíces características no tiene un punto de acumulación finito. Además, cabe notar que solo hay un número finito de raíces características con parte real positiva. Existen tres tipos de puntos que nos interesan a la hora de estudiar una función f compleja. Por un lado, los ceros de f, es decir, los puntos en los que se anula. Por otro lado, los puntos críticos. En ellos, se tiene que f0(a) = 0. Por último, nos interesan las singularidades de f. Detallemos estas últimas. 55 56 APÉNDICE A. RESULTADOS DE ANÁLISIS COMPLEJO Definición A.3. Sea Ω⊂Cun abierto, a∈Ωyf: Ω\{a} −→ Cholomorfa. Decimos que ftiene una singularidad eliminable en asi existe ˜ f: Ω −→ Cholomorfa tal que ˜ f|Ω\{a}=f. singularidad polar en asi existe n > 0tal que la función g: Ω\{a} −→ C, definida como g(z) = (z−a)nf(z), tiene una singularidad eliminable en ade valor no nulo. singularidad esencial en asi no es ni eliminable ni polar. El interés del estudio de las singularidades está en saber cómo se comporta fcuando se acerca a las mismas. Teorema A.4. (Teorema de Riemann) La función ftiene una singularidad eliminable en asi, y solamente si, existe un disco de radio rcentrado en a,Dr(a)⊂Ω, tal que sup{|f(z)|: 0 <|z−a|< r}<∞. Es decir, si la función está acotada lo suficientemente cerca de la singularidad, entonces esta es eliminable y fse extiende a una función holomorfa en Ω. Teorema A.5. La función ftiene una singularidad polar en asi, y solamente si, l´ım z→a|f(z)|= ∞. Por lo tanto, viendo estos dos resultados, se tiene que, alrededor de una singularidad esencial no tenemos caracterizado el comportamiento de f. Veamos el siguiente teorema. Teorema A.6. (Teorema de Casorati-Weierstrass) Si la función ftiene una singularidad esencial en a, entonces para todo ε > 0, con Dε(a)⊂ Ω, se tiene que f(Dε(a)\{a})es denso en C. Teorema A.7. (Teorema fuerte de Picard) Bajo las hipótesis del Teorema de Casorati-Weierstrass, la función falcanza cada valor complejo infinitas veces con, como máximo, una excepción. Es decir, para cualquier entorno de a, existe, como mucho, un punto del plano complejo que no esté en la imagen por fde dicho entorno. Por lo tanto, si ese punto existe y lo conocemos, sabemos que falcanza todos los demás un número infinito de veces. Como última observación, podemos decir que, así como una función fcon singularidad eliminable se extiende a otra función holomorfa ˜ f: Ω −→ C, en el caso de las singularidades polares, ftambién se extiende a una función holomorfa pero que ahora no tiene como codominio C, si no la recta proyectiva compleja. Es decir, ˜ f: Ω −→ P1(C). Por lo tanto, esto distingue las singularidades eliminables y polares de las esenciales. Apéndice B Código MATLAB utilizado En este apéndice mostraremos el código que se ha usado para las simulaciones en en los modelos SEIR, SIR con retardo y SIR con medidas de contención. SEIR 1c l e a r a l l 2c l c 3c l o s e a l l 4 5% i n t e r v a l o temporal 6tspan =[0 1 5 0 ] ; 7 8% parametros del modelo 9beta =0.315; 10 alpha =0.05; 11 k=0.2; 12 mu=0.01; 13 %comprobacion de R0 14 r0=k∗beta /(( alpha+mu) ∗(k+mu) ) 15 16 % vector de datos i n i c i a l e s 17 y0=[1 −0.0000036−0.00001489 0.00001489 0. 00 0 00 3 6] ; 18 19 % so lu ci on del sistema 57 64 BIBLIOGRAFÍA [13] Brauer, F., van den Driessche, P., Wu, J.; Mathematical Epidemiology, Springer, 2008. [14] Sen, D., Sen, D.; Use of a modified SIRD model to analyze COVID-19 data, Industrial and Engineering Chemistry Research, vol. 60, number 11, pp. 4251–4260, 2021. [15] Y. Li, M., Muldowney, J.; Global stability for the SEIR model in epidemiology, Mathematical Biosciences, vol. 125, issue 2, 1995. Pages 155-164. [16] Ma, W., Song, M., Takeuchi, Y.; Global Stability of an SIR Epidemic Model with Time Delay, Elseived Ltd, 2004. [17] McClouskey, C.C.; Complete global stability for an SIR epidemic model with delay — Distributed or discrete Nonlinear Analysis: Real World Applications, vol. 11, Issue 1, pages 55-59, 2010. [18] McClouskey, C.C.; An SIR Epidemic Model with Time Delay and General Nonlinear Incidence Rate, Abstract and Applied Analysis, vol. 2014, 2014 [19] Lina, Q., Zhaob, S., et al.; A conceptual model for the coronavirus disease 2019 (COVID-19) outbreak in Wuhan, China with individual reaction and governmental action, International Journal of Infectious Diseases 93 pages 211-216, 2020. [20] Datos de salud pública del Gobierno de España, disponibles en https://bit.ly/3hfdZRp. [21] Monitoraggio della situazione nazionale, https://bit.ly/3xhcJ5X. [22] Houcque, D.; Applications of MATLAB: Ordinary Differential Equations (ODE) , disponible en https://www.mccormick.northwestern.edu/docs/efirst/ode.pdf [23] Shampine, L. F., Thompson, S.; Solving DDEs in MATLAB, Elsevier Science B.V., 2001. [24] Stein, E., Shakarchi, R.; Complex Analysis, Princeton University Press, 2003.