Matemáticas y epidemias: modelos basados en ecuaciones diferenciales

Universidad de Sevilla Facultad de F´ısica Trabajo Fin de Grado: Matem´aticas y epidemias: modelos basados en ecuaciones diferenciales Autor: Pablo Navarro Torrero Tutor: Enrique Fern´andez Cara 3 de junio de 2022 ´ Indice general 1. Introducci´on a la modelizaci´on de epidemias 9 1.1. Historia y actualidad .............................. 9 1.2. La reproductividad R0............................. 11 1.3. Modelos epidemiol´ogicos ............................ 12 1.3.1. SIR .................................... 13 1.3.2. SIRS y vacunas ............................. 15 1.3.3. SEIR, SEIRS y vacunas ........................ 16 1.3.4. SEIRQ y m´as all´a ............................ 19 2. Estudio num´erico y simulaciones 21 2.1. Resultados para el modelo SIR ......................... 21 2.1.1. Confinamiento en el modelo SIR .................... 23 2.2. Resultados para el modelo SIRS ........................ 25 2.2.1. Confinamiento en el modelo SIRS ................... 26 2.2.2. Sensaci´on p´ublica de riesgo en el modelo SIRS ............ 29 2.2.3. Vacunas: modelo SIRS* ........................ 31 2.3. Resultados para el modelo SEIR ........................ 32 2.3.1. Confinamiento en el modelo SEIR ................... 33 2.3.2. Sensaci´on p´ublica de riesgo en el modelo SEIR ............ 34 2.4. Resultados para el modelo SEIRS ....................... 36 2.4.1. Confinamiento en el modelo SEIRS .................. 37 2.4.2. Sensaci´on p´ublica de riesgo en el modelo SEIRS ........... 40 2.4.3. Vacunas: modelo SEIRS* ........................ 41 2.5. Resultados para el modelo SEIRQ ....................... 42 2.5.1. Modificando el modelo SEIRQ ..................... 45 3 ´ Indice general 3. An´alisis y resoluci´on de problemas inversos 49 3.1. Procedimiento .................................. 49 3.2. Problema inverso: modelo SIR ......................... 51 3.3. Problema inverso: modelo SIRS ........................ 53 4 ´ Indice de figuras 1.1. Diagrama del modelo (SIR) .......................... 14 1.2. Diagrama del modelo (SIRS) .......................... 15 1.3. Diagrama del modelo (SIRS)∗......................... 16 1.4. Diagrama del modelo (SEIR) ......................... 17 1.5. Diagrama del modelo SEIRS .......................... 17 1.6. Diagrama del modelo SEIRS∗......................... 18 1.7. Diagrama del modelo SEIRQ ......................... 19 2.1. Evoluci´on de la epidemia seg´un el modelo SIR. ................ 22 2.2. Variaci´on del par´ametro βen el tiempo. ................... 23 2.3. Evoluci´on de la epidemia seg´un el modelo SIR con β=β(t), dada por (2.1). L´ıneas continuas: βvariable; l´ıneas de puntos: β= 0,125. .......... 24 2.4. Evoluci´on de la epidemia seg´un el modelo SIRS con tmax = 360 ....... 25 2.5. Evoluci´on de la epidemia seg´un el modelo SIRS con tmax = 900 ....... 26 2.6. Evoluci´on de la epidemia seg´un el modelo SIRS con β(t) y tmax = 360. L´ıneas continuas: βvariable; l´ıneas de puntos: β= 0,125. .......... 27 2.7. Variaci´on del par´ametro β(t) con el tiempo. ................. 28 2.8. Evoluci´on de la epidemia seg´un el modelo SIRS con β(t) y tmax = 900. L´ıneas continuas: βvariable; l´ıneas de puntos: β= 0,125. .......... 28 2.9. Evoluci´on de la epidemia seg´un el modelo SIRS con percepci´on del riesgo por parte de la poblaci´on. L´ıneas continuas: κ= 0; l´ıneas discontinuas: κ= 100; l´ıneas de puntos: κ= 200. ...................... 30 2.10. Comparaci´on de la evoluci´on de la epidemia seg´un el modelo SIRS con y sin vacunaci´on. L´ıneas continuas, con vacunaci´on; l´ıneas discontinuas, sin vacunaci´on. ................................... 31 5 ´ Indice de figuras 2.11. Evoluci´on de la epidemia seg´un el modelo SEIR ............... 32 2.12. Variaci´on del par´ametro βen el tiempo. ................... 33 2.13. Evoluci´on de la epidemia seg´un el modelo SEIR con β=β(t), dada por (2.6). L´ıneas continuas: βvariable; l´ıneas discontinuas: β= 0,125. ..... 34 2.14. Evoluci´on de la epidemia seg´un el modelo SEIR con percepci´on del riesgo por parte de la poblaci´on. L´ıneas continuas: κ= 0; l´ıneas discontinuas: κ= 100; l´ıneas de puntos: κ= 200. ...................... 35 2.15. Evoluci´on de la epidemia seg´un el modelo SEIRS con tmax = 360 ...... 36 2.16. Evoluci´on de la epidemia seg´un el modelo SEIRS con tmax = 1200 ..... 37 2.17. Variaci´on del par´ametro β(t) con el tiempo. ................. 38 2.18. Evoluci´on de la epidemia seg´un el modelo SEIRS con β(t) y tmax = 1200. L´ıneas continuas: βvariable; l´ıneas de puntos: β= 0,125. .......... 39 2.19. Evoluci´on de la epidemia seg´un el modelo SEIRS con percepci´on del riesgo por parte de la poblaci´on. L´ıneas continuas: κ= 0; l´ıneas discontinuas: κ= 100; l´ıneas de puntos: κ= 200. ...................... 40 2.20. Comparaci´on de la evoluci´on de la epidemia seg´un el modelo SEIRS con y sin vacunaci´on. L´ıneas continuas, con vacunaci´on; l´ıneas discontinuas, sin vacunaci´on. ................................... 41 2.21. Evoluci´on de infecciosos y expuestos seg´un el modelo SEIRS con y sin vacunaci´on. L´ıneas continuas, con vacunaci´on; l´ıneas discontinuas, sin vacunaci´on. 42 2.22. Evoluci´on de la epidemia seg´un el modelo SEIRQ con tmax = 540 ..... 43 2.23. Evoluci´on de la epidemia seg´un el modelo SEIRQ con p∈ {0,050; 0,058}. 44 2.24. Evoluci´on de IyEseg´un el modelo SEIRQ con p∈ {0,050; 0,058}. . . . 44 2.25. Diagrama del modelo SEIRQS∗........................ 46 2.26. Evoluci´on de la epidemia seg´un el modelo SEIRQS*. L´ıneas continuas: ρ= 0 y p= 0,01; l´ıneas discontinuas: ρ= 0 y p= 0,05; l´ıneas de puntos: ρ= 0,02 y p= 0,01. .............................. 47 2.27. Evoluci´on de infecciosos, expuestos y confinados seg´un el modelo SEIRQS*. 47 3.1. Resoluci´on num´erica modelo SIR ....................... 51 3.2. Resoluci´on num´erica problema inverso modelo SIR .............. 52 3.3. Resoluci´on num´erica modelo SIRS ....................... 53 3.4. Resoluci´on num´erica problema inverso modelo SIRS ............. 54 6 Introducci´on El principal objetivo de este trabajo es comprender los distintos modelos de epidemias, qu´e par´ametros los caracterizan, qu´e resultados se espera obtener, sus aplicaciones y sus limitaciones. El primer cap´ıtulo est´a dedicado a una presentaci´on de dichos modelos, incluyendo aspectos hist´oricos. En el segundo cap´ıtulo, mediante c´alculos num´ericos se simulan algunos de los modelos estudiados y se analizan los resultados obtenidos. Finalmente, en el tercer cap´ıtulo, se estudia la resoluci´on de problemas inversos aplicado a distintos modelos. 8 Cap´ıtulo 1 Introducci´on a la modelizaci´on de epidemias 1.1. Historia y actualidad Debido a la situaci´on actual en la que vivimos, pandemia de COVID-19, hoy en d´ıa se habla m´as que nunca, en todos los ´ambitos y sobre todo en los medios de comunicaci´on, sobre modelos matem´aticos, picos de contagiados, “aplanar” la curva, confinamientos, vacunaci´on, etc. Todo esto son conceptos relacionados con el estudio y modelizaci´on de epidemias, y se explicar´an a lo largo del trabajo. Para empezar, debemos entender que los seres humanos a lo largo de la historia han sufrido numerosas epidemias debidas a muy diversas enfermedades y agentes infecciosos. En la Tabla 1.1 se pueden ver las principales pandemias que asolaron a la humanidad junto a su fecha, localizaci´on, n´umero de fallecidos y causa. Desde las primeras epidemias, el ser humano siempre ha intentado entender la causa y el funcionamiento de ´estas. Ya en la antig¨uedad grandes cient´ıficos y fil´osofos escribieron sobre ellas: Plinio el viejo, Hip´ocrates, Arist´oteles, etc. El propio Arist´oteles sosten´ıa que las epidemias estaban relacionadas con eventos astron´omicos. Sin embargo, tendremos que avanzar hasta ´epocas m´as recientes para conseguir descifrar los mecanismos que hay detr´as. En el a˜no 1899, el microbi´ologo neerland´es Martinus Willem Beijem lleg´o a la conclusi´on de que exist´ıan unos entes m´as peque˜nos que las bacterias, los virus [2]. Aunque hasta el a˜no 1931, con la invenci´on del microscopio electr´onico, 9 Cap´ıtulo 1. Introducci´on a la modelizaci´on de epidemias Figura 1.3: Diagrama del modelo (SIRS)∗ En este caso al a˜nadir el efecto de la vacunaci´on las ecuaciones que se obtienen son las siguientes: (SIRS)∗            S′=−βI NS+µN −µS −ρS I′=βI NS−γI −µI R′=γI −µR +ρS (1.3) Tras unos sencillos c´alculos se obtienen los siguientes puntos de equilibrio: Equilibrio libre de dolencia E0= (N, 0,0). Si Rµ>1 equilibrio end´emico: E1,µ∗=N/Rµ,µ(Rµ−1) −ρN β,γ(Rµ−1) + ρN β Al igual que antes, el coeficiente de reproductividad es Rµ=β γ+µ. 1.3.3. SEIR, SEIRS y vacunas En los modelos precedentes, hemos visto c´omo aumenta la complejidad al a˜nadir nuevos par´ametros que relacionan los distintos compartimentos. En los siguientes modelos, (SEIR y variantes), aumentamos en una dimensi´on la complejidad del modelo a˜nadiendo un nuevo compartimento, constituidos por los individuos expuestos; ya que, durante los 16 Cap´ıtulo 1. Introducci´on a la modelizaci´on de epidemias procesos epid´emicos es usual que un individuo contaminado sufra un periodo de latencia durante el cual no contamina a ning´un otro. La poblaci´on de individuos expuestos se relaciona con el resto de poblaciones a trav´es del par´ametro σ, que se define como el inverso del tiempo medio durante el cual un individuo se mantiene expuesto (periodo de latencia). En la Figura 1.4 se resume la din´amica del modelo SEIR, donde ahora N=S+E+I+R. Figura 1.4: Diagrama del modelo (SEIR) Las EDOs que describen este modelo son: (SEIR)                  S′=−βI NS E′=βI NS−σE I′=σE −γI R′=γI (1.4) Al igual que en el modelo SIR, encontramos que R0=β γ. Los puntos de equilibrio son los siguientes: El equilibrio libre de dolencia E0= (N, 0,0,0). Si β > γ, el equilibrio end´emico E1,σ = (N/R0,0,0, N(1 −1/R0)). Tal y como vimos en la secci´on anterior, teniendo en cuenta la posibilidad de p´erdida de inmunidad, podemos pasar del modelo SEIR al modelos SEIRS; v´ease la Figura 1.5. Figura 1.5: Diagrama del modelo SEIRS 17 Cap´ıtulo 1. Introducci´on a la modelizaci´on de epidemias (SEIRS)                  S′=−βI NS+µN −µS E′=βI NS−σE −µE I′=σE −γI −µI R′=γI −µR (1.5) En este modelo, los puntos de equilibrio no tienen expresiones tan sencillas como en los modelos anteriores. Son los siguientes: El equilibrio libre de dolencia E0= (N, 0,0,0). El equilibrio end´emico: E1,µ,σ = (N/Rµσ,µ γσ(γ+µ)ΛN, µ γΛN, ΛN) donde Rµ,σ =β µ σ(γ+µ+σ) + γ Λ = 1−1 Rµσ µ γσ (γ+µ) + µ γ+ 1 Por analog´ıa con la secci´on anterior, podemos a˜nadir tambi´en la posibilidad de incluir efectos de vacunaci´on. As´ı, volver´ıamos a incluir el par´ametro ρ, como queda descrito en la Figura 1.6. Figura 1.6: Diagrama del modelo SEIRS∗ 18 Cap´ıtulo 1. Introducci´on a la modelizaci´on de epidemias (SEIRS∗)                  S′=−βI NS+µN −µS −ρS E′=βI NS−σE −µE I′=σE −γI −µI R′=γI −µR +ρS (1.6) 1.3.4. SEIRQ y m´as all´a Como hemos podido ver en las secciones previas, partiendo del modelo SIR, podemos ir aumentando gradualmente la complejidad de forma que permita simular con mayor precisi´on la realidad. Sin embargo, este incremento de precisi´on tiene un coste y este es el aumento de par´ametros y poblaciones (compartimentos), lo cual resta aplicabilidad (requiere de un mayor n´umero de datos que permita caracterizar correctamente los par´ametros del modelo). Presentamos a continuaci´on un ´ultimo modelo de epidemias, conocido como SEIRQ, donde se a˜naden los efectos producidos por el confinamiento. Bas´andonos en el modelo SEIR, a˜nadimos una nueva variable Q=Q(t) que representa el n´umero de individuos que se encuentran en cuarentena, aislados del resto de la poblaci´on. Tambi´en a˜nadimos dos nuevos par´ametros, pyλ, que se interpretan como la proporci´on de susceptibles que entran en cuarentena y el inverso del tiempo medio que dura la cuarentena de un individuo, respectivamente. La Figura 1.7 detalla c´omo se relacionan los compartimentos en este caso. Figura 1.7: Diagrama del modelo SEIRQ Para este modelo tendremos ahora que N=S+E+I+R+Q. Y las ecuaciones que 19 Cap´ıtulo 1. Introducci´on a la modelizaci´on de epidemias rigen este modelo son las siguientes: (SEIRQ)                        S′=−βI NS−pS +λQ E′=βI NS−σE I′=σE −γI R′=γI Q′=pS −λQ (1.7) Los puntos de equilibrio son algo m´as f´aciles de calcular, con R0=β γ: El equilibrio libre de dolencia E0= (0,0,0, N, 0). Si R0>1, el equilibrio end´emico: E1,p,λ =N/R0,0,0,1−1 R0 (1 + p λ)N, p λR0 N Para concluir esta secci´on, haremos algunos comentarios. En primer lugar, cabe destacar que el modelo SIR es tan solo un punto de partida o, como se conoce en la comunidad cient´ıfica anglosajona, un toy model, es decir, un modelo inicial que progresivamente se va mejorando. Como hemos podido ver a lo largo de estas secciones, existen muchas posibles variantes de este modelo que tienen en cuenta unos u otros efectos y sus implicaciones; sin embargo, a´un existen muchos otros efectos que podr´ıamos tener en cuenta, como la estructura por edades, las posibles patolog´ıas previas, la presencia de supercontagiadores, la efectividad de las vacunas y otros. Con esto quiero decir que, al menos en este trabajo, en vez de concentrarnos en conseguir el mejor modelo posible que tenga en cuenta la mayor cantidad de efectos, nos centraremos en modelos simples como los que hemos analizado hasta ahora, a˜nadiendo peque˜nas modificaciones e intentando afinar los par´ametros, para obtener as´ı los resultados m´as cercanos a la realidad. 20 Cap´ıtulo 2 Estudio num´erico y simulaciones En este cap´ıtulo, mediante el uso del software matem´atico MATLAB [6], procedemos a simular num´ericamente distintos modelos de tipo SIR. Consideraremos los modelos estudiados en las secciones del Cap´ıtulo 1 y algunos m´as. Antes de proceder con los resultados, vamos a introducir brevemente c´omo se utiliza MATLAB. Para resolver num´ericamente los sistemas de ecuaciones diferenciales, propios de cada modelo, se utiliza la funci´on ode45(odefun, [t0,tf], y0), cuyos par´ametros se definen como sigue: odefun: es un manejador de la funci´on que eval´ua el segundo miembro de la ecuaci´on, f(t, y). Puede ser el nombre de una funci´on an´onima dependiente de dos variables, siendo tla primera de ellas e yla segunda, o tambi´en una M-funci´on, en cuyo caso se escribe @odefun. [t0,tf]:es el intervalo en el que se quiere resolver la ecuaci´on. y0: es el valor del dato inicial. En nuestro caso, al ser un sistema de EDOs, los par´ametros en vez de ser escalares son vectores. Para informaci´on m´as detallada sobre la funci´on ode45, puede consultarse [7]. A continuaci´on, se presentan los resultados de las experiencias num´ericas realizadas. 2.1. Resultados para el modelo SIR Procedemos a simular el modelo SIR, resolviendo num´ericamente las Ecuaciones (1.1), con los siguientes valores para los par´ametros: la tasa de infecci´on, β= 0,125; tasa de 21 Cap´ıtulo 2. Estudio num´erico y simulaciones recuperaci´on, γ= 0,05. Por otra parte, los datos iniciales ser´an: S(0) = 1000, I(0) = 5 y R(0) = 0. De tal modo que N(0) = S(0) + I(0) + R(0) = 1005. A menos que se indique lo contrario, los valores de los par´ametros que vayamos definiendo y los datos iniciales se mantendr´an fijos en todas las simulaciones. Los resultados obtenidos para el modelo SIR se muestran en la Figura 2.1. Figura 2.1: Evoluci´on de la epidemia seg´un el modelo SIR. Hemos tomado un tiempo de simulaci´on igual a 360 d´ıas. Tal y como podemos apreciar, si analizamos las curvas, tras un primer momento en el que hay un gran aumento en el n´umero de infecciosos, se llega al m´aximo de contagios y, a partir de ese momento, las curvas tienden a valores constantes. Finalmente, la epidemia cesa, ya que, el n´umero de infecciosos cae a cero. Cabe destacar que, en esta simulaci´on, apenas un 10 % de la poblaci´on susceptible no se contagi´o de la enfermedad. En las siguientes simulaciones veremos que los resultados pueden variar. 22 Cap´ıtulo 2. Estudio num´erico y simulaciones 2.1.1. Confinamiento en el modelo SIR A continuaci´on, vamos a mostrar c´omo podemos ver el efecto de los confinamientos en el modelo SIR, sin necesidad de irnos a un modelo m´as complejo como el SEIRQ. Ser´a suficiente permitir a βque cambie con el tiempo. Para ello, supondremos que, en el inicio de la epidemia, tendremos βconstante e igual aβ1= 0,125. Tras un cierto tiempo, en el cual las autoridades se han percatado de la epidemia y sus riesgos para la poblaci´on, se procede a confinar a la poblaci´on. Esto, da lugar a una disminuci´on de la tasa de infecci´on β, que llegar´a a un valor m´ınimo β2= 0,065. Con el tiempo, la situaci´on mejora y las restricciones se van levantando. Nuevamente, el valor de βaumenta, volviendo a llegar a su valor original β1. Para m´as informaci´on sobre estos valores y una justificaci´on m´as completa, se puede consultar [8]. La Ecuaci´on (2.1) y la Figura 2.2 muestran el valor de β(t): Figura 2.2: Variaci´on del par´ametro βen el tiempo. 23 Cap´ıtulo 2. Estudio num´erico y simulaciones β(t) =                          β1si t≤ti β1+β2−β1 τβ (t−ti) si t∈[ti, ti+τβ] β2si t∈[ti+τβ, tf] β2+β1−β2 τβ (t−tf) si t∈[tf, tf+τβ] β1si t≥tf+τβ (2.1) donde, ti= 50, τβ= 50 y tf= 200. Hemos supuesto un comportamiento lineal en las transiciones entre β1yβ2. Con estos valores de β, los resultados obtenidos en la simulaci´on aparecen en la Figura 2.3. Figura 2.3: Evoluci´on de la epidemia seg´un el modelo SIR con β=β(t), dada por (2.1). L´ıneas continuas: βvariable; l´ıneas de puntos: β= 0,125. Las curvas con puntos se refieren a los resultados obtenidos con βfijo, v´ease Figura 2.1. De manera que podamos apreciar las diferencias. Como se puede observar, el m´aximo de contagios se adelanta unos d´ıas, el n´umero m´aximo de infecciosos es menor y el n´umero 24 Cap´ıtulo 2. Estudio num´erico y simulaciones final de susceptibles es mayor. Lo que indica que un menor porcentaje de la poblaci´on sufrir´a la enfermedad. 2.2. Resultados para el modelo SIRS En este caso, vamos a resolver las EDOs que rigen el modelo SIRS, v´ease el sistema (1.2). Este modelo es m´as interesante que el anterior, ya que tiene en cuenta la posibilidad de que los individuos pierdan la inmunidad. Para las simulaciones, fijemos µ= 0,0025. Recordemos que 1/µ se interpreta como el tiempo medio en el que desaparece la inmunidad en un individuo infeccioso o recuperado. Los resultados de la simulaci´on se muestran en las Figuras 2.4 y2.5. Figura 2.4: Evoluci´on de la epidemia seg´un el modelo SIRS con tmax = 360 Como vemos, al final de la simulaci´on, el n´umero de susceptibles empieza a crecer, mientras que los recuperados disminuyen. Veamos qu´e sucede si continuamos con la simulaci´on m´as all´a del tiempo final inicialmente fijado. 25 Cap´ıtulo 2. Estudio num´erico y simulaciones 2.3. Resultados para el modelo SEIR En esta secci´on, vamos a mostrar los resultados obtenidos al resolver num´ericamente el modelo SEIR, v´ease (1.4). En este modelo, se a˜nade un nuevo compartimento, E, que hace referencia a la poblaci´on expuesta. Para las siguientes simulaciones, fijamos σ= 0,1. Recordemos que 1/σ se interpreta como el tiempo medio durante el cual un individuo se mantiene expuesto. La poblaci´on total (constante) es ahora: N=S(t) + E(t) + I(t) + R(t) para todo tyE(0) = 0. Los resultados obtenidos para el modelo SEIR se muestran en la Figura 2.11. Figura 2.11: Evoluci´on de la epidemia seg´un el modelo SEIR Comparando con el modelo SIR, v´ease Figura 2.1, vemos c´omo las poblaciones de susceptibles y recuperados tienden a los mismos valores. Esto era de esperar, ya que los sistemas que rigen estos modelos, v´eanse (1.1)y(1.4), describen igual comportamiento para estas poblaciones. Sin embargo, se observa que el m´aximo de contagios es menor en el modelo SEIR y 32 Cap´ıtulo 2. Estudio num´erico y simulaciones tiene lugar m´as tarde. A diferencia del modelo SIRS, no tienen lugar m´as brotes, ya que no existe ning´un mecanismo para que los recuperados vuelvan a ser susceptibles. De forma que la epidemia cesa pasado un cierto tiempo. 2.3.1. Confinamiento en el modelo SEIR Al igual que en secciones anteriores, podemos preguntarnos qu´e ocurre cuando βvar´ıa con el tiempo. Como el m´aximo de contagios se da ahora m´as tarde, conviene modificar los valores de β(t). La funci´on βse muestra en (2.6) y en la Figura 2.12. Figura 2.12: Variaci´on del par´ametro βen el tiempo. M´as precisamente, tomaremos β(t) =                          β1si t≤ti β1+β2−β1 τβ (t−ti) si t∈[ti, ti+τβ] β2si t∈[ti+τβ, tf] β2+β1−β2 τβ (t−tf) si t∈[tf, tf+τβ] β1si t≥tf+τβ (2.6) donde ti= 100, τβ= 50 y tf= 250. 33 Cap´ıtulo 2. Estudio num´erico y simulaciones Con estos valores de β, los resultados obtenidos aparecen en la Figura 2.13. Figura 2.13: Evoluci´on de la epidemia seg´un el modelo SEIR con β=β(t), dada por (2.6). L´ıneas continuas: βvariable; l´ıneas discontinuas: β= 0,125. Al igual que ocurr´ıa en el modelo SIR, el m´aximo de contagios disminuye y tiene lugar antes, v´ease la Figura 2.3 para comparar. De nuevo se observa que, un menor porcentaje de la poblaci´on sufre la enfermedad. 2.3.2. Sensaci´on p´ublica de riesgo en el modelo SEIR A continuaci´on, veamos c´omo influye la sensaci´on p´ublica de riesgo en el modelo SEIR. La funci´on Ψβ(t, D) est´a dada por (2.4), donde hemos asumido las mismas aproximaciones. 34 Cap´ıtulo 2. Estudio num´erico y simulaciones El sistema de ecuaciones diferenciales que caracteriza este modelo es el siguiente:                  S′=−Ψβ I NS E′= Ψβ I NS−σE I′=σE −γI R′=γI (2.7) En las simulaciones, c= 0,01 y κ={0,100,200}. Los resultados obtenidos se muestran en la Figura 2.14. Figura 2.14: Evoluci´on de la epidemia seg´un el modelo SEIR con percepci´on del riesgo por parte de la poblaci´on. L´ıneas continuas: κ= 0; l´ıneas discontinuas: κ= 100; l´ıneas de puntos: κ= 200. Como cab´ıa esperar, un comportamiento m´as responsable por parte de la poblaci´on (con mayor κ) se ve reflejado en la evoluci´on de la epidemia. Se hace evidente una clara reducci´on en el n´umero de infecciosos y expuestos. 35 Cap´ıtulo 2. Estudio num´erico y simulaciones 2.4. Resultados para el modelo SEIRS Al igual que hicimos en la Secci´on 2.2, vamos a analizar c´omo se modifica el modelo SEIR, al tener en cuenta la posibilidad de que los individuos pierdan la inmunidad. EL correspondiente modelo ser´a l´ogicamente denominado SEIRS. Resolviendo las Ecuaciones (1.5) con µ= 0,0025, se obtienen los resultados mostrados en la Figura 2.15. Figura 2.15: Evoluci´on de la epidemia seg´un el modelo SEIRS con tmax = 360 Veamos c´omo se comportan las poblaciones cuando simulamos para tiempos mayores. Los resultados se muestran en la Figura 2.16. 36 Cap´ıtulo 2. Estudio num´erico y simulaciones Figura 2.16: Evoluci´on de la epidemia seg´un el modelo SEIRS con tmax = 1200 Como vemos, en los infecciosos se produce un segundo m´aximo alrededor de t= 800 d´ıas. A diferencia del modelo SEIR, la epidemia no cesa y se convierte en una endemia, ya que, en este caso, el coeficiente de reproductividad es Rµ,σ =β µ σ(γ+µ+σ) + γ>1. Comparando con la Figura 2.5, podemos ver que los m´aximos alcanzados por I=I(t) son algo menores y tienen lugar antes. Cabe destacar, claro est´a, que tenemos una nueva poblaci´on (los expuestos) que constituyen el compartimento previo al de los infecciosos. Esto explica la disminuci´on que observamos. 2.4.1. Confinamiento en el modelo SEIRS Vamos a permitir nuevamente que βpueda variar en el tiempo. Para ello, nos fijamos en la Figura 2.16 y analizamos los tiempos en los que ocurren los mayores n´umeros de 37 Cap´ıtulo 2. Estudio num´erico y simulaciones contagiados. Proponemos los siguientes valores de β(t) en la Ecuaci´on (2.8): β(t) =                                                    β1si t≤ti β1+β2−β1 τβ (t−ti) si t∈[ti, ti+τβ] β2si t∈[ti+τβ, tf] β2+β1−β2 τβ (t−tf) si t∈[tf, tf+τβ] β1si t∈[tf+τβ, ti2] β1+β2−β1 τβ (t−ti2) si t∈[ti2, ti2+τβ] β2si t∈[ti2+τβ, tf2] β2+β1−β2 τβ (t−tf2) si t∈[tf2, tf2+τβ] β1si t≥tf2+τβ (2.8) donde ti= 100, tf= 300, ti2= 650 y tf2= 850; v´ease la Figura 2.17. Figura 2.17: Variaci´on del par´ametro β(t) con el tiempo. Los resultados num´ericos obtenidos se muestran en la Figura 2.18. 38 Cap´ıtulo 2. Estudio num´erico y simulaciones Figura 2.18: Evoluci´on de la epidemia seg´un el modelo SEIRS con β(t) y tmax = 1200. L´ıneas continuas: βvariable; l´ıneas de puntos: β= 0,125. Analizando la Figura 2.18, se aprecia que los m´aximos se adelantan. El primero, es de menor valor que el obtenido con βfijo; sin embargo, tanto el segundo como el tercero son m´as acentuados. Por lo tanto, podemos concluir que, pese a que los confinamientos parec´ıan estar justificados y tener un orden temporal l´ogico, los resultados no son totalmente satisfactorios. Podr´ıa parecer que el uso del confinamiento como medida para combatir epidemias no es ´util, no obstante, hay muchos factores que no estamos teniendo en cuenta: el tiempo durante el cual se impone, el tiempo de adaptaci´on de la poblaci´on, lo estrictas que sean las medidas, etc. Todos estos aspectos modifican el valor de βy por tanto caracterizan la evoluci´on de la epidemia. Vemos por lo tanto que medidas como el confinamiento deben ser estudiadas cuidadosamente antes de ser impuestas, pues pese a parecer la opci´on m´as l´ogica, no siempre proporcionan los mejores resultados. 39 Cap´ıtulo 2. Estudio num´erico y simulaciones 2.4.2. Sensaci´on p´ublica de riesgo en el modelo SEIRS De nuevo, veamos como influye la sensaci´on p´ublica de riesgo en el modelo SEIRS. La funci´on Ψβ(t, D) que usaremos, sera la dada por (2.4), de nuevo con las mismas aproximaciones. El sistema de ecuaciones diferenciales que caracteriza este modelo es el siguiente:                  S′=−Ψβ I NS+µN −µS E′= Ψβ I NS−σE −µE I′=σE −γI −µI R′=γI −µR (2.9) En las simulaciones, c= 0,01 y κ={0,100,200}. Los resultados obtenidos se muestran en la Figura 2.19. Figura 2.19: Evoluci´on de la epidemia seg´un el modelo SEIRS con percepci´on del riesgo por parte de la poblaci´on. L´ıneas continuas: κ= 0; l´ıneas discontinuas: κ= 100; l´ıneas de puntos: κ= 200. 40 Cap´ıtulo 2. Estudio num´erico y simulaciones 2.4.3. Vacunas: modelo SEIRS* Veamos el efecto que tiene la vacunaci´on en el modelo SEIRS*. Recordamos que ρ es, por definici´on, el n´umero de individuos vacunados por unidad de tiempo. En nuestras simulaciones, ρ= 0,02. Resolviendo num´ericamente el sistema (1.6), obtenemos los resultados que se muestran en la Figura 2.20. Figura 2.20: Comparaci´on de la evoluci´on de la epidemia seg´un el modelo SEIRS con y sin vacunaci´on. L´ıneas continuas, con vacunaci´on; l´ıneas discontinuas, sin vacunaci´on. Como vemos, al aplicar vacunas, el n´umero de infecciosos y de expuestos disminuye radicalmente. Se puede observar la evoluci´on de estas dos poblaciones en la Figura 2.21. 41 Cap´ıtulo 2. Estudio num´erico y simulaciones Analizando los resultados, se observa que una mayor preocupaci´on de la poblaci´on a la hora de hacer cuarentena, junto con una vacuna eficaz, puede cambiar totalmente el curso de una epidemia. Adem´as, se puede ver claramente c´omo el n´umero de infecciosos y expuestos cae r´apidamente con estas medidas. Incluso se consigue evitar un segundo repunte en los casos, que s´ı se observa en las l´ıneas continuas. Se pueden comparar estas gr´aficas con las Figuras 2.19 y2.20 donde, a partir del modelo SEIRS, a˜nad´ıamos el efecto de la percepci´on del riesgo y el uso de vacunas. Para finalizar este cap´ıtulo, tal y c´omo dec´ıamos al final del primero, se˜nalemos que los modelos se pueden complicar considerablemente, a˜nadiendo uno u otros efectos. En este caso, hemos partido del modelo SEIRQ y hemos a˜nadido par´ametros que relacionan los diferentes compartimentos entre s´ı. Con ello conseguimos poder analizar los efectos de distintas medidas y las actitudes de la poblaci´on mediante la variaci´on de dichos par´ametros. 48 Cap´ıtulo 3 An´alisis y resoluci´on de problemas inversos En este cap´ıtulo vamos a considerar algunos problemas inversos. Ser´an resueltos num´ericamente haciendo uso del software MATLAB [6]. Existen muchos problemas inversos posibles en el contexto de los sistemas descritos en los cap´ıtulos precedentes.Los que vamos a resolver aqu´ı son de la siguiente forma: dados ciertos valores de tiempo (t) y conocidos los correspondientes valores de I, determinar las constantes caracter´ısticas del sistema. Veamos un ejemplo aplicado al modelo SIR: conocidos los n´umeros de infecciosos a lo largo del tiempo y las condiciones iniciales (S(0), I(0), R(0)), el problema inverso consiste en obtener los valores de los par´ametros βyγque mejor se ajustan a los datos. Este ajuste calcular´a usando una t´ecnica de m´ınimos cuadrados no lineales. El an´alisis y resoluci´on de estos problemas inversos es de gran utilidad, ya que permite caracterizar la infecci´on en los estadios m´as tempranos de una epidemia. Proporciona informaci´on sobre la transmisibilidad del pat´ogeno, la tasa de infecci´on, la tasa de recuperaci´on, etc. Todo esto es de vital importancia para estudiar en mayor profundidad una epidemia y as´ı comprenderla mejor y actuar en consecuencia de forma que se minimicen los da˜nos. 3.1. Procedimiento En esta secci´on vamos a hacer un breve repaso al c´odigo y funciones MATLAB que vamos a usar para resolver los problemas inversos. Vamos a empezar explicando la funci´on 49 Cap´ıtulo 3. An´alisis y resoluci´on de problemas inversos fmincon. Esta funci´on se usa para resolver problemas de minimizaci´on con restricciones. Nuestro problema inverso es de este tipo: ajuste de m´ınimos cuadrados no lineales. En general, esta funci´on se ejecuta de la siguiente manera: fmincon(funJ, x0, C, d, A, b, ci, cs, resnolin, options). Los argumentos se definen como sigue: funJ: es un manejador de la funci´on a minimizar J(x). Como en otros casos, puede ser el nombre de una funci´on an´onima o un manejador de una M-funci´on. x0: es una aproximaci´on inicial de la soluci´on. C, d, A, b: se usan para definir restricciones lineales de igualdad y desigualdad. ci, cs: son, respectivamente, las cotas inferior y superior ci ≤x≤cs. resnolin: es un manejador de la funci´on que define las restricciones no lineales, tanto de igualdad como de desigualdad. options: es una estructura que permite modificar algunos de los par´ametros internos de la funci´on fmincon. Para m´as informaci´on sobre problemas de minimizaci´on y la funci´on fmincon, puede consultarse [7]. Para la resoluci´on de problemas inversos, solo vamos a utilizar algunos de estos par´ametros. Definimos la funci´on funJ, que eval´ua la suma de los cuadrados de las distancias de los datos observados a los valores correspondientes a la elecci´on de los par´ametros. Estos valores se consiguen a partir de las soluciones al modelo SIR (o aqu´el con el que trabajemos) fijados unos valores de los par´ametros βyγ. Proponemos unos valores iniciales x0 = [β0, γ0] con los que empezar a trabajar y les imponemos unas cotas ci ycs. Y, por ´ultimo, utilizamos el argumento options, para definir el algoritmo que resuelve nuestro problema (“active-set”). Con el fin de mostrar c´omo se resuelve un problema inverso, vamos a trabajar de la siguiente manera. Tal y como hicimos en el Cap´ıtulo 2, vamos a resolver num´ericamente distintos modelos de epidemias, conocidos los par´ametros y las condiciones iniciales. Despu´es, vamos a quedarnos solo con los datos referentes a las condiciones iniciales y a los n´umeros de infecciosos a lo largo del tiempo, cada 30 d´ıas. Con estos datos, propondremos unos valores iniciales de los par´ametros y, con la herramienta MATLAB, resolveremos el 50 Cap´ıtulo 3. An´alisis y resoluci´on de problemas inversos problema inverso. Finalmente, comprobaremos que ha sido resulto el problema de manera correcta comparando los par´ametros originales con los obtenidos mediante minimizaci´on. 3.2. Problema inverso: modelo SIR Tal y como hicimos en la Secci´on 2.1, resolvemos num´ericamente el modelo SIR, v´ease (1.1), con las siguientes condiciones iniciales: S(0) = 1000, I(0) = 5, R(0) = 0. Y damos los siguientes valores a los par´ametros: la tasa de infecci´on β= 0,125; tasa de recuperaci´on γ= 0,05. Los resultados se muestran en la Figura 3.1. Figura 3.1: Resoluci´on num´erica modelo SIR Para resolver el problema inverso, los datos que tenemos sobre el n´umero de infecciosos aparecen en la Tabla 3.1. t (d´ıas) 30 60 90 120 150 180 210 240 270 300 330 360 I43 193 208 97 36 12 4 1 0 0 0 0 Tabla 3.1: N´umero de infecciosos cada 30 d´ıas en el modelo SIR. 51 Cap´ıtulo 3. An´alisis y resoluci´on de problemas inversos Las condiciones iniciales elegidas son: S(0) = 1000, I(0) = 5, R(0) = 0. Proponemos como aproximaci´on inicial β0= 0,03 y γ0= 0,1, con cotas inferior y superior ci = 0 y cs = 2. En la Figura 3.2 se muestra el resultado obtenido al resolver el problema inverso. Figura 3.2: Resoluci´on num´erica problema inverso modelo SIR En la Tabla 3.2 se muestran los valores obtenidos para los par´ametros βyγ. Propuesto Final Original Error Relativo ( %) β0,03 0,124999 0,125 0,0008 γ0,1 0,050001 0,05 0,002 Residuo 0,000096 Tabla 3.2: Resultados para el problema inverso correspondiente al modelo SIR Como vemos, tanto si comparamos las Figuras 3.1 y3.2, como los valores de βyγobtenidos, los resultados son muy buenos. Encontramos pr´acticamente una total coincidencia 52 Cap´ıtulo 3. An´alisis y resoluci´on de problemas inversos entre los valores observados y calculados de los par´ametros. El error relativo obtenido es m´ınimo. El residuo que se muestra en la parte inferior de la Tabla 3.2, es el valor de la suma de los cuadrados de las distancias en cada uno de los instantes de tiempo elegidos y podemos ver tiene un valor muy bajo. Si bien estamos trabajando con un modelo sencillo, este resultado nos muestra c´omo con pocos datos podemos obtener una gran cantidad de informaci´on. 3.3. Problema inverso: modelo SIRS Vamos a analizar otro problema inverso. En este caso resolveremos num´ericamente el modelo SIRS, v´ease (1.2). Dadas las mismas condiciones iniciales y los mismos valores para los par´ametros βyγ, a˜nadimos el par´ametro µ= 0,0025, que refleja la p´erdida de inmunidad. Los resultados se muestran en la Figura 3.3. Figura 3.3: Resoluci´on num´erica modelo SIRS Para resolver el problema inverso, los datos que tenemos referentes al n´umero de 53 Cap´ıtulo 3. An´alisis y resoluci´on de problemas inversos infecciosos son los que se muestran en la Tabla 3.3 t (d´ıas) 30 60 90 120 150 180 210 240 270 300 330 360 I40 179 213 114 52 25 14 9 6 5 5 6 Tabla 3.3: N´umero de infecciosos cada 30 d´ıas en el modelo SIRS. Sabiendo que las condiciones iniciales son: S(0) = 1000, I(0) = 5, R(0) = 0, proponemos como aproximaci´on inicial β0= 0,03, γ0= 0,1 y µ0= 0,01 con cotas ci = 0 y cs = 2. En la Figura 3.4 se muestra el resultado obtenido al resolver el problema inverso. Figura 3.4: Resoluci´on num´erica problema inverso modelo SIRS En la Tabla 3.4 se muestran los valores obtenidos para los par´ametros β,γyµ. 54 Cap´ıtulo 3. An´alisis y resoluci´on de problemas inversos Propuesto Final Original Error Relativo ( %) β0,03 0,124999 0,125 0,0008 γ0,1 0,049999 0,05 0,002 µ0,01 0,002499 0,0025 0,04 Residuo 0,000008 Tabla 3.4: Resultados problema inverso modelo SIRS Al igual que en la secci´on previa con el modelo SIR, el resultado obtenido es muy satisfactorio. Con estos ejemplos de resoluci´on de problemas inversos hemos mostrado un m´etodo con el que resolver problemas de este tipo. Aplicando esta estrategia, podemos obtener informaci´on de sistemas complejos en los cuales la cantidad de datos disponibles es limitado. Este es el caso de las epidemias. Durante una epidemia, no siempre es f´acil obtener datos, sobre todo en los primeros momentos, cuando a´un no se conocen bien las caracter´ısticas de la enfermedad que se est´a propagando. No obstante, hemos visto que con pocos datos podemos llegar a obtener bastante informaci´on, algo que m´as tarde puede ser usado para generar mejores predicciones y entender con mayor profundidad el origen, la evoluci´on y las consecuencias. Por todo esto, es importante coleccionar la mayor cantidad de datos con la mejor calidad posible. 55 Bibliograf´ıa [1] Principales pandemias en la historia de la humanidad. Carlos Casta˜neda Gullot y Gerardo Ramos Serpa 21/06/2020. Cuadro 1. http://www.revpediatria.s ld.cu/index.php/ped/article/view/1183 [2] Las matem´ aticas de la pandemia. Manuel de Le´on y Antonio G´omez Corral. CSIC-Catarata, 2020. [3] Cuarentena: origen del concepto, qu´ e significa y cu´ al es su implicancia como medida sanitaria. ICIM facultad de medicina del desarrollo 13/04/2020. ht tps://medicina.udd.cl/icim/2020/04/13/cuarentena-origen-del-concepto-q ue-significa-y-cual-es-su-implicancia-como-medida-sanitaria/ [4] The smallpox vaccine: the dispute between Bernoulli and d’Alembert and the calculus of probabilities. Colombo, C., Diamanti, M. Lett Mat Int 2, 185–192 (2015). https://doi.org/10.1007/s40329-015-0073-5 [5] Matem´ aticas y epidemias. Nicolas Baca¨er ©Cassini, Paris, 2021. [6] MATLAB. Software matem´atico. https://es.mathworks.com/products/matlab .html [7] Apuntes de MATLAB orientados a m´ etodos num´ ericos elementales. Rosa Echevarr´ıa. Universidad de Sevilla, 2022 https://personal.us.es/echevarria/do cumentos/APUNTESMATLAB.pdf [8] Modelando una epidemia I: Confinamientos. Renato Alvarez Nodarse y Fernando Carre˜no Navas. Blog del Instituto de Matem´aticas de la Universidad de Sevilla, 19 junio, 2021 https://institucional.us.es/blogimus/2021/06/modelando-una -epidemia-i-confinamientos/ 56 Bibliograf´ıa [9] Modelando una epidemia II: efecto del comportamiento de la poblaci´ on. Renato Alvarez Nodarse y Fernando Carre˜no Navas. Blog del Instituto de Matem´aticas de la Universidad de Sevilla, 29 junio, 2021 https://institucional.us.e s/blogimus/2021/07/modelando-una-epidemia-ii-efecto-del-comportamiento -de-la-poblacion/ 57