Dinámica de poblaciones: sobre la base del modelo de Lotka-Volterra

Facultad de Física D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA POPULATION DYNAMICS: ON THE BASIS OF THE MODEL OF LOTKA-VOLTERRA Trabajo Fin de Grado GRADO EN FÍSICA Alumno: Juan Escamilla Martín Tutora: Prof.ª Dra. M.ª Carmen Lemos Fernández Departamento de Física de la Materia Condensada Curso 2019-2020 A mi abuela Dionisia Cañeque Cogollor. Gracias por estar siempre a mi lado. Sin ti no habría sido posible haber llegado hasta aquí. A mi familia y amigos, por su cariño y apoyo. A mi buen amigo y compañero José Antonio Díaz Fraile, por su valiosa e impagable ayuda durante gran parte de mi etapa universitaria. Y mi más sincero agradecimiento a María del Carmen Lemos Fernández, por tutorizar este Trabajo Fin de Grado y por su total disponibilidad y trato cercano. El alumno 5 Índice de contenidos Resumen ...................................................................................................................... 7 Introducción ................................................................................................................. 8 1 Modelo de Malthus .............................................................................................. 10 2 Modelo de Verhulst ............................................................................................. 13 3 Modelo depredador-presa de Lotka y Volterra .................................................... 15 3.1 Solución analítica en torno al punto de silla ................................................. 18 3.2 Solución analítica en torno al punto vórtice ................................................. 19 3.3 Solución numérica general ............................................................................ 24 4 Modelo epidemiológico SIR de McKendrick y Kermack ...................................... 31 4.1 Ajuste del modelo SIR a los datos de la pandemia del COVID-19 en China ..... 34 4.1.1 Modelo SIR clásico: Susceptibles, Infectados y Eliminados .................... 34 4.1.2 Modelo SIR extendido: Susceptibles, Infectados, Fallecidos y Altas ....... 37 Conclusiones .............................................................................................................. 41 Referencias ................................................................................................................. 44 Apéndices ................................................................................................................... 46 Apéndice A: Demostraciones analíticas .................................................................. 46 A.1 Resolución de la ecuación (1.4) ..................................................................... 46 A.2 Demostración de la ecuación (2.2) ................................................................ 46 A.3 Demostración de la ecuación (3.7) ................................................................ 47 A.4 Resolución de la ecuación (3.17) ................................................................... 48 A.5 Resolución de la ecuación (3.18) ................................................................... 49 A.6 Demostración de las ecuaciones (3.19), (3.20), (3.21), (3.22) y (3.23) ........... 49 A.7 Demostración de la ecuación (3.26) .............................................................. 51 A.8 Demostración de la expresión (3.28) ............................................................. 51 A.9 Función W de Lambert en ℝ ........................................................................ 52 A.10 Integración de la ecuación (3.31) en un entorno próximo al punto vórtice . 53 6 Apéndice B: Rutinas en MATLAB ® R2017a ........................................................... 56 B.1 Figura 1.1 ..................................................................................................... 56 B.2 Figura 1.2 ..................................................................................................... 56 B.3 Figura 2.1 ..................................................................................................... 56 B.4 Figura 3.1 ..................................................................................................... 57 B.5 Figura 3.2 ..................................................................................................... 57 B.6 Figura 3.3 ..................................................................................................... 58 B.7 Figura 3.4 ..................................................................................................... 60 B.8 Figura 3.5 ..................................................................................................... 63 B.9 Figura 4.1 ..................................................................................................... 64 B.10 Figuras 4.2 y 4.3 ......................................................................................... 65 B.11 Figura I ....................................................................................................... 69 B.12 Figura II ...................................................................................................... 69 D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 7 Resumen El presente Trabajo Fin de Grado es una introducción al modelo depredadorpresa, también conocido como modelo de Lotka-Volterra, del cual se hace un estudio numérico y analítico alrededor de sus puntos de equilibrio. Como aplicación, se ajusta este modelo a los datos registrados en China sobre la actual pandemia del COVID-19 (modelo epidemiológico SIR). Todos los cálculos, simulaciones y figuras recogido s en esta memoria han sido realizados por el alumno con el software MATLAB ® (versión R2017a). Como se verá más adelante, dichas simulaciones explican cualitativamente bien el comportamiento periódico y desfasado en el tiempo de las poblaciones presa y depredadora, debido a su mutua interacción, así como el comportamiento de las poblaciones de individuos infectados, fallecidos y dados de alta en China a causa de la pandemia del COVID-19. Abstract The present dissertation is an introduction to the model predatorprey, also known as the Lotka-Volterra model, of which a numeric and analytic study is ma de around its balance points. This model is adjusted to the data registered in China of the current COVID-19 pandemic (epidemiologic model SIR) . All the calculations, simulations and figures collected in this work have been processed by the student with th e software MATLAB ® (version R2017a). As it will be shown, these simulations explain well qualitatively the periodic and the out of phase behaviour of those prey and predator populations, due to their mutual interaction, as well as the behaviour of the population of the infected, deceased and discharged in China as a result of the COVID-19 pandemic. D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 8 Introducción La dinámica de poblaciones es un área de la matemática que estudia cómo evoluciona el número de constituyentes de un conjunto a lo largo del tiempo; por ejemplo, los individuos de una especie biológica en unas determinadas condiciones. Además de tener grandes aplicaciones en el ámbito de la biología, cabe señalar que esta disciplina es muy utilizada también en otras materias como la economía y la demografía, entre otras. Para el estudio de muchos fenómenos, la dinámica de poblaciones se vale de modelos matemáticos diferenciales, que son aquellos que se describen mediante ecuaciones diferenciales (1) . En la literatura podemos encontrar muchos de estos modelos, pero en esta memoria nos vamos a centrar en hacer una introducción al modelo depredador-presa (o también llamado modelo de Lotka-Volterra) formulado por los matemáticos A. J. Lotka y V. Volterra durante la primera mitad de la década de 1920. A LFRED J AMES L OTKA V ITO V OLTERRA (1880-1949) (1860-1940) (1) Tarasiévich, 2004. (*) Las imágenes que aparecen en esta página han sido extraídas de la primera referencia bibliográfica (Bacaër, 2011). D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 9 Este modelo es muy utilizado en el campo de la ecología para estudiar la interacción entre dos especies biológicas cuando una se alimenta de la otra, presentando ambas una evolución periódica en el tiempo. El primero en proponer las ecuaciones que describen el modelo depredador-presa fue Lotka en 1920, aunque en un contexto bastante diferente. Él obtuvo estas ecuaciones en sus estudios sobre reacciones químicas en los que observó que algunas de ellas presentaban un comportamiento oscilatorio en el tiempo. Pero, no fue hasta 1925 cuando Lotka definió la interacción depredador-presa de dos especies que conviven (2) . Un año más tarde, en 1926, Volterra llegaría de forma independiente a la misma conclusión que Lotka tras analizar los datos sobre la pesca en el norte del Mar Adriático que su yerno, Umberto D’Ancona (1896-1964), había recogido durante los años de la Primera Guerra Mundial (3) . El modelo de Lotka-Volterra también se puede emplear para predecir y explicar el comportamiento cualitativo de epidemias como la que ocurrió entre 2003 y 2004 en Hong-Kong (China) a causa del coronavirus SARS-CoV (4) del cual proviene la mutación SARS-CoV-2 (más conocido como COVID-19), que ha causado la pandemia que actualmente estamos padeciendo y que analizaremos en la sección cuarta de esta memoria. Antes de abordar el modelo depredador-presa y de ahondarnos en su formalismo matemático, hablaremos brevemente de otros dos modelos claves en la historia de la dinámica de poblaciones: el modelo basado en la concepción social de T. R. Malthus, de 1798; y el modelo propuesto por P. F. Verhulst en 1838. Como veremos, el modelo de Lotka-Volterra se construye utilizando algunas de las premisas de estos dos modelos. (2) Lotka, 1925; Bacaër, 2011. (3) Volterra, 1926. (4) Small, 2012. D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 16   individuos decrecerán hasta su extinción, según predice el modelo de Malthus. Después, Volterra consideró a las dos especies homogéneamente mezcladas y, siguiendo un razonamiento análogo al de Verhulst, estableció que la variación temporal del número de individuos de ambas especies debía ser proporcional al número de contactos posibles que pueden darse entre los individuos de la especie presa 8 9 y los individuos de la especie depredadora 8  , el cual viene dado por el producto  9   . Esto le llevó a proponer el siguiente sistema de ecuaciones diferenciales: ⎩ ⎪ ⎨ ⎪ ⎧   9  =  9  9 − > 9  9   =  9   9 − > 9        = > 9  9   −     = −      − > 9  9  (3.1) que ponen de manifiesto la disminución en el número de individuos  9 de la especie presa 8 9 cuando el número de individuos   de la especie depredadora 8  aumenta. > 9 es un parámetro positivo que da cuenta de la tasa per cápita de decrecimiento de la especie presa 8 9 debido a su interacción con la especie depredadora 8  . Asimismo, > 9 también es un parámetro positivo que da cuenta de la tasa per cápita de crecimiento de la especie depredadora 8  debido a su interacción con la especie presa 8 9 . A priori, no hay razón para pensar que el decrecimiento de la especie 8 9 suponga un crecimiento instantáneo de la especie 8  y viceversa, por tanto, en general > 9 ≠> 9(14) . Un año antes, en 1925, el matemático austriaco Alfred James Lotka (1880-1949) también trató este problema y llegó a la misma conclusión que Volterra (15) . Aunque en 1920 Lotka ya había propuesto el sistema de ecuaciones (3.1) en sus estudios sobre reacciones químicas (16) , Volterra propuso este modelo depredador-presa sin tener conocimiento de los trabajos de Lotka; así lo atestiguan las cartas que ambos se (14) Volterra, 1926. (15) Lotka, 1925. (16) Lotka, 1920a, 1920b. D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 17 intercambiaron a finales de 1926 (17) . Es por ello que las ecuaciones (3.1) reciben comúnmente el nombre de ecuaciones de Lotka-Volterra. Estas ecuaciones están acopladas entre sí y no son lineales, por tanto, no podemos conocer las expresiones analíticas generales de las funciones  9  y   . No obstante, sí podemos obtener una aproximación de dichas funciones en las cercanías de los puntos de equilibrio (o puntos críticos de primera especie) del sistema de ecuaciones (3.1), lo cual proporciona información relevante acerca del comportamiento general de  9  y   . Los puntos de equilibrio son aquellos que cumplen que   9  @ A B - A B C = 0 ;     @ A E - A E C = 0 Así, los puntos de equilibrio del sistema de ecuaciones (3.1) son 0,0 y G 'E HEB , 'B HBE I. Si hacemos uso del desarrollo en serie de Taylor de una función de dos variables (=(),J en un entorno muy próximo al punto ),J=)*,J* ( = (  ) , J  | , - , * K - K * + L(  ) , J  L) . , - , * K - K *  ) − ) *  + L(  ) , J  LJ . , - , * K - K *  J − J *  + /  )  , J   en el sistema de ecuaciones (3.1) e ignoramos los términos de orden mayor que 1, dicho sistema toma la forma siguiente alrededor del punto de equilibrio  9 ,  = 9M , M : ⎩ ⎪ ⎨ ⎪ ⎧   9  =   9 − > 9   M    9 −  9 M  − > 9  9 M    −   M      = > 9   M   9 −  9 M  −    − > 9  9 M     −   M  que, reagrupando términos, se puede expresar en forma matricial: ⎝ ⎜ ⎛   9      ⎠ ⎟ ⎞ =  9 − > 9  M − > 9  9M > 9  M > 9  9M −    9   + 2 > 9  9M  M −  9  9M    M −2> 9  9M  M  (3.2) Se obtiene así un sistema de dos ecuaciones diferenciales ordinarias de primer orden, acopladas, lineales y no homogéneas, que se puede resolver analíticamente. (17) Israel y Millán Gasca, 2002. D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 18 3.1 Solución analítica en torno al punto de silla En un entorno muy próximo al punto de equilibrio 0,0, el sistema de ecuaciones (3.2) se reduce a dos ecuaciones diferenciales ordinarias de primer orden, lineales y homogéneas:   9  =  9  9 (3.3)     = −     (3.4) cuyas soluciones analíticas son (18)  9    =  9   ' B  (3.5)      =     T ' E  (3.6) donde los tamaños de la especie presa 8 9 y de la especie depredadora 8  en el instante inicial se denotan como  9 y   , respectivamente. De acuerdo con las ecuaciones (3.5) y (3.6), las trayectorias en el espacio de las fases son hipérbolas cuyas asíntotas son  9 =0 y   =0 (ver figura 3.2) y la ecuación que las describe es (19)   =    U  9  9  V T ' E ' B (3.7) Cuando las trayectorias en el espacio de las fases tienen un comportamiento hiperbólico en torno a un punto de equilibrio, se dice que éste es un punto de silla. Este tipo de punto es inestable, puesto que, las trayectorias divergen alrededor de él (20) . (18) La resolución de las ecuaciones (3.3) y (3.4) es análoga a la de la ecuación (1.4) (ver apéndice A.1). (19) La demostración de la ecuación (3.7) se encuentra en el apéndice A.3. (20) Samanta y Gómez Aíza, 2014. D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 19 Figura 3.2. Espacio de las fases en el que se representan las trayectorias que describe el sistema de ecuaciones (3.1) en un entorno muy próximo al punto de equilibrio 0,0. Las trayectorias que ocupan el primer cuadrante (en color azul) son las que tienen sentido físico, puesto que, el número de individuos de cada población es positivo. Se observa que el número de individuos  9 de la especie presa 8 9 crece mientras que el número de individuos   de la especie depredadora 8  disminuye. En esta representación se han tomado los siguientes valores hipotéticos de los parámetros del sistema de ecuaciones (3.1):  9 =0,15;   =0,1. 3.2 Solución analítica en torno al punto vórtice En un entorno muy próximo al punto de equilibrio G 'E HEB , 'B HBE I, el sistema de ecuaciones (3.2) toma la siguiente forma: ⎝ ⎜ ⎛   9     ⎠ ⎟ ⎞ = ⎝ ⎜ ⎛ 0 −   > 9 > 9  9 > 9 > 9 0 ⎠ ⎟ ⎞ W X X X X X Y X X X X X Z [   9   + ⎝ ⎜ ⎛  9   > 9 −  9   > 9 ⎠ ⎟ ⎞ (3.8) D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 20 Para resolver las ecuaciones (3.8) es necesario desacoplarlas transformando la matriz de coeficientes [ en una matriz diagonal. Cuando el producto de una matriz \ por un vector columna ]^ devuelve dicho vector multiplicado por un número _, se dice que _ es un valor propio o autovalor de \ y que ]^ es un vector propio o autovector de \: \ ] ^ = _ ] ^ (3.9) El sistema de ecuaciones algebraicas (3.9) también se puede expresar como  \ − _ ℐ  ] ^ = 0 a ^ (3.10) que tiene infinitas soluciones no triviales b]^≠0c si se cumple que det  \ − _ ℐ  = 0 (3.11) donde ℐ es la matriz identidad. La ecuación (3.11) se denomina ecuación secular (o característica) y permite calcular los autovalores de una matriz \. Asimismo, el sistema de ecuaciones (3.10) permite calcular los autovectores asociados a cada autovalor de una matriz \. Por tanto, para calcular los autovalores de [ hay que resolver la siguiente ecuación: det ⎝ ⎜ ⎛ − _ −   > 9 > 9  9 > 9 > 9 − _ ⎠ ⎟ ⎞ =_  + 9   =0 que es cuadrática en _ y sus soluciones son _ 9 =gh 9   ; _  =−gh 9   donde g=√−1, que es la unidad imaginaria. Seguidamente, vamos a calcular los autovectores asociados a cada autovalor de [. Para el autovalor _ 9 , el sistema de ecuaciones (3.10) toma la forma siguiente: D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 21 ⎝ ⎜ ⎛ − g h  9   −   > 9 > 9  9 > 9 > 9 − g h  9   ⎠ ⎟ ⎞  ] 99 ] 9 =G 0 0 I (3.12) donde ] 99 y ] 9 son las componentes de un autovector ]^ 9 asociado al autovalor _ 9 . Un posible autovector que satisface el sistema de ecuaciones (3.12) es ]^ 9 =U > 9 h  −g > 9 h 9 V De forma análoga, los autovectores asociados al autovalor _  cumplen que: ⎝ ⎜ ⎛ g h  9   −   > 9 > 9  9 > 9 > 9 g h  9   ⎠ ⎟ ⎞  ] 9 ]  =G 0 0 I (3.13) donde ] 9 y ]  son las componentes de un autovector ]^  asociado al autovalor _  . Un posible autovector que satisface el sistema de ecuaciones (3.13) es ]^  =U> 9 h  g > 9 h 9 V Con los autovectores ]^ 9 y ]^  procedemos a construir la siguiente matriz: j=b]^ 9 ,]^  c=U > 9 h  > 9 h  −g > 9 h 9 g > 9 h 9 V que recibe el nombre de matriz de cambio de base porque, para resolver el sistema de ecuaciones (3.8) es necesario desacoplarlo haciendo el siguiente cambio de base:   9    = j U  k 9  k  V (3.14) Si multiplicamos el sistema de ecuaciones (3.8) por la inversa de la matriz de cambio de base j T9 y sustituimos en él la ecuación (3.14), obtenemos lo siguiente: j T9 j ⎝ ⎜ ⎛   k 9    k   ⎠ ⎟ ⎞ =j T9 ⎝ ⎜ ⎛ 0 −   > 9 > 9  9 > 9 > 9 0 ⎠ ⎟ ⎞ jU  k 9  k  V+j T9 ⎝ ⎜ ⎛  9   > 9 −  9   > 9 ⎠ ⎟ ⎞ (3.15) D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 22 Llegados a este punto, es importante notar que los autovectores ]^ 9 y ]^  deben ser linealmente independientes para que exista j T9 . Como sabemos, la inversa de una matriz \ se define como el cociente entre su adjunto traspuesto y su determinante: \ T 9 = b Adj  \  c  det  \  (3.16) Haciendo uso de este teorema, la matriz inversa de j es: j T9 =12 ⎝ ⎜ ⎛ 1 > 9 √  g > 9 √ 9 1 > 9 √  −g > 9 √ 9 ⎠ ⎟ ⎞ Sustituyendo en el sistema de ecuaciones (3.15), se obtienen dos ecuaciones diferenciales ordinarias de primer orden, lineales y no homogéneas:   k 9  = g h  9    k 9 + 1 2 > 9 > 9 b  9 h   − g   h  9 c (3.17)   k   = − g h  9    k  + 1 2 > 9 > 9 b  9 h   + g   h  9 c (3.18) En los apéndices A.4 y A.5 se resuelven estas ecuaciones y se demuestra que, a partir de sus soluciones, se puede obtener la solución analítica del sistema de ecuaciones (3.8) que viene dada por (21) :  9    = n cos b h  9    + r c +   > 9 (3.19)      = s sin b h  9    + r c +  9 > 9 (3.20) donde la frecuencia angular (o pulsación) es u=√ 9   y n = v   9  −   > 9   + > 9    > 9   9     −  9 > 9   (3.21) (21) La demostración de las ecuaciones (3.19), (3.20), (3.21), (3.22) y (3.23) se encuentra en el apéndice A.6. D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 23 Figura 3.3. (a) Trayectorias en el espacio de las fases descritas por el sistema de ecuaciones (3.1). (b) Pequeñas fluctuaciones en el tiempo en el ciclo w 9 : el número de individuos  9 de la especie presa 8 9 fluctúa alrededor de su valor medio 'E HEB y el número de individuos   de la especie depredadora 8  fluctúa alrededor de su valor medio 'B HBE . Se han tomado los siguientes valores hipotéticos para los parámetros del sistema de ecuaciones (3.1):  9 =0,15; > 9 =3·10 Tz ;   =0,1; > 9 =2·10 Tz . s = > 9 > 9 v  9   n (3.22) r = arctan } > 9 > 9 v    9 G    −  9 > 9 I G  9  −   > 9 I ~ (3.23) donde  9 y   son los valores de  9 y   en el instante inicial (=0). Las funciones (3.19) y (3.20) son periódicas (ver figura 3.3.(b)) y su período vale • = 2 € u = 2 € √  9   (3.24) Como el desfase angular entre ellas es ∆‚= ƒ  radianes, el desfase temporal es ∆  = ∆ ‚ u = € 2 √  9   (3.25) D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 24 De las ecuaciones (3.19) y (3.20) se deduce que las trayectorias en el espacio de las fases son elipses (ver figura 3.3.(a)) descritas por la siguiente ecuación (22) : 1 n    9 −   > 9   + 1 s     −  9 > 9   = 1 (3.26) Cuando las trayectorias en el espacio de las fases son curvas cerradas alrededor del punto de equilibrio, se dice que es un punto vórtice o centro. Este tipo de punto es estable (23) . 3.3 Solución numérica general Del estudio analítico en un entorno muy próximo a los puntos de equilibrio del sistema de ecuaciones (3.1), hemos deducido que las trayectorias de interés en el espacio de las fases son ciclos recorridos en sentido antihorario que se encuentran restringidos en el primer cuadrante (24) . En esta sección vamos a ocuparnos de estudiar el comportamiento general de las funciones  9  y   . Como se adelantó en la introducción de la sección 3, no es posible obtener expresiones analíticas generales de las funciones  9  y    y, por ello, es necesario recurrir a métodos numéricos para la integración del sistema de ecuaciones (3.1) cuando la condición inicial  9 ,   está significativamente lejos de los puntos de equilibrio. En la figura 3.4.(a) podemos observar que, cuanto más se alejan las coordenadas  9 ,   del punto de equilibrio G 'E HEB , 'B HBE I, las trayectorias en el espacio de las fases continúan siendo curvas cerradas y, por tanto, las evoluciones temporales de las especies 8 9 (presa) y 8  (depredadora) siguen manteniendo un carácter periódico, aunque, cada vez menos sinusoidal. Esto hace que las trayectorias en el espacio de las fases vayan dejando de ser elípticas. (22) La demostración de la ecuación (3.26) se encuentra en el apéndice A.7. (23) Samanta y Gómez Aíza, 2014. (24) Lotka, 1920b, 1925; Volterra, 1926. D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 25 Figura 3.4. (a) Trayectorias en el espacio de las fases descritas por el sistema de ecuaciones (3.1). (b) Fluctuaciones en el tiempo en los ciclos w 9 y w „ : el número de individuos  9 de la especie presa 8 9 fluctúa alrededor de su valor medio 'E HEB y el número de individuos   de la especie depredadora 8  fluctúa alrededor de su valor medio 'B HBE . Para la integración del sistema de ecuaciones (3.1) se ha utilizado el método de Runge- Kutta de cuarto orden (algoritmo ode45 de MATLAB ® ) con los siguientes valores hipotéticos para los parámetros:  9 =0,15; > 9 =3·10 Tz ;   =0,1; > 9 =2·10 Tz . D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 32 Basándose en estas premisas, McKendrick y Kermack propusieron el siguiente sistema de ecuaciones diferenciales (33) : ⎩ ⎪ ⎪ ⎨ ⎪ ⎪ ⎧ 8  = − >8• • =>8•−• –  = • (4.1.a) (4.1.b) (4.1.c) donde el tamaño de la población total afectada por la enfermedad () viene dado por  = 8    + •    + –    Las ecuaciones (4.1) reciben el nombre de modelo epidemiológico SIR (34) , porque permite visualizar esquemáticamente el flujo de la población durante la epidemia: Cabe mencionar que las ecuaciones (4.1.a) y (4.1.b) constituyen un caso particular del modelo depredador-presa de Lotka y Volterra: la población susceptible de contraer la enfermedad es análoga a la especie presa  9 ≡8 y la población infectada es análoga a la especie depredadora   ≡• (ver ecuaciones (3.1) en la página 16). En el modelo epidemiológico SIR, el coeficiente asociado al crecimiento de la población susceptible de contraer la enfermedad ( 9 ) es cero, ya que, por hipótesis, el número de estos individuos únicamente puede disminuir. Además, como este decrecimiento supone un aumento instantáneo del número de individuos infectados, se cumple que la tasa per cápita de decrecimiento de la población susceptible debido a su interacción con la población infectada (> 9 ) es igual a la tasa per cápita de crecimiento de la población infectada debido a su interacción con la población susceptible (> 9 ) (35) . (33) McKendrick y Kermack, 1927. (34) Small, 2012. (35) Volterra, 1926. Susceptibles (Susceptible) Infectados (Infected) Eliminados (Removed) D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 33 Figura 4.1. Trayectorias en el espacio de las fases I-S descritas por el sistema de ecuaciones (4.1). Para integrar estas ecuaciones se ha utilizado el método de Runge- Kutta de cuarto orden (algoritmo ode45 de MATLAB ® ), suponiendo un total de cien mil individuos afectados por la epidemia =10   y los siguientes valores hipotéticos para los parámetros: =0,6; >=2·10 T . La recta =8  +•  indica el número de individuos susceptibles de contraer la enfermedad y de individuos infectados al comenzar la epidemia, cuando todavía no hay individuos fallecidos ni recuperados –  =0. En la figura 4.1 podemos observar que la integración del sistema de ecuaciones (4.1) conlleva a trayectorias en el espacio de las fases I-S que no son curvas cerradas y, por tanto, las funciones 8, • y – no son periódicas. Como ya sabemos, dado que las ecuaciones (4.1) no son lineales y están acopladas entre sí, solo podemos conocer el comportamiento de estas funciones mediante métodos numéricos de integración. D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 34 A continuación, vamos a comprobar el alcance de este modelo ajustando los datos que la Universidad de Johns Hopkins ha publicado sobre la pandemia causada por el coronavirus SARS-CoV-2 (más conocido como COVID-19) en China (36) . 4.1 Ajuste del modelo SIR a los datos de la pandemia del COVID-19 en China 4.1.1 Modelo SIR clásico: Susceptibles, Infectados y Eliminados Empecemos por definir en MATLAB ® los siguientes vectores: t_datos, que contiene los días de duración de la pandemia del COVID-19 en China; I_datos, que recoge el número total de individuos infectados que se registraron cada día; y R_datos, que recoge la suma del número total de individuos fallecidos y del número total de individuos recuperados que se registraron cada día. En el comienzo de una epidemia, el efecto del término −• en la ecuación (4.1.b) es prácticamente inapreciable frente al término >8•. Esto nos permite hacer una primera estimación del parámetro > mediante diferencias finitas: delta0 = (I_datos(2)-I_datos(1))/(S0*I_datos(1)*(t_datos(2)-t_datos(1))); suponiendo que, en el instante inicial, la población susceptible de contraer la enfermedad es del orden de cien mil individuos 8  =10  , puesto que, el dato exacto lo desconocemos. Por otro lado, cuando el número de individuos infectados es máximo se cumple que • @ ™-™š›œ• =0 de donde podemos inferir una primera estimación del parámetro :   ≈>  8 ∗  siendo  ∗ el instante de tiempo en el que se registra el número máximo de individuos infectados. (36) Johns Hopkins University, 2020. D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 35 El número de individuos susceptibles de contraer la enfermedad en el instante  ∗ se puede estimar integrando la ecuación (4.1.a): 8 ∗ ≈8   THŸ ™ ‰ ¡∗ Ÿ Asimismo, •  ∗  podemos calcularla utilizando el método de integración trapezoidal, lo cual permite ejecutar los siguientes comandos: k = find(I_datos==max(I_datos)); r0 = delta0*S0*exp(-delta0*trapz(I_datos(1:k))); A continuación, vamos a definir una función que tenga como argumentos de entrada el vector w0 = [S0,delta0,r0]; así como los vectores t_datos , I_datos y R_datos . Esta función debe integrar las ecuaciones (4.1) mediante el método de Runge-Kutta de cuarto orden (algoritmo ode45 de MATLAB ® ) tomando como condición inicial el vector [w0(1),I_datos(1),R_datos(1)] y, a su vez, devolver el siguiente valor: ¢=£b• ¤ −• ¤‰¥¦§ c  +b– ¤ −– ¤‰¥¦§ c  ¨ ¤- donde © es el número de días de duración de la epidemia y los parámetros • ¤ y – ¤ denotan, respectivamente, los individuos infectados y eliminados de la interacción susceptible-infectado estimados por el modelo SIR en el día j-ésimo: function Epsilon = epsilon(w0,t_datos,I_datos,R_datos) dNdt = @(t,N) [-w0(2)*N(1)*N(2);-N(2)*(w0(3)- w0(2)*N(1));w0(3)*N(2)]; [~,N] = ode45(dNdt,t_datos,[w0(1),I_datos(1),R_datos(1)],... odeset('RelTol',1e-5,'AbsTol',1e-6)); difI = N(:,2)-I_datos(1:length(t_datos)); difR = N(:,3)-R_datos(1:length(t_datos)); Epsilon = difI'*difI+difR'*difR; end D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 36 Finalmente, mediante el algoritmo fminsearch, vamos a optimizar la función definida anteriormente de manera que el valor de ¢ que devuelva sea mínimo: w = fminsearch(@epsilon,w0,[],t_datos,I_datos,R_datos); Así, obtenemos un vector w con los valores optimizados de los parámetros 8  , > y  que nos permite integrar el sistema de ecuaciones (4.1) tomando como condición inicial el vector [w(1),I_datos(1),R_datos(1)]: dNdt = @(t,N) [-w(2)*N(1)*N(2);-N(2)*(w(3)- w(2)*N(1));w(3)*N(2)]; [t,N] = ode45(dNdt,[0,300],[w(1),I_datos(1),R_datos(1)],... odeset('RelTol',1e-5,'AbsTol',1e-6)); De este modo, podemos comparar gráficamente el modelo epidemiológico SIR con el ajuste de éste a los datos experimentales: Figura 4.2. Ajuste del modelo epidemiológico SIR a los datos que la Universidad de Johns Hopkins ha publicado sobre la pandemia causada por el coronavirus COVID-19 en China. Para integrar el sistema de ecuaciones (4.1) se ha utilizado el método de Runge-Kutta de cuarto orden (algoritmo ode45 de MATLAB ® ) con los siguientes valores optimizados de los parámetros: 8  =85518; >=3,7·10 T¬ ; =0,047. D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 37 4.1.2 Modelo SIR extendido: Susceptibles, Infectados, Fallecidos y Altas Según las características que tenga la enfermedad epidémica en cuestión, se puede dividir la población total afectada incorporando tantas clases de individuos como sean necesarias, por ejemplo: aquellos que presentan síntomas de estar infectados, pero, que no son transmisores de la enfermedad; o los llamados asintomáticos, quienes aparentemente se encuentran sanos, pero, son portadores del agente infeccioso siendo capaces de transmitirlo (37) . Dado que la Universidad de Johns Hopkins da a conocer los datos relacionados con el número de individuos infectados, fallecidos y recuperados, es interesante afinar el modelo epidemiológico SIR modificando el sistema de ecuaciones (4.1) de la forma siguiente: ⎩ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎧ 8  = − >8• • =>8•− ® •− • • ¯ = ® •   =  • • (4.2) donde ¯ es el número de individuos fallecidos y  es el número de individuos dados de alta. Para ajustar los datos al sistema de ecuaciones (4.2) empezaremos por definir los vectores t_datos, I_datos, y estos otros: F_datos y A_datos que recogen, respectivamente, el número total de individuos fallecidos y dados de alta que se registraron cada día. Seguidamente, procederemos a destacar los cambios que hay que (37) Sánchez Garduño et al. (2002); Small, 2012. D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 38 realizar en la rutina de MATLAB ® descrita en la sección 4.1.1, la cual se encuentra íntegramente en el apéndice B.10. En primer lugar, puesto que = ® + • , vamos a considerar como primera aproximación de los parámetros  ® y  • que  ®Ÿ = •Ÿ ≈  2 donde   es una primera estimación del parámetro , que se obtuvo en la sección 4.1.1. Evidentemente, se trata de una burda aproximación, pero, concuerda en el orden de magnitud. Esto nos permite redefinir el vector w0 del siguiente modo: w0 = [S0,delta0,r0/2,r0/2]; Como sabemos, este vector es necesario para construir la siguiente función: function Epsilon_Mod = epsilon_mod(w0,t_datos,I_datos,F_datos,A_datos) dNdt = @(t,N) [-w0(2)*N(1)*N(2);-N(2)*(w0(3)+w0(4)- w0(2)*N(1));w0(3)*N(2);w0(4)*N(2)]; [~,N] = ode45(dNdt,t_datos,[w0(1),I_datos(1),F_datos(1),A_datos(1)],... odeset('RelTol',1e-5,'AbsTol',1e-6)); difI = N(:,2)-I_datos(1:length(t_datos)); difF = N(:,3)-F_datos(1:length(t_datos)); difA = N(:,4)-A_datos(1:length(t_datos)); Epsilon_Mod = difI'*difI+difF'*difF+difA'*difA; end la cual devuelve el valor ¢=£b• ¤ −• ¤‰¥¦§ c  +b¯ ¤ −¯ ¤‰¥¦§ c  +b ¤ − ¤‰¥¦§ c  ¨ ¤- Igualmente, vamos a minimizar el valor ¢ haciendo uso del algoritmo fminsearch: w = fminsearch(@epsilon_mod,w0,[],t_datos,I_datos,F_datos,A_datos); D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 39 Así, obtenemos un vector w con los valores optimizados de los parámetros 8  , >,  ® y  • que nos permite integrar el sistema de ecuaciones (4.2) tomando como condición inicial el vector [w(1),I_datos(1),F_datos(1),A_datos(1)]: dNdt = @(t,N) [-w(2)*N(1)*N(2);-N(2)*(w(3)+w(4)- w(2)*N(1));w(3)*N(2);w(4)*N(2)]; [t,N] = ode45(dNdt,[0,300],[w(1),I_datos(1),F_datos(1),A_datos(1)],... odeset('RelTol',1e-5,'AbsTol',1e-6)); De este modo, podemos comparar gráficamente el modelo extendido SIR con el ajuste de éste a los datos experimentales: Figura 4.3. Ajuste del modelo descrito por el sistema de ecuaciones (4.2) a los datos que la Universidad de Johns Hopkins ha publicado sobre la pandemia causada por el coronavirus COVID-19 en China. Para integrar esas ecuaciones se ha utilizado el método de Runge-Kutta de cuarto orden (algoritmo ode45 de MATLAB ® ) con los siguientes valores optimizados de los parámetros: 8  =85967; >=3,7·10 T¬ ;  ® =0,002,  • =0,045. D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 40 Los ajustes de los modelos SIR clásico y SIR extendido nos permiten estimar el instante de tiempo  ∗ en el que se produce el máximo en la curva de individuos infectados por COVID-19 en China desde el comienzo de la pandemia. En la siguiente tabla se muestran dichas estimaciones y sus desviaciones relativas (en %) al valor real de  ∗ :  '‘¥± ∗ (días)  ²™³ M± ¥ œ §´M¦ ∗ (días) |  '‘¥± ∗ −  ²™³ M± ¥ œ §´M¦ ∗ |  '‘¥± ∗ (%)  ²™³ ‘,‘¨‰´‰¦ ∗ (días) |  '‘¥± ∗ −  ²™³ ‘,‘¨‰´‰¦ ∗ |  '‘¥± ∗ (%) 26 25,91 0,3 25,98 0,08 Tabla 4.1. Valor real del instante de tiempo  ∗ en el que se produce el máximo de la curva de individuos infectados por COVID-19 en China desde el comienzo de la pandemia. También se presentan las estimaciones devueltas por los modelos SIR clásico y SIR extendido, acompañadas de sus desviaciones relativas (en %) al valor real de  ∗ . Como hemos podido comprobar, la versión clásica del modelo epidemiológico SIR determina con una precisión del 99,7 % que el máximo de la curva de individuos infectados en China por COVID-19 se produjo, aproximadamente, veintiséis días después del comienzo de la pandemia; mientras que, el modelo SIR extendido determina dicho máximo con una precisión del 99,9 %. D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 41 Conclusiones En esta memoria hemos analizado el modelo depredador-presa, definido por los matemáticos A. J. Lotka y V. Volterra entre los años 1925 y 1926. Mediante métodos numéricos y analíticos de integración, hemos podido verificar la ley establecida por Vito Volterra sobre la periodicidad de las fluctuaciones de las especies presa y depredadora, la cual afirma que, en general, el período de estas oscilaciones depende de los parámetros  9 ,   , > 9 y > 9 —todos ellos positivos— así como del valor que toman los individuos de ambas especies en el instante inicial; mientras que para fluctuaciones muy pequeñas solo depende de los parámetros  9 y   . Recordemos que los parámetros  9 y   indican, respectivamente, las tasas per cápita netas (nacimientos menos muertes) de las especies presa 8 9 y depredadora 8  ; > 9 es la tasa per cápita de decrecimiento de la especie presa 8 9 debido a su interacción con la especie depredadora 8  ; y > 9 es la tasa per cápita de crecimiento de la especie depredadora 8  debido a su interacción con la especie presa 8 9 . Por otra parte, también hemos podido comprobar la segunda ley enunciada por Volterra, la cual establece que, independientemente del valor inicial que tomen los individuos de ambas especies, el promedio del número de individuos de la especie presa es  9M = 'E HEB y el promedio del número de individuos de la especie depredadora es  M = 'B HBE . Llegados a este punto, no es difícil convencerse de la tercera y última ley expuesta por este insigne matemático. Podríamos decir que es la que da una respuesta a la cuestión planteada por U. D’Ancona en 1926 justificando que, ante una intensa injerencia del ser humano en la asociación de dos especies que conviven, el coeficiente  9 , que da cuenta del crecimiento de la especie más inofensiva, disminuye; mientras que el coeficiente   , asociado al decrecimiento de la especie más voraz, aumenta. En estas D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 48   =  U 9  9 V T'E 'B A.4 Resolución de la ecuación (3.17) Vamos a resolver la ecuación (3.17) por el método de variación de parámetros. Si ignoramos el término no homogéneo, nos queda una ecuación diferencial análoga a la ecuación (1.4): k 9º  =gh 9   k 9º cuya solución es  k 9 º    =  k 9   ´ √ ' B ' E  (a.4) Supongamos una solución particular de la ecuación (3.17) de la forma  k 9 »    = (     ´ √ ' B ' E  (a.5) Si introducimos la ecuación (a.5) en la ecuación (3.17), obtenemos que (  = 9 √  −g   √ 9 2> 9 > 9  T´√'B'E e integrando (=  √ 9 +g  9 √  2> 9 > 9 √ 9    T´√'B'E +… donde … es una constante arbitraria que, sin pérdida de generalidad, podemos tomar como cero. Por tanto, la solución particular (a.5) queda de la forma  k 9 »    =   √  9 + g  9 √   2 > 9 > 9 √  9   (a.6) Finalmente, la solución de la ecuación (3.17) podemos escribirla como la suma de las ecuaciones (a.4) y (a.6):  k 9    =  k 9   ´ √ ' B ' E  +   √  9 + g  9 √   2 > 9 > 9 √  9   (a.7) D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 49 A.5 Resolución de la ecuación (3.18) Al igual que hemos hecho en el apartado anterior, vamos a resolver la ecuación (3.18) por el método de variación de parámetros. Si ignoramos el término no homogéneo, nos queda una ecuación diferencial análoga a la ecuación (1.4): k º  =−gh 9   k º cuya solución es  k  º    =  k    T ´ √ ' B ' E  (a.8) Supongamos una solución particular de la ecuación (3.18) de la forma  k  »    = (     T ´ √ ' B ' E  (a.9) Si introducimos la ecuación (a.9) en la ecuación (3.18), obtenemos que (  = 9 √  +g   √ 9 2> 9 > 9  ´√'B'E e integrando (=  √ 9 −g  9 √  2> 9 > 9 √ 9    ´√'B'E +… donde … es una constante arbitraria que, sin pérdida de generalidad, podemos tomar como cero. Por tanto, la solución particular (a.9) queda de la forma  k  »    =   √  9 − g  9 √   2 > 9 > 9 √  9   (a.10) Finalmente, la solución de la ecuación (3.18) podemos escribirla como la suma de las ecuaciones (a.8) y (a.10):  k     =  k    T ´ √ ' B ' E  +   √  9 − g  9 √   2 > 9 > 9 √  9   (a.11) A.6 Demostración de las ecuaciones (3.19), (3.20), (3.21), (3.22) y (3.23) A partir de las ecuaciones (3.14), (a.7) y (a.11), se obtienen las soluciones  9  y    del sistema de ecuaciones (3.8): D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 50  9    = > 9 h   b  k 9   ´ √ ' B ' E  +  k    T ´ √ ' B ' E  c +   > 9 (a.12)      = − g > 9 h  9 b  k 9   ´ √ ' B ' E  −  k    T ´ √ ' B ' E  c +  9 > 9 (a.13) Si hacemos  9 =0= 9 y   =0=  , obtenemos un sistema de dos ecuaciones algebraicas con dos incógnitas (k 9 y k  ):  9 => 9 h  bk 9 +k  c+  > 9   =−g > 9 h 9 bk 9 −k  c+ 9 > 9 cuya solución es  k 9  = 1 2 > 9 √     9  −   > 9  + g 2 > 9 √  9     −  9 > 9  (a.14)  k   = 1 2 > 9 √     9  −   > 9  − g 2 > 9 √  9     −  9 > 9  (a.15) Si sustituimos las expresiones (a.14) y (a.15) en (a.12) y (a.13), obtenemos que  9    =   9  −   > 9  cos b h  9    c − > 9 > 9 v    9     −  9 > 9  sin b h  9    c +   > 9 (a.16)      =     −  9 > 9  cos b h  9    c + > 9 > 9 v  9     9  −   > 9  sin b h  9    c +  9 > 9 (a.17) No obstante, las expresiones (a.16) y (a.17) se pueden presentar de un modo más conveniente haciendo uso de los teoremas de la suma de dos ángulos:  9 =ncosr cosh 9   −nsinr sinh 9   +  > 9   =ssinr cosh 9   +scosr sinh 9   + 9 > 9 Inmediatamente podemos hacer las siguientes identificaciones: n sin r = > 9 > 9 v    9     −  9 > 9  (a.18) n cos r =  9  −   > 9 (a.19) s sin r =    −  9 > 9 (a.20) D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 51 s cos r = > 9 > 9 v  9     9  −   > 9  (a.21) Si elevamos al cuadrado las ecuaciones (a.18) y (a.19) y seguidamente las sumamos, obtenemos el valor del parámetro n: n=v 9 −  > 9   +> 9    > 9   9   − 9 > 9   Análogamente, si elevamos al cuadrado las ecuaciones (a.20) y (a.21) y seguidamente las sumamos, obtenemos el valor del parámetro s: s=v> 9   9 > 9     9 −  > 9   +  − 9 > 9   => 9 > 9 v 9   n Por otra parte, si dividimos la ecuación (a.18) entre la ecuación (a.19) o la ecuación (a.20) entre la ecuación (a.21), obtenemos el valor del parámetro r: r=arctan}> 9 > 9 v   9 G  − 9 > 9 I G 9 −  > 9 I~ A.7 Demostración de la ecuación (3.26) Comencemos por despejar el término cos√ 9   +r de la expresión (3.18) y el término sin√ 9   +r de la expresión (3.20): cos b h  9    + r c = 1 n   9    −   > 9  (a.22) sin b h  9    + r c = 1 s       −  9 > 9  (a.23) Si elevamos al cuadrado las ecuaciones (a.22) y (a.23) y seguidamente las sumamos, obtenemos: 1 n   9 −  > 9   +1 s    − 9 > 9   =1 A.8 Demostración de la expresión (3.28) Partamos de la ecuación diferencial (3.27). Si separamos sus variables de la forma    9 −> 9  9 =− 9   −> 9   D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 52 podemos integrarla trivialmente obteniendo ¼ = > 9   + > 9  9 −  9 ln   −   ln  9 (a.24) A continuación, vamos a restar a la expresión (a.24) su valor en el punto crítico G 'E HEB , 'B HBE I: w = > 9   + > 9  9 −  9 ln   −   ln  9 −  9 −   +  9 ln  9 > 9 +   ln   > 9 (a.25) siendo w=¼ 9 ,  −¼G 'E HEB , 'B HBE I. Por último, haciendo uso de las propiedades (a.1) y (a.3) y teniendo en cuenta que ln=1, la expresión (a.25) queda de la forma w = > 9  9 + > 9   + ln U  > 9   9    T ' B  > 9     9  T ' E V A.9 Función ½ de Lambert en ℝ La función Š de Lambert o función producto-logaritmo se define como la función inversa de ()=) , , es decir, )=Š¾()¿. Por ello, resulta de interés estudiar las características de (: Dominio: −∞,+∞ Comportamiento en ∞: lim ,→T7 ()= lim ,→T7 ) , = lim ,→67 G−)  , I=0⟹Asíntota en )=0 lim ,→67 ()= lim ,→67 ) , =+∞ Punto crítico de primera especie: () ) . ,-,C =1+) , | ,-,C =0⟹) M =−1⇒−1,−1   () )  . ,-T9 =2+) , | ,-T9 =1>0⟹Mínimo relativo Punto crítico de segunda especie:   () )  . ,-,C =2+) , | ,-,C =0⇒) M =−2⇒−2,−2   ⇒Cambio de curvatura D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 53 Figura II. Representación de la función ()=) , . Como se puede observar en la figura II, a cada valor de (∈G− 9‘ ,0I le corresponden dos valores de )∈−∞,0. En otras palabras, la función ( no es inyectiva. Esto obliga a descomponer la función de Lambert en las siguientes ramas: • A cada valor )∈−∞,−1 le corresponde un único valor de (∈G− 9‘ ,0I. La rama de la función de Lambert que cumple esta condición de denota como Š T9 . • A cada valor )∈−1,+∞ le corresponde un único valor de (∈G− 9‘ ,+∞I. La rama de la función de Lambert que cumple esta condición de denota como Š  . A.10 Integración de la ecuación (3.31) en un entorno próximo al punto vórtice De la ecuación (3.29) se infiere que > 9    9  T H EB ' E  A B T ' E H EB  = ‡ 9 ' E (a.26) > 9  9    T H BE ' B  A E T ' B H BE  = G ‡ … I T 9 ' B (a.27) D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 54 Si hacemos las siguientes definiciones © 9 =  9 −   > 9 (a.28) ©  =   −  9 > 9 (a.29) donde © 9 y ©  representan las fluctuaciones de  9 y   respecto al punto de equilibrio G 'E HEB , 'B HBE I, las ecuaciones (a.26) y (a.27) toman la siguiente forma:  1 + > 9   © 9   T H EB ' E ¨ B = ‡ 9 ' E (a.30)  1 + > 9  9 ©    T H BE ' B ¨ E = G ‡ … I T 9 ' B (a.31) Haciendo uso de los términos de primer orden Â=1 de los siguientes desarrollos en serie de potencias 1+n)  T¥, =1+n)  £−1 à  Â+1!n) ÃT9 7 Ã-9 n) Å =1+£GsÂIn)−1 à 7 Ã-9 las ecuaciones (a.30) y (a.31) se reducen a 12  > 9   © 9   = 1    1 − ‡  (a.32) 12  > 9  9 ©    = 1  9 G ‡ … − 1 I (a.33) en un entorno muy próximo al punto de equilibrio G 'E HEB , 'B HBE I. Trivialmente, se puede despejar © 9 y ©  de las ecuaciones (a.32) y (a.33): © 9 = ± 1 > 9 h 2    1 − ‡  (a.34) ©  = ± 1 > 9 Ç 2  9 G ‡ … − 1 I (a.35) Teniendo en cuenta las ecuaciones (a.28), (a.29), (a.34) y (a.35), podemos integrar la ecuación (3.31) en los procesos →Ž, Ž→…, …→• y •→: D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 55 ∆  • → ’ = 1 > 9 > 9 ‹ ‡ ‡ © 9 ©  “ 9 . ¨BÈ ¨ E É  = − 1 2 √  9   ‹ ‡ ‡ Ç  1 − ‡  G ‡ … − 1 I “ 9 (a.36) ∆  ’ → “ = 1 > 9 > 9 ‹ ‡ ‡ © 9 ©  9 “ . ¨BÈ ¨ E È  = 1 2 √  9   ‹ ‡ ‡ Ç  1 − ‡  G ‡ … − 1 I 9 “ (a.37) ∆  “ → ” = 1 > 9 > 9 ‹ ‡ ‡ © 9 ©  “ 9 . ¨BÉ ¨ E È  = − 1 2 √  9   ‹ ‡ ‡ Ç  1 − ‡  G ‡ … − 1 I “ 9 (a.38) ∆  ” → • = 1 > 9 > 9 ‹ ‡ ‡ © 9 ©  9 “ . ¨BÉ ¨ E É  = 1 2 √  9   ‹ ‡ ‡ Ç  1 − ‡  G ‡ … − 1 I 9 “ (a.39) Evidentemente, los límites de integración son los mismos que en las integrales (3.36), (3.37), (3.38) y (3.39). Si nos fijamos bien, las integrales (a.36), (a.37), (a.38) y (a.39), que indican los intervalos de tiempo de los diferentes procesos de una trayectoria de fase muy próxima al punto vórtice (ver figura 3.4.(a)), son idénticas: ∆  • → ’ = ∆  ’ → “ = ∆  “ → ” = ∆  ” → • = 1 2 √  9   ‹ ‡ ‡ Ç  1 − ‡  G ‡ … − 1 I 9 “ (a.40) La integral (a.40) se puede resolver analíticamente haciendo el cambio de variable Ê=v‡…−1 1−‡ obteniendo de una forma alternativa el resultado de la ecuación (3.25): ∆  • → ’ = ∆  ’ → “ = ∆  “ → ” = ∆  ” → • = 1 √  9   ‹ Ê 1 + Ê  7  = € 2 √  9   De acuerdo con la ecuación (3.40), el período viene dado por •= 2€ √ 9   coincidiendo con el resultado de la ecuación (3.24). D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 56 Apéndice B: Rutinas en MATLAB ® R2017a B.1 Figura 1.1 syms t; fplot(exp((log(2)/25)*t),'LineWidth',3); hold on; fplot(1+(1/25)*t,'LineWidth',3); plot([0,25,25],[2,2,0],'--k'); plot(25,2,'o','MarkerFaceColor',1/255*[255,147,0],'MarkerEdgeColor','k',... 'MarkerSize',10); axis([0,40,0.8,2.5]); xlabel('$$t\thinspace(a\tilde{n}os)$$','Interpreter','LaTeX'); ylabel('$$Cantidad$$ ','Rotation',0,... 'HorizontalAlignment','right','Interpreter','LaTeX'); legend({'$$\frac{N(t)}{N_0}$$','$$\frac{A(t)}{A_0}$$'},'Location','north',... 'Orientation','Horizontal','Interpreter','LaTeX'); set(gca,'layer','top'); set(gca,'FontSize',20,'TickLabelInterpreter','LaTeX'); B.2 Figura 1.2 syms x; fplot(exp(x),'LineWidth',3); hold on; fplot(exp(-x),'LineWidth',3); yticks(1:10); axis([0,5,0,10]); xlabel('$$rt$$','Interpreter','LaTeX'); ylabel('$$\frac{N}{N_0}$$ ','Rotation',0,'HorizontalAlignment','right',... 'Interpreter','LaTeX'); legend({'$$r>0$$','$$r<0$$'},'Interpreter','LaTeX'); set(gca,'layer','top'); set(gca,'FontSize',20,'TickLabelInterpreter','LaTeX'); B.3 Figura 2.1 syms x; fplot(exp(x)/((1/0.05)-(1-exp(x))),'LineWidth',3); hold on; fplot(exp(x)/((1/1.8)-(1-exp(x))),'LineWidth',3); plot([0,10],[1,1],'--k'); axis([0,10,0,2]); xlabel('$$rt$$','Interpreter','LaTeX'); ylabel('$$\frac{N}{\left(\frac{r}{\gamma}\right)}$$ ','Rotation',0,... 'HorizontalAlignment','right','Interpreter','LaTeX'); D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 57 legend({'$$N_0<\frac{r}{\gamma}$$','$$N_0>\frac{r}{\gamma}$$',... '$$N\left(t\right)=\frac{r}{\gamma}$$'},'Orientation','horizontal',... 'Location','north','Interpreter','LaTeX'); set(gca,'layer','top'); set(gca,'FontSize',20,'TickLabelInterpreter','LaTeX'); B.4 Figura 3.1 Selacios_Rijeka = [11.76 21.48 23.08 21 36.47 25.86 15.99 15.74 14.73 ... 10.52 13.24]; MenosVoraces_Rijeka = [54.86 22.25 22.58 20.25 27.19 24.86 35.57 31.99 ... 34.76 51.9 56.1]; plot(1914:1924,Selacios_Rijeka,'LineWidth',3); hold on; plot(1914:1924,MenosVoraces_Rijeka,'LineWidth',3); xlabel('$$A\tilde{n}o$$','Interpreter','LaTeX'); ylabel('$$Capturas \left(\%\right)$$ ','Rotation',0,... 'HorizontalAlignment','right','Interpreter','LaTeX'); legend({'Selacios','Especies menos voraces'},'Location','southoutside',... 'Orientation','horizontal','Interpreter','LaTeX'); set(gca,'layer','top'); set(gca,'FontSize',20,'TickLabelInterpreter','LaTeX'); xlim([1913,1925]); B.5 Figura 3.2 plot([0 0],[-6 6],'k--'); hold on; plot([-6 6],[0,0],'k--'); plot(0,0,'o','MarkerFaceColor',1/255*[255,147,0],'MarkerEdgeColor','k',... 'MarkerSize',10); r1 = 0.15; r2 = 0.1; N0 = [1,2,3]; N1 = linspace(1e-20,6); x = N1(1:end-1); u = diff(N1); n = [0 1 1 0 ; 0 0 1 1]; color = (1/255)*[0,114,189;190,190,190;190,190,190;190,190,190]; for i = 1:4 for j = 1:length(N0) N2 = N0(j)*(N1/N0(j)).^(-r2/r1); y = N2(1:end-1); v = diff(N2); plot((-1)^n(1,i)*N1,(-1)^n(2,i)*N2,'Color',color(i,:),'LineWidth',3); for k = 1:25:length(N1) D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 64 f4 = @(x) 1./(r1*r2*x.*(1+lambertw(0,(-(x.^(1/r2))/exp(1)))).*... (1+lambertw(0,(-((C(i)./x).^(1/r1))/exp(1))))); T4(i) = integral(f4,C(i),1,'RelTol',1e-1,'AbsTol',1e-1); end T = T1+T2+T3+T4; subplot(1,2,1); semilogx(V,T/((2*pi)/sqrt(r1*r2)),'LineWidth',3); hold on; fplot(1,'k--'); axis([V(1),1e2,0.5,4]); xlabel('$$V$$','Interpreter','LaTeX'); ylabel('$$\frac{T}{\left(\frac{2\pi}{\sqrt{r_1r_2}}\right)}$$ ',... 'Rotation',0,'HorizontalAlignment','right','Interpreter','LaTeX'); set(gca,'FontSize',20,'TickLabelInterpreter','LaTeX'); plot2 = subplot(1,2,2); plot2.Position(1)=0.6; semilogx(V,T1/(pi/(2*sqrt(r1*r2))),'LineWidth',3); hold on; semilogx(V,T2/(pi/(2*sqrt(r1*r2))),'LineWidth',3); semilogx(V,T3/(pi/(2*sqrt(r1*r2))),'LineWidth',3); semilogx(V,T4/(pi/(2*sqrt(r1*r2))),'LineWidth',3); fplot(1,'k--'); fplot(1,'k--'); axis([V(1),1e2,0,4]); xlabel('$$V$$','Interpreter','LaTeX'); ylabel('$$\frac{\Delta t}{\left(\frac{\pi}{2\sqrt{r_1r_2}}\right)}$$ ',... 'Rotation',0,'HorizontalAlignment','right','Interpreter','LaTeX'); legend({'$$A\rightarrow B$$','$$B\rightarrow C$$','$$C\rightarrow D$$',... '$$D\rightarrow A$$'},'Location','northwest','Interpreter','LaTeX'); set(gca,'FontSize',20,'TickLabelInterpreter','LaTeX'); B.9 Figura 4.1 d = 2e-5 ; r = 0.6; dSIdt = @(t,SI) [-d*SI(1)*SI(2);-SI(2)*(r-d*SI(1))]; N = 1e5; S0 = (0.1:0.1:0.9)*1e5; I0 = N - S0; for i = 1:length(S0) [t,SI] = ode45(dSIdt,[0,300],[S0(i),I0(i)],odeset('RelTol',1e-12,'AbsTol',1e-12)); x = SI(1:end-1,1); u = diff(SI(:,1)); y = SI(1:end-1,2); v = diff(SI(:,2)); D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 65 plot(SI(:,1),SI(:,2),'color',(1/255)*[0,114,189],'LineWidth',3); hold on; for k = 200:6000:length(SI) h = annotation('arrow','headStyle','cback2',... 'color',(1/255)*[0,114,189],'HeadLength',15,'HeadWidth',15); set(h,'parent',gca); set(h,'position',[x(k),y(k),u(k),v(k)]); end end f = plot([0,N],[N,0],'k','LineWidth',3); axis equal; axis([0,N,0,N]); xlabel('$$S$$','Interpreter','LaTeX'); ylabel('$$I$$ ','Rotation',0,'HorizontalAlignment','right',... 'Interpreter','LaTeX'); legend(f,{'$$N=S_0+I_0$$'},'Orientation','horizontal',... 'Interpreter','LaTeX'); set(gca,'layer','top'); set(gca,'FontSize',20,'TickLabelInterpreter','LaTeX'); B.10 Figuras 4.2 y 4.3 I_datos = [503 595 858 1325 1970 2737 5277 5834 7835 9375 11357 15806 ... 18677 22373 25762 28477 31393 33413 35705 37424 38638 38560 52309 ... 56860 57452 57992 58108 58002 56541 54825 54608 51859 51390 49631 ... 47413 45365 42924 39809 37199 34898 32368 29864 27402 25230 23702 ... 22159 20335 18933 17567 16116 14859 13569 12124 10783 9906 9030 ... 8106 7372 6731 6189 5799 5410 5030 4603 4310 3881 3600 3236 2967 ... 2764 2640 2545 2425 2267 2062 2020 1973 1905 1865 1810 1794 1835 ... 1829 1761 1699 1656 1572 1537 1479 1436 1418 1371 1265 1154 1098 ... 998 907 881 833 796 749 736 643 537 461 404 361 346 226 206 176 159 ... 141 131 120 113 110 115 115 115 115 111 111 114 112 107 101 97 103 ... 104 119]'; F_datos = [17 18 26 42 56 82 131 133 171 213 259 361 425 491 563 633 718 ... D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 66 805 905 1012 1112 1117 1369 1521 1663 1766 1864 2003 2116 2238 2238 ... 2443 2445 2595 2665 2717 2746 2790 2837 2872 2914 2947 2983 3015 ... 3044 3072 3100 3123 3139 3161 3172 3180 3193 3203 3217 3230 3241 ... 3249 3253 3259 3274 3274 3281 3285 3291 3296 3299 3304 3308 3309 ... 3316 3322 3326 3330 3333 3335 3335 3337 3339 3340 3343 3343 3345 ... 3345 3346 3346 4636 4636 4636 4636 4636 4636 4636 4636 4636 4637 ... 4637 4637 4637 4637 4637 4637 4637 4637 4637 4637 4637 4637 4637 ... 4637 4637 4637 4637 4637 4637 4638 4638 4638 4638 4638 4638 4638 ... 4638 4638 4638 4638 4638 4638 4638 4638 4638]'; A_datos = [28 30 36 39 49 58 101 120 135 214 275 463 614 843 1115 1477 ... 1999 2596 3219 3918 4636 5082 6217 7977 9298 10755 12462 14206 ... 15962 18014 18704 22699 23187 25015 27676 30084 32930 36329 39320 ... 42162 44854 47450 50001 52292 53944 55539 57388 58804 60181 61644 ... 62901 64196 65660 67017 67910 68798 69755 70535 71266 71857 72362 ... 72814 73280 73773 74181 74720 75100 75582 75923 76206 76405 76565 ... 76760 76946 77207 77310 77410 77567 77679 77791 77877 77956 78039 ... 78200 78311 78401 77552 77614 77690 77745 77799 77861 77983 78109 ... 78175 78277 78374 78422 78474 78523 78573 78586 78684 78792 78870 ... 78929 78977 78993 79127 79167 79198 79222 79246 79261 79281 79293 ... 79306 79310 79310 79310 79310 79332 79335 79343 79352 79358 79367 ... 79371 79382 79386 79389]'; R_datos = F_datos + A_datos; t_datos = (0:length(I_datos)-1)'; S0 = 1e5; delta0 = (I_datos(2)-I_datos(1))/(S0*I_datos(1)*(t_datos(2)-t_datos(1))); k = find(I_datos==max(I_datos)) ; r0 = delta0*S0*exp(- delta0*trapz(I_datos(1:k))); figure(1); % Figura 4.2. w0 = [S0,delta0,r0]; D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 67 w = fminsearch(@epsilon,w0,[],t_datos,I_datos,R_datos); dNdt = @(t,N) [-w(2)*N(1)*N(2);-N(2)*(w(3)-w(2)*N(1));w(3)*N(2)]; [t,N] = ode45(dNdt,[0,300],[w(1),I_datos(1),R_datos(1)],... odeset('RelTol',1e-5,'AbsTol',1e-6)); plot(t,N(:,1),'Color',(1/255)*[255,192,0],'LineWidth',3); hold on; plot(t,N(:,2),'r','LineWidth',3); plot(t_datos,I_datos,'r.','MarkerSize',10); plot(t,N(:,3),'Color',1/255*[68,114,196],'LineWidth',3); plot(t_datos,R_datos,'.','Color',1/255*[68,114,196],'MarkerSize',10); axis([0,length(t_datos)-1,0,9e4]); set(gca,'layer','top'); xticks(0:10:length(t_datos)-1); yticks(0:1e4:9e4); xticklabels({'22 ene','1 feb','11 feb','21 feb','2 mar','12 mar',... '22 mar','1 abr','11 abr','21 abr','1 may','11 may','21 may','31 may'}); xtickangle(45); xlabel('$$D\acute{\imath}a$$','Interpreter','LaTeX'); ylabel('$$Poblaci\acute{o}n$$ ','Rotation',0,... 'HorizontalAlignment','right','Interpreter','LaTeX'); legend({'$$S(t)$$','$$I(t)$$','$$I_{datos}$$','$$R(t)$$','$$R_{datos}$$'},... 'Location','east','Interpreter','LaTeX'); set(gca,'FontSize',20,'TickLabelInterpreter','LaTeX'); m = find(N(:,2)==max(N(:,2))); tMaxI_SIR_clasico = t(m); figure(2); % Figura 4.3. w0 = [S0,delta0,r0/2,r0/2]; w = fminsearch(@epsilon_mod,w0,[],t_datos,I_datos,F_datos,A_datos); dNdt = @(t,N) [-w(2)*N(1)*N(2);-N(2)*(w(3)+w(4)- w(2)*N(1));w(3)*N(2);w(4)*N(2)]; [t,N] = ode45(dNdt,[0,300],[w(1),I_datos(1),F_datos(1),A_datos(1)],... odeset('RelTol',1e-5,'AbsTol',1e-6)); D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 68 plot(t,N(:,1),'Color',(1/255)*[255,192,0],'LineWidth',3); hold on; plot(t,N(:,2),'r','LineWidth',3); plot(t_datos,I_datos,'r.','MarkerSize',10); plot(t,N(:,3),'k','LineWidth',3); plot(t_datos,F_datos,'k.','MarkerSize',10); plot(t,N(:,4),'Color',(1/255)*[34,177,76],'LineWidth',3); plot(t_datos,A_datos,'.','Color',(1/255)*[34,177,76],'MarkerSize',10); axis([0,length(t_datos)-1,0,9e4]); set(gca,'layer','top'); xticks(0:10:length(t_datos)); yticks(0:1e4:9e4); xticklabels({'22 ene','1 feb','11 feb','21 feb','2 mar','12 mar',... '22 mar','1 abr','11 abr','21 abr','1 may','11 may','21 may','31 may'}); xtickangle(45); xlabel('$$D\acute{\imath}a$$','Interpreter','LaTeX'); ylabel('$$Poblaci\acute{o}n$$ ','Rotation',0,... 'HorizontalAlignment','right','Interpreter','LaTeX'); legend({'$$S(t)$$','$$I(t)$$','$$I_{datos}$$','$$F(t)$$','$$F_{datos}$$',... '$$A(t)$$','$$A_{datos}$$'},'Location','east','Interpreter','LaTeX'); set(gca,'FontSize',20,'TickLabelInterpreter','LaTeX'); m = find(N(:,2)==max(N(:,2))); tMaxI_SIR_extendido = t(m); % Función para el modelo SIR: function Epsilon = epsilon(w0,t_datos,I_datos,R_datos) dNdt = @(t,N) [-w0(2)*N(1)*N(2);-N(2)*(w0(3)-w0(2)*N(1));w0(3)*N(2)]; [~,N] = ode45(dNdt,t_datos,[w0(1),I_datos(1),R_datos(1)],... odeset('RelTol',1e-5,'AbsTol',1e-6)); difI = N(:,2)-I_datos(1:length(t_datos)); difR = N(:,3)-R_datos(1:length(t_datos)); Epsilon = difI'*difI+difR'*difR; end % Función para la variante del modelo SIR: D INÁMICA DE POBLACIONES : SOBRE LA BASE DEL MODELO DE L OTKA -V OLTERRA 69 function Epsilon_Mod = epsilon_mod(w0,t_datos,I_datos,F_datos,A_datos) dNdt = @(t,N) [-w0(2)*N(1)*N(2);-N(2)*(w0(3)+w0(4)- w0(2)*N(1));w0(3)*N(2);w0(4)*N(2)]; [~,N] = ode45(dNdt,t_datos,[w0(1),I_datos(1),F_datos(1),A_datos(1)],... odeset('RelTol',1e-5,'AbsTol',1e-6)); difI = N(:,2)-I_datos(1:length(t_datos)); difF = N(:,3)-F_datos(1:length(t_datos)); difA = N(:,4)-A_datos(1:length(t_datos)); Epsilon_Mod = difI'*difI+difF'*difF+difA'*difA; end B.11 Figura I [x,y] = meshgrid(0.1:0.1:0.8,(1e-3)*(1:0.1:5)); z = x./y; p = subplot(1,1,1); plot3(x,y,z); surf(x,y,z); p.Position(1) = 0.15; xlabel('$$r_i$$','Interpreter','LaTeX'); ylabel('$$\delta_{ij}$$','Interpreter','LaTeX'); zlabel('$$N_j^c=\frac{r_i}{\delta_{ij}}$$ ','Rotation',0,... 'HorizontalAlignment','right','Interpreter','LaTeX'); set(gca,'layer','top'); set(gca,'FontSize',20,'TickLabelInterpreter','LaTeX'); B.12 Figura II syms x; fplot(x*exp(x),[-8,-1],'LineWidth',3); hold on; fplot(x*exp(x),[-1,1],'LineWidth',3); plot([0,0],[-0.4,0.4],'k--'); plot([-8,1],[0,0],'k--'); axis([-8,1,-0.4,0.4]); xticks(-8:1:1); yticks(-0.5:0.1:0.5); xlabel('$$x$$','Interpreter','LaTeX'); ylabel('$$f(x)$$ ','Rotation',0,'HorizontalAlignment','right',... 'Interpreter','LaTeX'); legend({'$$W_{-1}$$','$$W_0$$'},'Location','northwest','Interpreter','LaTeX'); set(gca,'FontSize',20,'TickLabelInterpreter','LaTeX');