´ nea Matema ´ tica 58 (2014) 77-110 Miscela

SMM

Modelos din´amicos de poblaciones simples y de sistemas depredador-presa Guruprasad Samanta Permanent address: Department of Mathematics, Bengal Engineering and Science University, Shibpur, Howrah-711103, India [email protected]

Ricardo G´omez A´ıza Instituto de Matem´aticas de la Universidad Nacional Aut´onoma de M´exico Circuito Exterior, Ciudad Universitaria C.P. 04510, M´exico D.F., M´exico [email protected]

Resumen En este art´ıculo describiremos el comportamiento din´amico de los modelos matem´aticos b´asicos de poblaciones simples y de sistemas de dos especies depredador-presa interactuando. El primer modelo matem´atico de estos sistemas fue concebido por Lotka y de manera independiente por Volterra alrededor de 1925. Este modelo es uno de los puntos de partida de la ecolog´ıa matem´atica. Tiene dos puntos (o estados) de equilibrio, uno es trivial y el otro es no trivial (o de coexistencia). El estado de equilibrio trivial siempre es inestable y el no-trivial es localmente estable pero no es localmente asint´oticamente estable, lo cual ocurre debido a que este modelo no posee un mecanismo para estabilizar (asint´oticamente) la posici´on de equilibrio de coexistencia. Discutiremos tambi´en los efectos de la competencia intraespec´ıfica entre la poblaci´on de presas as´ı como diversas versiones y posibles modificaciones de este modelo.

78

´ G. SAMANTA Y R. GOMEZ A´IZA

N (t)

N (t) = N0 ert r>0

N0

t

0

Figura 1. El cl´erigo ingl´es Thomas Robert Malthus (1766–1834) propuso el modelo exponencial de crecimiento de poblaciones alrededor del a˜ no 1798. La gr´ afica muestra una curva de crecimiento exponencial.

1. Introducci´ on 1.1 El modelo de crecimiento de poblaci´ on de Malthus

Thomas Robert Malthus (1766–1834) fue uno de los primeros investi´ propuso, alrededor gadores en estudiar la din´amica de poblaciones. El de 1798, un modelo matem´atico de crecimiento de poblaciones basado en la idea de que ((la tasa per c´apita de crecimiento de una poblaci´on es directamente proporcional a su tama˜ no)). Su modelo, aunque simple, ha sido la base de muchos modelos de crecimiento de poblaciones biol´ogicas. El modelo de Malthus se puede describir en t´erminos del tiempo como variable continua por la ecuaci´on diferencial 1 dN =b−d=r (1) N dt donde N = N (t) denota la densidad o tama˜ no de la poblaci´on de una especie al tiempo t y los par´ametros b > 0 y d > 0 son las tasas per c´apita de nacimiento y muerte respectivamente, las cuales suponemos que son constantes, de forma que r ∈ R es la tasa per c´apita neta de crecimiento de la poblaci´on. De la ecuaci´on (1) obtenemos Zt

t=0

dN = N

Zt

rdt

=⇒

N (t) = N0 ert ,

(2)

t=0

donde N0 = N (0) es el tama˜ no inicial de la poblaci´on (al tiempo t = 0). Dada la forma de la soluci´on, el modelo de Malthus es conocido como ((crecimiento exponencial)) (ver figura 1). El comportamiento de la poblaci´on ser´a distinto seg´ un los valores de las tasas per c´apita de nacimiento y muerte: 1. Si la tasa per c´apita de nacimiento es mayor que la tasa per c´apita de muerte, i. e. b > d, entonces r > 0 y por lo tanto N (t) → ∞ si

79

´ MODELOS DINAMICOS DE POBLACIONES SIMPLES . . .

N (t) N (t) = N0

N0 Kert N0 ert + K N0

r > 0, K > 0

K=

N (t) = K

r c

N0 0

t

Figura 2. Pierre Francois Verhulst (1804-1849) modific´ o el modelo de Malthus en 1938, introduciendo un t´ermino inhibidor. La gr´ afica muestra las curvas de crecimiento del modelo log´ıstico.

t → ∞, es decir, el tama˜ no de la poblaci´on crecer´a arbitrariamente al pasar el tiempo (en forma exponencial). 2. Si las tasas per c´apita de nacimiento y muerte coinciden, i. e. b = d, entonces r = 0 y por lo tanto N (t) = N (0) en todo momento t > 0, es decir, el tama˜ no de la poblaci´on permanecer´a constante a lo largo del tiempo. 3. Si la tasa per c´apita de nacimiento es menor que la tasa per c´apita de muerte, i. e. b < d, entonces r < 0 y por lo tanto N (t) → 0 si t → ∞, es decir, el tama˜ no de la poblaci´on ir´a disminuyendo arbitrariamente al pasar el tiempo, de forma que en este caso la poblaci´on tiende hacia la extinci´on. Los casos 2 y 3 son realistas, sin embargo el caso 1 no lo es ya que al crecer la poblaci´on debe aumentar la demanda de recursos como comida, agua, etc., sin embargo, en la pr´actica, los recursos son limitados y no pueden aumentar arbitrariamente y a la par con un crecimiento arbitrario de la poblaci´on. 1.2 El modelo log´ıstico de crecimiento de poblaciones

Cuarenta a˜ nos m´as tarde, en 1838, el matem´atico belga Pierre Francois Verhulst (1804-1849) modific´o el caso 1 al introducir un t´ermino inhibidor en la ecuaci´on (1), el cual es proporcional al tama˜ no de la poblaci´on N (t), para as´ı tomar en cuenta la competencia entre los miembros de la poblaci´on por los recursos (agua, comida, etc.), los cuales no estar´an disponibles en cantidades ilimitadas. El modelo de Verhulst se puede describir en t´erminos del tiempo como variable continua por la ecuaci´on diferencial 1 dN = r − cN (r > 0, c > 0) (3) N dt

80

´ G. SAMANTA Y R. GOMEZ A´IZA

donde nuevamente N = N (t) denota la densidad o tama˜ no de la poblaci´on de una especie al tiempo t y el t´ermino cN corresponde a la competencia intraespec´ıfica de la poblaci´on (c > 0 es una constante). Como en el modelo de Malthus, el modelo de Verhulst incluye una tasa per c´apita de crecimiento r > 0 (positiva porque se enmarca en el caso 1 descrito anteriormente) que representa la tasa per c´apita de crecimiento neta a la que la poblaci´on crecer´ıa en ausencia de competencia intraespec´ıfica (i. e. si c = 0). El modelo de Verhulst se conoce tambi´en como el ((modelo log´ıstico de crecimiento de poblaciones)). Separando variables e integrando, obtenemos, despu´es de algunas simplificaciones, la soluci´on de la ecuaci´on (3): N (t) =

N0 Kert = N0 ert + K − N0

1+



K

 , K −rt −1 e N0

r donde K = , de forma que N (t) → K si t → ∞ (ya que r > 0). c La ecuaci´on log´ıstica y su soluci´on ocurren en diversos campos. La funci´on log´ıstica (i. e. N (t)) se conoce tambi´en como funci´on sigmoide y su gr´afica es una curva sigmoidea conocida como curva S (ver figura 2). La constante K se conoce como la capacidad de carga del h´abitat de la especie. Desde el punto de vista biol´ogico, la caracter´ıstica que falta en el modelo de Malthus es el concepto de capacidad de carga. Al crecer el tama˜ no de una poblaci´on disminuye la capacidad del medio ambiente de albergarla debido a que la disponibilidad per c´apita de la comida disminuye, los desperdicios aumentan y se acumulan, y en consecuencia las tasas de nacimiento tienden a decrecer mientras que las tasas de muerte tienden a aumentar. La capacidad de carga es el nivel de poblaci´on en el que las tasas de nacimiento y muerte coinciden, resultando en un tama˜ no poblacional estable a lo largo del tiempo. El modelo log´ıstico de crecimiento de poblaci´on (3) tambi´en puede escribirse como   N dN = rN 1 − (r > 0, K > 0). (4) dt K

La ley log´ıstica describe muy bien el crecimiento de una colonia de bacterias en un medio nutriente. Esta ley tambi´en ha sido utilizada para ajustar poblaciones humanas. Las siguientes suposiciones son intr´ınsecas al modelo log´ıstico: • La poblaci´on est´a uniformemente distribuida en el h´abitat. Entonces la densidad de poblaci´on (i. e. el tama˜ no de la poblaci´on por unidad de a´rea) no var´ıa con la posici´on pero s´ı con el tiempo. • No hay migraci´on. La poblaci´on cambia u ´nicamente debido al nacimiento y muerte de los individuos.

´ MODELOS DINAMICOS DE POBLACIONES SIMPLES . . .

81

Observemos tambi´en lo siguiente: • Si N0 6= K, entonces N (t) no puede cruzar la l´ınea recta N = K. • Si N0 = K = rc , entonces la soluci´on es la funci´on N (t) = N0 = K para toda t (t´omese N = K en (3) o en (4)). • Desde el punto de vista biol´ogico, poblaciones negativas no tienen sentido. Definici´ on 1.1. Consideremos una ecuaci´on diferencial aut´onoma, ordinaria de primer orden dN = f (N ). (5) dt 1. Una soluci´on constante N (t) = C de (5) se llama punto de equilibrio (o estacionario). Un punto de equilibrio corresponde a dN = dt f (N ) = 0. 2. Un punto de equilibrio C es estable si N (t) → C si t → ∞, para toda condici´on inicial N (0). 3. Si un punto de equilibrio no es estable, entonces es inestable. El sistema descrito por la ecuaci´on diferencial (3) (o (4)) tiene dos puntos de equilibrio: • N = 0 es un punto de equilibrio inestable. • N = K = rc es un punto de equilibrio estable.

2. El modelo depredador-presa de Lotka-Volterra En general una poblaci´on no vive aislada y las especies tienen ((enemigos)): algunas se alimentan de otras; hay depredadores y presas. Uno de los modelos que incorporan interacciones entre depredadores y presas fue propuesto en 1925 por el biof´ısico americano Alfred Lotka y en forma independiente por el matem´atico italiano Vito Volterra. El modelo de Lotka-Volterra describe las interacciones entre dos especies en un ecosistema: una poblaci´on que consiste de presas (e. g. conejos) y una poblaci´on de que consiste de depredadores (e. g. zorros) y supone las siguientes hip´otesis: 1. Los conejos tienen una cantidad ilimitada de comida y por lo tanto, en ausencia de zorros, la densidad de conejos N1 = N1 (t) al tiempo 1 t crecer´a exponencialmente (ley de Malthus), es decir, dN = r1 N1 , dt con r1 > 0. 2. En presencia de zorros, N1 (t) es decreciente debido al consumo de los zorros a una tasa proporcional a N1 N2 , donde N2 = N2 (t) 1 = denota la densidad de zorros al tiempo t. Por lo tanto dN dt r1 N1 − b1 N1 N2 , donde r1 > 0 y b1 > 0 son constantes. Aqu´ı b1 N1 representa la densidad de especies presa consumidas por especies

82

´ G. SAMANTA Y R. GOMEZ A´IZA

depredador en una unidad de tiempo y se conoce como la respuesta funcional del depredador (o tambi´en como la respuesta funcional de Volterra o la respuesta Holling tipo I ) sobre la presa. 3. En ausencia de conejos, la densidad de zorros N2 (t) decrecer´a ex2 ponencialmente a cero (extinci´on), i. e. dN = −r2 N2 , donde r2 > 0 dt es una constante. Esto es porque si la poblaci´on de zorros se queda aislada, es decir sin presas para su alimento (los conejos), entonces morir´an m´as r´apido de lo que podr´an reproducirse, y experimentar´an un decrecimiento poblacional exponencial. 4. En presencia de conejos, la densidad de zorros N2 (t) aumentar´a con una tasa de crecimiento proporcional a N1 N2 , es decir, proporcional al n´ umero de encuentros entre zorros y conejos. Junto con el 2 inciso anterior 3 obtenemos dN = −r2 N2 + b2 N1 N2 , con r2 > 0 dt b2 y b2 > 0 constantes. Aqu´ı b1 es llamado el factor de conversi´on y claramente 0 < bb12 < 1, de forma que 0 < b2 < b1 . Resumiendo, obtenemos el siguiente modelo de Lotka-Volterra para sistemas de poblaciones depredador-presa:   r1 dN1 = r1 N1 − b1 N1 N2 = b1 N1 − N2 dt b1 (6)   dN2 r2 = −r2 N2 + b2 N1 N2 = b2 N2 − + N1 dt b2

El sistema anterior est´a conformado por un par de ecuaciones diferenciales cuadr´aticas y no posee una soluci´on anal´ıtica general. Sin embargo, dados los par´ametros r1 , r2 , b1 y b2 as´ı como las condiciones iniciales N1 (0) = N10 y N2 (0) = N20 , podemos realizar una integraci´on num´erica para aproximar la soluci´on. El modelo de Lotka-Volterra tambi´en es un modelo de retroalimentaci´on ya que la poblaci´on de presas tiene un efecto positivo en el tama˜ no de la poblaci´on de depredadores, mientras que esta u ´ltima tiene un efecto negativo (inhibidor) en el tama˜ no de la poblaci´on de presas. 2.1 Trayectorias

Del sistema (6) obtenemos dN2 N2 (b2 N1 − r2 ) = . dN1 N1 (r1 − b1 N2 ) Hacemos separaci´on de variables e integramos para as´ı obtener Zt Zt 1 1 (r1 − b1 N2 )dN2 = (b2 N1 − r2 )dN1 . N2 N1 t=0

t=0

(7)

83

´ MODELOS DINAMICOS DE POBLACIONES SIMPLES . . .

N2

N1 =

r2 b2

II N˙ 1 < 0, N˙ 2 < 0

I

N˙ 1 < 0, N˙ 2 > 0

A2

P0

A3

P

r2 r1 , b2 b1

⇥

N2 =

r1 b1

A1

A4 III N˙ 1 > 0, N˙ 2 < 0

IV

N˙ 1 > 0, N˙ 2 > 0

N1

0

Figura 3. Alfred Lotka (1880-1949, parte superior) y Vito Volterra (18601940, parte inferior) desarrollaron en 1925, de manera independiente, el modelo depredador-presa. En la gr´ afica se ilustran las trayectorias.

Entonces r1 log

N1 N2 − b1 (N2 − N20 ) = b2 (N1 − N10 ) − r2 log N20 N10

y as´ı N2r1 N1r2 b1 (N2 −N20 ) b2 (N1 −N10 ) e , r1 r2 = e N20 N10 donde N1 (0) = N10 y N2 (0) = N20 son las poblaciones iniciales de presas y depredadores respectivamente. A trav´es de cada punto P0 ≡ (N10 , N20 ) del primer cuadrante, obtenemos una curva de la familia e−b2 N1 N1r2 = Φeb1 N2 N2−r1 , donde Φ =

r1 r2 N20 N10 eb1 N20 +b2 N10

(8)

es una constante. Estas curvas en el plano N1 N2

son las trayectorias del modelo de Lotka-Volterra (6) y son las que se ilustran en la figura 3. El plano N1 N2 se llama plano fase. Para describir la forma de las trayectorias consideremos los signos de las derivadas 2 1 y N˙ 2 = dN en las cuatro regiones I, II, III y IV en las que N˙ 1 = dN dt dt el primer cuadrante se divide al considerar las l´ıneas rectas N1 = rb22 y N2 = rb11 . Los signos de N˙ 1 y N˙ 2 en estas cuatro regiones los obtenemos de (6). Supongamos primero que el punto inicial P0 ≡ (N10 , N20 ) se encuentra en la regi´on I. En esta regi´on N1 es decreciente y N2 es creciente, de forma que el punto se mover´a en direcci´on contraria a las manecillas del reloj al transcurrir el tiempo hasta que alcance el punto

84

´ G. SAMANTA Y R. GOMEZ A´IZA

dN2 A2 , donde N˙ 2 = dN = 0, por lo que la tangente a la trayectoria en A2 1 es paralela al eje N1 . Una vez en la regi´on II, tanto N1 como N2 son decrecientes y el punto contin´ ua movi´endose en direcci´on contraria a dN2 las manecillas del reloj hasta que alcanza el punto A3 donde N˙ 1 = dN 1 est´a indefinida y por lo tanto la tangente a la trayectoria en A3 es paralela al eje N2 . Similarmente, en la regi´on III, N1 es creciente mientras que N2 es decreciente hasta el punto A4 , y en la regi´on IV tanto N1 como N2 son crecientes hasta el punto A1 . Todas las trayectorias son entonces descritas por un movimiento en direcci´on contraria a las manecillas del reloj y se acumulan alrededor de los ejes ya que se acercan a ellos pero nunca son incidentes. Podemos entonces confirmar que todas las soluciones de (6) van a ser ´orbitas peri´odicas cerradas alrededor del punto P ∗ ≡ rb22 , rb11 .

Observaci´ on 2.1. Cambiar las tasas per c´apita de nacimiento y muerte tiene u ´nicamente el efecto de cambiar el periodo de la oscilaci´on, i. e. ninguna de las dos poblaciones ser´a dominante y tampoco hay posibilidad de que ninguna de las dos se extinga, de forma que tanto las presas como los depredadores coexisten juntos.

Definici´ on 2.1. Consideremos un sistema aut´onomo de ecuaciones diferenciales de primer orden dN1 = f1 (N1 , N2 ) dt (9) dN2 = f2 (N1 , N2 ) dt 1. Un punto de equilibrio ocurre cuando N1 (t) y N2 (t) son constantes, 1 2 de forma que dN = 0 y dN = 0. En un punto de equilibrio decimos dt dt que el sistema est´a en reposo. Los puntos de equilibrio se conocen tambi´en como puntos estacionarios, o puntos singulares, o puntos cr´ıticos o puntos de reposo. 2. Un punto de equilibrio (N1∗ , N2∗ ) es localmente asint´oticamente estable si la respuesta a una peque˜ na perturbaci´on tiende a cero cuando el tiempo tiende a infinito, o m´as precisamente, si l´ım N1 (t) = N1∗

t→∞

y

l´ım N2 (t) = N2∗

t→∞

siempre y cuando (N1 (t), N2 (t)) comience en valores iniciales P0 ≡ (N1 (0), N2 (0)) ≡ (N10 , N20 ) suficientemente cercanos al punto de equilibrio (N1∗ , N2∗ ). Los puntos de equilibrio localmente asint´oticamente estables se conocen tambi´en como pozos o atractores. Al conjunto de todos los valores iniciales (N10 , N20 ) cuyas trayectorias convergen a un atractor se le llama el dominio de atracci´on.

´ MODELOS DINAMICOS DE POBLACIONES SIMPLES . . .

85

3. Un punto de equilibrio (N1∗ , N2∗ ) es estable, pero no localmente asint´oticamente estable, si la respuesta a una peque˜ na perturbaci´on es peque˜ na pero no tiende a cero cuando el tiempo tiende a infinito. En este caso, el punto de equilibrio se llama centro. En este caso la soluci´on regresar´a a la condici´on inicial en forma peri´odica, de tal manera que su trayectoria es una curva cerrada alrededor del punto de equilibrio, por lo que es llamada ciclo. 4. Un punto de equilibrio es inestable si condiciones iniciales (cerca del punto de equilibrio) producen soluciones que se alejan del punto de equilibrio al transcurrir el tiempo. Cerca de un punto de equilibrio inestable, peque˜ nas perturbaciones a condiciones iniciales tienden a afectar dr´asticamente las soluciones. Los punto de equilibrio inestables se conocen tambi´en como fuentes o repulsores. Observaci´ on 2.2. El sistema (6) tiene dos puntos de equilibrio, a saber, 0 ≡ (0, 0)

∗

P ≡

y



r2 r1 , b 2 b1



.

De la discusi´on anterior deducimos que el punto de equilibrio 0 del sistema (6) es inestable y que el punto de equilibrio P ∗ es estable (pero no localmente asint´oticamente) y es un centro.

3. Ecuaci´ on secular (o caracter´ıstica) y estabilidad local En esta secci´on abordaremos el problema de determinar por medio de la ecuaci´on secular la estabilidad local de un sistema. Consideremos nuevamente el modelo aut´onomo (9) y supongamos que tiene un punto de equilibrio en (N 1 , N 2 ), de forma que f1 (N 1 , N 2 ) = 0

y

f2 (N 1 , N 2 ) = 0.

(10)

Para analizar la estabilidad lineal del modelo (9) en (N 1 , N 2 ) sustituimos N1 (t) = N 1 +u1 (t) y N2 (t) = N 2 +u2 (t) en (9), donde 0 < ui (t)  1 (con i = 1, 2), considerando as´ı una peque˜ na perturbaci´on. Usando la

86

´ G. SAMANTA Y R. GOMEZ A´IZA

expansi´on en serie de Taylor obtenemos
du1 = f1 N 1 + u1 , N 2 + u2 dt = f1 N 1 , N 2

1 + 2!






  ∂ ∂ + u2 f1 (N1 , N2 )|(N 1 ,N 2 ) + u1 ∂N1 ∂N2

∂ ∂ u1 + u2 ∂N1 ∂N2

2

f1 (N1 , N2 )|(N 1 ,N 2 ) + . . .

y
du2 = f2 N 1 + u1 , N 2 + u2 dt = f2 N 1 , N 2

1 + 2!






  ∂ ∂ + u2 f2 (N1 , N2 )|(N 1 ,N 2 ) + u1 ∂N1 ∂N2

∂ ∂ + u2 u1 ∂N1 ∂N2

2

f2 (N1 , N2 )|(N 1 ,N 2 ) + . . . .

Usando (10) y linealizando (ya que u1 (t) y u2 (t) son muy peque˜ nas), obtenemos du1 ∂f1 ∂f1 = u1 + u2 dt ∂N 1 ∂N 2 donde

∂fi ∂N j

du2 ∂f2 ∂f2 = u1 + u2 dt ∂N 1 ∂N 2

y

(con i, j = 1, 2) denota los valores de

∂fi ∂Nj

(11)

en el punto de

equilibrio (N 1 , N 2 ). Consideraremos la soluci´on en la forma u1 = A1 eλt

y

u2 = A2 eλt .

(12)

Observamos que u1 y u2 no se anulan simult´aneamente ya que hemos realmente generado una perturbaci´on. Sustituyendo u1 y u2 en (12) en (11) y simplificando llegamos a   ∂f1 ∂f1 A1 − A2 = 0 λ− ∂N 1 ∂N 2   ∂f2 ∂f2 − A1 + λ − A2 = 0, ∂N 1 ∂N 2

87

´ MODELOS DINAMICOS DE POBLACIONES SIMPLES . . .

P

P

(N 1 , N 2 )

P

(a) Nodo estable

P

(N 1 , N 2 )

(b) Nodo inestable

Figura 4. Cuando las ra´ıces de la ecuaci´ on caracter´ıstica son reales y tienen el mismo signo (caso 1), entonces el punto de equilibrio es estable o inestable, ya sea que el signo de las ra´ıces sea negativo o positivo respectivamente.

o bien, en forma matricial,  ∂f1 ∂f1 − λ−  ∂N 1 ∂N 2    ∂f2 ∂f2 λ− − ∂N 1 ∂N 2



    0  A1   =  .   A 0 2

Ya que u1 y u2 no se anulan simult´aneamente, lo mismo ocurre con A1 y A2 , de forma que ∂f1 ∂f1 − λ− ∂N 1 ∂N 2 =0 ∂f ∂f 2 2 − λ− ∂N 1 ∂N 2 y por lo tanto     ∂f1 ∂f2 ∂f1 ∂f2 ∂f1 ∂f2 2 λ − + λ+ − = 0. ∂N 1 ∂N 2 ∂N 1 ∂N 2 ∂N 2 ∂N 1

(13)

Esta ecuaci´on se conoce como la ecuaci´on secular (o caracter´ıstica). Es una ecuaci´on cuadr´atica en λ y por lo tanto tiene dos soluciones λ1 y λ2 . La soluci´on general la podemos entonces escribir como u1 (t) = a1 eλ1 t + a2 eλ2 t (14) u2 (t) = a3 eλ1 t + a4 eλ2 t

88

´ G. SAMANTA Y R. GOMEZ A´IZA

P

P

(a) Punto silla

N

(N 1 , N 2 )

(N 1 N 2 )

N

(b) Centro

Figura 5. (a) Si las ra´ıces de la ecuaci´ on caracter´ıstica son reales y tienen signos contrarios (caso 2), entonces el punto de equilibrio es un punto silla. (b) Si las ra´ıces de la ecuaci´ on caracter´ıstica son imaginarias (caso 4), entonces el punto de equilibrio es un centro.

b (con i = 1, . . . , 4). Usando N1 (t) = N 1 + u1 (t) y donde ai ∈ ER N2 (t) = N 2 + u2 (t) junto con (14), obtenemos los siguientes casos:

Caso 1: Nodo. Las ra´ıces λ1 y λ2 de la ecuaci´on caracter´ıstica (13) son reales y tienen el mismo signo. En este caso el punto de equilibrio (N 1 , N 2 ) se clasifica como nodo (o punto nodal ). Hay dos subcasos: • Estable. Si λ1 , λ2 < 0, entonces el punto de equilibrio (N 1 , N 2 ) es un nodo estable. En este caso u1 (t) → 0 y u2 (t) → 0 cuando t → ∞ y por lo tanto N1 (t) → N 1 y N2 (t) → N 2 cuando t → ∞ (ver figura 4.a). • Inestable. Si λ1 , λ2 > 0, entonces el punto de equilibrio (N 1 , N 2 ) es un nodo inestable. En este caso u1 (t) → ∞ y u2 (t) → ∞ cuando t → ∞ y por lo tanto N1 (t) → ∞ y N2 (t) → ∞ cuando t → ∞ (ver figura 4.b). Caso 2: Punto silla. Las ra´ıces λ1 y λ2 de la ecuaci´on caracter´ıstica (13) son reales y de signos diferentes. En este caso el punto de equilibrio (N 1 , N 2 ) se clasifica como punto silla. Ahora u1 (t) → ∞ y u2 (t) → ∞ cuando t → ∞ y por lo tanto N1 (t) → ∞ y N2 (t) → ∞ cuando t → ∞. De forma que (N 1 , N 2 ) es inestable (ver figura 5.a). Caso 3: Foco. Las ra´ıces λ1 y λ2 de la ecuaci´on caracter´ıstica (13) son n´ umeros complejos conjugados con parte real no nula. En este caso el punto de equilibrio (N 1 , N 2 ) se clasifica como foco (o punto focal ). Cerca del punto focal la las trayectorias que corresponden a las soluciones se comportan como espirales.

89

´ MODELOS DINAMICOS DE POBLACIONES SIMPLES . . .

N

N

(N 1 N 2 )

(N 1 N 2 )

N

N

(a) Foco estable

(b) Foco inestable

Figura 6. Cuando las ra´ıces de la ecuaci´ on caracter´ıstica son n´ umeros complejos conjugados con parte real no nula (Caso 3), entonces el punto de equilibrio es un foco. (a) Si las partes reales de las ra´ıces son negativas, entonces el foco es estable. (b) Si las partes reales de las ra´ıces son positivas, entonces el foco es inestable.

• Estable. Si las partes reales de λ1 y λ2 son ambas negativas, entonces el punto focal (N 1 , N 2 ) es estable y las trayectorias de las soluciones son espirales que tienden hacia el punto de equilibrio (N 1 , N 2 ) cuando t → ∞ (ver figura 6.a). • Inestable. Si las partes reales de λ1 y λ2 son ambas positivas, entonces el punto focal (N 1 , N 2 ) es inestable y las trayectorias de las soluciones son espirales que tienden a alejarse del punto de equilibrio (N 1 , N 2 ) cuando t → ∞ (ver figura 6.b). Caso 4: Centro. Las ra´ıces λ1 y λ2 de la ecuaci´on caracter´ıstica (13) son n´ umeros imaginarios. En este caso las trayectorias de las soluciones alrededor del punto de equilibrio (N 1 , N 2 ) son curvas cerradas que corresponden a soluciones peri´odicas de la ecuaci´on diferencial lineal (11). En este caso el punto de equilibrio (N 1 , N 2 ) se clasifica como centro (o punto v´ortice) y se ilustra en la figura 5.b. Observemos que el centro es un punto de equilibrio estable, pero no localmente asint´oticamente estable.

4. Estabilidad local (o lineal) del modelo de LotkaVolterra Sabemos que el modelo de Lotka-Volterra  (6) tiene dos puntos de equi r r librio, a saber 0 ≡ (0, 0) y P ∗ ≡ b22 , b11 . Sean f1 (N1 , N2 ) = r1 N1 −b1 N1 N2

y

f2 (N1 , N2 ) = −r2 N2 +b2 N1 N2 (15)

90

´ G. SAMANTA Y R. GOMEZ A´IZA

donde ri , bi > 0 (con i = 1, 2) son constantes positivas. Entonces ∂f1 = r1 − b1 N2 ∂N1

∂f1 = −b1 N1 ∂N2

∂f2 = b2 N2 ∂N1

∂f2 = −r2 + b2 N1 . ∂N2

(16)

Analicemos ahora los dos puntos de equilibrio susituyendolos en (16). I. Punto de equilibrio (0, 0). En este caso tenemos ∂f1 ∂f1 ∂f2 ∂f2 = r1 , = 0, = 0, = −r2 . (17) ∂N1 (0,0) ∂N2 (0,0) ∂N1 (0,0) ∂N2 (0,0)

Por lo tanto la ecuaci´on caracter´ıstica es (usando (13)) λ2 − (r1 − r2 )λ − r1 r2 = 0

=⇒ =⇒

(λ − r1 )(λ + r2 ) = 0 λ = r1 , −r2 .

Observamos entonces que las ra´ıces de la ecuaci´on caracter´ıstica son reales y tienen signos opuestos (caso 2), y por lo tanto el punto de equilibrio (0, 0) es inestable, es un  punto  silla. II. Punto de equilibrio P ∗ ≡ rb22 , rb21 . En este caso tenemos

∂f1 b1 r2 ∂f2 b2 r1 ∂f2 ∂f1 = 0, =− , = , = 0. (18) ∂N1 P ∗ ∂N2 P ∗ b2 ∂N1 P ∗ b1 ∂N2 P ∗

Por lo tanto la ecuaci´on caracter´ıstica es (usando (13)) √ √ λ2 + r1 r2 = 0 =⇒ λ = ±i r1 r2 (donde i = −1).

Entonces el punto de equilibrio P ∗ es estable pero no asint´oticamente estable, es un centro (caso 4). Concluimos de la discusi´on anterior que el punto de equilibrio ((interno)), o de coexistencia, del modelo de Lotka-Volterra (6), es decir P ∗ ≡  r2 r1 , , no es asint´oticamente estable. La ausencia de estabilidad asint´otib2 b1 ca de este punto de equilibrio indica que el sistema de Lotka-Volterra no posee un mecanismo para mantener un estado de coexistencia estable. Desde el punto de vista ecol´ogico, la causa detr´as de ´esto es la ausencia del concepto de capacidad de carga del h´abitat para la especie de presas (es decir, la ausencia de competencia intraespec´ıfica dentro en el h´abitat para la especie de presas). Desde el  punto  de vista de la r2 r1 teor´ıa de la estabilidad, el estado estable no trivial b2 , b1 es un estado de equilibrio neutral.

91

´ MODELOS DINAMICOS DE POBLACIONES SIMPLES . . .

5. Estabilizaci´ on de sistemas depredador-presa introduciendo competencia intraespec´ıfica en la especie de presas Introducir competencia intraespec´ıfica en la especie de presas dentro del modelo de Lotka-Volterra nos lleva al siguiente sistema: dN1 = r1 N1 − βN1 N2 − γN12 dt

(19)

dN2 = −r2 N2 + κβN1 N2 dt donde el t´ermino γN12 corresponde a la competencia intraespec´ıfica en el h´abitat para la especie de presas. Es entonces claro que en ausencia de depredadores, el valor l´ımite de la poblaci´on de la especie de presas e1 = r1 . Veremos ahora que el sistema (19) tiene un u es N ´nico punto γ r 2 de equilibrio interno N ∗ ≡ (N1∗ , N2∗ ), donde N1∗ = κβ2 y N2∗ = r1 κβ−γr . κβ 2 r1 r2 ∗ e El suponer la condici´on natural N1 > N , es decir, > , nos lleva 1

γ

κβ

a que N2∗ > 0 y por lo tanto existir´a el punto de equilibrio (N1∗ , N2∗ ) de coexistencia (o interno). (19) tiene otros dos puntos de  El sistema  r1 equilibrio, a saber (0, 0) y γ , 0 , mismos que tambi´en analizaremos a continuaci´on. El modelo (19) es conocido como el modelo log´ıstico de Lotka-Volterra de un sistema depredador-presa, o tambi´en como modelo amortiguado de Lotka-Volterra de un sistema depredador-presa.

5.1 An´ alisis de estabilidad lineal del modelo log´ıstico de LotkaVolterra

Para el an´alisis de la estabilidad lineal consideremos g1 (N1 , N2 ) = r1 N1 −βN1 N2 −γN12

y g2 (N1 , N2 ) = −r2 N2 +κβN1 N2 ,

de forma que ∂g1 = r1 − βN2 − 2γN1 , ∂N1 ∂g2 = κβN2 , ∂N1

∂g1 = −βN1 , ∂N2 ∂g2 = −r2 + κβN1 . ∂N2

I. Punto de equilibrio (0, 0). En este caso ∂g1 ∂g1 ∂g2 = r1 , = 0, = 0, ∂N1 (0,0) ∂N2 (0,0) ∂N1 (0,0)

∂g2 = −r2 , ∂N2 (0,0)

92

´ G. SAMANTA Y R. GOMEZ A´IZA

de forma que la ecuaci´on caracter´ıstica es (usando la ecuaci´on (13)) λ2 − (r1 − r2 )λ − r1 r2 = 0

=⇒

λ = r1 , −r2 .

As´ı, el punto de equilibrio (0, 0) es inestable, es un punto silla (caso 2).  r1 e II. Punto de equilibrio N ≡ γ , 0 . En este caso ∂g1 ∂g1 r1 β ∂g2 ∂g2 κβr1 − r2 γ = −r1 , =− , = 0, = , ∂N1 Ne ∂N2 Ne γ ∂N1 Ne ∂N2 Ne γ

de forma que la ecuaci´on caracter´ıstica es (usando la ecuaci´on (13))   κβr1 − r2 γ r1 (κβr1 − r2 γ) 2 λ − −r1 + λ− =0 γ γ y por lo tanto λ = −r1 (< 0),

κβr1 − r2 γ (> 0). γ

  As´ı, el punto de equilibrio rγ1 , 0 es inestable, es un punto silla (caso 2). III. Punto de equilibrio N ∗ ≡ (N1∗ , N2∗ ). En este caso ∂g1 ∂g1 ∗ = −γN1 , = −βN1∗ , ∂N1 N ∗ ∂N2 N ∗ ∂g2 ∂g2 ∗ = κβN2 , = 0, ∂N1 N ∗ ∂N2 N ∗ de forma que la ecuaci´on caracter´ıstica es (usando la ecuaci´on (13)) λ2 + γN1∗ λ + κβ 2 N1∗ N2∗ = 0

y as´ı √  1 −γN1∗ ± ∆ , donde ∆ = γ 2 (N1∗ )2 − 4κβ 2 N1∗ N2∗ . 2 Por lo tanto, o bien las dos ra´ıces de la ecuaci´on secular son reales y negativas (caso 1, nodo estable), o son complejas conjugadas con parte real negativa (caso 3, foco estable). Por lo tanto el punto de equilibrio (N1∗ , N2∗ ) es localmente asint´oticamente estable. λ=

5.2 Estabilidad global

A´ un queda la siguiente pregunta: ¿toda trayectoria que inicia dentro del primer cuadrante positivo eventualmente llega al estado de equilibrio de coexistencia N ∗ ≡ (N1∗ , N2∗ )? Una respuesta afirmativa implicar´ıa la estabilidad asint´otica global de N ∗ . Para buscar la respuesta a esta pregunta enunciemos el siguiente:

´ MODELOS DINAMICOS DE POBLACIONES SIMPLES . . .

93

Teorema 5.1 (Estabilidad de Lyapunov). Consideremos el siguiente sistema de ecuaciones diferenciales: dxi = fi (x1 , x2 , . . . , xn ) , i = 1, 2, . . . , n dt

(20)

y supongamos que tiene un punto de equilibrio en (x1 , x2 , . . . , xn ). Supongamos adem´as que existe una funci´on diferenciable ν(x1 , x2 , . . . , xn ) que satisfaga las siguientes condiciones: 1. ν(x1 , x2 , . . . , xn ) tiene un m´ınimo estricto en (x1 , x2 , . . . , xn ), es decir, ν > 0 y ν = 0 para xi = xi , i = 1, 2, . . . , n. 2. La derivada de ν a lo largo de las curvas integrales de (20) satisface dν X ∂ν dxi X ∂ν fi ≤ 0 = = dt ∂x dt ∂x i i i=1 i=1 n

n

y fuera de una vecindad arbitrariamente peque˜ na de (x1 , x2 , . . . , dν xn ) ocurre que dt < 0. Entonces el punto de equilibrio (x1 , x2 , . . . , xn ) es globalmente asint´oticamente estable. La funci´on ν(x1 , x2 , . . . , xn ) se conoce como la funci´on de Lyapunov. Ahora enunciemos y probemos una desigualdad que nos ser´a de utilidad. Proposici´ on 5.2. z − log(z) − 1 ≥ 0 para toda z > 0 y la igualdad ocurre si y solo si z = 1. Demostraci´on. Sea w = z − log(z) − 1, de forma que dw = 1 − z1 . dz dw Si 0 < z < 1, entonces dz < 0, por lo tanto w es decreciente y en consecuencia w = z − log z − 1 > 1 − log 1 − 1 = 0. Si z > 1, entonces dw > 0, por lo tanto w es creciente y en consecuencia w = z − log(z) − dz 1 > 1 − log(1) − 1 = 0. En virtud del teorema de Lyapunov y la desigualdad de la proposici´on 5.2, definimos la funci´on    N1 N1 N2∗ N2 N2 L(N1 , N2 ) = − log ∗ − 1 + − log ∗ − 1 , N1∗ N1 κ N2∗ N2 (21) de forma que L(N1 , N2 ) ≥ 0 en el cuadrante positivo del plano N1 N2 y L(N1 , N2 ) = 0 u ´nicamente en (N1∗ , N2∗ ). Tambi´en tenemos que, por N1∗



94

´ G. SAMANTA Y R. GOMEZ A´IZA

(19), dL ∂L dN1 ∂L dN2 = + dt ∂N1 dt ∂N2 dt  
N1∗ r1 N1 − βN1 N2 − γN12 = 1− N1 1 + κ y por lo tanto



N∗ 1− 2 N2



(−r2 N2 + κβN1 N2 )

N2 dL = −γN12 + (κβN1∗ − r2 ) dt κ r2 + (r1 + γN1∗ − βN2∗ ) N1 − r1 N1∗ + N2∗ . (22) κ Ya que r2 r1 − βN2∗ − γN1∗ = 0, −r2 + κβN1∗ = 0, −r1 N1∗ + N2∗ = −γ(N1∗ )2 , κ de (22) obtenemos dL = −γ (N1 − N1∗ )2 ≤ 0 dt y fuera de una vecindad arbitrariamente peque˜ na de (N1∗ , N2∗ ) tenemos dL < 0. Entonces la funci´on L(N1 , N2 ) dada por (21) es la funci´on de dt Lyapunov y el estado de equilibrio de coexistencia (N1∗ , N2∗ ) es globalmente asint´oticamente estable.

6. Respuesta funcional de los depredadores en las presas De la primera ecuaci´on en (19) observamos que el n´ umero de presas ((devoradas)) por los depredadores por unidad de tiempo es βN1 . Aqu´ı φ(N1 ) = βN1 se conoce como la respuesta funcional de los depredadores sobre las presas, o tambi´en respuesta funcional de Volterra, o bien respuesta funcional de Holling tipo I. As´ı, la respuesta funcional de Holling tipo I supone que los depredadores est´an siempre hambrientos y que un depredador comer´a m´as presas si estas u ´ltimas se encuentran disponibles. Una respuesta funcional de depredadores realista tomar´a en cuenta el que los depredadores queden satisfechos. Despu´es de todo, ¿cu´antos conejos se puede comer un zorro? Al aumentar la poblaci´on de presas,

´ MODELOS DINAMICOS DE POBLACIONES SIMPLES . . .

(a) Buscando la presa

95

(b) Manipulando la presa

Figura 7. La respuesta funcional del depredador sobre la presa est´ a determinada por dos acciones: (a) b´ usqueda y (b) manipulaci´ on de la presa.

el n´ umero de presas devoradas por un depredador aumenta en una tasa decreciente. Para obtener una respuesta matem´atica razonable, hagamos el siguiente an´alisis: • La depredaci´on involucra dos actividades: 1. Buscar la presa (ver figura 7.a). 2. Manipular la presa (ver figura 7.b). El tiempo de manipulaci´on se refiere al tiempo que le toma al zorro perseguir, capturar y comer al conejo. • El tiempo total disponible para la depredaci´on es T . • N1 es el n´ umero total de presas (i. e. el tama˜ no de la poblaci´on de presas). • V es el n´ umero de presas capturadas por un depredador por unidad de tiempo. • Th es el tiempo que toma un depredador en manipular a una presa. • T − V Th es el tiempo de b´ usqueda de presas que utiliza un depredador. • V es proporcional al n´ umero total de presas y al tiempo total de b´ usqueda, por lo que V = `N1 (T −V Th ), donde ` es una constante, y entonces `T N1 V = . 1 + `Th N1 As´ı, una respuesta funcional de los depredadores sobre las presas m´as realista es βN1 ψ(N1 ) = , α + N1 donde α, β > 0 son constantes positivas. La funci´on ψ se conoce como respuesta funcional de Holling tipo II y su gr´afica se ilustra en la figura 8. Ahora ψ(N1 ) → β cuando N1 → ∞ y entonces el consumo no puede ser arbitrariamente grande (en contraste con el tipo I). Observemos que

96

´ G. SAMANTA Y R. GOMEZ A´IZA

N1

0

Figura 8. Gr´ afica de la respuesta funcional de Holling tipo II.

ψ(α) = β2 . A α se le llama la constante de saturaci´on media. Notemos que si el tiempo total disponible para depredaci´on se usa u ´nicamente en b´ usqueda (i. e. Th = 0), entonces obtenemos la respuesta tipo I. Si incorporamos la respuesta funcional de Holling tipo II en el modelo log´ıstico de Lotka-Volterra, obtenemos   dN1 N1 βN1 N2 = rN1 1 − − dt K α + N1 (23) dN2 cβN1 N2 = −dN2 + . dt α + N1 Este modelo se conoce como sistema depredador-presa dependiente de presa.   1 1 Si en lugar de utilizar ψ(N1 ) usamos ψ N = αNβN como resN2 2 +N1 puesta funcional de los depredadores sobre las presas dentro del modelo log´ıstico de Lotka-Volterra, entonces obtenemos   dN1 N1 βN1 N2 = rN1 1 − − dt K αN2 + N1 (24) dN2 cβN1 N2 = −dN2 + , dt αN2 + N1 1 2 donde dN = dN = 0 en (N1 , N2 ) = (0, 0). Este modelo se conoce como dt dt sistema depredador–presa dependiente de raz´on. La respuesta funcional de Holling tipo III es una generalizaci´on de la tipo II, y es de la forma

ψ(N1 ) =

βN1k , α + N1k

la cual no se deduce f´acilmente de separar las acciones de b´ usqueda y manipulaci´on de presas. Su motivaci´on supone que existe un proceso de

´ MODELOS DINAMICOS DE POBLACIONES SIMPLES . . .

97

aprendizaje en los depredadores, dando como resultado un crecimiento de la tasa de descubrimiento de presas al incrementar la ocurrencia de encuentros previos de los depredadores con las presas, ya que a mayor densidad de presas, mayor el n´ umero de encuentros previos. Las respuestas funcionales de Holling tipos I, II y III son todas mon´otonas en el primer cuadrante. Por lo tanto, la tasa de consumo de presas por depredador se incrementa al aumentar la poblaci´on de presas. Sin embargo, se ha observado y hay evidencia experimental que indica que esta situaci´on no necesariamente ocurre, por ejemplo en el caso de defensa grupal en din´amica de poblaciones, t´ermino que se utiliza para describir el fen´omeno en el que el acto de depredaci´on disminuye o incluso desaparece debido a la creciente habilidad de las presas de defenderse u ocultarse cuando se encuentran agrupadas en grandes cantidades. Por ejemplo, es m´as complicado para un depredador de insectos identificar a un individuo particular en un enjambre. Para modelar este efecto, J. F. Andrews (1968) sugiri´o la funci´on ψ(N1 ) =

βN1 α + γN1 + N12

llamada respuesta funcional de Holling tipo IV.

7. Ciclos l´ımite y bifurcaci´ on de Hopf Consideremos el sistema de ecuaciones diferenciales dx = y − x(x2 + y 2 − 1) dt dy = −x − y(x2 + y 2 − 1) dt

(25)

con  un par´ametro escalar. Este sistema se deduce al hacer x = r cos θ y y = r sen θ, de forma que x˙ = r˙ cos θ − rθ˙ sen θ y y˙ = r˙ sen θ + rθ˙ cos θ. Sustituyendo en (25) obtenemos r˙ cos θ − rθ˙ sen θ =

r sen θ − r(r2 − 1) cos θ

(26)

r˙ sen θ + rθ˙ cos θ = −r cos θ − r(r2 − 1) sen θ

(27)

y observamos que (26) × cos θ + (27) × sen θ

⇒

(26) × (− sen θ) + (27) × cos θ

⇒

r˙ = −r(r2 − 1) rθ˙ = −r,

98

´ G. SAMANTA Y R. GOMEZ A´IZA

de forma que (25) en coordenadas polares es r˙ = −r(r2 − 1) θ˙ = −1.

(28)

Las funciones r(t) = 1 y θ(t) = −t son soluciones de (28) y por lo tanto el c´ırculo r = 1 es una ´orbita cerrada. Suponiendo que  > 0, si 0 < r < 1, entonces r˙ > 0, por lo tanto (0, 0) es inestable. Por otro lado, si r > 1, entonces r˙ < 0, de forma que la o´rbita cerrada r = 1, es decir x2 + y 2 = 1, es estable. Por lo tanto toda trayectoria que resuelva el sistema (25) tiende al c´ırculo x2 + y 2 = 1, el cual se denomina ciclo l´ımite del sistema (25). Definici´ on 7.1. Una trayectoria cerrada de un sistema din´amico es una ´orbita. Diremos que el movimiento a lo largo de una o´rbita es peri´odico. Un ciclo l´ımite (o atractor ) es una ´orbita tal que cualquier trayectoria que empiece en un punto cercano a ella converge a ella cuando t → ∞. Un ciclo origen (o repulsor ) es una o´rbita cuyas trayectorias vecinas divergen de ella. De acuerdo a la teor´ıa general de sistemas din´amicos, una o´rbita que no forme parte de una familia de ´orbitas conc´entricas debe de ser o bien un ciclo l´ımite o un ciclo origen. Un ciclo origen de un sistema din´amico invertible es un ciclo l´ımite del sistema din´amico a tiempo reversible, y viceversa. Como instancia notemos que si  < 0, entonces el c´ırculo unitario x2 + y 2 = 1 se vuelve un ciclo origen para el sistema (25). La ecuaci´on secular del sistema (25) es     ∂f1 ∂f2 ∂f1 ∂f2 ∂f1 ∂f2 2 λ− + − λ+ = 0, ∂x ∂y ∂x ∂y ∂y ∂x donde f1 (x, y) = y − x(x2 + y 2 − 1) por lo que ∂f1 = −(3x2 + y 2 − 1) ∂x ∂f1 = 1 − 2xy ∂y

y

f2 (x, y) = −x − y(x2 + y 2 − 1), ∂f2 = −(x2 + 3y 2 − 1) ∂y ∂f2 = −1 − 2xy. ∂x

Consideremos la posici´on de equilibrio (0, 0). En este punto la ecuaci´on secular es λ2 − 2λ + 2 + 1 = 0, o lo que es lo mismo, (λ − )2 = −1, de forma que λ =  ± i. Por lo tanto, para  = 0, el punto (0, 0) es un centro. Por otro lado, si  > 0, entonces (0, 0) es inestable y si  < 0, entonces (0, 0) es estable. En este caso decimos que una bifurcaci´on de Hopf ocurre en  = 0.

´ MODELOS DINAMICOS DE POBLACIONES SIMPLES . . .

99

Una bifurcaci´on de Hopf est´a caracterizada por el cambio en la estabilidad de un punto de equilibrio (o punto cr´ıtico) junto con la creaci´on de un ciclo l´ımite.

8. Bifurcaci´ on de Hopf en el sistema depredadorpresa dependiente de presa Podemos escribir al sistema (23) en la forma equivalente   dX X βXY = rX 1 − − , X(0) > 0, dt K α+X (29) dY cβXY = −dY + , Y (0) > 0, dt α+X donde X(t) y Y (t) denotan la densidad al tiempo t de presas y depredadores respectivamente, r, K y d son constantes positivas que representan la tasa per c´apita de crecimiento intr´ınseca de presas, la capacidad de carga de las presas, y la tasa per c´apita de muerte de los depredadores respectivamente, y α, β y c son constantes positivas que representan la constante de saturaci´on media de captura, la tasa de captura y la tasa de conversi´on respectivamente. 8.1 Equilibrio

El sistema (29) posee tres puntos de equilibrio: E0 = (0, 0) es trivial y E1 = (K, 0) es el u ´nico equilibrio axial. El punto de equilibrio en el interior es E ∗ = (X ∗ , Y ∗ ), donde X∗ =

αd cβ − d

y

Y∗ =

rcα(Kcβ − Kd − αd) , K(cβ − d)2

y cuya existencia ocurre cuando se satisface la condici´on K(cβ − d) > αd.

(30)

Esta desigualdad implica que el beneficio m´aximo para los depredadores a partir de la interacci´on con las presas (cβ) debe ser mayor que 1 + (α/K) veces la tasa per c´apita de muerte de los depredadores (d) para asegurar la existencia del punto de equilibrio interno E ∗ = (X ∗ , Y ∗ ), donde K y α son la capacidad de carga de las presas y la constante media de saturaci´on respectivamente. 8.2 An´ alisis de estabilidad y bifurcaci´ on

Procedamos ahora al an´alisis de la estabilidad de las ecuaciones diferenciales (29) que gobiernan la evoluci´on del sistema. Las matrices

100

´ G. SAMANTA Y R. GOMEZ A´IZA

variacionales que corresponden a E0 y E1 est´an dadas por   r 0 V (E0 ) = 0 −d ! βK −r − α+K y V (E1 ) = . 0 K(cβ−d)−αd α+K Claramente E0 es inestable y si E ∗ = (X ∗ , Y ∗ ) existe, entonces E1 es inestable. Ahora consideremos la estabilidad del m´as interesante de los puntos de equilibrio positivos (equilibrio interno), E ∗ = (X ∗ , Y ∗ ). La matriz variacional en E ∗ = (X ∗ , Y ∗ ) est´a dada por ! αβY ∗ βX ∗ 2r ∗ r − X − − ∗ 2 ∗ K (α+X ) α+X V (E ∗ ) = . (31) cαβY ∗ 0 ∗ 2 (α+X ) Teorema 8.1. Si existe E ∗ = (X ∗ , Y ∗ ) con K(cβ − d) < α(cβ + d), entonces E ∗ es localmente asint´oticamente estable. Demostraci´on. La ecuaci´on caracter´ıstica para la matriz variacional V (E ∗ ) (dada por (31)) es λ2 − A1 λ + A2 = 0,

donde

r ∗ βX ∗ Y ∗ cαβ 2 X ∗ Y ∗ X + y A = > 0. 2 K (α + X ∗ )2 (α + X ∗ )3 Claramente, si A1 < 0, entonces E ∗ = (X ∗ , Y ∗ ) es localmente asint´oticamente estable. Ya que A1 < 0, r βY ∗ − + 0, (α + X ∗ )3 ˆ entonces la ecuaci´on caracter´ıstica es 3. si E ∗ existe y E = E, λ2 + [det V (E ∗ )]K=K ∗ = 0, cuyas ra´ıces son puramente imaginarias, y d r(cβ − d) 4. [trV (E ∗ )]K=K ∗ = − 6= 0. dK cαβ Entonces se satisfacen todas las condiciones del teorema de bifurcaci´on de Hopf (Murray 2004), por lo que existen soluciones peri´odicas cerca de E ∗ = (X ∗ , Y ∗ ).

103

´ MODELOS DINAMICOS DE POBLACIONES SIMPLES . . .

2.5

2

2

1.5

1.5

1

1

0.5

0.5

Y

2.5

Y X

0

0

0.2

0.4

0.6

0.8

1 X

1.2

1.4

(a)

1.6

1.8

2

0

0

10

20

30

40

50 t

60

70

80

90

100

(b)

Figura 9. En este ejemplo X(0) = 2, Y (0) = 1, r = 1.1, K = 1.2, β = 2.1, α = 1.2, d = 1.3 y c = 2.2. (a) Retrato de fase del sistema (29), el cual muestra que E ∗ = (X ∗ , Y ∗ ) = (0.47, 0.53) es localmente asint´ oticamente estable. (b) La curva s´ olida ilustra la poblaci´ on de presas y la curva punteada ilustra a la poblaci´ on de depredadores. Ambas poblaciones convergen a sus respectivos valores del estado de equilibrio en tiempo finito.

8.3 Simulaci´ on num´ erica

Un estudio anal´ıtico siempre ser´a incompleto sin una verificaci´on num´erica de los resultados. Presentaremos ahora simulaciones por computadora de algunas soluciones del sistema (29). Escojamos la condici´on inicial X(0) = 2 y Y (0) = 1, y los valores de los par´ametros r = 1.1, K = 1.2, β = 2.1, α = 1.2, d = 1.3 y c = 2.2. Entonces la condici´on del teorema 8.1 se satisfacen y E ∗ = (X ∗ , Y ∗ ) = (0.47, 0.53) es localmente asint´oticamente estable. El retrato de fase se ilustra en la figura 9.a. En este caso, las poblaciones de presas y depredadores tienden a sus valores de equilibrio X ∗ and Y ∗ respectivamente, al transcurrir el tiempo (ver figura 9.b). Si ahora incrementamos el valor de K y dejamos fijos el resto de los par´ametros, entonces, por el teorema 8.3, hay un valor cr´ıtico K = K ∗ = 2.14 tal que E ∗ = (X ∗ , Y ∗ ) pierde su estabilidad cuando K pasa a trav´es de este valor cr´ıtico. Para K = 2.58 > K ∗ , observamos que E ∗ = (0.4699, 0.7154) es inestable y hay una o´rbita peri´odica cerca de E ∗ (ver figura 10.a). Las oscilaciones de las poblaciones de presas y depredadores al transcurrir el tiempo se ilustran en la figura 10.b.

104

´ G. SAMANTA Y R. GOMEZ A´IZA

3

3

X

Y

Y 2.5

2.5

2

2

1.5

1.5

1

1

0.5

0.5

0

0

0.2

0.4

0.6

0.8

1 X

1.2

1.4

1.6

1.8

2

0

0

50

...............

100

150

t

(a)

(b)

Figura 10. Aqu´ı todos los par´ ametros son como en la figura 9 excepto por K = 2.5 > K ∗ . (a) Retrato de fase del sistema (29) que muestra una o ´rbita peri´ odica cerca de E ∗ = (0.4699, 0.7154). (b) Oscilaciones de las poblaciones de presas y depredadores al transcurrir el tiempo (la curva s´ olida representa la densidad de poblaci´ on de presas y la curva punteada la densidad de poblaci´ on de depredadores).

9. Bifurcaci´ on de Hopf en el sistema depredadorpresa dependiente de raz´ on Escribamos al sistema (24) de la secci´on §6 como dX dt dY dt

cXY , mY + X f XY = −dY + , mY + X

= X(a − bX) −

(32)

donde X(t) y Y (t) denotan las densidades de poblaci´on de presas y depredadores al tiempo t respectivamente. Aqu´ı, a/b > 0 es la capacidad de carga de las presas, d > 0 es la tasa per c´apita de muerte de los depredadores y a, c, m y f son constantes positivas que representan la tasa de crecimiento intr´ınseco de las presas, la tasa de captura, la constante de saturaci´on media de captura y la tasa de conversi´on respectivamente. El sistema (32) no est´a definido en el origen (0, 0), por lo que ahora supondremos que X(0) > 0,

Y (0) > 0,

y

dX dY = = 0 en (X, Y ) = (0, 0). dt dt

Ahora bien, se puede demostrar que el sistema (32) es continuo y satisface la condici´on de Lipschitz en el primer cuadrante cerrado del plano XY . Estamos interesados en el comportamiento din´amico del sistema (32) en el interior de este primer cuadrante.

´ MODELOS DINAMICOS DE POBLACIONES SIMPLES . . .

105

El sistema (32) siempre tiene los puntos de equilibrio en la frontera E0 = (0, 0) y E1 = (a/b, 0). En esta secci´on estamos interesados en el comportamiento del punto de equilibrio interno E ∗ = (X ∗ , Y ∗ ), el cual existe y es u ´nico si y solo si una de las siguientes dos condiciones se satisface: cd 1. Si c > ma, entonces d < f < . c − ma 2. Si c ≤ ma, entonces f > d. En ambos casos, X ∗ y Y ∗ est´an dadas por X∗ =

f (am − c) + cd , bmf

Y∗ =

(f − d){f (am − c) + cd} . bdf m2

El sistema (32) no puede ser linealizado directamente en E0 = (0, 0). Un cambio dt = (mY + X)dτ en la escala del tiempo transforma el modelo (32) en un sistema polinomial equivalente y entonces la estabilidad de E0 se puede estudiar. As´ı, siguiendo a Xiao y Ruan (2001), observamos que si f ≥ a + d y c − am − dm ≥ 0, entonces E0 es un atractor global del sistema (32). La matriz variacional V (E1 ) en E1 est´a dada por   −a −c V (E1 ) = . 0 f −d Por lo tanto E1 es localmente asint´oticamente estable (resp. inestable) si f < d (resp. f > d). Tambi´en, nuevamente siguiendo a Xiao y Ruan (2001), es sencillo ver que si f ≤ d y c − am − dm < 0, entonces E1 es un atractor global, y si f ≤ d con c − am − dm ≥ 0, entonces E1 es un atractor. Para ver la estabilidad de E ∗ = (X ∗ , Y ∗ ) tenemos el siguiente resultado. Teorema 9.1. Supongamos que E ∗ existe y sea ∆ = (c−am−dm)f 2 + (mf −c)d2 . Entonces E ∗ es localmente asint´oticamente estable si ∆ > 0, y en otro caso es inestable. Demostraci´on. La matriz variacional en E ∗ = (X ∗ , Y ∗ ) es ! cX ∗ Y ∗ cX ∗2 ∗ −bX + − ∗ ∗ 2 ∗ ∗ 2 (mY +X ) (mY +X ) V (E ∗ ) = . f mY ∗2 f mX ∗ Y ∗ − ∗ ∗ 2 (mY +X ) (mY ∗ +X ∗ )2 Es sencillo ver que la traza de V (E ∗ ) es X ∗Y ∗ (mY ∗ + X ∗ )2 2 (c − am − dm)f + (mf − c)d2 = , mf 2

tr(V (E ∗ )) = −bX ∗ + (c − f m)

106

´ G. SAMANTA Y R. GOMEZ A´IZA

y su determinante es det(V (E ∗ )) =

bf mX ∗2 Y ∗ > 0. (mY ∗ + X ∗ )2

La ecuaci´on caracter´ıstica de V (E ∗ ) es λ2 + P λ + Q = 0, donde P = −tr(V (E ∗ )) y Q = det(V (E ∗ )). Ya que Q = det V (E ∗ ) > 0, es claro que E ∗ es localmente asint´oticamente estable o bien inestable, ya sea que P > 0 o P < 0. El siguiente resultado proporciona un criterio para la existencia de una bifurcaci´on de Hopf cerca de E ∗ = (X ∗ Y ∗ ). Teorema 9.2. Supongamos que E ∗ existe con c − m(a + d) > 0 y que p 2 −md + m2 d4 + 4cd2 (c − am − dm) f∗ = . 2(c − am − dm) Entonces una bifurcaci´on de Hopf ocurre en f = f ∗ si 2c − mf ∗ 6= 0. Demostraci´on. Observamos que cada una de las siguientes condiciones se satisface: 1. tr(V (E ∗ ))f =f ∗ = 0, 2. det(V (E ∗ ))f =f ∗ > 0, 3. si E ∗ existe con c = c∗ , entonces la ecuaci´on caracter´ıstica es λ2 + det(V (E ∗ ))f =f ∗ = 0, cuyas ra´ıces son puramente imaginarias, y d (2c − mf ∗ )d2 4. tr(V (E ∗ ))f =f ∗ = 6= 0. df mf 3 Por lo tanto se satisfacen las condiciones del teorema de bifurcaci´on de Hopf (Murray 2004). Observaci´ on 9.1. Cuando ocurre una bifurcaci´on de Hopf, existen o´rbitas peri´odicas de corta amplitud cerca de E ∗ = (X ∗ , Y ∗ ). 9.1 Simulaci´ on num´ erica

Ahora presentamos una simulaci´on computacional de algunas soluciones del sistema determinista (32). Adem´as de verificar los resultados anal´ıticos que hemos deducido, las soluciones num´ericas son muy importantes desde el punto de vista pr´actico. Escogemos X(0) = 0.6 y Y (0) = 0.4, y los par´ametro del sistema (32) los escogemos como a = 1.1, b = 0.7, c = 2.1, m = 1, d = 0.5 y f = 0.79. Entonces E ∗ = (X ∗ , Y ∗ ) = (0.4702, 0.2727) y ∆ = −0.0154 < 0. El teorema 9.1 implica que E ∗ es localmente asint´oticamente estable. El

107

´ MODELOS DINAMICOS DE POBLACIONES SIMPLES . . .

1

0.45

0.9 0.4

0.8 0.35

0.7

0.3 Y

0.6

0.5

0.25

X

0.4 0.2

Y

0.3 0.15

0.1 0.1

0.2

0.2

0.3

0.4

0.1

X

250 t

(a)

(b)

0.5

0.6

0.7

0.8

0.9

1

0

50

100

150

200

300

350

400

450

500

Figura 11. Aqu´ı a = 1.1, b = 0.7, c = 2.1, m = 1, d = 0.5, f = 0.79 y (X(0), Y (0)) = (0.6, 0.4). (a) Retrato de fase del sistema (32). Claramente es una espiral estable convergente a E ∗ = (X ∗ , Y ∗ ) = (0.4702, 0.2727). (b) La curva s´ olida ilustra la poblaci´ on de presas y la punteada la poblaci´ on de depredadores. Ambas poblaciones convergen a los valores de estado de equilibrio en tiempo finito. 1.5 0.4

X Y ............

0.35

0.3

1

Y

0.25

0.2

0.15

0.5 0.1

0.05

0

0

0

X

250 t

(a)

(b)

0.5

1

1.5

0

50

100

150

200

300

350

400

450

500

Figura 12. (a) Retrato fase del sistema (32) que ilustra una o ´rbita peri´ odica cerca de E ∗ = (X ∗ , Y ∗ ) = (0.4233, 0.2624) cuando los par´ ametros son como en la figura 11 excepto por f = 0.81 > f ∗ . (b) Oscilaciones de las poblaciones de presas y depredadores al transcurrir el tiempo. (La curva s´ olida representa la densidad de poblaci´ on de las presas y la punteada la de los depredadores.)

retrato fase de esta soluci´on se ilustra en la figura 11.a. Claramente la soluci´on es una espiral estable que converge a E ∗ . En la figura 11.b se ilustran las poblaciones de presas y depredadores que convergen a sus valores de estados estables X ∗ y Y ∗ respectivamente, en tiempo finito. Si gradualmente aumentamos el valor de la tasa de conversi´on f (manteniendo fijos el resto de los par´ametros), se observa que el comportamiento del sistema cambia cuando f pasa a trav´es del valor de bifurcaci´on f ∗ = 0.8048 (el cual se obtiene del teorema 9.2). Para

108

´ G. SAMANTA Y R. GOMEZ A´IZA

f = 0.81 > f ∗ , vemos que ∆ = 0.0056 > 0, y por lo tanto, por el teorema 9.1, E ∗ = (X ∗ , Y ∗ ) = (0.4233, 0.2624) es inestable. El retrato de fase correspondiente es una ´orbita peri´odica alrededor de E ∗ = (X ∗ , Y ∗ ) (ver figura 12.a). En este caso, las oscilaciones de las poblaciones de presas y depredadores se ilustran en la figura 12.b. As´ı pues, al usar la tasa de conversi´on f como control, es posible romper el comportamiento estable del sistema y llevarlo a un estado inestable. Con este control tambi´en es posible mantener las poblaciones dentro de ciertos niveles deseados.

10. Otras modificaciones 10.1 Efectos de retraso del tiempo

En la discusi´on anterior hemos supuesto que la tasa de cambio del tai ma˜ no de la poblaci´on dN depende u ´nicamente del tama˜ no instant´aneo dt de la poblaci´on Ni (t). Sin embargo, nos encontramos en la situaci´on en i la que dN tambi´en depende de los tama˜ nos de la poblaci´on en instantes dt de tiempo anteriores, i. e. de Ni (t − τ ) con τ ≥ 0. La respuesta de la tasa de crecimiento puede retrasarse debido a diversas razones: 1. Periodo de maduraci´on (e. g. el periodo requerido para que una larva se convierta en adulto). 2. Periodo de gestaci´on (e. g. el periodo necesario para que un depredador digiera una presa). 3. Etc´etera. As´ı, podemos seguir modificando el modelo (24) y obtener   dN1 βN1 N2 N1 = rN1 1 − − dt K αN2 + N1   dN2 cβN1 (t − τ ) = −d + N2 , dt αN2 (t − τ ) + N1 (t − τ )

1 2 con dN = dN = 0 para (N1 , N2 ) = (0, 0) y con condiciones iniciales dt dt N1 (θ) = φ1 (θ) y N2 (θ) = φ2 (θ) con φi (θ) ≥ 0 (con i = 1, 2), y para toda θ ∈ [−τ, 0], donde φi (θ) son funciones continuas en θ ∈ [−τ, 0]. Para tener significado biol´ogico, suponemos tambi´en que φi (0) > 0 (con i = 1, 2).

10.2 Ruido

Finalmente, podemos suponer que las fluctuaciones en el medio ambiente se manifestar´an principalmente como fluctuaciones tanto en el

´ MODELOS DINAMICOS DE POBLACIONES SIMPLES . . .

109

coeficiente de crecimiento intr´ınseco de las presas, r, como en la mortalidad de los depredadores, d, ya que estos son los par´ametros principales sujetos al acoplamiento de depredador-presa en estos medios ambientes. Entonces el comportamiento de (23) en un ambiente aleatorio se enmarcar´a en el contexto del siguiente modelo:   dN1 r βN1 N2 = r + η1 (t) − N1 N1 − dt K α + N1 cβN1 N2 dN2 = (−d + η2 (t)) N2 + , dt α + N1 donde los t´erminos de perturbaci´on η1 (t) y η2 (t) son ruidos coloreados independientes (o procesos Ornstein-Uhlenbeck), que son ruidos realistas. La esperanza matem´atica y la funci´on de correlaci´on de los procesos ηj (t) est´an dados por hηj (t)i = 0, hηj (t1 )ηj (t2 )i = j δj exp(−δj |t1 − t2 |) (con j = 1, 2), donde j (> 0) y δj−1 (> 0) son respectivamente la intensidad y el tiempo de correlaci´on del ruido ηj (t) y h·i representa el promedio en conjunto del proceso p estoc´astico. −1 Si δj → 0, entonces ηj (t) → 2j ξj (t), donde ξj (t) son variables aleatorias independientes con distribuci´on normal est´andar (esperanza cero y varianza uno), i. e. ruido blanco caracterizado por hξj (t)i = 0 y hξj (t1 )ξj (t2 )i = δ(t1 − t2 ) (con j = 1, 2), donde δ(t) es la funci´on delta de Dirac. Agradecimientos. El Dr. G. Samanta agradece a la TWAS-UNESCO y a la Universidad Nacional Aut´onoma de M´exico (UNAM) por el apoyo financiero. Agradece tambi´en al Prof. Javier Bracho Carpizo y al Prof. Marcelo Aguilar del Instituto de Matem´aticas de la UNAM por su apoyo y est´ımulo. Los autores agradecen tambi´en a los ´arbitros an´onimos y a la Coordinadora Editorial de Miscelanea Matem´atica, la Dra. Ana Meda Guardiola, por su cuidadosa lectura y por sus valiosos comentarios y sugerencias, mismos que nos ayudaron para mejorar significativamente la presentaci´on de este trabajo.

Bibliograf´ıa [1] J. F. Andrews, ((A mathematical model for the continuous culture of microorganisms utilizing inhibitory substrates)), Biotechnol. Bioeng., vol. 10, 1968, 707–723. ´ [2] E. J. Avila Vales, ((Coexistencia para un sistema con dos especies en competencia)), Miscelanea Matem´ atica, vol. 30, 2000, 17–25. [3] W. E. Boyce y R. DiPrima, Elementary differential equations and boundary value problems, Wiley, 2005. [4] M. Braun, Differential equations and their applications, Springer Verlag, Berlin, 1983. [5] L. Esteva Peralta y M. Falconi Maga˜ na, Biolog´ıa matem´ atica. un enfoque desde los sistemas din´ amicos, Las Prensas de Ciencias, 2012.

110

´ G. SAMANTA Y R. GOMEZ A´IZA

´ [6] A. G. Estrella Gonz´ alez, G. E. Garc´ıa Almeida y E. J. Avila Vales, ((Estabilidad local de ecuaciones diferenciales ordinarias con retardo y aplicaciones)), Miscelanea Matem´ atica, vol. 51, 2010, 73–92. [7] J. L. Guti´errez y F. S´ anchez-Gardu˜ no, Matem´ aticas para las ciencias naturales, Sociedad Matem´ atica Mexicana, 1998, N´ umero 11 de series de textos de la Sociedad Matem´ atica Mexicana. [8] J. Hale, Ordinry differential equations, Wiley-Interscience, New York, 1969. [9] C. Holling, ((The components of predation as revealed by a study of small–mammal predation of the european sawfly)), Can. Entomol., vol. 91, 1959, 293–320. [10] , ((Some characteristics of simple types of predation an parasitism)), Can. Entomol., vol. 91, 1959, 385–398. [11] W. Horsthemke y R. Lefever, Noise induced transitions, Springer-Verlag, Berlin, 1984. [12] M. Kot, Elements of mathematical ecology, Cambridge University Press, Cambridge, 2001. [13] A. Lotka, Elements of physical biology, The Williams and Wilkins Co., Baltimore, 1925. [14] A. Maiti y G. Samanta, ((Deterministic and stochastic analysis of a prey–dependent predator-prey system)), International Journal of Mathematical Education in Science and Technology, vol. 36, 2005, 65–83. [15] , ((Deterministic and stochastic analysis of a ratio–dependent prey-predator system)), International Journal of Systems Science, vol. 37, 2006, 817–826. [16] T. Malthus, An essay on the principle of population, as it affects the future improvement of society, with remarks on the speculations of mr. godwin, m. condorcet and other writers, J. Johnson, London, 1798. Reprint, University of Michigan Press, U.S.A., 1959. [17] R. May, Stability and complexity in model ecosystems, Princeton University Press, Princeton, 1973. , ((Stability in randomly fluctuating versus deterministic environment)), Am. [18] Nat., vol. 107, 1973, 621–650. [19] J. Murray, Mathematical biology i, Springer–Verlag, 2004. [20] M. Nu˜ nez L´ opez y J. X. Velasco Hern´ andez, ((Competencia y superinfecci´ on en plagas y enfermedades)), Miscelanea Matem´ atica, vol. 49, 2009, 63–82. [21] E. Pielou, Mathematical ecology, A Wiley–Interscience Publication, John–Wiley and Sons, New York, 1977. [22] F. S´ anchez Gardu˜ no, V. Castellanos Vargas, I. Quilant´ an Ortega y G. C. Vel´ azquez L´ opez, ((Matem´ aticas en la distribuci´ on espacial de poblaciones)), Miscelanea Matem´ atica, vol. 48, 2009, 75–101. [23] F. S´ anchez Gardu˜ no, P. Miramontes y J. L. Guti´errez, Cl´ asicos de biolog´ıa matem´ atica, Siglo XXI–CEICH, UNAM, 2002. [24] J. M. Smith, Models in ecology, Cambridge University Press, Cambridge, 1975. [25] Y. M. Svirezhev y D. O. Logofet, Stability of biological communities, Mir Publisher, Moscow, 1983. [26] J. X. Velasco Hern´ andez, ((Sobre enfermedades infecciosas)), Miscelanea Matem´ atica, vol. 29, 1999, 51–72. , ((Modelos matem´ aticos en epidemiolog´ıa: enfoques y alcances)), Miscelanea [27] Matem´ atica, vol. 44, 2007, 11–27. [28] P. F. Verhulst, ((Notice sur la loi que la population pursuit dans son accroissement)), Corr. Math. Phys., vol. 10, 1938, 113–121. [29] V. Volterra, ((Variazioni e fluttuazioni del numers d’individui in specie animali conviventi)), Mem. Acad. Lineii Roma 2, 1926, 31–113. [30] D. Xiao y S. Ruan, ((Global dynamics of a ratio-dependent predator–prey system)), Journal of Mathematical Biology, vol. 43, 2001, 268–290.