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 = r1N11 − 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