Academia.eduAcademia.edu

Libro ecología matemática.pdf

Ecología Matemática Principios y aplicaciones Fernando R. Momo Ángel F. Capurro Ecología Matemática Principios y aplicaciones Ediciones Cooperativas es un emprendimiento cooperativo de docentes de la Facultad de Ciencias Económicas de la Universidad de Buenos Aires para difundir sus trabajos e investigaciones Ni ng un a pa Momo, Fernando R. Ecología matemática : principios y aplicaciones / Fernando R. Momo y Angel F. Capurro - 1a ed. - Buenos Aires : Ediciones Cooperativas, 2006. 114 p. ; 21x14 cm. ISBN-10: 987-1246-21-8 ISBN-13: 978-987-1246-21-2 1. Ecología Matemática. I. Capurro, Angel F. II. Título  2006 Momo, Fernando R. Derechos exclusivos 1º edición, Julio 2006  2006 Ediciones Cooperativas Tucumán 3227 (1189) Buenos Aires – Argentina (54 011) 4864 5520 / (15) 4937 6915 http://www.edicionescoop.org.ar [email protected] Hecho el depósito que establece la ley 11.723 Impreso y encuadernado por: Imprenta Dorrego. Dorrego 1102, Cap. Fed. 1ª. ed. Tirada: 200 ejemplares. Se terminó de imprimir en Julio de 2006. Editorial asoci ada a: IMPRESO EN ARGENTINA – PRINTED IN ARGENTINE Dedicatoria Para mi amigo Ángel, in memoriam memoriam Conocí a Ángel Capurro cursando la Licenciatura en Ciencias Biológicas en la Facultad de Ciencias Exactas y Naturales de la Universidad de Buenos Aires. Muy poco después, ya recibidos, compartimos gran parte de nuestra formación de postgrado en Argentina y en Italia, donde asistimos en 1990 al Third Autumn Course on Mathematical Ecology en el ICTP en Trieste. Allí consolidamos nuestra amistad y colaboración profesional. Cuando en 1994 le propuse dictar juntos un curso de Ecología Matemática sólo tardó veinticuatro horas en enviarme su propuesta de programa que sería casi el programa definitivo y que fue sin duda la semilla de este libro que habíamos comenzado a escribir cuando un accidente automovilístico le quitó la vida en diciembre de 2000. Fue mi compañero de ruta, el Codirector de mi Programa de Investigación en Ecología Matemática que desarrollé en la Universidad Nacional de Luján, el director y orientador de muchísimos jóvenes investigadores, un científico brillante, entrañable amigo, esposo y padre. Por mi heterodoxa formación científica, muchas veces sentí que tenía pocos interlocutores y, en algunos temas, definitivamente Ángel era el único. Todas y cada una de las veces en que nos encontramos recuerdo haberme sentido bien, porque me contagiaba alegría, intensidad, ganas de vivir, algo de esa inagotable capacidad de trabajo que lo hacía diferente a todos y mejor que la mayoría de nosotros. Compartimos unas pocas pasiones, mayoritariamente académicas, y muchos pasos en muchas calles de varios lugares del mundo. Allí donde estaba, yo me sentía en casa, me sentía abriga- do, bien recibido; porque todos los lugares parecían su casa cuando él estaba allí y siempre tuvo el corazón abierto. Conocí pocas personas tan generosas, tan abiertas y tan chinchudas. ¡Pero con qué elegancia y buen humor sabía salir adelante, siempre al frente, siempre creando y abriendo el juego! Siempre lo admiré y siempre lo quise entrañablemente. Sé que lo sabía. Ante lo inabarcable de la muerte todos somos insignificantes, pero ante la vida, sólo algunos son como Ángel: un torbellino que intenta bebérsela, entenderla y mejorarla, todo al mismo tiempo. Este país, pobre de ideas y de ganas, lo extraña mucho; la ciencia lo extraña también. Pero sus amigos lo extrañamos mucho más, porque vislumbramos al magnífico ser humano que fue y que va a perdurar en todo lo que sembró. La concreción de este libro es un homenaje y un recuerdo de los que a él le hubiesen gustado: concretos y tangibles; un producto creativo que espero sea útil a muchos colegas. Amigo Ángel, ¡salud! Open Door, 5 de mayo de 2005 Introducción ntroducción ¿Por qué un libro de ecología matemática? Cualquier ecólogo que acostumbre leer las revistas científicas más difundidas de nuestra disciplina (digamos Ecology, American Naturalist, Oikos, Oecología, Ecosystems) se encontrará con multitud de ecuaciones y discusiones acerca de modelos matemáticos, con diagramas de fases, comentarios acerca de condiciones de estabilidad, coexistencia o exclusión. Incluso cualquier estudiante que explore uno de los textos clásicos o no tanto de la ecología (Margalef, Begon, Odum, Smith) se encontrará con que la mayoría de las teorías que sustentan la ciencia de la ecología desde el siglo diecinueve se expresan en forma de un par de ecuaciones diferenciales o se sintetizan con un gráfico. En este contexto no parece necesario justificar la existencia de este libro. Hay buenos textos de ecología matemática (incluímos un breve comentario sobre algunos en el final de este libro), incluso hay algunos en castellano, sobre todo de origen mexicano. Sin embargo no hay hasta la fecha un libro en español (aunque este está escrito en “argentino bonaerense”) que incluya una síntesis de las principales herramientas del ecólogo matemático a un nivel de interfase entre la aproximación típica del ecólogo (más intuitiva, más biológica) y la del matemático (más rigurosa, más abstracta). Después de dictar el curso de Ecología Matemática, Principios y Aplicaciones en universidades de toda la Argentina y Latinoamérica a un heterodoxo conjunto de profesionales: biólogos, agrónomos, matemáticos, físicos, químicos, ingenieros, que descubrieron la alegría de poder trabajar juntos en problemas ecológicos de interés común, aportando cada uno lo suyo, creemos que este libro se ubica precisamente en ese punto de interés común en el que se produce la interacción más fértil. Hay que agradecer a mucha gente. Por empezar al Instituto de Ciencias de la Universidad Nacional de General Sarmiento que ha apoyado y hecho posible esta publicación dándole su respaldo académico. También a todos los estudiantes que hicieron el curso 9 Fernando R. Momo - Ángel F. Capurro en sus diferentes épocas ya que siempre aportaron ideas nuevas y creativas. Hemos usado algunas de ellas como ejemplos en los capítulos y sus autores están allí debidamente mencionados en particular. También es necesario agradecer en particular a algunas personas por sus aportes esenciales: Juan Corley, Marcela Castelo, Leonardo Saravia, Emma Ferrero, Gustavo Ferreira, Irene Schloss, María Vernet, Néstor Claus, Pepe Ure, Matías Eöry, Carolina Eöry, Marisol Esusy, Julia Iribarren, María Eugenia Díaz, Carlos Coviella, Ariadna Hammann, Pedro Díaz. Por supuesto, ninguno de ellos es responsable por los errores que pueda contener el libro. También hay que agradecer a algunos maestros: Leonardo Malacalza, Jorge Rabinovich, Carlos Verona, Thomas Vevlen, Joel Cohen, Enrique Bucher, Fernando Kravetz, Eduardo Rapoport, Alfredo Salibián, Graciela Canziani, Alfredo Novelli, Alberto Jech, Simon Levin, Carlos Castillo Cháves, Joan Armengol, Ramón Margalef. Todos ellos fueron ejemplo, ayuda, inspiración. Por supuesto que tampoco tienen la culpa de nada en lo que a este libro respecta. El libro está en manos (y a merced) de sus lectores y potenciales críticos. Esperamos que sirva de ayuda para quienes quieran trabajar en una de las áreas más excitantes de la ecología actual: el punto en que los conocimientos se vuelven teoría cuantitativa y modelos generalizables, la zona que acerca nuestra ciencia a su forma más moderna y rigurosa, la parte más imaginativa y divertida de la ecología. Al menos para nosotros. Los autores 10 Capítulo 1 Los modelos matemáticos y su razón de ser Pareciera ser que los modelos matemáticos en ecología necesitan ser justificados de alguna manera, como si construirlos fuese una actividad censurable. Sin embargo, los modelos matemáticos se han ganado su lugar en la teoría ecológica como generadores de hipótesis, instrumentos de predicción y herramientas de descripción cuantitativa. Casi no hay proyecto reinvestigación en ecología que no pretenda culminar con algún tipo de “modelo matemático del proceso considerado” y hoy no es posible leer ninguna revista de investigación ecológica de primer nivel sin conocer aunque sea los rudimentos de las ecuaciones diferenciales y el álgebra de matrices. Desde los comienzos de la ecología como ciencia y aún antes, los modelos matemáticos de los procesos ecológicos fueron parte de los fundamentos de esta ciencia. Para muestra bastan algunos botones: los modelos de crecimiento poblacional de Malthus (1798) y Verhulst (1838, 1845 y 1847), los modelos de competencia y de depredación de Lotka (1925) y Volterra (1926), los trabajos sobre competencia, experimentales y matemáticos en una combinación magistral, de Gause (1934), los trabajos de Robert MacArthur entre 1950 y 1975, los de Hutchinson, especialmente a partir de 1957 desarrollando la moderna teoría del nicho ecológico. Si repasamos algunos de esos trabajos clásicos o si intentamos simplemente dar una forma más útil para la predicción a nuestras hipótesis experimentales o de campo, los modelos matemáticos aparecen como la herramienta por excelencia. Pero no sólo son una herramienta. Según Caswell, los modelos matemáticos son a las hipótesis teóricas generales lo que los diseños experimentales y los experimentos en sí son para las hipótesis particulares. Es decir, formular un concepto teórico en forma de modelo permite poner a prueba determinadas predicciones que derivan de esos con- 11 Fernando R. Momo - Ángel F. Capurro ceptos y explorar posibles consecuencias o comportamientos no previstos o no inmediatamente intuidos. El desarrollo de modelos permite también muchas veces explicar discrepancias entre los datos experimentales y las predicciones e identificar variables críticas para medir u observar en los experimentos o en los trabajos de campo. Comentemos un ejemplo simple1 que comienza con una observación casi trivial: cuando se representa la biomasa del perifiton (conjunto de algas que crecen sobre sustratos sumergidos) en función del tiempo los puntos se ajustan bien a un modelo de crecimiento logístico (curva sigmoidea) pero los valores altos de biomasa oscilan mucho alrededor de la curva teórica que predice el modelo. Esto podría ser algo trivial: simplemente esas fluctuaciones podrían deberse al error experimental o a variaciones espaciales alrededor de un valor medio que es, en definitiva lo que el modelo predice. Sin embargo hay algunos indicios de que estas fluctuaciones podrían ser indicativas de algo más interesante. Por ejemplo el hecho de que la magnitud de las mismas aumenta cuando el agua empieza a moverse (por ejemplo cuando se incrementa la velocidad de corriente en un río); sin embargo, cuando esa velocidad es muy alta, la biomasa deja de fluctuar. ¿Podría haber una relación causal entre la velocidad del agua y las fluctuaciones? ¿cuál podría ser esa relación exactamente? En el trabajo mencionado se planteó la siguiente idea: el perifiton crece siguiendo una dinámica logística pero a esa dinámica hay que agregarle un término que represente el efecto del movimiento del agua. Ese término podría ser positivo a velocidades bajas (porque hay células que caen y colonizan el sustrato) e incluso aumentar un poco con la velocidad (por ejemplo porque el movimiento del agua resuspende más células y aumenta la colonización). Sin embargo, cuando la velocidad del agua es mayor, el arrastre de células empieza a ser mayor que el aporte. 1 Momo, F.R. 1995. A new model for periphyiton growth in running waters. Hydrobiologia, 299: 215-218. Trabajo basado en observaciones de campo y sesudos comentarios de mi colega y compañero Adonis Giorgi. 12 Ecología Matemática - Principios y aplicaciones Hasta ahí no hay nada que incluya fluctuaciones sino que tendríamos sólo equilibrios diferentes. Pero, ¿qué pasa si la velocidad del agua no es constante? ¿y qué pasa cuando el arrastre se vuelve mayor que el crecimiento intrínseco del perifiton? Expresemos esto en forma de ecuaciones y gráficos: Velocidad de crecimiento de la biomasa = Crecimiento intrínseco logístico + Colonización o desprendimiento Que en ecuaciones queda de la siguiente forma: dB = r (B − U )(K − B ) + C dt Esta ecuación debe ser leída con cuidado. En primer lugar nos indica que la biomasa necesita alcanzar un cierto valor umbral para empezar a tener crecimiento positivo. Ese valor umbral es U. También nos dice que tiene un valor máximo para el cuál alcanza crecimiento cero (llamado en la jerga ecológica “capacidad de carga”); ese valor es K. Por último, vemos que al crecimiento intrínseco de la población (producido por su propia actividad biológica y el balance entre natalidad y mortalidad) se le agrega un término independiente C que representa una colonización si es mayor que cero o un desprendimiento de células por el arrastre del agua si es negativo. Supusimos que ese término varía en forma de parábola con la velocidad de corriente según el siguiente gráfico: 13 Fernando R. Momo - Ángel F. Capurro C versus velocidad de corriente C (mg de Clorofila a/día) 4 3 2 1 0 -1 0 2 4 6 8 10 12 -2 -3 -4 v (m/s) Este gráfico tiene valores más o menos arbitrarios pero ilustra el punto que se busca, es decir, a velocidad de corriente nula hay algún grado de colonización por las células que caen desde el fitoplancton; cuando aumenta la velocidad del agua a valores pequeños, la intensidad de colonización aumenta porque hay resuspensión de perifiton de las vecindades, pero si sigue aumentando el arrastre empieza a ser mayor y, por encima de una cierta velocidad crítica, el arrastre es mayor que el aporte y el término C pasa a ser negativo. La ecuación de este gráfico es una simple parábola: C = a1v 2 + a2 v + a3 con a1 < 0, a2 > 0 y a3 > 0. ¿Cuáles serían ahora las condiciones de equilibrio del sistema? Decimos que el sistema está en equilibrio cuando la biomasa permanece constante, es decir cuando no crece ni disminuye, es decir, cuando la velocidad de crecimiento es igual a cero. En otras palabras, cuando el primer término de la ecuación de velocidad (crecimiento intrínseco) se iguala con el segundo (colonización/ arrastre) en módulo pero con signo diferente. 14 Ecología Matemática - Principios y aplicaciones Permítanos el lector hacer un pequeño salto en la deducción ya que lo que nos interesa ahora no es desarrollar el álgebra del asunto sino explorar las consecuencias conceptuales de usar el modelo. Entonces digamos que igualando a cero la derivada y haciendo los malabarismos correspondientes de “pasaje de términos” llegamos a que, para que haya equilibrio, es necesario que se cumpla la siguiente igualdad: − rB + r (U + K ) = rUK − C B expresión que podemos analizar a partir del gráfico de las dos ecuaciones2 que se muestra a continuación, llamando al primer término F(B) y al segundo G(B): F(B) y G(B) en función de B 35 F(B), G(B) 30 25 F(B) 20 G(B) 15 10 5 0 0 5 10 15 20 25 Biomasa (mg Clorofila a) Es evidente aquí que los puntos en que se cruzan las dos funciones representan biomasas de equilibrio. El de la izquierda es un equilibrio inestable y el de la derecha uno estable. ¿Qué sucede 2 Esta elegante forma de analizar el modelo fue sugerida al autor del manuscrito por la Dra. Graciela Canziani y permitió hacer más comprensible para un público no matemático el razonamiento. 15 Fernando R. Momo - Ángel F. Capurro cuando la velocidad de corriente varía y con ella varía el término C? Pues bien, hacer C negativo es equivalente a “subir” la curva G(B) lo cual produce que los dos equilibrios, el estable y el inestable, tengan valores más cercanos. Si la velocidad del agua aumenta mucho y con ella C se hace negativo y de módulo muy grande, puede pasar que la curva G(B) quede completamente por encima de la curva F(B) y entonces no haya equilibrio estable posible. El perifiton no se desarrolla. Ahora bien, si C está variando (siguiendo las variaciones de velocidad de corriente), el equilibrio estable fluctúa todo el tiempo. ¡Como lo que se observa en los datos! Y si el valor de velocidad de corriente es muy alto, podemos pasar a la situación de inexistencia de equilibrio estable, lo cual biológicamente sucede cuando el agua produce el desprendimiento de un gran pedazo o “costra” de perifiton generando una brusca caída de la biomasa observada. El modelo presenta lo que técnicamente se llama “catástrofe de pliegue”3. El aumento del valor del parámetro r, la tasa intrínseca de crecimiento, también produce una mayor inestabilidad del sistema. Para nuestra agradable sorpresa, al colocar en el modelo los valores reales de velocidad de corriente, las predicciones del modelo representan muy bien las observaciones ¡fluctuación incluida! Moraleja: al construir el modelo matemático con supuestos biológicos razonables y explorar sus consecuencias y predicciones, logramos transformar lo que parecía “ruido” en nueva información. Y logramos predecir no sólo valores medios sino fluctuaciones. A lo largo de este libro veremos que esta es una situación que se repite. Los modelos matemáticos permiten descubrir dinámicas que a primera vista no eran obvias. Y puesto que la ciencia es una aventura de descubrimiento, ¿no es esta una razón importante para aprender algo más de modelos matemáticos? 3 En el sentido de René TOM y su teoría de catástrofes. 16 Capítulo 2 Modelos de crecimiento poblacional Sin duda, una de las aplicaciones más útiles de los modelos matemáticos en ecología la constituyen los modelos que representan o intentan representar las variaciones de densidad de las poblaciones. Se trata de expresar mediante una ecuación la densidad de la población (número de individuos o biomasa por unidad de espacio) en función del tiempo dada una población inicial o en función de una densidad anterior de la misma población. ¿Qué características debe tener un modelo matemático del crecimiento poblacional para que lo consideremos útil? En principio tiene que cumplir con algunas restricciones propias del fenómeno biológico que representan; por ejemplo, si la población es nula (vale cero) el crecimiento tiene que ser nulo también. También debe evitar las poblaciones negativas que pueden tener sentido en la ecuación pero no en el mundo real. Esta última restricción se soluciona poniendo condiciones de borde o eligiendo el tipo adecuado de funciones (por ejemplo, para una densidad poblacional que sea función decreciente de una variable se puede elegir una aproximación exponencial, que evita valores negativos, en vez de una lineal). Además necesitamos que el modelo represente razonablemente bien el comportamiento de los datos. Veremos que esto es posible aún con modelos muy sencillos y poco “realistas”. Un ejemplo de esto es el clásico modelo de crecimiento exponencial o modelo de Malthus, que es el más simple posible de los modelos continuos de crecimiento poblacional. En este modelo simplemente se afirma que la velocidad de crecimiento de la población es proporcional a su densidad. Esto, puesto en términos algebraicos, es lo siguiente: dN = r⋅N dt (2.1) 17 Fernando R. Momo - Ángel F. Capurro Como se ve, cumple con nuestras sencillas restricciones: la velocidad de crecimiento es cero cuando N es cero. Dada una población inicial mayor que cero, no hay posibilidad de que la población sea nunca negativa, aún cuando la constante r sea menor que cero. En este modelo, r representa la tasa de crecimiento per cápita de la población y tiene unidades de T-1. Se la llama generalmente tasa intrínseca de crecimiento y es fácil ver que: r= 1 dN ⋅ N dt (2.2) Podemos comprender que el modelo asume implícitamente que el efecto de la densidad poblacional sobre la velocidad de crecimiento es instantáneo. La tasa de crecimiento r del modelo exponencial es una constante. La ecuación diferencial se resuelve y se obtiene: N (t ) = N 0 ⋅ e rt (2.3) A pesar de su pasmosa sencillez, el modelo exponencial se aproxima sorprendentemente bien a la dinámica real de muchas poblaciones, sobre todo en períodos cortos. Crecen exponencialmente (o casi) la mayoría de las poblaciones consideradas “plaga” cuando los recursos de alimento y espacio son suficientes (imaginemos por ejemplo el crecimiento de una población de moscas en un criadero de pollos durante el verano, o el crecimiento de una población de bacterias en un recipiente con caldo nutritivo), pero el modelo exponencial se usa para aproximar el crecimiento de otras poblaciones donde el carácter de plaga no se aplica (o es mucho menos evidente) como por ejemplo la población humana. Observemos que un parámetro demográfico usado frecuentemente es el de “tiempo de duplicación de la población”, es decir el tiempo que tardaría tal o cual población en crecer hasta una densidad del doble que la actual. Ese tiempo es análogo a una vida media en el caso en que la población tenga un r menor que cero. Veamos los formalismos; si quiero calcular el tiempo de duplicación, cualquiera sea el valor inicial de la población N0 al cabo de un tiempo igual al tiempo de duplicación TD deberá ser el doble de N0; entonces: 18 Ecología Matemática - Principios y aplicaciones N TD = 2 ⋅ N 0 = N 0 ⋅ e rTD ∴ ln 2 = r ⋅ TD TD = (2.4) ln 2 r En el caso que mencionamos de la población humana, este valor se utiliza para comparar el crecimiento entre países y eventualmente en la planificación de gobierno. Cuando la población que se está estudiando es de una especie en peligro de extinción (caso en que r suele ser menor que cero), el TD es la vida media de la población, es decir el tiempo que tardaría en reducirse a la mitad la población actual. En este caso, lo que se planifica son las medidas de conservación o el modo de revertir la tendencia. El crecimiento dependiente de la densidad y el modelo logístico El siguiente paso en el camino de formular un modelo algo más realista, es considerar que la población no puede seguir creciendo hasta el infinito cada vez más rápido (que es en esencia lo que el modelo exponencial representa). Si esto no sucede, debe haber un “freno” al crecimiento poblacional y, si hemos de escribir un modelo matemático de ese freno, éste debe ser algún tipo de función de la propia densidad poblacional. En términos algo más precisos, diremos que la tasa de crecimiento r ya no será una constante sino que dependerá de la densidad poblacional. El modelo clásico para esta situación es el modelo logístico de crecimiento o modelo de Verhulst en el cual, la tasa de crecimiento r se considera una función lineal decreciente de la densidad poblacional: r = −m ⋅ N + c (2.5) donde m y c son constantes que tienen un sentido ecológico como veremos. Reemplazando esta fórmula en la ecuación 2.1 obtenemos 19 Fernando R. Momo - Ángel F. Capurro dN = (− m ⋅ N + c ) ⋅ N = c ⋅ N − m ⋅ N 2 dt (2.6) Como vemos, hemos pasado de una ecuación lineal para la velocidad de crecimiento a una ecuación cuadrática. Junto con el cambio de potencia aparece ahora otra novedad: la población no crece infinitamente. Es más, tiene puntos de equilibrio, es decir, valores de N para los cuales la velocidad de crecimiento es nula. Uno de esos puntos, por supuesto (esa fue nuestra condición para todo modelo “razonable” de crecimiento poblacional), es el cero: si la población es nula, la velocidad de crecimiento también lo es. Pero ahora hay otro punto de equilibrio. Veamos: 0 = c⋅ N −m⋅N2 ⇒ N ⋅ (c − m ⋅ N ) = 0 (2.7) ∴ N 1* = 0 N 2* = c m Es bueno aclarar aquí algo de la notación; el asterisco (*) denota que el valor de la variable corresponde a un valor de equilibrio (en el mundo de los modelos continuos, un punto para el cual la derivada de la variable respecto del tiempo es nula). Al valor de equilibrio distinto de cero c/m (en la jerga de los modelos, el equilibrio “no trivial”) le llamaremos K y asumiremos que representa el máximo tamamño que la población puede alcanzar en equilibrio. El nombre que recibe ese valor en ecología es el de “capacidad de carga”. A la constante c, que representa el valor máximo posible de la tasa intrínseca de crecimiento, le llamaremos por eso r0. Reemplazando ahora en la ecuación 2.5 obtenemos: r dN K −N  N = r0 N − 0 N 2 = r0 N ⋅ 1 −  = r0 N ⋅   dt K  K   K (2.8) Esta última es la formulación más conocida del modelo logístico que se interpreta biológicamente diciendo que la velocidad de crecimiento de la población es a la vez proporcional a su tamaño y a “cuanto le falta” para alcanzar su tamaño máximo K. 20 Ecología Matemática - Principios y aplicaciones Esta ecuación se puede integrar de diferentes maneras y una de las posibles expresiones de N en función del tiempo es: Nt = K K − N 0 − r0t 1+ e N0 (2.9) ¿Qué características tiene el punto de equilibrio no trivial N2*=K? El criterio de estabilidad Intuitivamente es sencillo entender que un punto de equilibrio es estable si la variable tiene tendencia a mantenerse en él. Haciendo esto un poco más formal (pero todavía poco) diríamos que si, partiendo del valor de equilibrio se aumenta o disminuye el valor de la variable (no olvidemos que se trata aquí de una densidad poblacional), el sistema cambiará en el sentido de compensar esa variación. Haciendo esto del todo formal, diremos que si construimos una nueva variable z = N – N*, el valor de equilibrio será estable si en su entorno el signo de la derivada es contrario al de la variable, o sea: z* = 0 estable ⇔ ∀ε dado z = z* + ε = ε → dz <0 dt (2.10) dado z = z * − ε = −ε → dz >0 dt Esto parece confuso y probablemente lo es, pero si lo expresamos en palabras es muy simple: cuando la población es mayor que el valor de equilibrio, tiende a bajar, cuando es menor, tiende a 21 Fernando R. Momo - Ángel F. Capurro subir; esto al menos en las cercanías del punto de equilibrio (esto es muy importante); es decir que el valor de equilibrio funciona como un atractor de la variable. La variable “artificial” z tiene una función derivada que es siempre de signo contrario al de z cuando su módulo es muy pequeño (z . 0). Es aún más fácil analizarlo gráficamente. Veamos el caso del modelo logístico. Aquí la función que relaciona la velocidad de crecimiento con el tamaño poblacional es una parábola con raíces en 0 y en K. Es fácil ver que si nos alejamos apenas del cero (poniendo un solo individuo) la velocidad de crecimiento es positiva y la población se aleja cada vez más del tamaño cero, hasta que se acerca al otro equilibrio, K, en ese caso, la derivada del tamaño respecto del tiempo es negativa para valores algo mayores que K y positiva para valores algo menores; entonces K es un atractor (punto de equilibrio estable). Velocidad de crecimiento vs Tamaño poblacional r0 = 0.8 K = 230 50 dN/dt 40 30 20 10 0 -10 0 -20 50 100 150 200 250 300 -30 N El criterio de examinar el gráfico y considerar estable cualquier punto de equilibrio en donde la función velocidad de crecimiento atraviese el eje del tamaño de arriba hacia abajo cuando nos movemos de izquierda a derecha, nos servirá para cualquier modelo de crecimiento. 22 Ecología Matemática - Principios y aplicaciones Modelos generalizados Volvamos a los modelos en sí. Si nos concentramos en la forma generalizada de todo modelo de crecimiento dN/dt = f(N) y en la condición de que f(0) = 0, vemos que los modelos estudiados hasta ahora (exponencial y logístico) son solamente dos posibles de una familia formada por una expansión en serie de la expresión generalizada, es decir dN = f ( N ) = a 0 + a1 N + a 2 N 2 + a 3 N 3 + ... + a k N n dt (2.11) Lo que hemos hecho en cada modelo es quedarnos con los términos de menor orden, el de primer orden para el modelo exponencial, los dos primeros para el logístico; pusimos como condición al empezar a formular los modelos que a0 = 0 siempre. Ni siquiera esto es necesario cuando consideramos poblaciones abiertas (con migración), pero es conveniente para desarrollar las herramientas matemáticas. ¿Tiene sentido biológico por ejemplo un modelo de crecimiento de orden 3? La respuesta es sí. No sólo eso, puede dar cuenta de un importante hecho biológico: el efecto Allee. Este efecto hace referencia al efecto de “rareza” poblacional, es decir, tiene en cuenta que la población puede necesitar una densidad mínima crítica para poder reproducirse (ya que así se aseguraría que existan encuentros suficientes de parejas). Esto quiere decir que cuando la densidad es menor que un valor umbral, la población tiene velocidad de crecimiento negativa. Gráficamente: 23 Fernando R. Momo - Ángel F. Capurro Velocidad de crecimiento vs Tamaño poblacional Efecto Allee con umbral en 75 y K en 200 15000 dN/dt 10000 5000 0 0 50 100 150 200 250 -5000 -10000 N Como vemos, si el tamaño poblacional es menor que 75 individuos (referidos al área unitaria que hayamos elegido), la población marcha hacia su extinción; cuando supera ese valor crítico, crece hasta estabilizarse en una capacidad de carga de 200 individuos. Nótese que ahora, el cero es un equilibrio estable (atractor) mientras que el valor 75 es un equilibrio inestable (repulsor), basta aumentar la población en un individuo sobre ese valor para que la misma crezca hasta K, del mismo modo, basta disminuir la población en un individuo a partir de 75 para que la población se extinga. Llamaremos a ese valor umbral crítico U. Entonces, la función que vemos en el gráfico queda representada como: dN = h ⋅ N ⋅ ( N − U ) ⋅ (K − N ) dt (2.12) que es una ecuación cúbica. Para el caso graficado h = 0.22, U = 75 y K = 200. Biológicamente, también puede interpretarse un modelo cúbico como este suponiendo que así como el término cuadrático representa la competencia intraespecífica (competencia con los congéneres) directa, el término cúbico representa la competencia 24 Ecología Matemática - Principios y aplicaciones difusa. No obstante, estas interpretaciones son dudosas y dependen del gusto del consumidor/usuario del modelo. Algunas variaciones sobre el modelo logístico y su dinámica Existe abundante literatura con modelos logísticos levemente modificados. No obstante, nuestro criterio general para analizar la presencia y estabilidad de puntos de equilibrio sigue funcionando. Veamos algunos ejemplos.   N θ  dN = r ⋅ N ⋅ 1 −    con θ diferente de 1 y mayor dt   K   1. que cero. Esta es una variante de la logística que permite explicar cambios en las tasas de competencia intraespecífica a medida que aumenta la densidad poblacional. Pertenece a Gilpin y Ayala (1973). El efecto que produce θ es el de tornar asimétrica la curva, desplazando el vértice hacia la derecha cuando es mayor que 1 (ver figura) y hacia la izquierda cuando es menor. Velocidad de crecimiento vs Tamaño poblacional r0 = 0.8 K = 230 π = 3.7 150 dN/dt 100 50 0 -50 0 50 100 150 -100 N 25 200 250 300 Fernando R. Momo - Ángel F. Capurro N dN = −α ⋅ N ⋅ ln K dt 2. de Swan (1977) parece un modelo algo absurdo para una población. Sin embargo fue desarrollado para representar el decrecimiento de un tipo particular de tumor (población celular) ante la quimioterapia. Nótese que mientras el tamaño de la población es mayor que el valor K, la velocidad de crecimiento es negativa (se asume α > 0); cuando N < K el logaritmo en la ecuación es menor que cero y la velocidad de crecimiento es positiva. Aquí estamos ante un estado de equilibrio K que es único y estable. El tumor no puede ser llevado a tamaño cero. El modelo no está definido para N = 0. Velocidad de crecimiento vs Tamaño poblacional K = 210 ฀ = 0.7 60 dN/dt 40 20 0 0 50 100 150 200 250 300 -20 -40 N 3. dN  a = r ⋅ N ⋅  − c − b ⋅ N  de Schoener. Tampoco tiedt  N ne definición en N = 0 ; sin embargo presenta otra vez un equilibrio estable distinto de cero. Sirve sólo para poblaciones abiertas (con inmigración). 26 Ecología Matemática - Principios y aplicaciones Velocidad de crecimiento vs Tamaño poblacional r = 0.8 a = 3300 c = 12 b = 0.4 3000 dN/dt 2000 1000 0 -1000 0 20 40 60 80 100 120 -2000 -3000 N Dejaremos para un capítulo posterior los modelos en los que se introduce un factor que representa la extracción de parte de la población. No obstante, los criterios de análisis siempre son los mismos. Los modelos discretos de crecimiento Cuando el tiempo no es considerado como una variable continua, se obtiene una familia de modelos diferente. En este caso la aplicación es inmediata a aquellas poblaciones con generaciones no superpuestas y episodios reproductivos anuales, aunque se pueden utilizar en otros casos. Este tipo de modelos no está definido por una ecuación diferencial, que expresa la velocidad de cambio de la variable en función de ella, sino por ecuaciones en diferencias, que expresan el valor que toma la variable en un momento como una función del valor que tomó en pasos de tiempo anteriores. Las ecuaciones son del tipo: N t = F ( N t −k ) (2.13) donde tanto t como k son números enteros. 27 Fernando R. Momo - Ángel F. Capurro El tamaño de los pasos de tiempo es arbitrario, es decir, el tiempo puede estar expresado en meses, años, semanas, etc., pero una vez elegido ese paso, sólo se anota el número entero que indica el número de paso en que estamos. Cuando el tamaño poblacional en el tiempo t depende sólo del tamaño poblacional en el tiempo t – 1 se dice que el modelo discreto es de primer orden, o de orden uno. El caso más simple de un modelo de primer orden es el modelo exponencial discreto: N t = φ ⋅ N t −1 (2.14) que se llama así porque puede reducirse a una fórmula cerrada análoga a la del modelo exponencial continuo. En efecto: N 2 = φ ⋅ N1 N 3 = φ ⋅ N 2 = φ ⋅ (φ ⋅ N1 ) = φ 2 ⋅ N1 N 4 = φ ⋅ N 3 = φ ⋅ φ ⋅ (φ ⋅ N1 ) = φ 3 ⋅ N1 (2.15) N k = φ k −1 ⋅ N1 siendo N1 el valor inicial de la población y φ una tasa discreta de crecimiento equivalente a er con r significando lo mismo que en (2.2) pero con un cociente incremental finito. No obstante, no todo es tan simple con los modelos discretos. Por ejemplo, la introducción de un término no lineal hace variar dramáticamente el comportamiento de la dinámica. Veamos el caso de la “logística discreta”, es decir, una posible imitación discreta del modelo logístico dada por: N t = a ⋅ N t −1 ⋅ (1 − N t −1 ) . (2.16) Este modelo está expresado en a forma llamada “canónica”, es decir, normalizado al intervalo 0-1 de manera que el único parámetro sea a. ¿Cómo definimos en un modelo así el valor de equilibrio? En este caso lo llamamos punto fijo y lo designaremos así cuando se cumpla la condición de que N* = F(N*), es decir que al aplicar la función que define el modelo al valor del punto fijo, el valor obtenido sea el mismo. Entonces el proceso iterativo de aplicar una y otra vez la ecuación no puede salirse de allí. En el caso del modelo logístico discreto enunciado, el valor del punto fijo es: 28 Ecología Matemática - Principios y aplicaciones N* = a −1 a (2.17) sólo definido con sentido biológico (poblaciones mayores o iguales que cero) para a > 1, lo que implica una población que crece. Este modelo no puede reducirse a una fórmula cerrada y sólo nos queda calcular los valores. Por suerte esto es muy fácil en la época de las planillas de cálculo y s interesante ver qué sucede a medida que aumentamos el valor de a. En todos los casos, los valores calculados corresponden a los puntos, la línea sólo es una ayuda visual ya que el modelo sólo permite calcular valores a intervalos regulares. No obstante podemos observar que el comportamiento del sistema pasa de parecerse mucho al de un modelo logístico continuo a presentar oscilaciones amortiguadas primero y luego oscilaciones no periódicas, aparentemente irregulares (caos). Todo con el mismo modelo cambiando sólo el valor de un parámetro. Esto nos abre un panorama completamente nuevo de análisis. La estabilidad o no del punto fijo depende ahora de los valores de los parámetros del modelo. Esto se estudiará más profundamente en el capítulo 6, pero es interesante preguntarse ¿cómo sé ahora si el punto fijo (“equilibrio”) es un atractor o no? Cuando se trabaja con modelos continuos, el criterio que usamos para el signo de la derivada puede generalizarse a sistemas con varias ecuaciones diferenciales y demostrar que un punto de equilibrio es estable cuando la matriz de Jacobi o “jacobiana” del sistema tiene un autovalor o valor propio dominante con parte real negativa (véase capítulo 3 para más detalles). Esto en el caso de una sola ecuación diferencial (una sola población en nuestro caso) es equivalente a que la derivada respecto del tamaño (N) de la función de velocidad de crecimiento sea negativa en el punto de equilibrio. 29 Fernando R. Momo - Ángel F. Capurro Nt vs. tiempo a = 1.2 0.18 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 0 5 10 15 20 25 30 Nt vs. tiempo a = 2.2 0.6 0.5 0.4 0.3 0.2 0.1 0 0 5 10 15 20 30 25 30 35 Ecología Matemática - Principios y aplicaciones Nt vs. tiempo a = 2.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 5 10 15 20 25 30 35 25 30 35 Nt vs. tiempo a = 3.7 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 5 10 15 20 Al pasar al dominio de lo discreto, existe un homeomorfismo entre ese criterio y otro. En el caso de modelos discretos, el módulo de la matriz jacobiana del sistema debe ser menor que uno para que el punto fijo sea un atractor (ver figura). En el caso en que tenemos una sola ecuación en diferencias, este criterio se transforma en que la derivada respecto de N de la función que representa Nt+1 vs. Nt 31 Fernando R. Momo - Ángel F. Capurro en el punto fijo tenga módulo menor que uno (esté comprendida entre –1 y 1). 32 Ecología Matemática - Principios y aplicaciones Apliquemos el criterio para nuestro simplísimo (pero no tanto) modelo discreto del crecimiento logístico. Dado que nuestro punto fijo está dado por el valor (a-1)/a (ecuación 2.17), tenemos que examinar qué valores toma la derivada de la función de Nt en ese punto. Veamos: F ( N ) = a ⋅ N ⋅ (1 − N ) dF ( N ) = a − 2aN dN dF ( N * )  a −1  = a − 2aN * = a − 2a  = a − 2a + 2 = 2 − a dN  a  (2.18) Entonces, para cumplir con nuestro criterio de estabilidad, a tiene que ser mayor que 1 y menor que 3; a partir de a > 3, el punto fijo ya no es un atractor y aparecen los comportamientos o dinámicas “raras” que incluyen ciclos límite (la población salta entre dos, cuatro, ocho o dieciséis tamaños en forma periódica) y caos. Gráficamente todavía se puede rizar el rizo (y esto es muy cómodo para los biólogos que tienen cierta pereza para calcular derivadas). Obsérvese que es muy sencillo representar F(N) en función de N y en el mismo gráfico representar la función identidad (G(N) = N). Entonces, el punto fijo N* es aquel donde F(N) se cruza con G(N), es decir, el punto en el cual, aplicar la función iterativa a un valor dado de N (en este caso N*) vuelve a dar el mismo valor. El gráfico, en el caso de nuestra versión del crecimiento logístico discreto luce como sigue: 33 Fernando R. Momo - Ángel F. Capurro 1.2 F(N); G(N) 1 0.8 F(N) 0.6 G(N) 0.4 0.2 N* 0 0 0.2 0.4 0.6 0.8 1 1.2 N La parábola representada en este gráfico nos permite “leer” el valor del siguiente N para cualquier valor inicial de Nt-1 (recordemos que normalizamos nuestro modelo para que los tamaños poblacionales queden comprendidos siempre en el intervalo [0; 1]). En rigor, estos gráficos representan en realidad mapas dinámicos del sistema y no tiene por qué cumplir estrictamente con las propiedades de lo que en matemática se denominan “funciones”. En este caso en particular sí las cumple porque es un modelo muy simple. Pero volvamos a nuestro criterio “chapucero”4 para determinar si el punto fijo es un atractor o no. Recordemos que la condición era que la derivada respecto de N de la función F(N) en dicho punto fuese de módulo menor que 1; y la derivada es la pendiente de la curva para ese punto, es decir, en el lugar donde se cruza con la función identidad G(N). Así que basta examinar la pendiente allí y si está comprendida entre 1 (caso en que el punto fijo es cero) y –1, caso en que comienzan los ciclos límite sabremos que el punto fijo de marras es un atractor (en el caso de modelos continuos hubiésemos dicho también “equilibrio estable”). Pero 4 Mi amigo y colega Jorge Velasco-Hernández me contó que en México le llaman “talacha” a la acción de resolver las cosas de modo desprolijo o poco ortodoxo. En Argentina se suele llamar a esto un modo chapucero de hacer las cosas o se dice simplemente que “lo atamos con alambre”. 34 Ecología Matemática - Principios y aplicaciones esto no es sólo válido para este caso en particular ¡sino para cualquier modelo que podamos imaginar por complicada que sea la función F(N)! Tenemos algo práctico. El lector inquieto puede ensayar estas propiedades iterando “a mano” la función en el mapa dinámico. En efecto, tome usted cualquier punto inicial en el eje de abscisas, suba verticalmente hasta la curva de F(N), allí podría leer en el eje de ordenadas el valor siguiente de la población y luego buscar ese valor en el eje de abscisas otra vez y repetir el proceso, pero podríamos ahorrar tiempo y trabajo moviéndonos horizontalmente hasta la curva identidad (función G(N)) y desde allí subir o bajar verticalmente otra vez hasta F(N). Si el punto fijo es un atractor, la repetición de este proceso nos llevará hasta el punto de cruce de las dos funciones. Este proceso es un poco engorroso pero es entretenido y didáctico. Recomendamos repetirlo por ejemplo con gráficos que este modelo (logístico discreto) realizados para distintos valores de a. El resultado es muy instructivo. En la siguiente figura se ve un recorrido posible empezando con N = 0.14. Se puede notar que se alcanza el punto fijo con unas pequeñas oscilaciones amortiguadas previas. Aquí a = 2.7 de manera que ese resultado coincide con lo que vimos en las simulaciones de N vs. T. 35 Fernando R. Momo - Ángel F. Capurro Para el lector más inquieto, proponemos analizar de esta simple forma el comportamiento de los dos modelos siguientes: 1.   N  N t = N t −1 1 + ρ ⋅ 1 − t −1  que es una versión K    no canónica. 2. N t = N t −1 ⋅ e  N  r ⋅ 1− t −1  K   que produce unas bonitas curvas logísticas para algunos valores de r. El ejercicio consiste aquí en encontrar el o los puntos fijos y deducir sus condiciones de estabilidad de acuerdo al valor de los parámetros. 36 Capítulo Capítulo 3 Depredación y competencia Los modelos de depredación y competencia forman parte de la batería de herramientas clásicas del ecólogo. Vito Volterra en Italia y Alfred Lotka en Estados Unidos fueron los precursores en este tema y crearon los modelos que, con diversas modificaciones y mejoras, seguimos usando hoy. El modelo de Volterra para depredación comienza suponiendo la existencia de dos poblaciones de animales, una de las cuales (el depredador) se alimenta de la otra (la presa). Se supone que las dos poblaciones están formadas por individuos idénticos, mezclados en forma más o menos homogénea en el espacio (veremos que esta suposición es importante a la hora de modelizar la interacción entre ambas poblaciones). Cada una de las poblaciones, en ausencia de la otra, tiene una dinámica “malthusiana”, es decir, exponencial: si no hay depredador, la población de la presa crece exponencialmente; si no hay presa, la población del depredador decrece exponencialmente. Podemos expresar esto en dos ecuaciones (nuevamente usaremos la convención de considerar a las constantes siempre positivas, para que el signo de los parámetros sea explícito): dN1 = r1N1 dt dN 2 = −r 2 N 2 dt (3.1) donde N1 representa la densidad de presas y N2 la densidad del depredador, mientras que r1 y r2 son los valores absolutos de las tasas intrínsecas de crecimiento. Ahora introduciremos la interacción, suponiendo encuentros casuales entre los depredadores y las presas, gobernados por la ley de acción de masas (como si fuesen moléculas que reaccionan entre sí en una mezcla homogénea); esto hace que la probabilidad 37 Fernando R. Momo - Ángel F. Capurro de un encuentro exitoso5 entre depredador y presa (p) multiplicada por las densidades de ambas poblaciones equivalga a la velocidad de depredación. dN1 = r1N1 − pN1N 2 dt dN 2 = −r2 N 2 + epN1N 2 dt (3.2) Estas dos ecuaciones forman ahora un sistema que representa la dinámica conjunta de las dos poblaciones. Lo que se resta a la presa, es decir la velocidad de depredación, se le suma al depredador pero multiplicado por un factor e que representa la eficiencia de transformación de presa comida en depredador. Cálculo del equilibrio y estudio de la dinámica Como siempre, una de las primeras cuestiones de interés es la búsqueda de puntos de equilibrio y el estudio de su dinámica. Para esto debemos igualar a cero ambas ecuaciones ya que nos interesa el caso en que las dos poblaciones están en equilibrio juntas. dN1 = r1N1 − pN1N 2 = 0 dt r ⇒ N 2* = 1 p dN 2 = −r2 N 2 + epN1N 2 = 0 dt r ⇒ N1* = 2 ep 5 (3.3) Entendemos por encuentro “exitoso” aquel en que el depredador efectivamente atrapa y se come a la presa, así que estamos mirando el éxito desde el punto de vista del depredador. 38 Ecología Matemática - Principios y aplicaciones donde las densidades con asterisco representan valores de equilibrio6. Llegando a este punto tenemos varias formas de continuar analizando el sistema. Una de ellas es en análisis gráfico. Si representamos la densidad del depredador en función de la densidad de la presa obtenemos lo que hemos llamado un “mapa” del sistema o diagrama de estado donde cada punto del plano representa una posible combinación de densidades de presa y depredador7. Todos los valores de depredador para los cuáles la densidad de presas es igual a N1*, constituyen una línea a la que llamaremos isolínea8 del depredador, que incluye un conjunto de estados para los cuales el crecimiento de la población del depredador es nulo (la derivada de N2 respecto del tiempo es cero). En nuestro mapa de estados esto se representa como una línea vertical gris. La población del depredador crece en cualquier estado que esté a la derecha de esa línea (porque hay más presas que las correspondientes al equilibrio) y decrece en cualquier estado representado por un punto a la izquierda de esa línea. Esto se representa por flechas verticales grises. Análogamente, los valores representados por una línea horizontal negra, constituyen la isolínea de la presa, es decir, el conjunto de estados para los cuales el crecimiento de la población de presa es nulo (la derivada de N1 respecto del tiempo vale cero). La población de la presa crece cuando hay menos depredadores (debajo de la isolínea) y decrece cuando hay más (encima de la isolínea); representamos este comportamiento con las flechas negras. 6 Por supuesto, no tiene interés biológico el caso en que ambas poblaciones tienen valor cero, aunque cumple con la condición de equilibrio. 7 Obviamente, por razones biológicas, sólo tienen sentido para nosotros las poblaciones con densidades mayores o iguales que cero; por lo tanto nos interesa sólo el primer cuadrante en este espacio definido por los dos ejes cartesianos. 8 También llamada “nulclina”. 39 Fernando R. Momo - Ángel F. Capurro N2 r1/p r2/ep N1 Si pensamos que cuanto más cerca está un punto de una de las isolíneas, menor es la velocidad con que crece la población a la que esa isolínea corresponde (por ejemplo, si estamos cerca de la isolínea de la presa, la “velocidad” horizontal del sistema es menor) y que, en cada punto posible del espacio de estados hay una velocidad dada por la resultante de la velocidad horizontal (la velocidad de crecimiento de la presa) y la velocidad vertical (la velocidad de crecimiento del depredador), podemos trazar a mano alzada las posibles órbitas del sistema a partir de cualquier punto. Sin embargo, nada nos asegura que estas órbitas sean realmente así; debemos rizar el rizo matemático para describirlas. La obtención de la fórmula de las órbitas Para estudiar cómo son realmente las órbitas, podemos solucionar el sistema de ecuaciones diferenciales. Para esto, el primer paso es reemplazar los valores de r1 y r2 obtenidos a partir de los valores de equilibrio de las poblaciones: 40 Ecología Matemática - Principios y aplicaciones r1 = pN 2* (3.4) r2 = epN1* Reemplazando estos valores en el sistema de ecuaciones diferenciales obtenemos ( ) dN1 = pN 2* N1 − pN1N 2 = − pN1 N 2 − N 2* dt dN 2 = −epN1* N 2 + epN1N 2 = epN 2 N1 − N1* dt ( ) (3.5) que es una bonita expresión porque las velocidades de cambio de las variables quedan en función de la diferencia entre su valor actual y su valor de equilibrio. Dividiendo miembro a miembro la segunda ecuación por la primera obtenemos dN 2 N 2 (N1 − N1* ) = −e dN1 N1(N 2 − N 2* ) ⇒ N 2 − N 2* ( )dN N 2 ( = −e N1 − N1* 2 )dN N 1 1 (3.6) Aplicando la propiedad distributiva, simplificando y planteando las integrales, tenemos ∫ dN 2 −N 2 ∫ * ⇒ ( dN 2 dN = −e ∫ dN1 + eN1* ∫ 1 N2 N1 ) 1 N 2 − N 2* ln N 2 + N1 − N1* ln N1 = C e (3.7) 41 Fernando R. Momo - Ángel F. Capurro donde C es una constante de integración. Esta ecuación representa una serie de órbitas elípticas, una para cada valor de C suponiendo que los demás parámetros no cambian. El sistema entonces oscila siguiendo alguna de esas órbitas, dependiendo de las condiciones iniciales. El período de oscilación es proporcional a 1 r1r2 . Estudio del equilibrio usando la matriz jacobiana Otra forma de estudiar el comportamiento del sistema, una vez obtenido un punto de equilibrio, es mediante la linealización del sistema en un entorno de ese punto. Esto quiere decir que construimos una nueva variable que es simplemente la diferencia entre el valor instantáneo de la variable primitiva y su valor de equilibrio (Ni – Ni*); esta nueva variable es cero en el punto de equilibrio y distinta de cero en todo otro punto. Para saber si el punto de equilibrio funciona como un atractor, es decir, si el sistema tiende a volver al equilibrio después de una perturbación, tenemos que estudiar el valor de la derivada respecto del tiempo de esta nueva variable en un entorno del cero. Esto se hace en un entorno pequeño para poder asumir un comportamiento cuasi lineal. La técnica que permite realizar esta tarea consiste en calcular una matriz, llamada Jacobiana. Llamaremos x a la primera variable del sistema (en nuestro caso N1) e y a la segunda; y llamaremos F(x,y) a la primera función del sistema (en nuestro caso dN1/dt) y G(x,y) a la segunda. Con esta nomenclatura, la matriz jacobiana se construye como:  dF  dx J=   dG  dx dF  dy    dG  dy  (3.8) ¿Cómo es entonces la matriz jacobiana de nuestro sistema de ecuaciones de presa-depredador? 42 Ecología Matemática - Principios y aplicaciones r1 − pN 2  J=   epN 2     − r2 + epN1  − pN1 (3.9) Como lo que estamos interesados en estudiar es el comportamiento del sistema cerca del punto de equilibrio, reemplazamos los valores de N1 y N2 en la matriz por sus valores en equilibrio (que obtuvimos antes). Está claro que el punto de cruce de las dos isolíneas representa un estado del sistema para el cual las derivadas de las dos variables respecto del tiempo son cero y el sistema en su conjunto está en equilibrio. Reemplazando y simplificando vemos que la matriz jacobiana en equilibrio es  0 * J =  er  1 − r2  e   0  (3.10) El comportamiento del sistema y la estabilidad del punto de equilibrio se pueden conocer ahora simplemente obteniendo los autovalores de esa matriz. Si uno de los autovalores es dominante (tiene módulo más grande que el resto) será el que determine el comportamiento a largo plazo del sistema. Cuando ese autovalor tiene parte real negativa, el punto de equilibrio es un atractor, si su parte real es positiva, el punto de equilibrio es un repulsor, si es cero, el equilibrio es neutro. Si la parte imaginaria de alguno de los autovalores es distinta de cero, el sistema presenta oscilaciones, que pueden ser amortiguadas, constantes o amplificadas, dependiendo del signo de la parte real del autovalor dominante. Una matriz cuadrada tiene tantos autovalores como filas o columnas tiene; esto nos proporciona una información adicional: de la misma manera en que el autovalor dominante determina el comportamiento asintótico, es decir, a largo plazo del sistema, los 43 Fernando R. Momo - Ángel F. Capurro autovalores subordinados pueden tener un gran peso en el comportamiento transitorio9 del sistema. Pero veamos qué nos dice la matriz jacobiana del modelo de Volterra. Para eso debemos obtener sus autovalores en el punto de equilibrio, igualando a cero el determinante de (J* - λI) donde λ representa los autovalores e I es la matriz identidad (con valor 1 en las diagonales y cero en los demás elementos). ( ) −λ det J * − λI = er1 − r2 e (3.11) −λ El determinante de una matriz de 2 x 2 es sencillo y obtenemos para λ lo siguiente λ2 + r1 r2 = 0 ⇒ λ = ± − r1 r2 ⇒ λ1 = i r1r2 (3.12) λ2 = −i r1 r2 Es decir que los dos autovalores tienen el mismo módulo (ninguno es dominante) y ambos son imaginarios puros. Esto significa que el sistema presenta oscilaciones de amplitud constante que depende de las condiciones iniciales. Sabemos que las oscilaciones son constantes (no se amortiguan ni se amplifican) porque las partes reales de los autovalores valen cero. 9 Usaremos la denominación “transitorio”, que, en castellano, parece más correcta que la habitual “trasiente”. 44 Ecología Matemática - Principios y aplicaciones Modificaciones al modelo básico de Volterra Una primera modificación al modelo básico, que fue propuesta por el mismo Volterra, es suponer que la presa, en ausencia del depredador, tiene crecimiento logístico, es decir, que su población se autorregula. Matemáticamente, esto consiste simplemente en agregar un término cuadrático a la ecuación de la presa y nuestro sistema queda expresado como  N  dN1 r N2 = r1N1 − 1 1 − pN1N 2 = r1N11 − 1  − pN1N 2 dt K1  K1  dN 2 = −r2 N 2 + epN1N 2 dt (3.13) Cuando calculamos las isolíneas, para buscar el punto de equilibrio si es que existe, vemos que la correspondiente al depredador no cambia respecto al modelo básico y sigue siendo N1* = r2 ep , pero la de la presa tiene ahora una pendiente negati- va en el mapa de estados y está expresada por la función N 2 = r1 p − N1 r1 pK 1 ; esto nos indica, como primer punto, que para que exista un equilibrio con las dos poblaciones no nulas, tiene que cumplirse la condición de que r2/ep < K1, porque si la presa tiene una capacidad de carga demasiado baja o el depredador tiene muy poca eficiencia, no tendremos presas suficientes para mantener una población viable de depredadores. La otra sorpresa la encontramos al estudiar la dinámica. Usted puede verificar fácilmente (¿?) que los autovalores de la matriz jacobiana en el equilibrio están dados por 45 Fernando R. Momo - Ángel F. Capurro λ1,2 = rr − 12 ± K 1ep  r1r2  r r2   − 4er1r2 + 4 1 2 K 1p  K 1ep  2 2 (3.14) Esta expresión parece algo complicada como para sacar alguna conclusión respecto a la dinámica del sistema, sin embargo, teniendo en cuenta la condición para que exista un punto de equilibrio con las dos poblaciones positivas (r2/ep < K1) se puede ver que η= r2 <1 epK 1  r  µ =  e − 2  pK 1   (3.15) de manera que podemos simplificar la expresión para los autovalores como λ1, 2 = − η ⋅ r1 ± (η ⋅ r1 )2 − 4 ⋅ r1 ⋅ r2 ⋅ µ 2 (3.16) Si examinamos esta expresión vemos que la parte real de los autovalores siempre es negativa y que la parte imaginaria puede ser distinta de cero dependiendo del valor de los parámetros, sobre todo de la relación entre r1 y r2. La conclusión es que el punto de equilibrio funciona como un atractor y que las poblaciones alcanzarán su valor de equilibrio bien en forma asintótica sin oscilaciones o en forma de oscilaciones amortiguadas, dependiendo del signo de la expresión bajo la raíz cuadrada. Otra modificación interesante del modelo base es aquella en que suponemos que el depredador se alimenta de un estadío avanzado de la presa que ya no es reproductivo y por lo tanto no influye en el crecimiento de la población de presas. Sería el caso en que un organismo (que pude ser el hombre) consume los estadíos post-reproductivos de un pez. También puede asimilarse al caso de 46 Ecología Matemática - Principios y aplicaciones los organismos comedores de hojarasca en el suelo o materia orgánica caída en un río. En todos los casos, si modelizamos la población de presas que entra en el sistema (el estadío postreproductivo o el stock de hojarasca) tenemos que considerar que, en ausencia del depredador, crece con una velocidad constante w y decrece con una tasa de descomposición, decaimiento o arrastre, según el caso, que llamaremos d. El sistema de ecuaciones queda entonces dN1 = w − dN1 − pN1N 2 dt (3.17) dN 2 = −r2 N 2 + epN1N 2 dt La isolínea de la “presa” queda entonces expresada por una hipérbola rectangular desplazada hacia abajo: N 2* = w d − pN1 p (3.18) En este caso, si existe el punto de intersección entre las dos isolíneas (la del depredador es igual que en los casos anteriores), éste representa un equilibrio estable, es decir, un atractor. Queda a cargo del lector analizar la matriz jacobiana. Respuestas funcionales y respuestas numéricas Se llama respuesta funcional a la función que representa como varía la cantidad de presas devoradas por un depredador en función de la oferta de presas. Se llama respuesta numérica del depredador a la función que representa cuantos depredadores nuevos se agregan a la población del mismo por presa devorada. En el modelo clásico de Volterra, la respuesta funcional y la respuesta numérica son lineales y esta última es simplemente el producto de la primera por una eficiencia. Esto se comprende mejor si escribimos las ecuaciones de un sistema depredador-presa de forma genérica como 47 Fernando R. Momo - Ángel F. Capurro dN1 = N1G1(N1 ) − N 2 H (N1 ) dt (3.19) dN 2 = N 2G2 (N 2 ) + N 2 J (N1 ) dt Aquí Gi(Ni) representa la tasa intrínseca de crecimiento de la población i, H(N1) representa la respuesta funcional, que en el modelo de Volterra es igual a pN1, y J(N1) representa la respuesta numérica que en el mismo modelo es epN1. Holling desarrolló un estudio clásico acerca de las respuesta funcionales considerando el uso del tiempo disponible por parte del depredador. Algunos animales responden en forma lineal a la disponibilidad de presas hasta llegar a un punto límite. Tal es el caso de organismos filtradores que comen cada vez más presas cuantas más hay hasta que se satura su sistema filtrador y un aumento ulterior en la oferta de presas no produce ningún efecto en la captura; esto es lo que se conoce como una respuesta funcional de tipo I. La cantidad de presas filtradas en un tiempo dado, antes de llegar al máximo, es HT (N1 ) = qN1T donde q es un parámetro que representa el volumen filtrado por unidad de tiempo. Cuando el tiempo total requerido para manipular las presas y comerlas es igual a T la captura per cápita se hace constante. Si suponemos que el depredador dispone de un tiempo total dedicado a nutrirse y que ese tiempo total lo reparte entre buscar presas y manipularlas para comérselas, tenemos que T = TS + TH donde TS representa el tiempo de búsqueda y TH el tiempo de manipulación. Si llamamos q al volumen de espacio explorado por unidad de tiempo y pA a la probabilidad de ataque exitoso a una presa tenemos que el número de presas encontradas durante el tiempo TS será qN1TS y el número de presas devoradas será D = p A qN1TS . Así las cosas, el tiempo que el depredador dedica a manipular esas presas será el número de presas devoradas multiplicado por el tiempo de manipulación que lleva cada una de ellas: TH = D ⋅ t H = p A qN 1TS t H (3.20) 48 Ecología Matemática - Principios y aplicaciones Dado que el tiempo total es la suma de los tiempos de manipulación y búsqueda como ya se dijo, si queremos expresar el número de presas devoradas por unidad de tiempo obtendremos E= p A qN1TS p A qN1 D = = T p A qN1TS t H + TS 1+ p A qN1t H (3.21) que es la respuesta funcional tipo II de Holling, en la cual el depredador se va “saturando” de presa paulatinamente. En el caso en que, para bajas densidades de presa, un aumento en el número de capturas haga que el depredador busque con más intensidad, produciéndose una aceleración inicial en la respuesta funcional, tenemos que pAq = E= aN1 1 + bN1 (3.22) aN12 1 + bN1 + aN12 t H que es la respuesta funcional de tipo III. Es interesante hacer notar que las respuestas funcionales no lineales pueden producir algunos afectos en la dinámica que se asemejan a los producidos por el efecto Allee (véase el capítulo 2). Si consideramos una presa que crece en forma logística y que es sometida a una depredación con una respuesta funcional tipo II, vemos que la velocidad neta de crecimiento de la presa puede calcularse como la diferencia entre su velocidad de crecimiento propia (la logística) y la velocidad de consumo por parte del depredador, representada por una familia de curvas dependiendo del número de depredadores. Esto permite la existencia de múltiples puntos de equilibrio que pueden ser estables o inestables según el caso. En este modelo, uno puede imaginar que la presa es, por ejemplo, una hierba y el depredador es un herbívoro; al ir aumentando el número de herbívoros (es decir, la presión de pastoreo), hay un desplaza- 49 Fernando R. Momo - Ángel F. Capurro miento de los puntos de equilibrio y se produce una catástrofe de pliegue10. Los modelos de competencia La competencia es la otra interacción binaria estudiada desde los orígenes de la ecología matemática; en este caso, hablaremos por prioridad del modelo de Lotka, aunque Volterra también desarrolló el tema. Si bien la competencia es, en rigor, una relación ternaria (están los dos competidores y el recurso) el tratamiento matemático se puede realizar obviando la población recurso y considerando sólo la retroalimentación positiva entre ambos competidores. Un primer planteo consiste en suponer que la presencia de la segunda población influye en la primera como si agregáramos más individuos de la misma, con una cierta equivalencia. Este es el caso en la competencia por interferencia, donde la presencia de individuos de la otra población hace que los de la primer tengan menor posibilidad de acceder al recurso. Expresamos esto sumando a los individuos de cada población un número que se obtiene multiplicando la abundancia de la otra población por un coeficiente y suponemos que ambas poblaciones tienen, en ausencia de la otra, un crecimiento logístico:  (N + αN 2 )  dN1  = r1N1 1 − 1 dt K 1    (N + βN1 )  dN 2  = r2 N 2 1 − 2 K2 dt   (3.23) A los factores de equivalencia, α y β, los llamaremos coeficientes de competencia. Para estudiar la dinámica de este modelo seguimos el mismo camino que utilizamos con el modelo de depredador y pre10 Un aumento pequeño en el número de herbívoros puede producir un descenso repentino en la biomasa de hierba; esa biomasa no se recupera inmediatamente al disminuir la presión de pastoreo. 50 Ecología Matemática - Principios y aplicaciones sa, es decir, calculamos lo que pasa cuando las dos derivadas son cero. Veremos que este modelo admite varios casos con resultados finales distintos y por esa razón es más fácil estudiarlo gráficamente (con el diagrama de fases) que utilizando el Jacobiano y demás herramientas analíticas relacionadas. Sin embargo veremos que otro método de análisis de estabilidad (las funciones de Lyapunov) funciona bien en este caso. Cuando igualamos a cero las dos derivadas, es obvio que existe un punto de equilibrio “trivial” cuando ambas poblaciones valen cero. En tal caso el problema biológico ya no existe, así que nos limitaremos al caso en que se eligen las soluciones diferentes de cero, es decir, cuando la porción entre paréntesis de las ecuaciones es igual a cero y, por lo tanto, las derivadas también. En ese caso la isolínea de la población 1, corresponde en un sistema de ejes N2 vs. N1 a: 1− N1 + αN 2 =0 K1 ∴ K 1 − N1 − αN 2 = 0 (3.24) ∴ N2 = K1 α − N1 α Esta última ecuación es una recta que corta al eje vertical (N2) en K1/α y al eje horizontal (N1) en K1. Este último valor representa el tamaño máximo de la población 1 que sólo se alcanza en equilibrio si la población competidora no tiene individuos. Haciendo un razonamiento similar, se obtiene la isolínea de la población 2, que también es una recta: 51 Fernando R. Momo - Ángel F. Capurro 1− N 2 + βN1 =0 K2 ∴ K 2 − N 2 − βN1 = 0 (3.25) ∴ N 2 = K 2 − βN 1 Es decir que la población 2 llega a su capacidad de carga en equilibrio cuando el número de individuos de la población 1 es cero. Para valores crecientes de N1, el número de individuos en equilibrio de la población 2 (N2) desciende en proporción β a N1. ¿Existe una situación en la cual ambas poblaciones estén en equilibrio y por lo tanto el sistema lo esté? Esto se comprende mejor gráficamente. Aquí representamos las cuatro posibles situaciones según los juegos de valores de los parámetros K1, K2, α y β. En todos los gráficos se utilizó K1 = K2 = 200 y los valores de α y β fueron los siguientes: α β Caso 1 0.8 1.2 Caso 2 1.3 0.75 Caso 3 1.2 1.3 Caso 4 0.8 0.75 Los casos 1 y 2 son relativamente triviales y corresponden a las situaciones en que una de las dos especies es mejor competidora que la otra y termina desplazándola (exclusión competitiva). El equilibrio del sistema es ese caso (descartado el mencionado doble cero) corresponde a una de las poblaciones en cero y la otra en su capacidad de carga (K). Cual de las dos poblaciones “gana” la competencia depende de los valores de los parámetros, especialmente α y β que determinan el “peso” que cada individuo de una especie representa para la otra en unidades de individuos propios. 52 Ecología Matemática - Principios y aplicaciones El caso 3 empieza a ponerse interesante porque existe un punto de equilibrio donde las dos poblaciones tienen un tamaño no nulo; este punto corresponde al punto de cruce de las isolíneas. Por definición, en ese punto ambas poblaciones tienen velocidad de crecimiento cero; ¿es este un equilibrio estable? Si intentamos contestar esa pregunta utilizando la linealización y calculando la matriz Jacobiana nos encontraremos con una serie de cálculos bastante engorrosos y poco claros. Sin embargo, es más fácil analizar directamente los gráficos teniendo en cuenta lo que las isolíneas representan (derivadas iguales a cero). Está claro que cualquier estado del sistema de competidores está representado por un punto en el gráfico N2 vs. N1, porque ese punto representa una tamaño para cada población. Ahora bien, si el punto se siúa por encima de la isolínea de N1, el sistema se moverá hacia la izquierda en el plano (en el sentido de disminuir el tamaño de la población 1; del mismo modo, si el punto se encuentra por encima de la isolínea de N2, el sistema se moverá hacia abajo, en el sentido de disminuir el tamaño de la población 2. El movimiento resultante (es decir el cambio en las dos poblaciones) es como la suma vectorial de los dos cambios. Lo simétrico sucede si el punto se encuentra por debajo de alguna de las isolíneas. En la siguiente figura se representa lo que pasa en un diagrama de fases como el del caso 3. Caso 1: N2 vs N1 300 250 isolínea de N2 200 isolínea de N1 150 100 50 0 0 50 100 53 150 200 250 Fernando R. Momo - Ángel F. Capurro Caso 2: N2 vs N1 250 200 isolínea de N2 150 isolínea de N1 100 50 0 0 100 200 300 Caso 3: N2 vs N1 250 200 isolínea de N2 150 isolínea de N1 100 50 0 0 50 100 150 200 250 Caso 4: N2 vs N1 300 250 isolínea de N2 200 isolínea de N1 150 100 50 0 0 100 200 54 300 Ecología Matemática - Principios y aplicaciones En negro se representa la isolínea de N1 y las flechas que indican si N1 sube o baja, en gris el equivalente para N2 . Es fácil concluir que el punto de equilibrio es estable ya que el sistema tiende a él desde cualquier punto del primer cuadrante (poblaciones iniciales distintas de cero y obviamente positivas). Dejamos para el lector inquieto realizar el mismo ejercicio con el caso 4. Podrá comprobar que el punto de cruce entre las dos isolíneas funciona como un repulsor y, dependiendo de los valores iniciales de las poblaciones, gana una de las especies o gana la otra, alcanzando su población la capacidad de carga que le corresponde en cada caso. N2 K1/a K2 K1 K2/b N1 Las funciones de Lyapunov Las funciones de Lyapunov otra forma de estudiar estabilidad de un sistema sin conocer la solución de las ecuaciones diferenciales. La función de Lyapunov puede ser entendida como una energía del sistema. Si hay un punto de equilibrio x* para la variable (punto que hace cero la derivada de x respecto del tiempo), la función de Lyapunov indicará estabilidad si es decreciente para 55 Fernando R. Momo - Ángel F. Capurro todo punto excepto x*. La función de Lyapunov V(x) debe ser tal que cumpla lo siguiente: V(x) > 0 y dV(x)/dt < 0 para todo x distinto de x* V(x*) = dV(x*)/dt = 0 Cuando un sistema tiene estabilidad de Lyapunov, cumple con la definición más fuerte de estabilidad (estabilidad uniforme asintótica “en grande”) siempre que se cumpla que V(x) > 0 para todo x distinto de x* y dV(x)/dt < 0 para todo x distinto de x*. Ejemplo: Evaluemos un sistema depredador presa dado por las siguientes ecuaciones: dN1/dt = N1 (10 - 2.2 N1 - 2 N2) - 1.25 dN2/dt = N2 (-0.3 + 0.2 N1 - 0.2 N2) sabiendo que la función de Lyapunov correspondiente es: V(N1, N2) = N1 - N1* - N1* ln (N1*/N1) + 10 [N2 - N2* - N2* ln (N2*/N2)] Haciendo los cálculos se ve que el punto de equilibrio de este sistema sí es estable lo cual es razonable dado que tanto la presa como el depredador tienen crecimiento logístico. Se puede demostrar que, en general, para todos los modelos de dinámica dados por una ecuación diferencial del tipo: dx/dt = x f(x) tipo del cual el modelo logístico es un caso particular, la función de Lyapunov correspondiente es: V(x) = x - x* - x* ln(x*/x) entonces es fácilmente demostrable (¡puede usted hacerlo!) que para que la derivada de la función de Lyapunov sea menor que cero 56 Ecología Matemática - Principios y aplicaciones para todo punto distinto de x* debe pasar que la función f(x) sea menor que cero a la derecha de x* y mayor que cero a la izquierda de x*. ¿Servirá la misma función de Lyapunov para un sistema en competencia? Desgraciadamente no hay una “receta” para encontrar funciones de Lyapunov en cualquier sistema. No obstante, como veremos en otra parte de este libro hay funciones que tienen sentido biológico y que se comportan como si fuesen funciones de Lyapunov matemáticamente hablando. Algunas de esas funciones son derivadas de la teoría termodinámica de sistemas abiertos. 57 Capítulo 4 Algunas aplicaciones de los modelos continuos Los modelos matemáticos continuos tienen múltiples aplicaciones a problemas ecológicos y no tan ecológicos en áreas tan diversas como la epidemiología, la economía, la conservación o la producción agropecuaria. En este capítulo intentaremos mostrar algunas de esas aplicaciones con ejemplos sencillos. El primer ejemplo es la aplicación a la extracción de un recurso biológico como una pesquería o la explotación maderera de un bosque. Inicialmente supondremos que el recurso que se explota es renovable, que tiene una biomasa (que es el objeto de nuestra extracción) y que esa biomasa tiene una dinámica. Una aproximación sencilla para esa dinámica es el modelo logístico. Es decir que la biomasa del recurso crece hasta alcanzar una cierta “capacidad de carga” caracterizada por un parámetro K. Así, la dinámica del recurso en ausencia de explotación puede ser descripta por la ecuación del modelo logístico usando como variable no el número de individuos sino la biomasa; por supuesto asumimos que la variables es una densidad es decir que la biomasa está expresada por unidad de espacio, área en el caso de un bosque, volumen de agua en el caso de un pez. La ecuación es, una vez más: dB  B = r ⋅ B ⋅ 1 −  dt  K (4.1) Ahora tenemos que agregarle a nuestro modelo la explotación del recurso, ¿cómo lo hacemos? Una forma simple, que da lugar a dinámicas muy interesantes, es suponer una velocidad de explotación constante h que se resta al crecimiento del recurso. No parece muy realista en el caso de una pesquería, pero podría funcionar en el caso de un bosque o de otro recurso fijo de fácil acceso. La ecuación que describe la dinámica del recurso incluyendo la explotación quedaría entonces expresada como: dB  B = r ⋅ B ⋅ 1 −  − h dt  K 59 (4.2) Fernando R. Momo - Ángel F. Capurro Utilizando las herramientas que introdujimos en los capítulos anteriores podemos analizar la existencia de puntos de equilibrio en este modelo y también la estabilidad de esos puntos. Veamos: dB B r  = r ⋅ B ⋅ 1 −  − h = rB − B 2 − h dt K  K r rB − B 2 − h = 0 ⇔ B = B * K ⇒ B* = − r ± r2 − 4h ⋅ r −2r K K = K ± K2 − 4hK r 2 Nótese que en el caso particular de que la velocidad de explotación h sea nula, la solución se simplifica y coincide con el caso ya conocido del modelo logístico que tiene dos puntos de equilibrio B*1 = 0 y B*2 = K. Cuando la explotación h es mayor que cero (es decir, se está extrayendo recurso) el punto de equilibrio inestable es mayor que cero y el punto de equilibrio estable es menor que K. No hemos analizado la estabilidad de esos puntos aunque afirmamos que el menor es inestable y el mayor es estable, ¿se atrevería el lector a hacer ese análisis? Vamos a suponer que sí y nosotros usaremos el camino del perezoso (en este caso el camino del biólogo), es decir, el análisis gráfico. Para hacerlo, mantendremos dos líneas separadas para representar la dinámica intrínseca del recurso y la explotación (esta última en color rojo); de esta manera quedará claro en qué situaciones la explotación excede al crecimiento del recurso y viceversa. El gráfico resultante es: 60 Ecología Matemática - Principios y aplicaciones Velocidades (Ton.ha-1.mes-1) Recurso explotado a velocidad constante K = 300 Ton/ha, r = 1.9 mes-1 160 crecimiento 140 explotación 120 100 80 60 40 20 0 0 B*1 50 100 150 200 B*2 250 300 350 Biomasa (Ton/ha) Como podemos ver, hay dos puntos de equilibrio (caracterizados porque las velocidades de crecimiento y de explotación del recurso son iguales y por lo tanto su crecimiento neto es nulo) B*1 y B*2. Sin embargo, su estabilidad es diferente. Si nos ubicamos a la izquierda de B*1, o sea, si tenemos una biomasa del recurso menor que B*1, la explotación supera al crecimiento intrínseco del recurso y su biomasa disminuirá irremediablemente hasta llegar a cero. Por el contrario, si nos ubicamos apenas a la derecha de B*1, veremos que la velocidad de crecimiento del recurso supera a la velocidad de explotación y por lo tanto la biomasa del recurso aumentará (velocidad neta de crecimiento positiva). Está claro que el punto B*1 es un repulsor y, por lo tanto, un punto de equilibrio inestable. Si sometemos al mismo análisis al punto B*2 nos encontramos con que es un atractor y representa un equilibrio estable. En síntesis, si el recurso es explotado a velocidad constante, podrá mantenerse en una biomasa estable (o mejor dicho “estacionaria”) menor que K pero tiene por así decirlo un “umbral de riesgo” representado por el equilibrio inestable B*1. Obsérvese que si la velocidad de explotación aumenta, la línea roja “sube” y corta la parábola en puntos más cercanos entre sí, es decir que los valo61 Fernando R. Momo - Ángel F. Capurro res de equilibrio estable e inestable se hacen más parecidos y cualquier fluctuación azarosa de la biomasa puede hacer que, estando en equilibrio estable, el sistema “salte” el umbral de riesgo y la biomasa empiece a disminuir. Incluso, podemos ver que si la velocidad de explotación alcanza un cierto valor crítico, ambos puntos colapsan en uno que pasa a ser un “punto de ensilladura” inestable. Ahora bien, sucede que eso pasa precisamente cuando la extracción es máxima en equilibrio, o sea, cuando la ganancia es máxima si el recurso tiene un valor comercial, o sea que ese riesgo no es despreciable. Obviamente, si la velocidad de explotación es mayor en cualquier biomasa que la velocidad intrínseca de crecimiento del recurso (la línea roja está por encima de la parábola) no existe ningún punto de equilibrio y la población explotada se extingue irremediablemente porque su velocidad neta de crecimiento es negativa para cualquier biomasa. El caso de la explotación proporcional a la biomasa En la explotación de una pesquería o de cualquier otro recurso biológico móvil al que hay que “atrapar”, la velocidad de explotación no suele ser constante. Cuando el recurso tiene menos biomasa es más difícil encontrarlo o más caro y la velocidad a la que se extrae baja. Optando nuevamente por el modelo más sencillo posible que incorpore esta característica, podemos suponer que h es ahora proporcional a B y que la constante de proporcionalidad está formada por la combinación de dos parámetros: el esfuerzo de explotación E y un coeficiente que depende de las características del recurso y de las artes que se utilicen para encontrarlo y atraparlo al cuál se denomina coeficiente de capturabilidad, simbolizado por la letra q. La ecuación (4.2) se transforma entonces en dB B  = r ⋅ B ⋅ 1 −  − q ⋅ E ⋅ B dt K  (4.3) Manipulando un poco esta ecuación obtenemos los valores de los puntos de equilibrio: 62 Ecología Matemática - Principios y aplicaciones dB r = (r − qE ) ⋅ B − B 2 = 0 dt K r   ⇒  r − qE − B  ⋅ B = 0 K   ∴ B1* = 0  qE  B2* = K ⋅ 1 −  r   Como puede verse, el equilibrio estable se ubica un poco por debajo de K, salvo que el esfuerzo de explotación E sea nulo (no se extrae el recurso) o el coeficiente de capturabilidad q sea nulo (el recurso no se puede atrapar y extraer con las artes disponibles). También se ve que el equilibrio B*2 es tanto menor que K cuanto más grandes sean el esfuerzo y la capturabilidad y cuanto más pequeña sea en proporción la tasa de crecimiento intrínseca (r) del recurso. Suena razonable, ¿no? Dejamos una vez más al lector inquieto las demostraciones matemáticas para concentrarnos en la representación gráfica (esto tiene una razón de ser como demostraremos más adelante). Manteniendo como líneas separadas el crecimiento intrínseco y la velocidad de extracción, llegamos al siguiente gráfico: Recurso explotado a velocidad variable Velocidades (Ton.ha-1.mes-1) K = 300 Ton/ha, r = 1.9 mes-1, q = 0.02, E = 5 160 140 120 100 80 60 40 20 0 crecimiento explotación 0 B*1 50 100 150 200 B*2 Biomasa (Ton/ha) 63 250 300 350 Fernando R. Momo - Ángel F. Capurro Ahora dejamos para el lector el siguiente ejercicio mental: ¿qué sucede con el gráfico cuando aumentamos los parámetros de proporcionalidad de la explotación, o sea q o E? Si se hace el ejercicio se observará que hay un cierto valor para la combinación de esos parámetros que hace que la línea de explotación corte a la de crecimiento en su punto más alto (el vértice de la parábola) y que entonces la velocidad de extracción del recurso (que se lee en la proyección del punto de intersección sobre el eje de ordenadas) es la máxima posible en equilibrio. Aquí ya no hay problemas de desestabilización. Existe un valor máximo de velocidad de extracción del recurso que lo mantiene en equilibrio estable con biomasa igual a la mitad de su capacidad de carga (K/2). Pueden deducirse fácilmente tanto el valor de esa velocidad de extracción, a la que llamaremos rendimiento sostenible máximo o RSM, como las condiciones que deben cumplir los parámetros para que el sistema tenga su punto de equilibrio allí: K r K rK rK rK RSM = hB* = r −   = − = 2 2 K 2 2 4 4 rK 4 = 1r qE = RSM = K K 2 2 2 2 ( ) Como se ve, este sencillo modelo es bastante versátil y puede ayudar a calcular cómo explotar más convenientemente un recurso. Para hacer eso, se deben introducir algunas cuestiones económicas relacionadas con la venta del recurso, su precio de mercado, etc. Pero esto excede nuestros modestos objetivos. Para el lector interesado recomendamos el libro de Colin Clark incluido en la bibliografía. Modelos epidemiológicos básicos El segundo tipo de ejemplos que examinaremos es el de un modelo epidemiológico que intenta capturar la dinámica de un sistema en el cual una población es infectada por algún tipo de agente patógeno y según el estado de salud de los individuos puede ser subdividida en subpoblaciones (por ejemplo en sanos, enfermos y 64 Ecología Matemática - Principios y aplicaciones recuperados). Los modelos continuos se utilizan en el caso en que la dinámica del agente patógeno es mucho más rápida que la del organismo enfermo y por lo tanto las transiciones entre estados (subpoblaciones) en los individuos de dicho organismo pueden considerarse instantáneas a los fines prácticos. Del mismo modo, la dinámica puede considerarse gobernada por la ley de acción de masas como veremos. Vamos a dar un ejemplo muy sencillo: una enfermedad en la que los individuos sanos se contagian de los enfermos (aquí los llamaremos “infecciosos” para usar la nomenclatura más común), los enfermos pueden recuperarse y volverse inmunes (a los que llamaremos simplemente “recuperados”). Si la enfermedad tiene un curso rápido, la población es estable y la enfermedad no es mortal, podremos expresar un modelo de su dinámica bastante simple. La representación “pictórica” de tal modelo es la siguiente: S βSI I νI R Aquí las letras mayúsculas representan los números de individuos en cada estado (S son susceptibles o sanos, I son los infecciosos o enfermos y R son los recuperados o inmunizados). Como se ve no aparecen ni la tasa de natalidad ni la de mortalidad de la población porque consideramos estable a la población mientras la enfermedad se propaga. Las flechas, que indican las velocidades de los flujos entre compartimientos responden a la ley de acción de masas. La primera, βSI, denota que el contagio de la enfermedad (el paso del estado sano al de enfermo) depende del producto del número de individuos sanos por el número de individuos enfermos por una tasa de contagio (β) que será mayor cuanto más contagiosa sea la enfermedad. Estamos evidentemente suponiendo que el contacto entre enfermos y sanos es azaroso y que se mezclan ambas subpoblaciones libremente o que, al menos, existe una cierta chance de que se encuentren que depende de cuantos son (ley de acción de masas). La segunda flecha representa la velocidad a la 65 Fernando R. Momo - Ángel F. Capurro que se recuperan y adquieren inmunidad los enfermos y por lo tanto ν es una tasa de recuperación o sea, la inversa del tiempo promedio que dura la enfermedad. No hace falta pensar ni calcular mucho para darse cuenta de que un sistema así tiende a un equilibrio en el que todos son finalmente inmunes, ¿o sí? Podríamos calcular el equilibrio igualando a cero las derivadas y recordando que la suma de los tres compartimientos o subpoblaciones permanece constante. O bien realizar algunas simulaciones. 1200 1000 S, I, R 800 600 S 400 I R 200 0 0 20 40 60 80 100 120 Tiempo (días) En la figura vemos cuál es la dinámica cuando β = 0.001 ind-1día-1 y ν = 0.02 día-1. Como vemos la subpoblación de sanos desciende rápidamente (a los diez días no queda casi ninguno) y la de enfermos alcanza su máximo al mismo tiempo. Ese máximo se visualiza como un pico “epidémico” de la enfermedad que poco a poco desaparece dando lugar a una población creciente de inmunizados. La población total permanece constante en 1003 individuos. Inicialmente sólo había 3 enfermos. Una dinámica así podría representar la dinámica de una enfermedad tal como el sarampión o la varicela introducida por unos pocos enfermos en una población totalmente susceptible (no hay 66 Ecología Matemática - Principios y aplicaciones inmunes ni por vacunación ni porque ya hayan sufrido la enfermedad). Un modelo como este, por los símbolos que utiliza, se llama modelo SIR y en general este tipo de modelos epidemiológicos básicos se conoce como modelos alfabéticos. Supongamos ahora que la inmunidad no es permanente y que se pierde a una tasa µ. ¿Variará la dinámica? El modelo sería ahora: µR S βSI I νI R Hemos agregado un flujo que, sin embargo, no hace variar sustancialmente la dinámica aunque sí cambia la posición del equilibrio. Estos modelos sirven además para explorar lo que podría pasar en diferentes circunstancias. Por ejemplo si introducimos una flecha (un flujo) que vaya directamente de S a R lo cuál equivale a vacunar. O lo que pasaría si la natalidad no está compensada por la mortalidad y todo el tiempo entran nuevos sanos a la población (suponemos que no hay transmisión “vertical” y por lo tanto los recién nacidos son sanos aunque sus padres no lo sean). Así podemos “jugar” y predecir lo que puede pasar con el curso de una enfermedad. Las herramientas matemáticas para hacerlo son siempre las mismas que venimos usando hasta ahora para analizar cualquier modelo dinámico continuo. Si se desea ver una dinámica más interesante, por ejemplo con presencia de epidemias y picos diferentes y con oscilaciones, lo más conveniente es agregar una natalidad no nula y repetir las simulaciones. El lector interesado puede encarar esa empresa. 67 Fernando R. Momo - Ángel F. Capurro Dinámica de un ecosistema marino planctónico Los ecosistemas marinos planctónicos han sido objeto de un estudio bastante extenso desde la ecología matemática. Probablemente el primer problema que se abordó con este enfoque fue la existencia, bastante evidente, de “parches” o manchones de plancton. En efecto, los organismos microscópicos que flotan en el agua del mar ven sujeto su destino al movimiento de las masas de agua que los arrastran y a sus propias tasas de natalidad y mortalidad. Sin embargo, en lugar de presentarse en forma infinitamente diluida o mezclada en el agua, sus poblaciones forman aglomeraciones separadas por aguas virtualmente “vacías” de ellos. Los primeros modelos para explicar ese patrón consideraron que los organismos se difundían en el agua en forma pasiva (como un soluto en una solución) pero existía además una producción de nuevos organismos en una forma no independiente de su densidad. Si el crecimiento de la población es lo suficientemente rápido como para compensar la dilución, el manchón se mantiene. Otro de los aspectos interesantes de discutir en estos sistemas es el hecho de que, además de diluirse, sedimentarse o reproducirse, los organismos del plancton se comen unos a otros. Concretamente, el zooplancton, o una fracción de él, se alimenta de fitoplancton y éste necesita de la luz del sol y los nutrientes del agua para fotosintetizar. Una de las preguntas clave en estos sistemas y en la ecología en general es ¿quién regula a quién? Es decir, ¿son los nutrientes y la luz quienes regulan todo a partir de un tope a la cantidad de alimento producida o son los depredadores los que limitan un crecimiento descontrolado de los autótrofos? A partir de algunos trabajos realizados en zonas costeras de la Antártida11 y como parte de una tesis de doctorado12 se desarrollaron algunos modelos matemáticos simples que ilustran el uso de las técnicas que hemos estudiado. Veamos: las zonas costeras de 11 Véase por ejemplo Schloss et al. 2002. Journal of Marine Systems, 39. O también Schloss y Ferreyra 2002. polar Biology, 25. 12 La tesista es Emma Ferrero, de la Universidad Nacional de Luján y a ella pertenece la elaboración de este modelo. 68 Ecología Matemática - Principios y aplicaciones aguas someras presentan en general muy baja biomasa y productividad biológica. Algunos de los factores que limitan la productividad fitoplanctónica son comunes a todos los ambientes antárticos; por ejemplo la intensa variación estacional de la luz, la presencia de hielos marinos, la presión de pastoreo o la intensa mezcla vertical. A estos factores se les suman entre otros la incidencia de fuertes vientos. Se puede plantear un modelo con tres compartimientos biológicos que representan la biomasa, tiene entradas y salidas que están influenciadas por las variables ambientales. Los tres compartimientos son fitoplancton, zooplancton y bacterias. Una vez planteadas las ecuaciones, se analiza la existencia de puntos de equilibrio y su estabilidad. El modelo matemático planteado está esquematizado en la siguiente figura, donde las flechas anchas representan flujos y las negras indican las influencias externas y sus efectos. El significado de cada símbolo y las unidades de las variables y parámetros se resumen en la tabla 4.1. Deshielo PASTOREO FOTOSINTESIS FITOPLANCTON ZOOPLANCTON MEZCLA VERTICAL MORTALIDAD HUNDIMIENTO MORTALIDAD Se plantea las ecuaciones diferenciales correspondientes a cada compartimiento y se obtiene el siguiente sistema: 69 Fernando R. Momo - Ángel F. Capurro α .I  dF  − R . F − ( c − b H m ) − F . Z q1 = F . Pm .tanh dt  Pm  dZ = − µ . z + q1 . e1. Z .F + q2 B.Z .e2 dt dB  B + d .F  = rB .B .1 −  − q2 . B. Z dt K   Donde F representa la densidad de las algas (por ejemplo en microgramos de clorofila-a por litro de agua), Z la densidad de los consumidores y B la densidad de las bacterias, siendo Pm la tasa máxima de fotosíntesis (por ejemplo en mg de carbono fijados por microgramo de clorofila-a), I la intensidad de la luz, R la tasa de respiración, H m la profundidad de la mezcla, q1 un factor de eficiencia de la depredación, q2 un factor de eficiencia de las bacterias, e1 el factor de eficiencia de asimilación del zooplancton, e2 el factor de eficiencia de asimilación del zooplancton, µ la tasa de mortalidad del zooplancton. La función c – b Hm representa la suspensión de algas (hundimiento), rB la tasa de crecimiento intrínseca de las bacterias, K máxima densidad de población de bacterias para la cual la tasa de crecimiento por capita es cero. Se supone que los nutrientes no son limitantes. En el modelo se considera que la tasa intrínseca de crecimiento del fitoplancton depende de su velocidad de fotosíntesis, expresada mediante una función tangente hiperbólica de la intensidad de luz . A ese crecimiento se le restan dos términos, el del hundimiento de células y el de depredación. El hundimiento se expresa como una función lineal de la profundidad de mezcla con pendiente negativa; esto permite que cuando la profundidad de mezcla supera un cierto valor, el hundimiento es “negativo”, es decir que la resuspensión supera a las pérdidas. La depredación se 70 Ecología Matemática - Principios y aplicaciones considera como una función lineal de las densidades poblacionales, siguiendo una ley de acción de masas. Como la velocidad de fotosíntesis es proporcional a la intensidad de la luz, pero las algas sufren una mezcla vertical por el arrastre del agua, es necesario integrar la luz en toda la columna de agua según la profundidad de la mezcla. El modelo para la intensidad de la luz es le siguiente: I promedio = 1 Zm ∫ Zm 0 I ( z ) dz = 1 Zm ∫ Zm 0 I 0 e −kz dz =I 0 e −kz Z mk 0 = Zm I0 1 − e −kZ m Z mk ( ) Donde el coeficiente de extinción de la luz, k, depende de N Np es decir, la concentración de partículas k = γ C p = γ p = γ Vm S .Zm es directamente proporcional al número de partículas que entraron por deshielo e inversamente proporcional a la superficie considerada y a la profundidad de la mezcla de la columna de agua, o sea, es mayor cuanto mayor es la concentración de partículas, por lo tanto, reemplazando k, se puede escribir I prom = 71 N I0S  − p  1 − e S  . Np   Fernando R. Momo - Ángel F. Capurro Tabla 4.1: Significado de los símbolos utilizados y las letras que simbolizan los diferentes parámetros del modelo. Símbolo F Z B Pm Significado Unidades Biomasa del fitoplancton Biomasa del zooplancton Biomasa del las bacterias Tasa máxima de fotosíntesis Intensidad de las luz Tasa de respiración Hm Profundidad de la mezcla -1 -4 Día 4,4 10 -1 -2 g C m día -2 -1 Mol m día -1 día 1,6 10 1 -4 1. 10 Metros -1 30 -1 qi Eficiencias del depredador Día (g de C) m ei Eficiencias de asimilación µ Tasa de crecimiento del zoo. Tasa de crecimiento de las bacterias rB 10 0,5 mg Cl-a m -2 -1 g Cl-a m día -2 C I R Valores -2 -2 -5 2. 10 0,15 -1 Día 0,07 Cuando se trabaja un poco con este modelo se obtienen diferentes dinámicas posibles para el sistema que dependen de los valores particulares que toman los parámetros. Sin embargo las conclusiones más fuertes son dos: por un lado, la persistencia del sistema como tal depende de que el fitoplancton no se hunda, es decir, que sea resuspendido por el movimiento del agua lo suficiente como para compensar las pérdidas. Por otro lado se puede ver que la eficiencia de depredación produce cambios dramáticos en el equilibrio del sistema. En síntesis, es un sistema controlado por dos vías: una ambiental (la resuspensión por mezcla) y otra del tipo top-down o de arriba hacia abajo (la depredación). Esta es una conclusión para nada trivial. Los sistemas planctónicos marinos en general y los antárticos en particular, están sometidos al efecto 72 Ecología Matemática - Principios y aplicaciones dramático de la radiación ultravioleta que se ha incrementado por el adelgazamiento de la capa de ozono estratosférico13. La evidencia experimental más reciente parece indicar que la radiación UV afecta en gran medida a los microdepredadores, sobre todo a los ciliados, y esto produce cambios en la composición y en la abundancia total de fitoplancton. En este contexto los modelos son de primera importancia porque permiten simular distintas situaciones y escenarios con incrementos probables de UV en el futuro. 13 De hecho este trabajo se realizó en el marco de un proyecto del Inter American Institute for Global Change research (IAI) sobre el problema del ozono y la radiación UV. Debemos agradecer al IAI su aporte monetario. 73 Capítulo 5 Modelos matriciales de poblaciones Cuando utilizamos modelos matemáticos para describir el crecimiento de las poblaciones o las interacciones entre ellas, asumimos que los individuos que las componen son iguales entre sí; de esta manera "olvidamos" intencionalmente, para no complicar excesivamente los modelos, la información acerca de las diferentes edades, estadíos, tamaños individuales o diferentes potenciales reproductivos de los individuos. Una de las formas de incluir parte de esa información en los modelos es el uso de modelos matriciales. En este tipo de modelos, la población es definida por un vector de estado cuyos elementos representan las abundancias de diferentes subgrupos de la población, por ejemplo edades o estadíos. Así, si en un momento dado la población que estamos estudiando está formada por 30 individuos juveniles, 57 adultos y 12 seniles, podemos representarla por un vector 30 r   n = 57  12  (5.1) donde los números son los elementos del vector y representan la abundancia de cada estadío. Ahora podemos incluir en nuestro modelo información acerca de las diferentes probabilidades de supervivencia y potenciales reproductivos de cada estadío. Nótese que estamos considerando un crecimiento discreto de la población y, entonces, esa información estará sintetizada en una matriz de transición cuyos elementos representan los aportes de cada elemento del vector (cada estadío) en un momento dado (por ejemplo un año) a cada uno de los estadíos en el momento siguiente (por ejemplo el año siguiente). A la matriz de transición le llamaremos A, de manera que podemos sintetizar nuestro modelo como: 75 Fernando R. Momo - Ángel F. Capurro r r nt +1 = A ⋅ nt (5.2) Por ejemplo podemos escribir el modelo de una población con tres estadíos:  n1   s11 n  =  s  2  21  n3  t +1  s 31 f2 s 22 s 32 f 3   n1  a 23  n2  s 33   n3  t (5.3) Vemos que los elementos de la primer fila de la matriz son los aportes de cada estadío al estadío inicial en el año siguiente, es decir la supervivencia del estadío 1 y las fecundidades de los estadíos 2 y 3. Los elementos de la subdiagonal representan los aportes de cada estadío al siguiente en el año siguiente. Los de la diagonal, los aportes de cada estadío a sí mismo en el año siguiente. Algunos de estos elementos o de los restantes pueden tener sentido biológico sólo cuando son nulos y eso depende de cómo estemos definiendo los estadíos y el modelo. Por ejemplo, si los estadíos son edades cronológicas absolutas, los elementos diagonales distintos de cero no tendrían sentido porque ningún organismo puede tener la misma edad cronológica dos años seguidos, tampoco el elemento a23 del ejemplo porque un organismo no puede ser más joven al año siguiente de lo que era este año; en cambio, si los estadíos son diferentes estados fenológicos de una planta, puede haber "vueltas atrás" o aportes de un elemento a sí mismo. Esto se comprende mejor si representamos el ciclo de vida en forma de grafo: 76 Ecología Matemática - Principios y aplicaciones f2 s33 s22 s11 1 s21 2 s32 3 a23 f3 Todo grafo orientado como este tiene asociada una matriz, cada flecha corresponde a un elemento no nulo de la matriz donde el elemento que aporta es el que corresponde a la columna y el que recibe a la fila. Si damos valores arbitrarios a las flechas del grafo anterior y escribimos la matriz resultante, obtenemos: 18  0  5  = 1    2  4   0 0 0 1 3 6 10 0 12 0  3  (5.4) ¿Por qué decimos que los modelos matriciales nos aportan más información que sus equivalentes simplificados? Para comprender esta afirmación vamos a empezar por un ejemplo clásico: la matriz de Bernardelli que es la matriz de transición del ejemplo numérico anterior. Ésta matriz fue originalmente (1941) planteada para un problema demográfico pero para nuestros fines inventaremos una interpretación más ecológica. Supondremos que representa el ciclo de vida de un insecto; también que representa sólo la dinámica de las hembras, tal y como se hace en las tablas de vida, para que tenga sentido el valor de fecundidad. Se consideran tres estadíos: larvas, pupas y adultos: Esto significa que en el valor de fecundidad estamos colapsando la fecundidad sensu stricto y la 77 Fernando R. Momo - Ángel F. Capurro supervivencia de los huevos a larvas; suponemos además que ningún estadío sobrevive como tal al año siguiente: las larvas pasan a pupas o mueren, las pupas pasan a adultos o mueren y las hembras adultas m,ueren después de reproducirse; esto asegura que todos los elementos diagonales tendrán valor 0. Por último, y obviamente, sólo las hembras adultas se reproducen y por lo tanto sólo la f3 no es nula. El grafo que representa este ciclo de vida es el siguiente: Si construímos la tabla de vida correspondiente y calculamos la tasa de crecimiento discreta R0 de esta población obtenemos: R0 = f 1 + s 21 ⋅ f 2 + s 21 ⋅ s 32 ⋅ f 3 = 0 + 1 ⋅ 0 + 1 ⋅ 1 ⋅ 6 = 1 2 3 2 (5.5) Como vemos, la población es estable a largo plazo, no crece ni decrece, ¿eso es todo lo que podemos decir? Observemos qué pasa si simulamos el comportamiento de la población utilizando nuestro modelo matricial; definimos un vector de estado inicial para la población, multiplicamos la matriz de transición por ese vector y obtenemos el vector de la población al año uno (es lo que está representado en la ecuación 4), multiplicamos la matriz por este vector y obtenemos el vector del año dos y así sucesivamente; ¿cuál es el resultado de este proceso? 78 Ecología Matemática - Principios y aplicaciones Tiempo 0 1 2 3 4 5 6 7 8 9 10 Larvas 10 18 24 10 18 24 10 18 24 10 18 Pupas 12 5 9 12 5 9 12 5 9 12 5 Adultos 3 4 1.67 3 4 1.67 3 4 1.67 3 4 Más allá de la presencia de valores fraccionarios que, en realidad, deberíamos interpretar como densidades, aparece aquí un comportamiento que la tabla de vida no permitía sospechar, hay oscilaciones periódicas en la abundancia de la población, ¿podríamos haber predicho esto a partir del análisis de la matriz? La respuesta es sí. El análisis básico de un modelo matricial Hay una manera simple de conocer cuál será la dinámica de la población representada en un modelo matricial: se trata de calcular los autovalores (o valores propios) de la matriz. Este análisis es semejante a los que realizamos anteriormente con modelos discretos sin estructura poblacional; si el autovalor dominante tiene módulo 1 (salvo en los casos especiales en que vale -1, i o -i), la población tendrá un comportamiento asintótico estable. La población puede tener una estructura estable de edades, a la que tenderá tras un número muy grande de iteraciones. Esto se produce si la matriz es primitiva y por lo tanto el primer valor propio es estrictamente dominante. La estructura estable de edades está dada por el autovector derecho (columna) del autovalor dominante (aclararemos esto más abajo con un ejemplo). Si el autovalor dominante tiene una parte imaginaria no nula, la población presentará un comportamiento oscilante. Si el autovalor dominante toma 79 Fernando R. Momo - Ángel F. Capurro valor -1, la población "salta" entre dos estados, si toma valor i o -i el salto se produce entre cuatro estados distintos, en los tres casos tenemos comportamientos periódicos. Pero no toda la información acerca de la dinámica está contenida en el valor propio dominante. Los otros autovalores, cuando tienen módulo menor que el dominante, gobiernan el comportamiento transitorio o temporal de la población (comportamiento "trasiente" si efectuamos una traducción algo bestial del inglés). Estos comportamientos temporales pueden ser muy importantes porque, a veces, el comportamiento asintótico se alcanza luego de un número de iteraciones tan grande que deja de tener interés biológico ya que lo que realmente se puede observar es el comportamiento transitorio. En el ejemplo de la matriz de Bernardelli aludido anteriormente, los tres autovalores tienen el mismo módulo (la matriz es “imprimitiva” de índice tres) y ninguno domina sobre los demás; como dos de ellos tienen parte imaginaria no nula (de hecho son complejos conjugados) la población oscila con período tres aunque en promedio sobre una serie muy larga de datos, la población es estable. Veamos cuál sería el resultado de la simulación si la fecundidad de la población fuese menor, por ejemplo 4. En el gráfico se observa que, si bien se mantiene el comportamiento oscilatorio, el tamaño promedio de la población declina, como lo indica su autovalor dominante con módulo menor que uno. También podemos ver que la proporción promedio de cada estadío sobre el total de la población tiende a estabilizarse y es la que predice el autovector derecho del autovalor dominante. 80 Ecología Matemática - Principios y aplicaciones Número de individuos Tamaño poblacional vs. tiempo 18 16 14 12 10 8 6 4 2 0 0 5 10 15 20 25 30 tiempo (generaciones) Más información disponible en la matriz Cada valor propio (autovalor) de la matriz tiene asociado dos autovectores, uno llamado “izquierdo” que es un vector fila, y otro llamado “derecho” que es un vector columna. El autovector izquierdo del autovalor dominante de una matriz poblacional indica el valor reproductivo de cada estadío en la población final, es decir, en qué proporción aporta cada estadío al futuro de la población; mientras que el autovector derecho indica la proporción asintótica de estadíos en la población si la hubiere, es decir, si la población alcanza una distribución estable de edades, ese autovector nos dice cuál será. También se puede saber cómo se modifica el autovalor (y por lo tanto la dinámica de la población) al variar cualquiera de los elementos de la matriz. Esto se conoce como análisis de sensibilidad de la matriz. Consiste técnicamente en medir el porcentaje de variación del módulo del autovalor en función del porcentaje de variación del valor de cada elemento suponiendo que los demás permanecen constantes. En términos más elegantes: 81 Fernando R. Momo - Ángel F. Capurro S ij = ∂ λ1 ∂aij (5.6) donde Sij indica “sensibilidad”. Señalemos que este tipo de análisis puede ser de gran importancia en el caso de estar pretendiendo conservar una población, o estudiando los efectos de un xenobiótico sobre su dinámica o en el caso de querer controlarla. Podemos usar este análisis para optimizar nuestra intervención y detectar cuáles son los estadíos más sensibles a lo largo del ciclo de vida de la especie. Llamando vi al iésimo elemento del autovector izquierdo (valor reproductivo) y wj al jotaésimo elemento del autovector derecho (el de distribución asintótica de edades); entonces: ∂λ1 ∝ v iw j , ∂aij (5.7) es decir que los autovectores del valor propio dominante contienen la información acerca de la sensibilidad de ese valor propio a cada elemento de la matriz de transición. En los ejemplos que hemos mostrado supusimos que la matriz de transición es una matriz constante, sin embargo, esa matriz puede tener coeficientes variables que dependen del tiempo o bien de la densidad de la población o de ambos. Esto da una gran riqueza de posibilidades aunque complica el tratamiento matemá– tico. 82 Capítulo 6 De la ecología teórica a la teoría ecológica En este capítulo discutiremos algunos temas teóricos que no tiene que ver estrictamente con la modelación matemática aunque la utilicen. La idea es ir un poco más allá y profundizar las cuestiones conceptuales que hacen a la ecología y sus temas conflictivos donde se necesita desarrollar teoría y donde la ecología matemática tiene cosas que decir. Las poblaciones humanas La población humana tiene características particulares que hacen que su papel y su efecto en los ecosistemas revistan una escala única. ¿Cuáles son esas características distintivas? Para contestar esto comencemos por repasar nuestro esquema “preescolar” del ecosistema (figura 6.1). Hemos distinguido tres grandes compartimientos biológicos, los productores, los consumidores y los descomponedores, y hemos esquematizado el flujo de la energía y el ciclo de la materia a través de esos compartimientos. Nadie dudaría en incluir a la población humana dentro del compartimiento de los consumidores, pero esto no nos aclararía demasiado. Si queremos destacar lo distintivo del Homo sapiens, hay que completar el esquema y mejorarlo. Está claro que los límites de los ecosistemas tienen cierto grado de arbitrariedad pero, sin embargo, no son totalmente arbitrarios. De hecho sabemos que hay ciclos de producción, consumo y descomposición que “naturalmente” se concentran en una cierta región a la que reconocemos como un ecosistema. También sabemos que entre distintos ecosistemas hay un cierto grado de intercambio de materia y energía, casi siempre menos intenso. Agreguemos estos últimos conceptos al esquema preescolar (figura 6.2). Para comprender la escala de impacto del ser humano tenemos que incluir algo más: la energía fósil. ¿Por qué? La población humana se distingue ecológicamente del resto de las pobla83 Fernando R. Momo - Ángel F. Capurro ciones del planeta por su enorme capacidad de transporte de materiales. Homo sapiens puede transportar grandes cantidades de materiales a grandes distancias; de esta manera “abre” los ciclos de la materia o deberíamos decir que los “estira” horizontalmente. Pero el transporte requiere energía y la que los humanos usamos para esta tarea no es energía de nuestro metabolismo sino un tipo particular de la energía llamada exosomática, es decir, energía que circula por afuera de los organismos vivos; es la energía fósil, es decir, energía fotosintética acumulada hace millones de años y transformada por el tiempo y las condiciones ambientales y geológicas en combustibles tales como carbón o petróleo. Figura 6.1: El esquema básico del ecosistema 84 Ecología Matemática - Principios y aplicaciones Figura 6.2: el mismo esquema pero con las importaciones y exportaciones de energía y materia. En algunos momentos de su historia, Homo sapiens ha utilizado otros tipos de energía exosomática no fósil, como por ejemplo la energía del viento que impulsó las naves que comerciaron, robaron o esquilmaron cubriendo distancias cada vez mayores hasta llegar a las del orden del planeta entero. No obstante, nuestra fuente actual de energía exosomática es principalmente fósil. De esta manera, cerramos un ciclo de la materia que se pierde en el tiempo. ¿Para qué usamos esa energía y ese transporte y cómo influye esto en la dinámica de los ecosistemas? Evidentemente, el uso de energía fósil no es la única novedad en la población humana; también se trata de la primera población que puede transmitir a su descendencia una enorme cantidad de información por vía no genética. En otras palabras, Homo sapiens tiene cultura; una forma particular de información y un cúmulo particular de actividades que no siempre persiguen fines inmediatos y evidentes de supervi85 Fernando R. Momo - Ángel F. Capurro vencia y reproducción (los dos “objetivos” biológicos básicos). La energía fósil no sólo se utiliza para obtener y distribuir grandes cantidades de comida y construir refugios y habitáculos diversos, sino también para mantener una enorme variedad y cantidad de actividades sociales de tipo cultural en el más amplio sentido de la palabra. Gracias a la energía fósil podemos alumbrarnos y leer, enseñar y aprender, conservar nuestros alimentos, fabricar diversos soportes de información, desde libros a computadoras, comunicarnos a grandes distancias, comerciar, explorar o construir rascacielos. Pero volvamos al uso de energía fósil en los procesos de producción y transporte de alimentos. Aquí, esa energía se usa para rotular la tierra, borrando así heterogeneidades y pudiendo cultivar unas pocas especies de plantas en grandes extensiones; también usamos energía fósil para producir fertilizantes y muchos de los pesticidas que nos permiten controlar ese proceso productivo. Usamos energía fósil para extraer agua y regar, para cosechar y transportar, para cocinar, moler, tostar, separar, fraccionar, envasar y todas las operaciones elementales del proceso de producción agropecuaria e industrialización de sus productos. Estas acciones tienen consecuencias. La población humana no está distribuida uniformemente en el espacio y desde el siglo dieciocho en adelante (aunque esta tendencia ya existía) se ha concentrado aceleradamente en algunos lugares que no son los mismos en que se producen los alimentos y otras materias primas. Es decir que las cosas se consumen lejos de donde se producen. Además, las grandes concentraciones humanas “producen” (usando ahora el término en su acepción común, no con su significado ecológico) enormes cantidades de desechos (materia “usada”); esos desechos generalmente se transportan activamente, como la basura, o pasivamente, como los desechos cloacales en los ríos, a otros lugares. Así, las poblaciones humanas siguen “abriendo” el ciclo de la materia (figura 6.3) y esto genera varios, sino todos, los problemas que reconocemos como problemas ambientales o ecológicos: contaminación, sobreexplotación, erosión, pérdida de diversidad ecológica, etc. 86 Ecología Matemática - Principios y aplicaciones Figura 3: varios ecosistemas interrelacionados merced a la acción humana. El hecho de que la producción de materias primas, la producción de bienes, el consumo de ambos y la acumulación de desechos sean procesos que tienen lugar en momentos y lugares diferentes, a veces muy lejanos tanto en tiempo como en espacio, tiene otras consecuencias ecológicas importantes. Una de ellas es que se pierde la percepción social e individual de las consecuencias ambientales de las demandas y de las actividades. Dicho de otra manera y con un ejemplo: dado que consumimos papel en un lugar muy lejano al de su producción y dado que esa producción es anterior a nuestra demanda del producto, no tenemos ninguna percepción de los posibles efectos ambientales de esa producción (tala de árboles, contaminación de ríos, etc.); tampoco tenemos por lo tanto ninguna capacidad ni posibilidad de regular nuestro consumo en función de sus consecuencias ambientales. La legislación intenta compensar esta falta de lazos de retroalimentación, sin embargo es 87 Fernando R. Momo - Ángel F. Capurro un instrumento tosco y lento para un conjunto de procesos sutiles y vertiginosos. La termodinámica de los sistemas ecológicos Varios ecólogos, pero especialmente Margalef, Jorgensen y Howard Odum, han insistido en la necesidad de desarrollar una teoría adecuada para enmarcar la ecología en una termodinámica de los sistemas abiertos y alejados del equilibrio. Los intentos de formalización de esta teoría han sido variados y realmente audaces. Veamos algunas de estas formalizaciones y cómo se relacionan con los modelos matemáticos14. Los ecosistemas, las poblaciones y, en general los sistemas biológicos son estructuras disipativas. Esto es, sistemas termodinámicos abiertos, alejados del equilibrio pero que mantienen un estado más o menos estacionario a costa de un intercambio permanente de materia y energía con el exterior. Según Prigogine15, un sistema de estas características tiene una tendencia natural a minimizar su producción de entropía por unidad de tiempo dS/dt. ¿Cómo lo hace? Sabemos por el teorema de fluctuacióndisipación que el sistema alcanza el estado estacionario disipando el exceso de entropía a través de fluctuaciones en sus "modos naturales" de fluctuación. En las poblaciones, un modo natural de responder a la disponibilidad de energía del medio, es decir, a la cantidad de energía (en forma de alimento) que fluye a través de la población, es variar la proporción de energía invertida en diferentes actividades (crecimiento, reproducción). Como la cantidad de energía que debe invertirse en el crecimiento de una hembra, normalmente es diferente de la que debe invertirse en el crecimiento de un macho, es natural suponer que los individuos de cada dos sexo tendrán diferentes probabilidades 14 Este trabajo fue realizado y publicado por uno de nosotros (Momo) y José Ure de la Universidad Nacional de General Sarmiento. 15 Prigogine, Y., 1967. Introduction to thermodynamics of irreversible processes. 88 Ecología Matemática - Principios y aplicaciones de sobrevivir. El resultado serán variaciones en la proporción de sexos según la energía disponible. Aquí es donde se plantea una analogía que puede ser fructífera: consideramos que la población de adultos es un macroestado con dos microestados elementales posibles: macho o hembra. Si no consideramos el valor energético, la probabilidad de ambos microestados es igual y vale 1/2. Considerando un pequeño desvío del valor medio (m) podemos estudiar su probabilidad: N1 = N/2 + m N2 = N/2 - m o sea que: N1 - N2 = 2 m Pero ese desvío depende de la energía del sistema según la siguiente relación: ENERGÍA = - m H (N/2 + m) Si asignamos a los dos microestados elementales energías diferentes (lo que implica probabilidades diferentes) y sabiendo que no hay independencia entre las elecciones individuales (la probabilidad de cada individuo de ser hembra o macho depende de las "elecciones" de los demás) obtenemos que la proporción entre microestados (hembras/machos = β) depende sólo de la diferencia energética entre ambos 2ε (equivalente a 2m en la formulación equiprobable). Entonces: β = exp (- 2 ε / T) donde el parámetro T se relaciona con la energía media U del sistema (población) mediante: dS / dU = 1 / T 89 Fernando R. Momo - Ángel F. Capurro Si estuviéramos hablando de partículas elementales, T sería la temperatura, pero aquí tenemos que darle otro nombre o no darle ninguno. Sin embargo, su significado biológico puede ser más jugoso. Veamos; podemos escribir: dS / dU = (dS/dt) / (dU/dt) Aquí, el denominador (derivada de la energía respecto del tiempo) es la producción bruta (PB) del sistema; y el numerador (derivada de la entropía respecto del tiempo) es velocidad de disipación, es decir, una medida de la ineficiencia del sistema, ya que mide su velocidad de producción de entropía. Cada sistema ecológico tenderá a minimizar el cociente dS/dt; en qué valor lo consigue depende de restricciones genéticas y fisiológicas y de la disponibilidad energética del medio (su producción). Sistemas inmersos en un ambiente rico (mayor dU/dt) podrán mantener estados estacionarios con valores mayores de disipación. El modelo se aplicó al estudio de poblaciones del anfípodo costero antártico Gondogeneia antarctica; parametrizando el costo energético en proporción a las biomasas corporales. Se encontró que, en este organismo, la proporción de sexos (hembras/machos) responde a la disponibilidad de alimento en forma de microalgas pero con un desfasaje de 98 días (3.26 meses). Dado que en el período observado las biomasas de machos y hembras varían, la eficiencia de la población en el uso de la energía no es constante; sin embargo, la biomasa de este anfípodo guarda una estrecha relación con la disponibilidad de algas (considerado el desfasaje). Calculamos el valor de 1/T para el período considerado a partir de los valores de producción secundaria, mayor en primavera y verano. Esta magnitud sería un indicador de la ineficiencia poblacional. Los máximos de esa ineficiencia se producen a fin del invierno (septiembre), cuando el estado nutricional de los organismos es peor; y en noviembre, después de la reproducción primaveral, cuando la mortalidad de hembras es muy grande. Por otra parte se observa un pico de valor mínimo donde 1/T es menor que cero, esto significa que el sistema está organizándose y dis90 Ecología Matemática - Principios y aplicaciones minuyendo su entropía interna; esto coincide con el momento de máxima inversión energética en crecimiento y producción de juveniles (octubre y mediados de noviembre). En las comunidades, consideradas como macroestados, los microestados elementales son las especies. En ese caso, la entropía de la comunidad viene dada por el conocido índice de ShannonWienner, cuya fórmula es equivalente a la definición de entropía de Boltzman: H = - Σ pi log pi Sabemos que existe una relación entre las variaciones de la entropía en el tiempo, el flujo de energía en el sistema y una "fuerza" de flujo o propiedad intensiva: dS/dt = X dQ/dt es decir que podemos estudiar los cambios de diversidad de las comunidades en el tiempo en relación con la producción bruta (dQ/dt) y obtener así valores para la fuerza de flujo X. Esta fuerza nos daría una medida de la sensibilidad que la comunidad manifiesta con respecto a la energía, es decir que debería estar relacionada con su inestabilidad intrínseca o sensibilidad y puede usarse para evaluar la fragilidad de los sistemas. Como se ve, aunque en forma precaria todavía, algunos temas teóricos de la ecología pueden abordarse mejor con el uso de modelos matemáticos que nos permiten formalizar y razonar más claramente sobre los sistemas que estamos intentando estudiar y sobre las consecuencias dinámicas de las suposiciones biológicas o físicas que hacemos. 91 Capítulo 7 Libro de cocina para hacer tu propio modelo matemático Ha llegado un momento crucial de este libro. Se supone que con lo que el amable lector ha devorado y aprendido hasta aquí, ya debería ser capaz de enfrentarse casi con cualquier modelo matemático de un proceso ecológico. Sin embargo, no hemos hablado casi del grave problema de hacer un modelo por nuestros propios medios cuando solo tenemos unos pocos datos, o alguna intuición sobre el problema, o a veces sólo un enunciado verbal de lo que está pasando. En este capítulo daremos algunas indicaciones simples para poder llevar a buen término la empresa de escribir el modelo matemático, analizarlo y, si es posible, publicarlo. Primer paso: delimitar claramente el problema ecológico El problema a encarar tendrá seguramente su origen en la observación, ya sea del resultado de experimentos o de relevamientos de campo. Seguramente esas observaciones sugieren algún tipo de hipótesis explicativa del tipo lo que está pasando parece ser...; o bien algún tipo de pregunta como ¿qué pasaría si...? Esa hipótesis o esa pregunta pueden ser expresadas en forma de un modelo. Veamos un ejemplo que es el que desarrollaremos en todo el capítulo16. En una serie de experimentos para evaluar el efecto de un aumento en la radiación ultravioleta (en adelante RUV) sobre las comunidades de algas del fitoplancton, se observó que, con las dosis más altas de RUV, las comunidades se simplificaban al punto de ser casi uniespecíficas. Pero además, el alga sobreviviente resultó ser un alga responsable de mareas tóxicas (tipo “marea roja”). Surgió inmediatamente la pregunta, ¿podría ser posible que un au16 El problema ecológico planteado aquí corresponde a trabajos realizados por Marcelo Hernando, de CONICET, bajo la dirección de Gustavo Ferreira, del Instituto Antártico Argentino. 93 Fernando R. Momo - Ángel F. Capurro mento de la RUV ambiental (producto de una disminución en el ozono estratosférico) aumente la probabilidad de floraciones tóxicas en aguas subantárticas? He allí el problema, enunciado en forma de pregunta o posibilidad. Es obvio que este tipo de preguntas no tienen una respuesta experimental inmediata y como, además, están formuladas en términos de una posibilidad, no estamos seguros de que valga la pena gastar tiempo y dinero en realizar una serie de experimentos a gran escala o durante más tiempo. ¿Podrán ayudarnos los modelos matemáticos? Si logramos pasar a la segunda etapa, probablemente la ecología matemática no proporcione parte de la respuesta. Segundo paso: identificar qué se puede modelizar y cuáles son los fenómenos relevantes Ahora tenemos que examinar qué cosas estamos pensando pero no estamos diciendo (lo que los científicos llamamos pomposamente “los supuestos”). Si tememos que la RUV favorezca de alguna manera a las algas tóxicas estamos asumiendo que éstas son menos dañadas, es decir, son más resistentes17, cuando les llega RUV. El origen de esa resistencia no nos preocupa todavía, basta que pueda expresarse en una forma sencilla matemáticamente. Como en las condiciones originales, antes de agregar RUV, estas algas tóxicas no estaban produciendo floraciones, supondremos que originalmente las algas no tóxicas les llevaban ventaja a sus compañeras venenosas. Como un ecólogo matemático es fundamentalmente una persona perezosa, lo primero que hacemos es echar mano de modelos preexistentes. Parece que estamos expresando en palabras un fenómeno de competencia con coeficientes que varían según la dosis de RUV que los organismos reciben. Esto puede trabajarse modificando levemente el modelo de Lotka (véase el capítulo 3 de 17 Otra posibilidad es que la RUV favorezca a estas algas indirectamente, por ejemplo mineralizando nutrientes que las algas aprovechan. No usaremos esta idea para no tener que agregar otra variable al modelo. 94 Ecología Matemática - Principios y aplicaciones este libro). Las variables a estudiar son dos: biomasa de algas tóxicas (N2) y biomasa de algas no tóxicas (N1). Nótese que hemos simplificado el sistema poniendo todas las especies tóxicas en un solo paquete y todas las “buenitas” en otro. Los parámetros que entrarán en el modelo son las tasas de crecimiento de las algas, los coeficientes de saturación (1/K), los coeficientes de competencia entre las algas y la mortalidad extra producida en cada una por la RUV. Aquí estamos suponiendo que la RUV no afecta los coeficientes de competencia en sí, sino que simplemente mata más o menos a cada tipo de alga. Ya estamos listos para dar el siguiente paso. Pero antes resumamos lo que acabamos de explicar: 1. Suponemos que hay dos “especies”, una buena y una tóxica. 2. Suponemos que ambas crecen según un modelo logístico y que compiten entre sí. En condiciones normales, la buena es mejor competidora que la otra pero no necesariamente la excluye. 3. Suponemos que la RUV introduce un factor de mortalidad o de frenado del crecimiento y que ese factor crece con la dosis de RUV linealmente (¡!). 4. Suponemos que la pendiente de esa relación lineal es mayor para las algas buenas (más sensibles) y menor para las tóxicas (resistentes). Hecha esta arbitrariedad, se procede a expresarlo como ecuaciones diferenciales: Tercer paso: escribir las ecuaciones del modelo Este paso supone tener en la cabeza una lista de las variables y parámetros que comentamos antes. Hagámoslo, a partir del modelo de competencia clásico: 95 Fernando R. Momo - Ángel F. Capurro dN1 = N1 (r10 − aN1 − b12 N 2 − s1UV ) dt dN 2 = N 2 (r20 − gN 2 − b21 N1 − s2UV ) dt donde: N1 es la abundancia de la buenita, r10 su potencial biótico (máxima tasa de crecimiento), a es su coeficiente de competencia intraespecífica o sea r10/K1, b12 es el efecto de competencia interspecífica que le producen las malitas, s1 es la sensibilidad que tiene a la RUV y UV es la dosis de RUV que recibe. Los coeficientes de la segunda ecuación son simétricos y equivalentes a los que mencionamos recién. Suponemos que b12 < b21 y que s1 > s2 y hacemos la primera trampa del proceso. Esta trampa la podemos hacer gracias a las computadoras y consiste en realizar algunas simulaciones numéricas y “jugar” con los coeficientes para ver qué pasa en distintas situaciones hipotéticas aunque siempre cuidando de no violar las suposiciones iniciales y las relaciones entre parámetros que acabamos de enunciar al principio del párrafo. El resultado es alentador: biomasa (mg de Cl-a/L) 0.6 Competencia entre algas con baja RUV 0.5 0.4 N1 0.3 N2 0.2 0.1 0 0 10 20 tiempo (días) 96 30 Ecología Matemática - Principios y aplicaciones Competencia entre algas con alta RUV biomasa (mg de Cl-a/L) 0.6 0.5 0.4 0.3 N1 0.2 N2 0.1 0 0 10 20 30 tiem po (días) Considere el lector que todos los parámetros se mantuvieron constantes y sólo se aumentó la cantidad o dosis de RUV que las algas reciben. Eso bastó para invertir el resultado de la interacción competitiva. ¿Podemos estar conformes con esto? Podemos, pero no deberíamos. Hemos aprendido que podemos extraer del modelo más información. Por ejemplo condiciones de equilibrio, estabilidad, sensibilidad, etc. Pasemos entonces al cuarto paso. Cuarto paso: estudiar la existencia, condiciones y tipos de equilibrio Recordemos nuestra caja de herramientas. Sabemos encontrar la expresión de los puntos de equilibrio del sistema. Sólo tenemos que igualar a cero las dos ecuaciones diferenciales y despejar los valores de N1 y N2 que verifican esas igualdades. Veamos: 0 = N1 (r10 − aN1 − b12 N 2 − s1UV ) 0 = N 2 (r20 − gN 2 − b21 N1 − s2UV ) 97 Fernando R. Momo - Ángel F. Capurro Hay un caso obvio que verifica las igualdades y es el caso en que ambas poblaciones son nulas. Ese es un caso que no nos interesa. Veamos la otra posibilidad: 0 = r10 − aN1 − b12 N 2 − s1UV 0 = r20 − gN 2 − b21 N1 − s2UV Esto determina dos isolíneas, similares a las que vimos en el capítulo tres con el modelo de competencia. Las dos rectas están definidas de la siguiente forma: N2A = r10 − s1UV a − N1 b12 b12 N 2B = r20 − s 2UV b21 − N1 g g Si hacemos variar los valores de los diferentes parámetros veremos que podemos generar las mismas situaciones que en un modelo clásico de competencia: que gane una u otra especie siempre, que coexistan en forma estable o que coexistan en forma inestable y gane una u otra según las condiciones iniciales. Sin embargo nos interesa especialmente el efecto del parámetro UV manteniendo los demás valores constantes. Cuando la cantidad de RUV recibida es pequeña, las algas no tóxicas triunfan en la competencia y eliminan a las tóxicas: 98 Ecología Matemática - Principios y aplicaciones Diagrama de fases del modelo caso de eliminación de las algas tóxicas 1.6 1.4 N2A N2B 1.2 N2 1 0.8 0.6 0.4 0.2 0 0 0.2 0.4 0.6 0.8 1 N1 Con RUV algo más alta, podemos tener coexistencia estable de ambos tipos de algas: Diagrama de fases del modelo caso de coexistencia estable 0.9 0.8 N2A N2B 0.7 0.6 N2 0.5 0.4 0.3 0.2 0.1 0 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 N1 Por último, si la radiación es más intensa, las algas tóxicas ganan la competencia siempre: 99 Fernando R. Momo - Ángel F. Capurro Diagrama de fases del modelo caso de floraciones tóxicas 0.6 N2A 0.5 N2B 0.4 N2 0.3 0.2 0.1 0 0 0.1 0.2 0.3 0.4 0.5 0.6 -0.1 N1 Hemos avanzado bastante. Con un modelo muy sencillo, casi diríamos caricaturesco, logramos demostrar que existe la posibilidad de que el aumento de la RUV genere floraciones de algas tóxicas simplemente porque éstas son menos afectadas por este tipo de radiación, cosa que podría suceder si son capaces de sintetizar alguna sustancia protectora de la radiación como las micosporinas u otro tipo de “pantallas” solares naturales. Quinto paso: hilando finito Lo que viene después es a gusto del consumidor pero básicamente consiste en mejorar el modelo para que represente mejor el fenómeno y, en cada mejora, ir analizando otra vez su dinámica y predicciones. Por ejemplo, se podría mejorar suponiendo que la mortalidad producida por la RUV no crece linealmente con la dosis sino de otra manera, ya sea acelerada o asintótica. Podríamos incluir alguna expresión de la fotosíntesis en la ecuación de crecimiento (el parámetro r ya no sería constante sino que dependería de la cantidad de luz visible recibida). Y así podríamos ir agregando elementos de realismo al sistema. ¿Vale la pena? Recordemos el viejo axioma que dice que de un modelo pedimos tres cosas: simplici100 Ecología Matemática - Principios y aplicaciones dad, realismo, exactitud. Pero no podemos maximizar las tres simultáneamente. Un modelo simple puede ser poco exacto o poco realista, pero tiene la ventaja de ser muy general. En el caso ejemplificado, y seguramente en muchos otros, conviene no exagerar con el realismo porque no sabemos tanto del sistema y porque estamos haciendo generalizaciones gruesas. Entonces el modelo cumple la función de ayudarnos a mejorar nuestras preguntas y, eventualmente, a rediseñar nuestros experimentos o métodos de observación. Otro tópico interesante es la comparación del modelo con los datos. Esto puede hacerse al mismo tiempo que un ajuste para determinar ciertos parámetros o bien usando parámetros medidos independientemente, ver si el modelo es capaz de predecir razonablemente los datos, ya sea sus valores concretos, su tendencia, la existencia o no de fluctuaciones, etc. Si tenemos más de un modelo o más de una versión, podemos comparar sus ajustes a los datos para decidirnos por uno de ellos. Además, nos puede interesar estudiar la sensibilidad del modelo a los parámetros. Esto es: ¿cuánto varía el resultado del modelo cuando varío un parámetro dado en un cierto porcentaje? Un buen modelo debe variar poco con los parámetros que no nos interesan o que son muy difíciles de medir (se dice que el modelo es robusto a estos parámetros) y ser sensible a los parámetros que más nos interesan. Como esto no se sabe a priori, también podemos usar el estudio de sensibilidad para decidir qué parámetros hay que medir con más exactitud si queremos que el resultado sea confiable y tenga un margen de error pequeño. A modo de despedida Aquí llegamos casi al final del libro (en realidad sólo falta la bibliografía recomendada y algunos ejercicios de aplicación) y esperamos que haya sido de provecho. Hemos contado sintéticamente los principios fundamentales para trabajar en ecología matemática. Esto puede servir tanto para hacer modelos en forma autónoma, como para discutir y trabajar en grupo con colegas de 101 Fernando R. Momo - Ángel F. Capurro otras disciplinas, como para leer con más provecho los manuscritos de otros ecólogos. El libro puede usarse como texto para un curso básico de ecología matemática o de teoría ecológica y creemos que su incorporación a cursos de grado en la formación de los ecólogos podría ser bastante provechosa. El lector dirá si hemos logrado nuestros objetivos. Permítasenos decir que la ecología matemática nos parece divertida, agradable, elegante e interesante. Es más de lo que se puede decir de algunas personas. 102 Algunos ejercicios de ecología matemática 1. Mostrar que la siguiente solución de la ecuación logística: -r t N (t) = No K / (No + (K - No) e tiene las siguientes propiedades: a) N tiende a K cuando t tiende a infinito b) El gráfico es cóncavo hacia arriba para No < N < K/2 c) El gráfico es cóncavo hacia abajo para K/2 < N < K d) El gráfico es cóncavo hacia arriba para No > K. 2. Considere el modelo dN /dt = r N (K - N) (N - M) con r > 0, 0 < M < K. Exprese la tasa de crecimiento g(N) como una polinomial en N y encuentre los coeficientes ao, a1 y a2. Muestre que N = 0, N = M y N = K son estados de equilibrio y diga cómo es su estabilidad. 3. Los modelos comúnmente usados en pesquerías son del tipo dN/dt = N f(N) con f(N) dado por: -b N a) Modelo de Ricker f(N) = r e b) Modelo de Beverton-Holt f(N) = r / (a + N) Analice el comportamiento de ambos modelos asumiendo a, b, y r mayores que cero. 4. Para modelos de poblaciones únicas, decir cuál o cuáles de las siguientes tasas de crecimiento densodependiente conducirían a una desaceleración del crecimiento cuando la población se incrementa. ¿Cuáles conducirían a poblaciones estables? f(N) = a) β / (1 + N) con β > 0 b) β - β N con β > 0 αN c) N - e con α > 0 d) log N 103 Fernando R. Momo - Ángel F. Capurro 5. Smith (1963) observó que en un cultivo de algas unicelulares, el g(N) para poblaciones de Daphnia magna decrece en forma no lineal cuando N se incrementa. De acuerdo a este hecho sugirió que dicha tasa de crecimiento depende de la tasa de consumo de comida. La expresión que propuso es: g(N) = r (T - F) / T donde F es la tasa de utilización cuando la población es N, y T es la tasa máxima cuando la población ha alcanzado un nivel de saturación. Asumió además que: F = C1 N + C2 dN/dt con C1 y C2 mayores que cero y dN/dt > 0. a) Explique este supuesto para F. b) La ecuación logística modificada queda expresada por: dN/dt = r N [(K - N) / (K + q N)] donde q = r C2 / C1 y K = T / C1. Represente entonces la expresión entre corchetes en función de N. c) ¿Qué comportamiento tiene, cualitativamente hablando, el crecimiento poblacional? 6. a) Analizar los puntos de equilibrio y su estabilidad para el caso de una población con crecimiento logístico que es explotada a una velocidad constante h. b) Analizar el caso en que la velocidad de explotación es proporcional a la biomasa (E.B) y, además, la población presenta efecto Allee. 7. Formular el modelo de depredador-presa para el caso de una población que se alimenta de materia orgánica en descomposición (ésta última considerada como presa). La materia orgánica tiene un aporte externo a velocidad constante y una tasa de descomposición constante. 8. Un trabajo de Livdahl y Willey (1991, Science 253) estudia la competencia entre larvas de mosquito de Aedes triseriatus (nativa) y Aedes albopictus (invasora y transmisora del dengue) en Texas, 104 Ecología Matemática - Principios y aplicaciones en dos ambientes diferentes: a) árboles huecos, b) cubiertas de auto abandonadas. Asumiendo un modelo de competencia de Lotka-Volterra, los parámetros para ambas especies se muestran en la tabla: En troncos huecos En cubiertas abandonadas Aedes triseriatus K1 α12 100 2.38 38 1.12 Aedes albopictus K2 α21 74 1.42 184 4.29 a) Analice los puntos de equilibrio y su estabilidad para ambos hábitats. b) ¿Cuáles son su conclusiones respecto al control de la invasión de A. albopictus?   N  N t +1 = N t 1 + λ 1 − t  K    9. La siguiente ecuación representa el crecimiento de un población: a) Halle el o los puntos de equilibrio del modelo. b) Analice qué sucede con la estabilidad de los puntos diferentes de cero dependiendo del valor de λ. 10. Para la ecuación de crecimiento discreta N r  1− t  K  t +1 t probar que N = 0 y N = K son puntos estables. Hallar el valor de Nt para el cual Nt+1 es máximo. ¿Cuánto vale el máximo? Analizar la dinámica del modelo según el valor que toma r. N =Ne 11. Construir la matriz de transición del siguiente grafo que representa un ciclo de vida: 105 Fernando R. Momo - Ángel F. Capurro Hallar los autovalores de la matriz y, en base a ellos, predecir ´la dinámica de la población. Si existe, hallar la estructura estable de edades de la población. 12. Suponga una especie con dos estadíos (Larva y Adulto); la supervivencia de larva a adulto es de 0.405; la fecundidad del adulto es 2. a) Dibuje el grafo del ciclo de vida y escriba la matriz correspondiente. b) ¿Cuál es la dinámica de la población? c) ¿Cómo variaría esa dinámica si los adultos sobrevivieran de un año a otro con una probabildad de 0.2? 13. Representar con un modelo discreto matricial la dinámica de una población biparental cuyos dos sexos tienen la misma probabilidad de nacimiento pero diferentes mortalidades. ¿Es posible hallar una proporción de sexos en equilibrio? 14. ¿Cómo se comporta una población con crecimiento logístico cuando es explotada en proporción directa a su densidad? ¿Qué efecto tiene en ese comportamiento la introducción de los costos y beneficios económicos de la explotación considerando la existencia de costos fijos no nulos? 15. Para el siguiente modelo de una enfermedad infecciosa, plantee las ecuaciones diferenciales correspondientes y describa el comportamiento del sistema: 106 Ecología Matemática - Principios y aplicaciones Dado que existe pérdida de la inmunidad, ¿de qué dependerá la utilidad o no de la vacunación y qué efecto tendría sobre la dinámica de la enfermedad? 16. ¿Cómo plantearía un modelo SIR en tiempo discreto? ¿Cuál es el estado final de equilibrio? 17. Elabore un modelo que represente la concentración de un contaminante en un río en función del tiempo considerando: a) la entrada de contaminante está dada por un efluente de caudal y concentración constantes, b) la velocidad del río es constante, c) el contaminante es biodegradable en el agua, d) una vez absorbido por las algas el contaminante no se degrada y puede ser acumulado también por el zooplancton y los peces a través de su alimentación. 107 Bibliografía básica La que sigue no es una lista exhaustiva pero incluye libros que sin duda vale la pena consultar y considerar a la hora de ponerse a estudiar la ecología matemática con cierto rigor. Pedimos disculpas por las omisiones. 1. Abraira Santos, V. y Pérez de Vargas Luque, A., 1996. Métodos multivariantes en bioestadística. Vol. Editorial Centro de Estudios Ramón Areces S.A., Madrid, 452 pp. 2. Caswell, H. 1989. Matrix Population Models. Sinauer Associated Inc. 3. COMAP, 1999. Las matemáticas en la vida cotidiana. Vol. Addison-Wesley Iberoamericana España, Madrid, 722 y apéndices pp. 4. DeAngelis, D. L. y L. J. Gross (Eds.) 1992. IndividualBased Models and Approaches in Ecology. Populations, Communities and Ecosystems. Chapman & Hall. 5. Esteva, L. y Falconi, M. compiladores. 2002. Biología Matemática. Un enfoque desde los sistemas dinámicos. UNAM. 6. Gillman, M. y Hails, R., 1997. An introduction to ecological modelling: Putting practice into theory. Methods in ecology Vol. Blackwell Science, Victoria, 202 pp. 7. González Guzmán, J. 1999. Ecología Matemática. Tomo I. Modelos de tiempo discreto de poblaciones sin estructura. Editorial de la Facultad de Ciencias Básicas y Matemáticas de la Universidad Católica de Valparaíso. 163 pp. 8. González Guzmán, J. 2001. Ecología Matemática. Tomo II. Modelos de tiempo discreto con estructura etárea, genética y espacial. Editorial de la Facultad de Ciencias Básicas y Matemáticas de la Universidad Católica de Valparaíso. 189 pp. 9. González Manteiga, M.T. 2003. Modelos matemáticos discretos en las ciencias de la naturaleza. Teoría y problemas. Editorial Díaz de Santos. Madrid. 223 pp. 109 Fernando R. Momo - Ángel F. Capurro 10. Haken, H., 1981. Secreto de los éxitos de la naturaleza. Sinergética: la doctrina de la acción de conjunto. Vol. ArgosVergara, Barcelona, 242 pp. 11. Hall, C. (ed.) 1995. Maximum power: the ideas and applications of H. T. Odum. University Press of Colorado. 12. Hallam, T. G. y S. A. Levin (Eds.) 1986. Mathematical Ecology. An Introduction. Springer-Verlag. Biomathematics Vol. 17. 13. Hannon, B. y Ruth, M., 1997. Modeling dynamic biological systems. Modeling Dynamic Systems Vol. Springer, New York, 338 pp. 14. Jeffers, J., 1991. Modelos en ecología. Vol. Oikos-Tau, Barcelona, 96 pp. 15. Kahn, P. B. 1990. Mathematical Methods for Scientists and Engineers. Linear and Nonlinear Systems. John Wiley & Sons. 16. Levin, S. A.; T. G. Hallam y L. J. Gross (eds.) 1988. Applied Mathematical Ecology. Springer-Verlag. Biomathematics Vol. 18. 491pp. 17. Newman, E. 1993. Applied Ecology. Blackwell Sicence. 18. Platt, T.; K. H. Mann y R. E. Ulanowicz (Eds.) 1981. Mathematcial Models in Biological Oceanography. UNESCO Press. 150pp. 19. Ruth, M. y Hannon, B., 1997. Modeling dynamic economic systems. Modeling Dynamic Systems Vol. Springer, New York, 339 pp. 20. Sánchez Garduño, F. Miramontes, P. y Gutiérrez Sánchez, J.L. Coordinadores. 2002. Clásicos de la biología matemática. Siglo Veintiuno editores. 21. Valderrama Bonnet, M., 1995. Modelos matemáticos en las ciencias experimentales. Ciencia y Técnica Vol. Ediciones Pirámide, Madrid, 222 pp. 22. Yodzis, P. 1989. Introduction to Theoretical Ecology. Harper & Row. 384 pp. 110 Tabla de contenidos Introducción: ¿Por qué un libro de ecología matemática?....... 9 Capítulo 1: los modelos matemáticos y su razón de ser.......... 11 Capítulo 2: modelos de crecimiento poblacional...................... 17 El crecimiento dependiente de la densidad y el modelo logístico...................................................................... 19 El criterio de estabilidad................................................................ 21 Modelos generalizados.................................................................. 23 Algunas variaciones sobre el modelo logístico y su dinámica...... 25 Los modelos discretos de crecimiento........................................... 28 Capítulo 3: depredación y competencia..................................... 37 Cálculo del equilibrio y estudio de la dinámica............................ 38 La obtención de la fórmula de las órbitas...................................... 40 Estudio del equilibrio usando la matriz jacobiana......................... 42 Modificaciones al modelo básico de Volterra............................... 45 Respuestas funcionales y respuestas numéricas............................ 49 Los modelos de competencia......................................................... 50 Las funciones de Lyapunov........................................................... 55 Capítulo 4: algunas aplicaciones de los modelos continuos..... 59 El caso de la explotación proporcional a la biomasa..................... 62 Modelos epidemiológicos básicos................................................. 64 Dinámica de un ecosistema marino planctónico........................... 68 Capítulo 5: modelos matriciales de poblaciones....................... 75 El análisis básico de un modelo matricial..................................... 79 Más información disponible en la matriz...................................... 81 Capítulo 6: de la ecología teórica a la teoría ecológica............ 83 Las poblaciones humanas.............................................................. 83 La termodinámica de los sistemas ecológicos............................... 88 Capítulo 7: libro de cocina para hacer tu propio modelo matemático...................................................................... 93 Primer paso: delimitar claramente el problema ecológico............ 93 Segundo paso: identificar qué se puede modelizar y cuáles son los fenómenos relevantes.......................................... 94 Tercer paso: escribir las ecuaciones del modelo............................ 95 Cuarto paso: estudiar la existencia, condiciones y tipos de equilibrio....................................................................... 97 Quinto paso: hilando finito.......................................................... 100 A modo de despedida................................................................... 101 Algunos ejercicios de ecología matemática............................. 103 Bibliografía básica..................................................................... 109