Sistema Presa-Depredador con coordenada espacial

Introducción

Es común que en Sistemas Dinámicos se estudien las ecuaciones de Lotka-Volterra, las cuales sirven para modelar cómo una presa y un depredador interactúan en un mismo ecosistema para así poder predecir su dinámica poblacional. Dicho sistema consta de dos ecuaciones diferenciales ordinarias de primer orden, acopladas, autónomas y no lineales: \begin{aligned}
\frac{d N}{d t} &= N (a-bP) \\
\frac{d P}{d t} &= P (cN-d)\end{aligned} donde $ a,b,c,d $ son constantes positivas y $ N, P $ denotan a la población de la presa y el depredador respectivamente. De igual forma, es importante resaltar que $ N $ y $ P $ son funciones que sólo dependen del tiempo. Si bien este modelo es útil para describir lo que pasa en la naturaleza, no se debe olvidar que las presas y los depredadores no sólo se mueven en el tiempo, sino también en el espacio. De aquí nace la motivación de este artículo, pues al agregar una coordenada espacial y reajustar ciertos términos, obtenemos un sistema de reacción-difusión con ecuaciones diferenciales parciales donde ahora interesa encontrar la formación de patrones espaciales. Esto sucede siempre que las condiciones de patrones de Turing se cumplan.¹

Se comenzará dando una breve introducción sobre qué son las condiciones de Turing y de dónde nacen. Después se planteará y explicará el nuevo sistema de ecuaciones que incluye a la coordenada espacial, para realizar su análisis y llegar a las condiciones que se deben cumplir para obtener dicha formación de patrones.

Patrones de Turing

La teoría planteada por Charles Darwin $(1809 – 1882)$ dice que el jaguar tiene un pelaje manchado pues la selección natural así lo prefirió para su supervivencia; sin embargo, esto no explica el proceso mediante el cual se forman dichos patrones. El matemático James D. Murray $ (1931) $ propuso que estos patrones se pueden explicar a partir de las ideas de Turing.

El matemático Alan Turing $(1912 – 1954)$ notó que si bien los procesos de reacciones químicas y difusión molecular, por separado, son procesos homogeneizadores, cuando se les combina de cierta forma pueden llegar a formar patrones estacionarios, ¹A. 2012. “Patrones de Turing En Sistemas Biológicos.” Universidad Autónoma Metropolitana. http://mat.izt.uam.mx/mcmai/documentos/tesis/Gen.10-O/Ledesma-DA-Tesis.pdf.. Esto lo llevó a plantear un modelo de reacción-difusión en un espacio unidimensional. Para entender mejor esta idea, se tomará el siguiente ejemplo: supóngase que se tiene un bosque en el cual comienza a haber pequeños incendios en distintas regiones, y los bomberos buscan acabar con ello. En esta situación se considera al fuego como activador y a los bomberos como inhibidor del proceso de reacción-difusión. Si el fuego se difunde más lento que los bomberos, se podrá llegar a la formación de patrones; de lo contrario se tendría un bosque completamente quemado y homogéneo. Así, para poder tener patrones, es necesario que la sustancia activadora se difunda más lento que la inhibidora.

Turing descubrió que las ecuaciones que describen al modelo reacción-difusión tienen soluciones en las que ambas sustancias están homogéneamente distribuidas en el espacio cuando no consideramos la difusión; sin embargo, esas soluciones se vuelven inestables cuando se incluye a esta última. Entonces, si aparecen desviaciones en la homogeneidad, estas irán aumentando para dar paso a la formación de patrones,²Sánchez, A. 2012. “Turing Y Sus Patrones: El Pionero de La Biología Matemática.” 2012. https://nadaesgratis.es/anxo-sanchez/turing-y-sus-patrones-el-pionero-de-la-biologia-matematica.. También demostró matemáticamente la existencia de dichos patrones, aunque bajo fuertes suposiciones, como son un dominio unidimensional, una reacción química sencilla y el uso de ecuaciones diferenciales de orden lineal.

Años después, Murray extendió el estudio de Turing a dos dimensiones y calculó bajo qué condiciones se desestabiliza la solución uniforme. Dedujo que un sistema de reacción-difusión tiene inestabilidades creadas por la difusión si el estado homogéneo se mantiene estable incluso al realizar pequeñas perturbaciones bajo la ausencia de la difusión, pero se vuelve inestable con pequeñas perturbaciones espaciales donde la difusión sí está presente. Al linealizar el modelo para estudiar las estabilidades, se llega a un eigenproblema que sólo depende del espacio, en el cual se buscará que su parte real sea mayor a cero para poder obtener la formación de patrones.

Los sistemas de reacción-difusión se pueden reescalar y desdimensionar para obtener su forma general: \begin{aligned}
u_{t} &= \gamma f(u,v) + \nabla^2u \\
v_{t} &= \gamma g(u,v)+ d \, \nabla^2 v.\end{aligned}

Donde $ d $ es la proporción de la difusión empleada en el sistema y $ \gamma $ representa la fuerza de los términos reactivos. Las condiciones de Turing para este sistema son ³Murray, J. D. 2011. Mathematical Biology: I. An Introduction. Interdisciplinary Applied Mathematics. Springer New York.: $$\begin{array}{cc}
f_{u}+g_{v}<0, & f_{u}g_{v}-f_{v}g_{u}>0 , \\
df_{u}+g_{v}>0 ,& (df_{u}+g_{v})^2-4d(f_{u}g_{v}-f_{v}g_{u})>0.
\end{array} \quad (1)$$Como el modelo a tratar en este artículo consta de ecuaciones integro-diferenciales, resulta muy difícil aplicar las condiciones de Turing de manera directa. Dado esto, para buscar la formación de patrones, se hará el mismo análisis que Murray hizo en su libro para obtener las mismas. Dicho análisis consta en hacerle pequeñas perturbaciones espaciales al sistema. Para esto es importante tener claro que un estado del sistema es estable cuando al perturbar levemente una solución estacionaria el sistema regresa a su estado original; mientras que si la perturbación saca al sistema del estado estacionario en el que se encontraba, es inestable ⁴A. 2012. “Patrones de Turing En Sistemas Biológicos.” Universidad Autónoma Metropolitana. http://mat.izt.uam.mx/mcmai/documentos/tesis/Gen.10-O/Ledesma-DA-Tesis.pdf..

Planteamiento del modelo

Como ya se ha mencionado previamente, lo que se busca es analizar el comportamiento de un sistema Lotka-Volterra al considerar no sólo el factor temporal, sino también el espacial. Para esto tenemos el siguiente sistema ⁵Núñez López, M., E. Brigatti, and M. Oliva. n.d. “Analysis of a Spatial Lotka-Volterra Model with a Finite Range Predator-Prey Interaction.” European Physical Journal B. :

\begin{array}{cc}\dfrac{\partial N(x, t)}{\partial t} = D_{N} \dfrac{\partial ^2N(x,t)}{\partial x^{2}} + rN(x,t) – \alpha N(x,t) \int_{x-L_{1}}^{x + L_{1}} P(s,t) ds,\\ \dfrac{\partial P(x,t)}{\partial t} = D_{P}\dfrac{\partial^{2}P(x,t)}{\partial x^{2}} – mP(x,t) + \beta P(x,t) \int_{x-L_{2}}^{x+L_{2}}N(s,t) ds.&(2)\end{array}

Las ecuaciones anteriores representan el comportamiento de los depredadores y de las presas. El primero consume al segundo a un ritmo $ \alpha $ mientras que se reproduce a una tasa $ \beta, $ $ r $ es el ritmo de crecimiento de la presa, y $ m $ es la tasa espontánea a la que mueren los depredadores. Por otro lado, $ D_{N} $ y $ D_{P} $ son los coeficientes de difusión de las presas y los depredadores, respectivamente, ⁶Núñez López, M., E. Brigatti, and M. Oliva. n.d. “Analysis of a Spatial Lotka-Volterra Model with a Finite Range Predator-Prey Interaction.” European Physical Journal B. .
Este sistema tiene dos soluciones estacionarias y espacialmente homogéneas que se obtienen de la siguiente forma: \begin{cases}
\displaystyle 0 = rN – \alpha N \int_{x-L_{1}}^{x + L_{1}} P ds\\
\displaystyle 0 = -mP + \beta P \int_{x-L{2}}^{x+L_{2}} N ds\\
\end{cases}

\begin{cases}
\displaystyle N\left(r – \alpha \int_{x-L_{1}}^{x + L_{1}} P ds\right) = 0 \\
\displaystyle P\left(\beta \int_{x-L_{2}}^{x + L_{2}} N ds -m\right) = 0
\end{cases} La primera solución es la trivial, cuando $ N(x,t) = 0 $ y $ P(x,t) = 0. $ La segunda se obtiene como sigue: \begin{cases}
\displaystyle \int_{x-L_{1}}^{x + L_{1}} P ds = \dfrac{r}{\alpha}\\
\displaystyle \int_{x-L_{2}}^{x + L_{2}} N ds = \dfrac{m}{\beta}.
\end{cases} Al resolver la integral, se llega al resultado \begin{cases}
\displaystyle P(x + L_{1} -x + L_{1}) = \dfrac{r}{\alpha}\\
\displaystyle N(x + L_{2} -x + L_{2}) = \dfrac{m}{\beta},\\
\end{cases} por lo que la segunda solución está dada por: \begin{aligned}
\implies \overline{P}(x,t) &\equiv \dfrac{r}{2\alpha L_{1}}\\
\overline{N}(x,t) &\equiv \dfrac{m}{2\beta L_{2}}.\\\end{aligned} Ahora, con el fin de conocer el comportamiento del sistema y la formación de patrones, se hará una perturbación armónica a cada una de las soluciones con el fin de hacer un análisis de estabilidad. Entonces se tomarán a consideración las soluciones perturbadas \begin{aligned}
N(x,t) &= \overline{N} + A_{N} \exp(\lambda t + ikx), \\
(x,t) &= \overline{P} + A_{P} \exp(\lambda t + ikx). \end{aligned}

Sustituyéndolas en el sistema original (2), se llega a lo siguiente: \begin{aligned}
A_{N}(\lambda + k^{2}D_{N}) + \dfrac{\alpha m}{\beta L_{2}}A_{P}\dfrac{\sin(kL_{1})}{k} &= 0, \\
A_{P}(\lambda + k^{2}D_{P}) – \dfrac{\beta r A_{N}}{\alpha L_{1}}\dfrac{\sin(kL_{2})}{k} &= 0.\end{aligned}

Un sistema lineal que se puede expresar de la siguiente forma:
$$
\begin{bmatrix}
\lambda + k^{2} D_{N} &\dfrac{\alpha m}{\beta}\dfrac{\sin(kL_{1})}{kL_{2}} \\
\dfrac{-\beta r}{\alpha}\dfrac{\sin(kL_{2})}{kL_{1}} & \lambda + k^{2}D_{P} \\
\end{bmatrix}
\begin{bmatrix}
A_{N} \\
A_{P} \\
\end{bmatrix}=\begin{bmatrix}
0 \\
0 \\
\end{bmatrix}
$$
,y que además tiene una solución distinta a la trivial cuando el determinante de la matriz es distinto de cero. Por lo tanto, se busca que $$ \det(A) = \lambda^{2} + (D_{N} + D_{P})k^{2}\lambda + D_{N}D_{P}k^{4} + rm \dfrac{\sin(kL_{1})\sin(kL_{2})}{k^{2}L_{1}L_{2}} = 0. $$ Resolviendo la ecuación cuadrática para $ \lambda $ usando la fórmula general se llega a $$
\lambda(k) = \dfrac{-k^{2}}{2}(D_{N}+D_{P}) \pm \sqrt{\dfrac{k^{4}}{4}(D_{N}-D_{P})^{2} – rm \dfrac{\sin(kL_{1})\sin(kL_{2})}{k^{2}L_{1}L_{2}}}, \quad (3)$$ la cual es una función de $ k $.

Por las condiciones de Turing se sabe que los patrones espaciales ocurren si $ \Re(\lambda(k)) > 0. $ Para que esto suceda se requiere que $ L_{1}\neq L_{2}. $ Esto muestra que la inestabilidad está dada por el rango de interacción y es independiente del proceso de difusión ⁷Núñez López, M., E. Brigatti, and M. Oliva. n.d. “Analysis of a Spatial Lotka-Volterra Model with a Finite Range Predator-Prey Interaction.” European Physical Journal B..

Para simplificar el análisis, se tomará el caso $ D_{N} = D_{P} = D, L_{2} = 2L_{1} $ y se introducirán las variables reescaladas $ \widehat{k} = kL_{1} $ y $ \widehat{\lambda} = \frac{\lambda L_{1}^{2}}{D}, $ por lo que la ecuación (3) se reduce a: \begin{aligned}
\widehat{\lambda}(\widehat{k}) &= -\widehat{k}^{2} + \frac{\sqrt{rm}L_{1}^{2}}{Dk}\sqrt{-\sin^{2}(\widehat{k})\cos({\widehat{k}})}\\
\iff \widehat{\lambda}(\widehat{k}) &= -\widehat{k}^{2} + \dfrac{\mu}{k}\sqrt{-\sin^{2}(\widehat{k})\cos({\widehat{k}})},\end{aligned} donde $ \mu = \dfrac{\sqrt{rm}L_{1}^{2}}{D}. $

Resultados

Para encontrar las restricciones de los valores de los parámetros que garanticen la formación de patrones se quiere encontrar el valor de $ \mu $ tal que el máximo de la función $ \widehat{\lambda}(\widehat{k}) $ sea cero. Con el fin de obtener dichos valores se creó un código en MATLAB, el cual dio los resultados $ \widehat{k}_{m} = 1.83 $ y $ \mu = 12.52 $. Por lo que las variables originales deben de cumplir las siguientes condiciones para que $ \Re(\lambda(k)) > 0: $ $$ k_{m} \approx \dfrac{1.83}{L_{1}} \; , \qquad \dfrac{\sqrt{rm}L_{1}^{2}}{D} \gtrsim 12.5232. $$

En la figura 1 se puede apreciar la función $ \widehat{\lambda}(\widehat{k}) $ para distintos valores del parámetro $ \mu $.

Simulación numérica

Se puede notar que para el caso $ \mu = 13.52 $, la formación de patrones se da cuando $ \widehat{k} $ toma valores entre $ 1.71 $ y $ 1.99 $.

Conclusiones

En este artículo se buscaron las condiciones que deben cumplir los distintos parámetros en un sistema Lotka-Volterra espacial para tener una configuración no-homogénea. Las condiciones de Turing nos dicen que las inestabilidades en un sistema de reacción-difusión típicamente son creadas por los términos difusivos, contrario a lo que uno podría pensar. Sin embargo, en el sistema que se analizó, se obtuvo que las inestabilidades del sistema son independientes del proceso de difusión y, más bien, dependen directamente del rango de interacción entre la presa y el depredador. Una explicación física y biológica de lo anterior se da gracias a los procesos de muerte y nacimiento. En este estudio se vio que para poder tener formación de patrones se debe cumplir que $ L_{1}\neq L_{2} $, donde $ L_{1} $ está relacionado con la distancia necesaria para la muerte de la presa $($en presencia de un depredador$)$ y $ L_{2} $ se relaciona con la distancia requerida para la reproducción del depredador $($ en presencia de la presa$)$. Como la muerte puede ocurrir en cualquier lugar, mientras que el nacimiento de un ser vivo se da de manera adyacente a otros seres vivos, obtenemos una correlación espacial, la cual deriva en la formación de patrones bajo ciertas condiciones reproductivas. El sistema presentado es un caso muy peculiar, pues justamente no requiere del término difusivo para poder llegar a la formación de patrones. Esto se debe a que se trabajó con ecuaciones integro-diferenciales que afectan de manera muy especial la forma en la que la presa y el depredador interactúan.