Modelo de Lotka-Volterra

Modelo de Lotka-Volterra Alfonso de Miguel July 1, 2014 Abstract En este trabajo se plantea un estudio básico de modelos matemáticos usados tı́picamente en el estudio de ecologı́a y dinámica de poblaciones. Daremos un repaso al archiconocido modelo de Lotka-Volterra desde el punto de vista de la teorı́a de sistemas dinámicos (determinista). A continuación se trabajará su deducción y análisis considerando los procesos elementales del sistema (procesos de salto) para deducir una ecuación maestra del sistema. Compararemos el modelo determinista y este último, basado en el análisis de las fluctuaciones internas. 1 Introducción. La incursión de la fı́sica y de la matemática en terrenos no convencionales es desde hace largo tiempo una barrera más 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étodos y conclusiones. La generalidad que ofrece la teorı́a de sistemas dinámicos, ası́ como la teorı́a de procesos estocásticos, nos permite modelar diversos fenómenos caracterı́sticos de cada campo de estudio y dotarlos de un cuerpo matemático y analı́tico, con el fin de lograr una comprensión más profunda y una mayor capacidad de predicción y análisis del fenómeno en cuestión. 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ámica de poblaciones en general. Dado su atractivo, el modelo aparece en numerosos textos y artı́culos, bien sea enfocando el análisis desde el punto de vista de la teorı́a de sistemas dinámicos (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ácter más general [2]. También existen, por supuesto, textos dedicados a un estudio estocástico 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ón desde esta última perspectiva. 1 Lo que aquı́ se pretende es realizar un estudio básico del modelo de Lotka-Volterra en dos dimensiones. En la primera sección se presenta el sistema dinámico y se repasan las caracterı́sticas básicas del modelo: construcción del modelo, significado, soluciones, puntos fijos y estabilidad. En la segunda sección entramos a considerar ya la cuestión estocástica: tomando el modelo como consecuencia de la ocurrencia de una serie de procesos elementales de variable discreta llegaremos a la deducción de su ecuación maestra, abordaremos su solución estacionaria y de ahı́ deduciremos las ecuaciones macroscópicas de Lotka-Volterra. Compararemos los resultados entre el modelo determinista y el modelo basado en las fluctuaciones intrı́nsecas. 2 Modelos de dinámica poblacional: Lotka-Volterra. 2.1 Construcción de modelos. Vamos a mostrar o, mejor dicho, a construir el modelo de Lotka-Volterra de depredadorpresa paso a paso para ver su razón de ser como modelo descriptivo (y muy simplificado) de la dinámica de poblaciones de un ecosistema 1 . Sea entonces un ecosistema donde habita una única especie y cuya población (densidad) en el hábitat designamos por la variable x. Parece razonable postular que la evolución en la población x de la especie presa en un tiempo t + dt será proporcional al tamaño de la población en t, x(t), estableciendo la igualdad un multiplicador r tal que dx = rx, (1) dt donde uno puede definir r como el número de nacimientos menos el número de muertes de miembros de la especie (crecimiento natural), esto es: r es la tasa neta de crecimiento per cápita. 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én de gran popularidad, ofrece una visión más 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érmino cuadrático negativo en x que actúa como regulador de la población, 1 Para ver una exposición histórica del desarrollo del modelo puede consultarse [4], [5] La capacidad de carga de una especie biológica en un ambiente es el tamaño máximo de población que el ambiente puede soportar indefinidamente en un periodo determinado, teniendo en cuenta el alimento, agua, hábitat, y otros elementos necesarios disponibles en ese ambiente. Más 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ón no es la única manera de describir la saturación o limitación al crecimiento de una población y dependiendo del problema concreto la modelización puede requerir una función 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ón o de otros factores. Ası́ que, desde una perspectiva más general, uno puede escribir la evolución en la población de una especie como dx = ρ(x)x. dt (3) Consideremos la situación donde tenemos una especie depredadora pero no hay presa. En ausencia de especies inferiores en la cadena de alimentaria se supone una extinción 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ón entre ambas es de presa-depredador, la evolución dinámica del número de miembros de estas especies ya no dependerá únicamente de la propia especie sino de cómo sea la interacción con la rival. De este modo, la presencia del depredador tendrá un efecto negativo en la tasa de crecimiento per cápita 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ón per cápita de la presa debido a la depredación y a la tasa de crecimiento del depredador debido al consumo per cápita de la presa como σ(x). Por lo tanto, de manera general, se tendrá este sistema de ecuaciones diferenciales ẋ = ρ(x) − φ(x, y), x ẏ = σ(x) − µ, y (5) (6) Se asume que tanto φ como σ son funciones suaves de sus argumentos. El coeficiente que acompaña a φ(x, y) es negativo, ya que la interacción presa-depredador sobre la presa se traduce en un perjuicio, en un decrecimiento del número de miembros de esta especie. Por contra, para el depredador esta interacción es beneficiosa en términos de su crecimiento y σ(x) va con signo positivo3 Reajustando los términos, 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é ser iguales aunque se trate de la misma interacción presa-depredador. De todos modos, en lo posterior se asumirán iguales. 3 dy = y[σ(x) − µ] = G(x, y). (8) dt que es una versión general de modelo de Lotka-Volterra bidimensional de depredadorpresa con crecimiento logı́stico de la presa. Uno puede bajo esta generalidad, y con intención de refinar el modelo, seguir aportando condiciones lógicas o esperables a satisfacer para las funciones φ(x, y) y σ(x) (véase por ejemplo [7]), pero opto ya por centrar la cuestión en un modelo concreto que será sobre el que luego haremos el trabajo estocástico. 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ápita (como en la ecuación (1)); ω1 es la constante de acoplo presa-depredador; nótese que φ(x, y) = αy, σ(x) = βx. En el caso de la tasa de desaparición per cápita de la presa, función φ(x, y), esta es proporcional a la población 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ápita de la presa. El producto αx representa la llamada respuesta funcional del depredador, esta describe cómo 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ápita de la presa, función σ(x), es proporcional a la población de presas en un tiempo t. La constante β recibe el nombre de eficiencia de conversión, esto es, la habilidad de los depredadores para convertir presa en tasa de crecimiento per cápita adicional para la población de depredadores. El producto βx se denomina también respuesta numérica de la población depredadora. Resumamos, para que queden claras las hipótesis 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ático (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ás realistas y concretas, véase [8]. En el caso de la ecuación (9), esta es lineal con x. 4 - El efecto de los depredadores es reducir el número 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ón de depredadores de forma proporcional a su propia población. Esto es una simplificación considerable de algún mecanismo complejo que convertirı́a carne de presas en hijos de depredadores. Evidentemente estamos en una situación muy simplista. Aparte de la no inclusión de términos logı́sticos o algún otro tipo de saturación 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ás, cuantas más presas hay, más comen. También se da por hecho que todos los depredadores tienen las mismas opciones de atrapar a sus presas. Hay modelos con difusión, 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ón de las nulclinas de x e y nos dan los puntos fijos del sistema, ası́ que tenemos: (x01 , y01 ) = (0, 0), (13)   µ r , . (14) (x02 , y02 ) = β α La primera solución implica la extinción total de ambas especies. La segunda significa una situación de coexistencia de ambas especies. Aparte uno tiene la solución (∞, 0) que acaece en ausencia de depredadores como condición 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ón del jacobiano   r − αy −αx , (15) J(x, y) =  β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 √ , = ±i rµ. σ J (17) β α Los autovalores son imaginarios puros, en consecuencia estamos ante un centro, este tipo de punto de equilibrio es estable pero no asintóticamente estable, se habla de estabilidad neutra o marginal. Las trayectorias (x(t), y(t)) que uno obtiene mediante esta análisis lineal implican soluciones periódicas entorno al punto fijo de coexistencia de especies. Por otro lado, podemos encontrar las ecuaciones de las órbitas del modelo presapredador. Tomando las ecuaciones (9) y (10), despejamos en ambas dt e igualamos, reajustando términos se llega a dy y(βx − µ) = dx x(r − αy) (18) las órbitas del sistema para x, y ≡ 0 son las curvas soluciones a esta ecuación 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 órbitas del sistema son la familia de curvas definidas por la ecuación anterior6 La función H(x, y) es constante a lo largo de las trayectorias (x(t), y(t)) solución del sistema de ecuaciones diferenciales de Lotka-Volterra (9-10). La ecuación (19) muestra que tenemos una integral o constante del movimiento. De hecho, el sistema de LotkaVolterra puede plantearse como un sistema dinámico Hamiltoniano [9],[10], siendo el propio H(x, y) el Hamiltoniano de este sistema. - Conservación de los promedios: Supóngase condiciones iniciales positivas x(0) = x0 > 0 e y(0) = y0 > 0. Si T es el periodo de una órbita cerrada que pasa por (x0 , y0 ), entonces ẋ = r − αy → ln x(t) − ln x(0) = x Z T 0 [r − αy(t)]dt → 0 = rT − α Z T y(t)dt, 0 usando la propiedad de periodicidad de las órbitas. 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ón de que estamos ante una solución periódica véase con más detalle []. 7 Haciendo lo propio para la ecuación de los depredadores ẏ = βx − µ → ln y(t) − ln y(0) = y Z T 0 [βx(t) − µ]dt → 0 = −µT + β Y análogamente hxiT = 1 T Z T x(t)dt = 0 µ . β Z T x(t)dt. 0 (21) El valor mediado de las variables de estado (x(t), y(t)) sobre el periodo de una órbita 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ón anterior, y cualquiera de sus variaciones mencionadas y referenciadas, son al fin y al cabo sistemas de ecuaciones diferenciales deterministas y su estudio matemático está basado en el uso de conceptos y resultados de la teorı́a de sistemas dinámicos convencional (ecuaciones diferenciales, teorı́a de bifurcaciones, sistemas Hamiltonianos, etc.). Los modelos tradicionales de ecologı́a matemática describen en forma diferencial la evolución de las densidades. Como se ha visto en la sección previa, en general, podemos 8 construir estos modelos sin necesidad de conocer al detalle la interacción entre los individuos que conforman la población. Se supone, no obstante, que cada término tiene su fundamento en las interacciones puntuales y que representan efectos macroscópicos que involucran a muchos individuos, siendo las ecuaciones de evolución suaves y continuas. Sin embargo, las poblaciones son discretas y las interacciones están individualizadas. Por medio de un análisis de todos estos procesos de interacción a escala microscópica deberı́a ser posible la obtención de las leyes macroscópicas que obtuvimos, surgiendo estas como un lı́mite macroscópico, un coarse-graining, de los fenómenos observados. Algo que nos enseña la teorı́a de fluctuaciones en sistemas dinámicos es esta posibilidad de obtener ecuaciones deterministas/macroscópicas a partir de una consideración estocástica del proceso o fenómeno 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ón uno puede describir el modelo de Lotka-Volterra en términos de reacciones de la siguiente manera: X + A → 2X, (22) X +Y (23) Y → 2Y, → ∅. (24) Donde cada reacción se produce a una tasa ki=0,1,2 . X e Y siguen designando respectivamente a la población de presa y depredador, pero ahora son variables discretas: representan en el número de miembros de dicha especie. La primera reacción representa el proceso de crecimiento, más concretamente de autorreplicación, de la especie presa a partir de una fuente de alimento A que se toma constante7 La segunda reacción representarı́a el proceso de depredación. La interacción de una especie presa y una depredadora dan lugar a la desaparición de una tal presa y al surgimiento de una depredadora8 Por último, la reacción más 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ásticas y su mera especificación en un instante de tiempo t concreto no aporta la misma información que en una situación determinista. Es por ello que debemos trabajar mediante la función de probabilidad del sistema: P (X, Y ; t) denotará entonces la probabilidad conjunto de que en un instante de 7 Tomar constante la alimentación de esta especie serı́a equivalente a la consideración de crecimiento explosivo de la especie à la Malthus. Por otro lado, puede verse aquı́ que no se plantea ningún tipo de reproducción sexual ni gestación sino que sencillamente el organismo X se autorreplica a partir de la ingesta de A. 8 Como ya mencionamos en la sección previa, esta hipótesis también 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ón maestra que buscamos es una ecuación que nos informa sobre la evolución 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ón 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ón era (X − 1, Y ) en el instante t y una presa se reprodujo (tuvo lugar la reacción primera, ecuación (22)) durante dt. Suceso elemental 3: La población era (X, Y + 1) en el instante t y un depredador murió (reacción tercera, ecuación (24)) durante dt. Suceso elemental 4: La población era (X + 1, Y − 1) en el instante t y un predador devoró a una presa (reacción segunda, ecuación (23)) durante dt. El hecho de que ocurra una reacción u otra es independiente. Además, 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ón o reproducción de la presa X, que es proporcional al número de presas X en el instante t y a la constante de reacción 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ón depredador-presa, que es proporcional al número de presas X y de depredadores Y en el instante t y a la constante de reacción k̄1 (la interacción está regida por la misma constante de acoplo desde el punto de vista de ambos). Y por último Y k̄2 dt, (28) es la probabilidad de muerte de un depredador, y es proporcional al número 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úmero de especies que habrı́a en el instante t considerado. 10 3.2 Ecuación Maestra y solución estacionaria. Sabiendo la forma de las probabilidades de transición, 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á conformado por la probabilidad de estar en un estado (m, n) por la probabilidad de transición de ir de un estado (m, n) a otro (m′ , n′ ). Desechando términos de orden O(dt2 ) y superiores, dividiendo por dt y tomando el lı́mite dt → 0 tenemos la ecuación de evolución para la probabilidad conjunta de X e Y , la ecuación 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ás compacto escribimos esta ecuación por medio de los operadores de traslación: ∂P (X, Y ; t) −1 = (ÊX −1)X k̄0 P (X, Y ; t)+(ÊY −1)k̄2 Y P (X, Y ; t)+(ÊX ÊY−1 −1)XY k̄1 P (X, Y ; t). ∂t (31) La bidimensionalidad del sistema complica enormemente, si no la hace imposible, la resolución de la ecuación maestra a tiempo t. Tampoco intentando la construcción de una función generatriz facilita el problema. La solución estacionaria verifica que ∂Pst (X, Y ) = 0. ∂t (32) Ası́: −1 0 = (ÊX − 1)X k̄0 P (X, Y ; t) + (ÊY − 1)k̄2 Y P (X, Y ; t) + (ÊX ÊY−1 − 1)XY k̄1 P (X, Y ; t). (33) Extraer de aquı́ la solución no parace nada claro tampoco pero recurro a un argumento de Haken [11] (recogido también 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ón maestra para enfatizar la consideración de los procesos elementales dados. Una forma más directa para construir las probabilidades de transición habrı́a sido seguir el procedimiento usado para las reacciones quı́micas [18], por medio de la combinatoria. 11 Sustituyéndose en la ecuación (33), esta se verifica. Pero, ¿es única esta solución? Aquı́ viene el segundo punto: todos los puntos (X, Y ) del sistema están conectados por al menos una trayectoria con cualquier otro punto (X ′ , Y ′ ) 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éase la figura 4. Ahora, este resultado puede parecer inesperado. Nuestro sistema determinista estudiado en la sección primera ofrecı́a dos puntos fijos: la solución de coexistencia y la extinción total. Aquı́ no hay ni rastro de la primera y la única perspectiva es ciertamente apocalı́ptica. ¿Cuál es la explicación a esta discrepancia entre el modelo estocástico y el determinista? La solución de existencia era un estado punto-silla (inestable) en el modelo determinista, pero la de coexistencia implicaba órbitas cerradas en el espacio de fases del sistema y el problema radica en su estabilidad. Ası́ es, obtuvimos una solución de estabilidad neutra, esto significa que en presencia de fluctuaciones, hecho inevitable en un sistema estocástico, las trayectorias cerradas van a verse perturbadas en mayor o menor medida y a la vista de la solución estacionaria obtenida, en t → ∞ la probabilidad de que haya una fluctuación que lleve al sistema a la extinción tiende a uno. Debido a la naturaleza discreta del sistema estocástico, las fluctuaciones en el número de especies pueden ser determinantes para la extinción si el número de especies se va haciendo pequeño. Además, este resultado significa que cualquier sistema con un número finito de partı́culas siempre alcanzará el estado de extinción, pero la escala temporal de extinción puede ser muy grande para sistemas razonablemente grandes [12]. Esta caracterı́stica está ausente en el llamado modelo de campo-medio, la descripción puramente determinista, ya que las densidades de población pueden ser arbitrariamente pequeñas sin que el destino de la población se vea afectado dramáticamente. El punto fijo de coexistencia de especies marginalmente estable se convierte en meta-estable en un sistema estocástico. Si bien las fluctuaciones acabarán extinguiendo al sistema, el estado de coexistencia puede ser de larga duración. En tiempo finito el sistema puede mostrar una dinámica de oscilaciones erráticas que aleje temporalmente al sistema de la 12 solución estacionaria, retrasándola considerablemente [14]. Según se afirma también en [14], simulaciones de Monte Carlo muestran que las densidades de población se aproximan a las densidades de coexistencia a través de oscilaciones amortiguadas a partir de la condición 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ón extrema del modelo obvia algunas cuestiones que puedan dar un vuelco a este final. Es el caso de la introducción 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ón de coexistencia entre especies. Más detalles sobre un tratamiento del modelp de Lotka-Volterra a partir de las fluctuaciones intrı́nsecas puede verse [15],[3]. 3.3 Derivación de las ecuaciones de Lotka-Volterra. Como conclusión al estudio estocástico interno de este modelo vamos a deducir las ecuaciones de Lotka-Volterra macroscópicas (deterministas) utilizadas en la sección primera. Vamos a trabajar dos métodos, el primero a partir de las ecuaciones de evolución de los momentos diferenciales y el segundo desde el punto de vista de la aproximación de Van Kampen. - Momentos diferenciales: Nos interesa obtener las ecuaciones de evolución del primer momento y de la varianza (en su versión intensiva): dni 1 = a1 (hni i) + a′′1 (hni i)σn2 i , dt 2 (35) dσn2 i (36) = a2 (hni i) + 2a′1 (hni i)σn2 i , dt donde se ha usado la notación n1 ≡ x, n2 ≡ y por simplicidad. Necesitamos conocer las magnitudes ap , a′′p (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 a′p (X) = (−1)p k̄0 A + (−1)p k̄1 Y → ap (x) = (−1)p k0 ax + (−1)p k1 xy, (38) → a′p (x) = (−1)p k0 a + (−1)p k1 y, (39) 13 a′′p (X) = 0 → a′′p (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, a′p (Y p p p p ) = (−1) k̄1 x + (−1) k̄2 → ap (y) = (−1) k1 x + (−1) k2 , a′′p (Y ) = 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ón entre la población 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ón de las fluctuaciones internas: dhxi = rhxi − αhxihyi, dt (49) dhyi = αhxihyi − µhyi. (50) dt Estas ecuaciones bajo esta aproximación de correlación nula es lo que se conoce también como ecuaciones o aproximación de campo-medio [15]. Vamos a comprobar que, aunque hayamos obtenido el modelo macroscópico, son un tanto problemáticas. 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ón (0, 0) obtenemos que la varianza es nula σx2 = 0, σy2 = 0. 14 (52) Las fluctuaciones en el punto de extinción son nulas lo cual caracteriza al punto (0, 0) como estable y absorbente. Es imposible escapar de él, no hay fluctuación que saque al sistema de este sumidero. Este resultado vuelve a estar en discrepancia con la descripción puramente determinista pero concuerda perfectamente con la predicción de la solución estacionaria de la ecuación maestra. Si ahora evaluamos para el punto de coexistencia (µ/α, r/α) ambos miembros se anulan idénticamente y no podemos obtener información sobre la varianza. Evidentemente esto no es bueno y es consecuencia de la aproximación de covarianza cero. Un análisis más riguroso exigirı́a una hipótesis 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ón de Van Kampen: Es muy instructivo aunque algo tedioso trabajar con un sistema bidimensional donde desarrollar la aproximación de Van Kampen. Es lo que aquı́ voy a hacer. No entraré en detalles sobre la aproximación (véase [18]). Lo que sı́ que es preciso remarcar es que aquı́ se va a utilizar simplemente como forma alternativa para la deducción de las ecuaciones macroscópicas sin entrar al análisis de la estabilidad, puesto que este método 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ópico. Siguiendo a van Kampen [19], en un sistema bidimensional suponemos que las variables relevantes son resultado de un término macroscópico del orden del volumen del sistema Ω y un término fluctuante, aleatorio, de orden Ω1/2 , ası́: X = Ωφ(t) + Ω1/2 ξ, (53) Y = Ωψ(t) + Ω1/2 η, (54) con la relación 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ópico de X e Y . Los operadores de traslación puede escribirse como k ∂ k2 ∂ 2 k + + ... , (56) ÊX ≈1+ √ Ω ∂ξ 2Ω ∂ξ 2 k2 ∂ 2 k ∂ + + ... (57) ÊYk ≈ 1 + √ Ω ∂η 2Ω ∂η 2 Entonces, uno parte de la ecuación maestra escrita en función de los operadores de traslación, ecuación (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 ∂P (X, Y ; t) ∂Π(ξ, η; t) dφ ∂Π(ξ, η; t) dψ ∂Π(ξ, η; t) = − Ω1/2 − Ω1/2 . ∂t ∂t dt ∂ξ dt ∂η 15 (58) El miembro derecho es más enrevesado y por no cargar el escrito con largos desarrollos expondré la situación 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) 1 1 (ÊX Êy−1 − 1) = Ω−1/2 ∂η − Ω−1/2 ∂ξ + Ω−1 ∂ηη − Ω−1 ∂ξ ∂η − Ω−1 ∂ξξ + O(Ω−3/2 ), (62) 2 2 Ahora hay que multiplicar estos resultados de acuerdo a como aparecen en la ecuación maestra, es decir −1 (ÊX − 1)k̄0 XP (X, Y ; t), (63) (ÊY − 1)k̄2 Y P (X, Y ; t), (64) (ÊX ÊY−1 − 1)k̄1 XY P (X, Y ; t). (65) Sustituyendo, por supuesto, además 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ón maestra [ecuación (58)]. El orden Ω0 nos da la ecuación de Fokker-Planck para la distribución bivariada de probabilidad de las fluctuaciones y el orden Ω1/2 nos da las ecuaciones macroscópicas/deterministas de Lotka-Volterra. Como no se va a hacer un análisis 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ópicas de Lotka-Volterra: dφ = k0 aφ − k1 φψ, dt (66) dψ = −k0 aφ + k1 φψ. (67) dt Donde φ es el valor macroscópico asociado a la presa y ψ es el valor macroscópico 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 las mismas conclusiones que en el apartado anterior. Las fluctuaciones entorno al punto de extinción son nulas y no podemos obtener información de los momentos de orden superior respecto al punto de coexistencia. 16 4 Resumen y conclusiones. En este trabajo se ha hecho un análisis básico del modelo de Lotka-Volterra, de amplio uso en la teorı́a de dinámica de poblaciones. En la primera parte hemos analizado cuestiones básicas propias de la teorı́a de sistemas dinámicos determinista. Hemos hecho énfasis en la edificación del modelo de Lotka-Volterra de presa-depredador a partir de algún modelo más simple y de ciertas hipótesis para clarificar su significado. Estudiamos las soluciones estacionarias del sistema, puntos fijos, y analizamos los tipos de soluciones. El modelo determinista presenta una solución de coexistencia de especies, matemáticamente 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ón de una (el depredador) o de ambas especies, extinción total. Esta solución es un punto-silla y por tanto inestable. A continuación se ha propuesto la construcción 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ón maestra para el conjunto de reacciones de Lotka-Volterra y se ha hallado la solución estacionaria. La situación del estado estacionario respecto al modelo determinista supone una serie discrepancia: la probabilidad es total para la solución de extinción total. Las causas son la consideración de variables discretas en el formalismo y la presencia de fluctuaciones importantes que acabarán tarde o temprano con la solución pseudo-estacionaria de estabilidad neutra, la coexistencia. Finalmente, se han deducido las ecuaciones macroscópicas de Lotka-Volterra depredador-presa a partir de los momentos diferenciales y de la expansión de van Kampen. El hecho de cuadrar las ecuaciones obtenidas con las deterministas exige asumir covarianza nula entre presa y depredador, una hipótesis injustificada que afecta posteriormente a la obtención 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ática de los sistemas biológicos. Universidad Nacional de Cuyo. Centro Atómico 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äuber, 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äuber 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ámicos Máster en Fı́sica de Sistemas Complejos, UNED [19] N.G. van Kampen Stochastic Processes in Physics and Chemistry Third edition, 2007 18