Modelo de Lotka-Volterra

Modelo de Lotka-Volterra Alfonso de Miguel July 1, 2014

Abstract En este trabajo se plantea un estudio b´asico de modelos matem´aticos usados t´ıpicamente en el estudio de ecolog´ıa y din´amica de poblaciones. Daremos un repaso al archiconocido modelo de Lotka-Volterra desde el punto de vista de la teor´ıa de sistemas din´ amicos (determinista). A continuaci´on se trabajar´a su deducci´on y an´ alisis considerando los procesos elementales del sistema (procesos de salto) para deducir una ecuaci´ on maestra del sistema. Compararemos el modelo determinista y este u ´ltimo, basado en el an´ alisis de las fluctuaciones internas.

1 Introducci´on. La incursi´ on de la f´ısica y de la matem´atica en terrenos no convencionales es desde hace largo tiempo una barrera m´ as que superada y el uso de herramientas t´ıpicas de estos campos tan fundamentales ha de ser siempre bienvenido por su capacidad de aportar nuevos puntos de vista, m´etodos y conclusiones. La generalidad que ofrece la teor´ıa de sistemas din´ amicos, as´ı como la teor´ıa de procesos estoc´asticos, nos permite modelar diversos fen´ omenos caracter´ısticos de cada campo de estudio y dotarlos de un cuerpo matem´atico y anal´ıtico, con el fin de lograr una comprensi´on m´as profunda y una mayor capacidad de predicci´ on y an´ alisis del fen´omeno en cuesti´on. En este trabajo me centro en el llamado modelo o ecuaciones de Lotka-Volterra [1], de amplio uso y estudio en ecolog´ıa, biolog´ıa y din´amica de poblaciones en general. Dado su atractivo, el modelo aparece en numerosos textos y art´ıculos, bien sea enfocando el an´alisis desde el punto de vista de la teor´ıa de sistemas din´amicos (estabilidad, bifurcaciones, etc.) o como posible modelo descriptivo de sistemas depredador-presa o de otro tipo de relaciones entre elementos de un ecosistema o sistema de car´acter m´as general [2]. Tambi´en existen, por supuesto, textos dedicados a un estudio estoc´astico del modelo, bien sea desde el punto de vista de las fluctuaciones internas del sistema o por causas de perturbaciones externas [3]. Siguiendo la filosof´ıa del curso lo propio aqu´ı es enfocar la cuesti´ on desde esta u ´ltima perspectiva.

1

Lo que aqu´ı se pretende es realizar un estudio b´asico del modelo de Lotka-Volterra en dos dimensiones. En la primera secci´on se presenta el sistema din´amico y se repasan las caracter´ısticas b´ asicas del modelo: construcci´on del modelo, significado, soluciones, puntos fijos y estabilidad. En la segunda secci´on entramos a considerar ya la cuesti´on estoc´astica: tomando el modelo como consecuencia de la ocurrencia de una serie de procesos elementales de variable discreta llegaremos a la deducci´on de su ecuaci´on maestra, abordaremos su soluci´ on estacionaria y de ah´ı deduciremos las ecuaciones macrosc´opicas de Lotka-Volterra. Compararemos los resultados entre el modelo determinista y el modelo basado en las fluctuaciones intr´ınsecas.

2 Modelos de din´amica poblacional: Lotka-Volterra. 2.1 Construcci´on de modelos. Vamos a mostrar o, mejor dicho, a construir el modelo de Lotka-Volterra de depredadorpresa paso a paso para ver su raz´ on de ser como modelo descriptivo (y muy simplificado) de la din´ amica de poblaciones de un ecosistema 1 . Sea entonces un ecosistema donde habita una u ´nica especie y cuya poblaci´on (densidad) en el h´ abitat designamos por la variable x. Parece razonable postular que la evoluci´on en la poblaci´ on x de la especie presa en un tiempo t + dt ser´a proporcional al tama˜ no de la poblaci´ on en t, x(t), estableciendo la igualdad un multiplicador r tal que dx = rx, (1) dt donde uno puede definir r como el n´ umero de nacimientos menos el n´ umero de muertes de miembros de la especie (crecimiento natural), esto es: r es la tasa neta de crecimiento per c´apita. Este modelo para r > 0 recibe el nombre de modelo de Malthus [6] y, consecuentemente, implica un crecimiento exponencial de la especie; lo cual parece a todas luces irreal. El modelo log´ıstico o modelo de Verlhust, tambi´en de gran popularidad, ofrece una visi´ on m´ as realista introduciendo el concepto o efecto de competencia entre miembros de una propia especie en lucha por los recursos limitados que el ecosistema ofrezca. El modelo log´ıstico se leer´ıa as´ı  dx x =r 1− x, (2) dt K donde K es la denominada capacidad de carga 2 . La novedad del modelo log´ıstico es que introduce un t´ermino cuadr´ atico negativo en x que act´ ua como regulador de la poblaci´on, 1

Para ver una exposici´ on hist´ orica del desarrollo del modelo puede consultarse [4], [5] La capacidad de carga de una especie biol´ ogica en un ambiente es el tama˜ no m´ aximo de poblaci´ on que el ambiente puede soportar indefinidamente en un periodo determinado, teniendo en cuenta el alimento, agua, h´ abitat, y otros elementos necesarios disponibles en ese ambiente. M´ as detalles: http: //es.wikipedia.org/wiki/Capacidad_de_carga 2

2

es decir, suprime el crecimiento exponencial incontrolado del modelo de Malthus. Por supuesto, esta correcci´ on no es la u ´nica manera de describir la saturaci´on o limitaci´on al crecimiento de una poblaci´ on y dependiendo del problema concreto la modelizaci´on puede requerir una funci´ on diferente [6]; es simplemente un punto de partida. Tampoco hay nada que limite a la tasa de crecimiento r como un factor constante, sino que puede depender de la propia poblaci´ on o de otros factores. As´ı que, desde una perspectiva m´as general, uno puede escribir la evoluci´on en la poblaci´on de una especie como dx = ρ(x)x. dt

(3)

Consideremos la situaci´ on donde tenemos una especie depredadora pero no hay presa. En ausencia de especies inferiores en la cadena de alimentaria se supone una extinci´on exponencial a un ritmo µ, tal que dy = −µy, (4) dt donde y es la densidad de miembros de la especie depredadora. Ahora introduzcamos ambas especies en el mismo ecosistema. Bajo la premisa de que la relaci´ on entre ambas es de presa-depredador, la evoluci´on din´amica del n´ umero de miembros de estas especies ya no depender´a u ´nicamente de la propia especie sino de c´omo sea la interacci´ on con la rival. De este modo, la presencia del depredador tendr´a un efecto negativo en la tasa de crecimiento per c´apita de la presa, mientras que el depredador basa su crecimiento neto en tanto se alimente de la presa. Denotaremos por φ(x, y) a la tasa de desaparici´ on per c´apita de la presa debido a la depredaci´on y a la tasa de crecimiento del depredador debido al consumo per c´apita de la presa como σ(x). Por lo tanto, de manera general, se tendr´a este sistema de ecuaciones diferenciales x˙ = ρ(x) − φ(x, y), x y˙ = σ(x) − µ, y

(5) (6)

Se asume que tanto φ como σ son funciones suaves de sus argumentos. El coeficiente que acompa˜ na a φ(x, y) es negativo, ya que la interacci´on presa-depredador sobre la presa se traduce en un perjuicio, en un decrecimiento del n´ umero de miembros de esta especie. Por contra, para el depredador esta interacci´on es beneficiosa en t´erminos de su crecimiento y σ(x) va con signo positivo3 Reajustando los t´erminos, podr´ıamos expresar definitivamente el modelo como dx = x[ρ(x) − φ(x, y)] = F (x, y), dt

(7)

3 De manera general estas funciones, φ(x, y) y σ(x) no tienen porqu´e ser iguales aunque se trate de la misma interacci´ on presa-depredador. De todos modos, en lo posterior se asumir´ an iguales.

3

dy = y[σ(x) − µ] = G(x, y). (8) dt que es una versi´ on general de modelo de Lotka-Volterra bidimensional de depredadorpresa con crecimiento log´ıstico de la presa. Uno puede bajo esta generalidad, y con intenci´on de refinar el modelo, seguir aportando condiciones l´ ogicas o esperables a satisfacer para las funciones φ(x, y) y σ(x) (v´ease por ejemplo [7]), pero opto ya por centrar la cuesti´ on en un modelo concreto que ser´a sobre el que luego haremos el trabajo estoc´astico.

2.2 Lotka-Volterra de presa-depredador. El modelo concreto de trabajo es este Lotka-Volterra de presa-depredador4 : dx = rx − αyx, dt

(9)

dy = βyx − µy (10) dt Donde r es la tasa de crecimiento natural per c´apita (como en la ecuaci´on (1)); ω1 es la constante de acoplo presa-depredador; n´otese que φ(x, y) = αy, σ(x) = βx. En el caso de la tasa de desaparici´on per c´apita de la presa, funci´on φ(x, y), esta es proporcional a la poblaci´ on de depredadores en un tiempo t; el factor α se suele denominar eficiencia de captura y es una medida del efecto del depredador en la tasa de crecimiento per c´ apita de la presa. El producto αx representa la llamada respuesta funcional del depredador, esta describe c´omo la tasa de captura de presas viene afectada por la abundancia de presas 5 . Por otro lado, la tasa de crecimiento del depredador debido al consumo per c´ apita de la presa, funci´on σ(x), es proporcional a la poblaci´on de presas en un tiempo t. La constante β recibe el nombre de eficiencia de conversi´ on, esto es, la habilidad de los depredadores para convertir presa en tasa de crecimiento per c´apita adicional para la poblaci´on de depredadores. El producto βx se denomina tambi´en respuesta num´ erica de la poblaci´on depredadora. Resumamos, para que queden claras las hip´otesis del modelo: - En ausencia de depredeadores, las presas tienen un comportamiento malthusiano. 4

Este modelo es el propuesto orginalmente por Volterra en 1926, para explicar las oscilaciones de los bancos de pesca en el mar Adri´ atico (y por Lotka, independientemente en 1920). 5 Existen distintos tipos de respuesta funcional con el fin de modelar el sistema de acuerdo a unas exigencias m´ as realistas y concretas, v´ease [8]. En el caso de la ecuaci´ on (9), esta es lineal con x.

4

- El efecto de los depredadores es reducir el n´ umero de presas de forma proporcional a ambas poblaciones. - En ausencia de presas, los depredadores se extinguen exponencialmente. - Las presas contribuyen a hacer crecer la poblaci´on de depredadores de forma proporcional a su propia poblaci´ on. Esto es una simplificaci´on considerable de alg´ un mecanismo complejo que convertir´ıa carne de presas en hijos de depredadores. Evidentemente estamos en una situaci´on muy simplista. Aparte de la no inclusi´on de t´erminos log´ısticos o alg´ un otro tipo de saturaci´on para modelar un crecimiento controlado de las poblaciones, se obvian otros detalles importantes como el hecho de que se asume que los depredadores son insaciables y, adem´as, cuantas m´as presas hay, m´as comen. Tambi´en se da por hecho que todos los depredadores tienen las mismas opciones de atrapar a sus presas. Hay modelos con difusi´on, donde se introduce una dependencia espacial, que intentan aplacar estas cuestiones. - Puntos fijos: Las nulclinas del sistema son aquellas curvas de pendiente nula, es decir 0 = x(r − αy) y 0 = y(βx − µ) → Nulclinas de x: y=

r α

y

x = 0.

(11)

Nulclinas de y:

µ y y = 0. (12) β La intersecci´ on de las nulclinas de x e y nos dan los puntos fijos del sistema, as´ı que tenemos: (x01 , y01 ) = (0, 0), (13)   µ r (x02 , y02 ) = , . (14) β α La primera soluci´ on implica la extinci´on total de ambas especies. La segunda significa una situaci´ on de coexistencia de ambas especies. Aparte uno tiene la soluci´on (∞, 0) que acaece en ausencia de depredadores como condici´on inicial, y(0) = 0. x=

- Trayectorias y estabilidad: Para estudiar la estabilidad del sistema es preciso linealizar el sistema de ecuaciones diferenciales en un entorno de los puntos fijos. Por medio de la evaluaci´ on del jacobiano   r − αy −αx , J(x, y) =  (15) βy βx − µ en los puntos fijos, extraemos los autovalores asociados para conocer la estabilidad de cada uno. Para (0, 0) el espectro de J es σ(J(0, 0)) = {r, −µ}, 5

(16)

Para el sistema presa-depredador r > 0 y µ > 0 y por lo tanto tenemos un autovalor positivo y otro negativo, estamos ante un punto-silla y por tanto un punto fijo inestable. Calculando los autovectores asociados se puede ver que las soluciones crecen exponencialmente para x(y) y decrecen exponencialmente para Para ( βµ , αr ) el espectro de J es    µ r √ σ J = ±i rµ. , (17) β α Los autovalores son imaginarios puros, en consecuencia estamos ante un centro, este tipo de punto de equilibrio es estable pero no asint´oticamente estable, se habla de estabilidad neutra o marginal. Las trayectorias (x(t), y(t)) que uno obtiene mediante esta an´alisis lineal implican soluciones peri´ odicas entorno al punto fijo de coexistencia de especies.

Por otro lado, podemos encontrar las ecuaciones de las o´rbitas del modelo presapredador. Tomando las ecuaciones (9) y (10), despejamos en ambas dt e igualamos, reajustando t´erminos se llega a dy y(βx − µ) = dx x(r − αy)

(18)

las ´orbitas del sistema para x, y ≡ 0 son las curvas soluciones a esta ecuaci´on diferencial de 1-orden.  Z  Z  dy dx µ r → (r − αy) = (βx − µ) → β− dx = − α dy + C y x x y βx − µ ln x + αy − r ln y = h = constante → 6

→

eβx eαy = H(x, y) = H = constante. x−µ y −r

(19)

Las ´orbitas del sistema son la familia de curvas definidas por la ecuaci´on anterior6 La funci´on H(x, y) es constante a lo largo de las trayectorias (x(t), y(t)) soluci´on del sistema de ecuaciones diferenciales de Lotka-Volterra (9-10). La ecuaci´on (19) muestra que tenemos una integral o constante del movimiento. De hecho, el sistema de LotkaVolterra puede plantearse como un sistema din´amico Hamiltoniano [9],[10], siendo el propio H(x, y) el Hamiltoniano de este sistema. - Conservaci´ on de los promedios: Sup´ongase condiciones iniciales positivas x(0) = x0 > 0 e y(0) = y0 > 0. Si T es el periodo de una ´orbita cerrada que pasa por (x0 , y0 ), entonces x˙ = r − αy → ln x(t) − ln x(0) = x

Z

T

Z [r − αy(t)]dt → 0 = rT − α

0

T

y(t)dt, 0

usando la propiedad de periodicidad de las ´orbitas. As´ı que el valor medio sobre un periodo T es Z 1 T r hyiT = y(t)dt = , (20) T 0 α 6

Para ver una demostraci´ on de que estamos ante una soluci´ on peri´ odica v´ease con m´ as detalle [].

7

Haciendo lo propio para la ecuaci´ on de los depredadores y˙ = βx − µ → ln y(t) − ln y(0) = y

Z

T

Z [βx(t) − µ]dt → 0 = −µT + β

0

T

x(t)dt. 0

Y an´alogamente hxiT =

1 T

Z

T

x(t)dt = 0

µ . β

(21)

El valor mediado de las variables de estado (x(t), y(t)) sobre el periodo de una ´orbita arbitraria es igual al valor del punto fijo de coexistencia.

3 Lotka-Volterra a partir de las fluctuaciones internas El modelo de Lotka-Volterra que hemos revisado en la secci´on anterior, y cualquiera de sus variaciones mencionadas y referenciadas, son al fin y al cabo sistemas de ecuaciones diferenciales deterministas y su estudio matem´atico est´a basado en el uso de conceptos y resultados de la teor´ıa de sistemas din´amicos convencional (ecuaciones diferenciales, teor´ıa de bifurcaciones, sistemas Hamiltonianos, etc.). Los modelos tradicionales de ecolog´ıa matem´atica describen en forma diferencial la evoluci´on de las densidades. Como se ha visto en la secci´on previa, en general, podemos

8

construir estos modelos sin necesidad de conocer al detalle la interacci´on entre los individuos que conforman la poblaci´ on. Se supone, no obstante, que cada t´ermino tiene su fundamento en las interacciones puntuales y que representan efectos macrosc´opicos que involucran a muchos individuos, siendo las ecuaciones de evoluci´on suaves y continuas. Sin embargo, las poblaciones son discretas y las interacciones est´an individualizadas. Por medio de un an´ alisis de todos estos procesos de interacci´on a escala microsc´opica deber´ıa ser posible la obtenci´ on de las leyes macrosc´opicas que obtuvimos, surgiendo estas como un l´ımite macrosc´ opico, un coarse-graining, de los fen´omenos observados. Algo que nos ense˜ na la teor´ıa de fluctuaciones en sistemas din´amicos es esta posibilidad de obtener ecuaciones deterministas/macrosc´opicas a partir de una consideraci´on estoc´astica del proceso o fen´ omeno de turno. Es lo que se denomina analizar el sistema desde el punto de vista de las fluctuaciones internas del sistema.

3.1 Reacciones elementales. Entremos en materia. Mediante una considerable simplificaci´on uno puede describir el modelo de Lotka-Volterra en t´erminos de reacciones de la siguiente manera: X + A → 2X,

(22)

→ 2Y,

(23)

X +Y Y

→ ∅.

(24)

Donde cada reacci´ on se produce a una tasa ki=0,1,2 . X e Y siguen designando respectivamente a la poblaci´ on de presa y depredador, pero ahora son variables discretas: representan en el n´ umero de miembros de dicha especie. La primera reacci´on representa el proceso de crecimiento, m´ as concretamente de autorreplicaci´on, de la especie presa a partir de una fuente de alimento A que se toma constante7 La segunda reacci´on representar´ıa el proceso de depredaci´ on. La interacci´on de una especie presa y una depredadora dan lugar a la desaparici´ on de una tal presa y al surgimiento de una depredadora8 Por u ´ltimo, la reacci´ on m´ as realista representa simple y llanamente la muerte de un predador. Estas reacciones presentadas ocurren con una cierta probabilidad y en consecuencia las variables de estado (X, Y ) ahora son variables estoc´asticas y su mera especificaci´on en un instante de tiempo t concreto no aporta la misma informaci´on que en una situaci´on determinista. Es por ello que debemos trabajar mediante la funci´on de probabilidad del sistema: P (X, Y ; t) denotar´ a entonces la probabilidad conjunto de que en un instante de 7 Tomar constante la alimentaci´ on de esta especie ser´ıa equivalente a la consideraci´ on de crecimiento explosivo de la especie a ` la Malthus. Por otro lado, puede verse aqu´ı que no se plantea ning´ un tipo de reproducci´ on sexual ni gestaci´ on sino que sencillamente el organismo X se autorreplica a partir de la ingesta de A. 8 Como ya mencionamos en la secci´ on previa, esta hip´ otesis tambi´en es muy brusca. Se presupone una completa eficiencia y obscuro mecanismo por el cual el depredador transforma a la presa consumida en un nuevo miembro Y de pleno de derecho.

9

tiempo t haya X miembros de la especie presa e Y miembros de la especie predadora. La ecuaci´on maestra que buscamos es una ecuaci´on que nos informa sobre la evoluci´on temporal de esta probabilidad; esta puede ser obtenida al considerar los procesos elementales que ocurran en un intervalo de tiempo (t, t + dt) que puedan contribuir a P (X, Y ; t + dt), que son Suceso elemental 1: La poblaci´ on era (X, Y ) en el instante t y no ocurren cambios durante dt. No nacen ni mueren ni presas ni depredadores. Suceso elemental 2: La poblaci´on era (X − 1, Y ) en el instante t y una presa se reprodujo (tuvo lugar la reacci´ on primera, ecuaci´on (22)) durante dt. Suceso elemental 3: La poblaci´on era (X, Y + 1) en el instante t y un depredador muri´o (reacci´ on tercera, ecuaci´ on (24)) durante dt. Suceso elemental 4: La poblaci´ on era (X + 1, Y − 1) en el instante t y un predador devor´o a una presa (reacci´ on segunda, ecuaci´on (23)) durante dt. El hecho de que ocurra una reacci´on u otra es independiente. Adem´as, siendo p la probabilidad de que un suceso ocurra, uno se refiere a 1 − p como la probabilidad de que un tal suceso no tome lugar. Con todo esto, la probabilidad asociada al primer suceso elemental es (1 − X k¯0 dt)(1 − XY k¯1 dt)(1 − Y k¯2 dt). (25) Siendo entonces X k¯0 dt,

(26)

la probabilidad de autorreplicaci´ on o reproducci´on de la presa X, que es proporcional al n´ umero de presas X en el instante t y a la constante de reacci´on k¯0 que vendr´ıa a ser un equivalente de la tasa de crecimiento natural de la especie. Mientras que XY k¯1 dt,

(27)

es la probabilidad de ocurrencia de la interacci´on depredador-presa, que es proporcional al n´ umero de presas X y de depredadores Y en el instante t y a la constante de reacci´on ¯ k1 (la interacci´ on est´ a regida por la misma constante de acoplo desde el punto de vista de ambos). Y por u ´ltimo Y k¯2 dt, (28) es la probabilidad de muerte de un depredador, y es proporcional al n´ umero de depredadores en el instante t y a la constante k¯2 , que medir´ıa la tasa de mortalidad del predador. Las probabilidades asociadas al resto de sucesos se construir´ıan de la misma forma a la que siguen las ecuaciones (26-28), teniendo en cuenta el n´ umero de especies que habr´ıa en el instante t considerado.

10

3.2 Ecuaci´on Maestra y soluci´on estacionaria. Sabiendo la forma de las probabilidades de transici´on, expresar´ıamos la probabilidad conjunta P (X, Y, t + dt) como la siguiente suma de contribuciones P (X, Y ; t + dt) = P (X, Y ; t)[1 − X k¯0 dt][1 − XY k¯1 dt][1 − Y k¯2 dt]+ +P (X − 1, Y ; t)(X − 1)k¯0 dt+ +P (X, Y + 1; t)(Y + 1)k¯2 + + P (X + 1, Y − 1)(Y − 1)(X + 1)k¯1 .

(29)

Cada sumando est´ a conformado por la probabilidad de estar en un estado (m, n) por la probabilidad de transici´ on de ir de un estado (m, n) a otro (m0 , n0 ). Desechando t´erminos de orden O(dt2 ) y superiores, dividiendo por dt y tomando el l´ımite dt → 0 tenemos la ecuaci´on de evoluci´ on para la probabilidad conjunta de X e Y , la ecuaci´ on maestra9 : ∂P (X, Y ; t) = −(X k¯0 + XY k¯1 + Y k¯2 )P (X, Y ; t)+ ∂t +(X − 1)k¯0 P (X − 1, Y ; t)+ +(Y + 1)k¯2 P (X, Y + 1; t)+ + (Y − 1)k¯1 (X + 1)P (X + 1, Y − 1; t).

(30)

Para un manejo m´ as compacto escribimos esta ecuaci´on por medio de los operadores de traslaci´on: ∂P (X, Y ; t) ˆ −1 −1)XY k¯1 P (X, Y ; t). ˆY −1)k¯2 Y P (X, Y ; t)+(E ˆX E ˆ −1 −1)X k¯0 P (X, Y ; t)+(E = (E Y X ∂t (31) La bidimensionalidad del sistema complica enormemente, si no la hace imposible, la resoluci´on de la ecuaci´ on maestra a tiempo t. Tampoco intentando la construcci´on de una funci´ on generatriz facilita el problema. La soluci´on estacionaria verifica que ∂Pst (X, Y ) = 0. ∂t

(32)

As´ı: ˆ −1 − 1)X k¯0 P (X, Y ; t) + (E ˆY − 1)k¯2 Y P (X, Y ; t) + (E ˆX E ˆ −1 − 1)XY k¯1 P (X, Y ; t). 0 = (E X Y (33) Extraer de aqu´ı la soluci´ on no parace nada claro tampoco pero recurro a un argumento de Haken [11] (recogido tambi´en en [12],[13]), basado en dos puntos: Primero suponer el ansantz Pst (X, Y ) = δX,0 δY,0 . (34) 9 Se ha elegido esta forma de construir la ecuaci´ on maestra para enfatizar la consideraci´ on de los procesos elementales dados. Una forma m´ as directa para construir las probabilidades de transici´ on habr´ıa sido seguir el procedimiento usado para las reacciones qu´ımicas [18], por medio de la combinatoria.

11

Sustituy´endose en la ecuaci´ on (33), esta se verifica. Pero, ¿es u ´nica esta soluci´on? Aqu´ı viene el segundo punto: todos los puntos (X, Y ) del sistema est´an conectados por al menos una trayectoria con cualquier otro punto (X 0 , Y 0 ) del sistema, lo cual permite siempre una posibilidad no nula de caer al punto (0, 0) y este estado es absorbente y estable; una vez se alcanza es imposible recuperarse, no se puede crear vida de la nada. V´ease la figura 4.

Ahora, este resultado puede parecer inesperado. Nuestro sistema determinista estudiado en la secci´ on primera ofrec´ıa dos puntos fijos: la soluci´on de coexistencia y la extinci´on total. Aqu´ı no hay ni rastro de la primera y la u ´nica perspectiva es ciertamente apocal´ıptica. ¿Cu´ al es la explicaci´ on a esta discrepancia entre el modelo estoc´astico y el determinista? La soluci´ on de existencia era un estado punto-silla (inestable) en el modelo determinista, pero la de coexistencia implicaba ´orbitas cerradas en el espacio de fases del sistema y el problema radica en su estabilidad. As´ı es, obtuvimos una soluci´on de estabilidad neutra, esto significa que en presencia de fluctuaciones, hecho inevitable en un sistema estoc´ astico, las trayectorias cerradas van a verse perturbadas en mayor o menor medida y a la vista de la soluci´on estacionaria obtenida, en t → ∞ la probabilidad de que haya una fluctuaci´ on que lleve al sistema a la extinci´on tiende a uno. Debido a la naturaleza discreta del sistema estoc´astico, las fluctuaciones en el n´ umero de especies pueden ser determinantes para la extinci´on si el n´ umero de especies se va haciendo peque˜ no. Adem´ as, este resultado significa que cualquier sistema con un n´ umero finito de part´ıculas siempre alcanzar´a el estado de extinci´on, pero la escala temporal de extinci´ on puede ser muy grande para sistemas razonablemente grandes [12]. Esta caracter´ıstica est´ a ausente en el llamado modelo de campo-medio, la descripci´on puramente determinista, ya que las densidades de poblaci´on pueden ser arbitrariamente peque˜ nas sin que el destino de la poblaci´on se vea afectado dram´aticamente. El punto fijo de coexistencia de especies marginalmente estable se convierte en meta-estable en un sistema estoc´ astico. Si bien las fluctuaciones acabar´an extinguiendo al sistema, el estado de coexistencia puede ser de larga duraci´on. En tiempo finito el sistema puede mostrar una din´ amica de oscilaciones err´aticas que aleje temporalmente al sistema de la 12

soluci´on estacionaria, retras´ andola considerablemente [14]. Seg´ un se afirma tambi´en en [14], simulaciones de Monte Carlo muestran que las densidades de poblaci´on se aproximan a las densidades de coexistencia a trav´es de oscilaciones amortiguadas a partir de la condici´ on inicial. As´ı las cosas, el destino fatal del ecosistema regido por el modelo de Lotka-Volterra parece inevitable. Y aunque esto puede ocurrir en la naturaleza, la simplificaci´on extrema del modelo obvia algunas cuestiones que puedan dar un vuelco a este final. Es el caso de la introducci´ on de espacialidad en el sistema, implicando este una posibilidad de refugio o acogedor de flujos migratorios para la supervivencia. Existen resultados [13],[16],[17] que demuestran que estas modificaciones permiten de nuevo una situaci´on de coexistencia entre especies. M´as detalles sobre un tratamiento del modelp de Lotka-Volterra a partir de las fluctuaciones intr´ınsecas puede verse [15],[3].

3.3 Derivaci´on de las ecuaciones de Lotka-Volterra. Como conclusi´ on al estudio estoc´ astico interno de este modelo vamos a deducir las ecuaciones de Lotka-Volterra macrosc´ opicas (deterministas) utilizadas en la secci´on primera. Vamos a trabajar dos m´etodos, el primero a partir de las ecuaciones de evoluci´on de los momentos diferenciales y el segundo desde el punto de vista de la aproximaci´on de Van Kampen. - Momentos diferenciales: Nos interesa obtener las ecuaciones de evoluci´on del primer momento y de la varianza (en su versi´on intensiva): dni 1 = a1 (hni i) + a001 (hni i)σn2 i , dt 2

(35)

dσn2 i = a2 (hni i) + 2a01 (hni i)σn2 i , (36) dt donde se ha usado la notaci´ on n1 ≡ x, n2 ≡ y por simplicidad. Necesitamos conocer las magnitudes ap , a00p (con p = 1, 2) para x e y. Se calculan a partir de ap (X) =

X

rρX ω(X + rρX |X), ap (Y ) =

ρ

X

rρY ω(Y + rρY |Y ),

(37)

ρ

el sumatorio se extiende a todas las reacciones. La variable discreta se relaciona con la continua mediante ap (X, Y ) = Ωap (x, y), donde Ω es el volumen del sistema. Estos son: ap (X) = (−1)p k¯0 AX + (−1)p k¯1 XY a0p (X) = (−1)p k¯0 A + (−1)p k¯1 Y

→ ap (x) = (−1)p k0 ax + (−1)p k1 xy,

(38)

→ a0p (x) = (−1)p k0 a + (−1)p k1 y,

(39)

13

a00p (X) = 0 → a00p (x) = 0.

(40)

ap (Y ) = (−1)p k¯1 xy + (−1)p k¯2 y → ap (y) = (−1)p k1 xy + (−1)p k2 y, a0p (Y

)=

(−1)p k¯1 x

+

(−1)p k¯2

a00p (Y

p

p

→ ap (y) = (−1) k1 x + (−1) k2 ,

) = 0 → ap (y) = 0.

(41) (42) (43)

Con k0 ≡ k¯0 Ω, k1 ≡ k¯1 Ω, k2 ≡ k¯2 .

(44)

Sustituyendo adecuadamente los resultados de las ecuaciones (38-43) obtenemos el sistema dhxi = k0 ahxi − k1 hxyi, (45) dt dhyi = k1 hxyi − k2 hyi, (46) dt dσx2 = k0 ahxi + k1 hxyi + 2(k0 a − k1 hyi)σx2 , dt

(47)

dσy2 = k1 hxyi + k2 hyi + 2(k1 hxi − k2 )σy2 . dt

(48)

Ahora tenemos que dar un paso ciertamente injustificado: suponer que no hay correlaci´on entre la poblaci´ on de presa y de depredadores, de manera que su covarianza sea nula, cov(x, y) = 0 y en consecuencia hxyi = hxihyi. Entonces, considerando que los valores medios se aproximan a los deterministas y redefiniendo r ≡ k0 a, α = β ≡ k1 , µ ≡ k2 , tenemos las ecuaciones de Lotka-Volterra presa-depredador obtenidas a partir de las consideraci´ on de las fluctuaciones internas: dhxi = rhxi − αhxihyi, dt

(49)

dhyi = αhxihyi − µhyi. (50) dt Estas ecuaciones bajo esta aproximaci´on de correlaci´on nula es lo que se conoce tambi´en como ecuaciones o aproximaci´ on de campo-medio [15]. Vamos a comprobar que, aunque hayamos obtenido el modelo macrosc´opico, son un tanto problem´aticas. En el estado estacionario encontramos los mismos puntos fijos que en el sistema determinista: (hxist , hyist ) = (0, 0), (hxist , hyist ) = (µ/α, r/α). (51) Sustituyendo los puntos fijos en las ecuaciones (47-48) y haciendo dσx2 /dt = dσy2 /dt = 0 obtenemos la varianza para cada punto estacionario. Para el punto de extinci´on (0, 0) obtenemos que la varianza es nula σx2 = 0, σy2 = 0. 14

(52)

Las fluctuaciones en el punto de extinci´on son nulas lo cual caracteriza al punto (0, 0) como estable y absorbente. Es imposible escapar de ´el, no hay fluctuaci´on que saque al sistema de este sumidero. Este resultado vuelve a estar en discrepancia con la descripci´on puramente determinista pero concuerda perfectamente con la predicci´on de la soluci´on estacionaria de la ecuaci´ on maestra. Si ahora evaluamos para el punto de coexistencia (µ/α, r/α) ambos miembros se anulan id´enticamente y no podemos obtener informaci´on sobre la varianza. Evidentemente esto no es bueno y es consecuencia de la aproximaci´on de covarianza cero. Un an´ alisis m´as riguroso exigir´ıa una hip´otesis menos restrictiva y nos llevar´ıa a considerar momentos de orden mayor y en consecuencia a ampliar y enrevesar el sistema de ecuaciones diferenciales. - Aproximaci´ on de Van Kampen: Es muy instructivo aunque algo tedioso trabajar con un sistema bidimensional donde desarrollar la aproximaci´on de Van Kampen. Es lo que aqu´ı voy a hacer. No entrar´e en detalles sobre la aproximaci´on (v´ease [18]). Lo que s´ı que es preciso remarcar es que aqu´ı se va a utilizar simplemente como forma alternativa para la deducci´ on de las ecuaciones macrosc´opicas sin entrar al an´alisis de la estabilidad, puesto que este m´etodo tiene una validez local: describe las fluctuaciones de la variable relevante alrededor de su valor medio. Es por tanto inadecuado para estudiar cuestiones de estabilidad relativa, como es el caso de nuestro sistema, donde tenemos dos puntos fijos de diferente estabilidad. Vamos a obtener de nuevo las ecuaciones de Lotka-Volterra en el orden macrosc´opico. Siguiendo a van Kampen [19], en un sistema bidimensional suponemos que las variables relevantes son resultado de un t´ermino macrosc´opico del orden del volumen del sistema Ω y un t´ermino fluctuante, aleatorio, de orden Ω1/2 , as´ı: X = Ωφ(t) + Ω1/2 ξ,

(53)

Y = Ωψ(t) + Ω1/2 η,

(54)

con la relaci´ on entre probabilidades: P (X, Y ; t) = Π(ξ, η; t).

(55)

Ahora Π(ξ, η; t) designa la probabilidad asociada a las nuevas variables aleatorias, las fluctuaciones respecto al valor macrosc´opico de X e Y . Los operadores de traslaci´on puede escribirse como k ∂ k2 ∂ 2 k ˆX E ≈1+ √ + + ... , (56) Ω ∂ξ 2Ω ∂ξ 2 2 2 ˆ k ≈ 1 + √k ∂ + k ∂ + ... E (57) Y Ω ∂η 2Ω ∂η 2 Entonces, uno parte de la ecuaci´ on maestra escrita en funci´on de los operadores de traslaci´on, ecuaci´ on (31), y reescribe sus miembros de acuerdo a las nuevas transformaciones establecidas. As´ı, el miembro izquierdo, el de ∂t P (X, Y ; t), pasar´ıa a ser

∂Π(ξ, η; t) dφ ∂Π(ξ, η; t) dψ ∂Π(ξ, η; t) ∂P (X, Y ; t) = − Ω1/2 − Ω1/2 . ∂t ∂t dt ∂ξ dt ∂η 15

(58)

El miembro derecho es m´ as enrevesado y por no cargar el escrito con largos desarrollos expondr´e la situaci´ on inicial y final. Los distintos miembros que tenemos para transformar son k¯1 XY = k¯1 (Ωφ + Ω1/2 ξ)(Ωψ + Ω1/2 η) = k1 (Ωφψ + Ω1/2 φη + Ω1/2 ψξ + ξη), k¯0 AX = k¯0 A(Ωφ + Ω1/2 ξ) = k0 a(Ωφ + Ω1/2 ξ), k¯2 Y = k¯2 (Ωψ + Ω

1/2

η) = k2 (Ωψ + Ω

1/2

η),

(59) (60) (61)

ˆX E ˆ −1 − 1) = Ω−1/2 ∂η − Ω−1/2 ∂ξ + 1 Ω−1 ∂ηη − Ω−1 ∂ξ ∂η − 1 Ω−1 ∂ξξ + O(Ω−3/2 ), (62) (E y 2 2 Ahora hay que multiplicar estos resultados de acuerdo a como aparecen en la ecuaci´on maestra, es decir ˆ −1 − 1)k¯0 XP (X, Y ; t), (E (63) X ˆY − 1)k¯2 Y P (X, Y ; t), (E

(64)

ˆX E ˆ −1 − 1)k¯1 XY P (X, Y ; t). (E Y

(65)

Sustituyendo, por supuesto, adem´as P (X, Y ; t) por Π(ξ, η; t). Una vez desarollados estos miembros (63-65) tenemos que cortar el desarrollo y mantener solo el orden Ω0 y Ω1/2 . Posteriormente hay que igualar por potencias de igual orden en Ω con el miembro izquierdo de la ecuaci´ on maestra [ecuaci´on (58)]. El orden Ω0 nos da la ecuaci´on de Fokker-Planck para la distribuci´ on bivariada de probabilidad de las fluctuaciones y el orden Ω1/2 nos da las ecuaciones macrosc´opicas/deterministas de Lotka-Volterra. Como no se va a hacer un an´ alisis posterior de las fluctuaciones, nos centramos en el orden Ω1/2 . As´ı tenemos: −

dφ ∂Π dψ ∂Π ∂Π ∂Π ∂Π ∂Π − = −k0 aφ + k2 ψ − k1 φψ + k1 φψ dt ∂ξ dt ∂η ∂ξ ∂η ∂η ∂ξ

Separando en dos ecuaciones, una para ∂ξ Π y otra para ∂η Π, llegamos definitivamente a las ecuaciones macrosc´ opicas de Lotka-Volterra: dφ = k0 aφ − k1 φψ, dt

(66)

dψ = −k0 aφ + k1 φψ. (67) dt Donde φ es el valor macrosc´ opico asociado a la presa y ψ es el valor macrosc´opico asociado al depredador. Evidentemente, los puntos fijos de este sistema de ecuaciones son exactamente los mismos que los obtenidos en los casos previos: sistema determinista y sistema deducido a partir de los momentos. Si uno calcula ahora los valores medios y momentos de orden superior, hξi, hηi, hξ 2 i, hη 2 i, ... llegar´a a las mismas conclusiones que en el apartado anterior. Las fluctuaciones entorno al punto de extinci´on son nulas y no podemos obtener informaci´ on de los momentos de orden superior respecto al punto de coexistencia. 16

4 Resumen y conclusiones. En este trabajo se ha hecho un an´alisis b´asico del modelo de Lotka-Volterra, de amplio uso en la teor´ıa de din´ amica de poblaciones. En la primera parte hemos analizado cuestiones b´ asicas propias de la teor´ıa de sistemas din´amicos determinista. Hemos hecho ´enfasis en la edificaci´ on del modelo de Lotka-Volterra de presa-depredador a partir de alg´ un modelo m´ as simple y de ciertas hip´otesis para clarificar su significado. Estudiamos las soluciones estacionarias del sistema, puntos fijos, y analizamos los tipos de soluciones. El modelo determinista presenta una soluci´on de coexistencia de especies, matem´aticamente caracterizada por un punto fijo de estabilidad neutra y dando lugar a trayectorias oscilantes para las poblaciones de presa y depredador. Las otras soluciones implican la desaparici´ on de una (el depredador) o de ambas especies, extinci´on total. Esta soluci´ on es un punto-silla y por tanto inestable. A continuaci´on se ha propuesto la construcci´ on del mismo modelo a partir del tratamiento de las fluctuaciones internas del sistema, es decir, considerando con una cierta probabilidad de ocurrencia los distintos sucesos elementales que constituyen las interacciones presa-depredador. Se ha obtenido la ecuaci´ on maestra para el conjunto de reacciones de Lotka-Volterra y se ha hallado la soluci´ on estacionaria. La situaci´on del estado estacionario respecto al modelo determinista supone una serie discrepancia: la probabilidad es total para la soluci´on de extinci´ on total. Las causas son la consideraci´on de variables discretas en el formalismo y la presencia de fluctuaciones importantes que acabar´an tarde o temprano con la soluci´ on pseudo-estacionaria de estabilidad neutra, la coexistencia. Finalmente, se han deducido las ecuaciones macrosc´opicas de Lotka-Volterra depredador-presa a partir de los momentos diferenciales y de la expansi´on de van Kampen. El hecho de cuadrar las ecuaciones obtenidas con las deterministas exige asumir covarianza nula entre presa y depredador, una hip´ otesis injustificada que afecta posteriormente a la obtenci´on de resultados relativos al punto fijo de coexistencia.

References [1] Wikipedia, http://en.wikipedia.org/wiki/Lotka%E2%80%93Volterra_equation [2] Guillermo Abramson, La matem´ atica de los sistemas biol´ ogicos. Universidad Nacional de Cuyo. Centro At´ omico Bariloche, Instituto Balseiro y CONICET [3] H.K. Leung, Stochastic Transients of Noisy Lotka-Volterra Model Chinese Journal of Physics. Vol 29, No. 6. [4] K. Sigmund, Kolmogorov and population dynamics Faculty for Mathematics, University of Vienna, and International Institute for Applied Systems Analysis, Austria [5] G. Belcredi, Modelo presa - depredador

17

[6] L.C.M Miranda, C.A.S. Lima, On the logistic modeling and forecasting of evolutionary processes: Application to human population dynamics Technological Forecasting & Social Change 77 (2010) 699-771 [7] http://www.ucl.ac.uk/~ucess29/page2/MATH3506Chapters/ /MATH3506Chapter4.pdf [8] J.H.P. Dawes, M.O. Souza A derivation of Holling’s type I, II and III functional responses in predator-prey systems. [9] C.M. Evans, G.L. Findley Analytic Solutions to the Lotka-Volterra Model For Sustained Chemical Oscillations. [10] S. Baigent Lotka-Volterra Dynamics - An introduction. 2010 [11] H. Haken Synergetics, An Introduction: Nonequilibrium Phase Transitions and SelfOrganization in Physics, Chemistry and Biology. Springer-Verlag 1977 [12] U. Dobramysi On the Relaxation Dynamics of Disordered Systems. Dissertation. Doctor of Philosophy in Physics [13] M. Mobilia, I. T. Georgiev, U. C. T¨auber, Phase Transitions and Spatio-Temporal Fluctuations in Stochastic Lattice Lotka-Volterra Models. Journal of Statistical Physics 128 447-483 (2007). [14] M. Parker, A. Kamenev Extinction in Lotka-Volterra model 2009. [15] R. Toral Nonlineal phenomena in biology 2012. [16] U.C. T¨ auber Population oscillations in spatial stochastic Lotka-Volterra models: A field-theoretic perturbational analysis [17] M.J. Keeling, H.B. Wilson, S.W. Pascala (200) Science 290 1758-1761 y G. Marion Stochastic modelling in mathematical biology: Predator-Prey systems and epidemic models [18] J. de la Rubia Fluctuaciones en Sistemas Din´ amicos M´aster en F´ısica de Sistemas Complejos, UNED [19] N.G. van Kampen Stochastic Processes in Physics and Chemistry Third edition, 2007

18