Sistemas dinámicos: El sistema de Lotka-Volterra

Resumen

Un sistema es un conjunto de "cosas" relacionadas entre sí, un sistema dinámico es aquel en el cual sus variables evolucionan a través del tiempo. Los sistemas dinámicos se clasifican en lineales o no lineales; continuos o discretos; y variantes o invariantes en el tiempo; mientras que la interacción entre sus variables se describe mediante modelos matemáticos formulados con base en un conjunto determinista de reglas que establecen cómo será el siguiente estado futuro a partir del actual, donde el valor de cada variable en un instante de tiempo determina el estado del sistema. Para conocer el comportamiento de un sistema a partir de un punto inicial es necesario resolverlo o integrarlo mediante iteraciones en pequeños intervalos de tiempo, lo que permite obtener un conjunto de puntos conocidos como solución si se ilustran con respecto al tiempo o trayectoria si se ilustran en un plano n-dimensional (espacio de estados), esto se realiza utilizando software de programación para simular o emular la evolución del sistema dinámico al calcular los valores de sus variables mediante el desarrollo de algoritmos computacionales.

En el modelizado de sistemas biológicos, la ley de acción de masas se aplica para describir que el número de encuentros entre los integrantes de dos poblaciones es proporcional al producto de los tamaños de ambas poblaciones en cualquier instante de tiempo t; un ejemplo típico son las ecuaciones de Lotka-Volterra, descubiertas de forma independiente en la década de 1920 por Alfred Lotka y Vito Volterra. Estas ecuaciones describen la dinámica del denominado sistema presa-depredador, el cual supone que la tasa de depredación sobre la presa es proporcional a la tasa a la que se encuentran los depredadores y la presa al igual que la tasa de crecimiento de los depredadores, la anterior se formula mediante un término no lineal xy, donde x es el número de presas y y es el número de depredadores en un instante de tiempo t, es decir, xy=x(t)y(t). El modelo matemático se presenta a continuación:

dx/dt = alpha*x - beta*x*y; (1)
dy/dt = delta*x*y - gamma*y; (2)

La descripción de las ecuaciones del sistema de Lotka-Volterra (1)-(2) se realiza a continuación. La Ecuación (1) modeliza la dinámica de la población presa en la cual el término +alpha*x define el crecimiento de esta población de manera exponencial y el término -beta*x*y describe la tasa de depredación como una proporción entre los encuentros de la población presa y los depredadores. La Ecuación (2) modeliza la evolución de los depredadores, el término +delta*x*y implica que el crecimiento de la población depredadora depende del encuentro con las presas, sin embargo, la tasa de crecimiento no necesariamente es igual a la tasa en que consumen a las presas, por eso se emplean parámetros distintos, el término - gamma*y representa la muerte natural o emigración de los depredadores como un decrecimiento exponencial en la ausencia de presas.

Contexto

En el contexto de sistemas dinámicos que describen sistemas biológicos o fisiológicos, el modelizado in silico es una extensión lógica de la experimentación in vitro controlada, es el resultado natural del gran aumento de la potencia computacional disponible a un costo que disminuye continuamente, combinando las ventajas de la experimentación in vivo e in vitro, sin someterse a las consideraciones éticas y la falta de control asociadas con los experimentos in vivo. A diferencia de los experimentos in vitro, que existen de forma aislada, los modelos in silico permiten incluir un conjunto prácticamente ilimitado de variables y parámetros, lo que hace que los resultados sean más aplicables en problemas del mundo real.

La práctica forma parte del conjunto de prácticas de la asignatura "Biología de Sistemas", la cual es parte de la especialidad de Ingeniería Biomédica en el Instituto Tecnológico de Tijuana. Para esto los estudiantes necesitan comprender conceptos básicos relacionados con leyes de crecimiento, sistemas dinámicos no lineales, estabilidad de puntos equilibrio y tipos de conjuntos compactos invariantes que se pueden presentar en un sistema de Ecuaciones Diferenciales Ordinarias de primer orden.

Objetivos de la Actividad

La actividad tiene como objetivo principal introducir a estudiantes de licenciatura y posgrado al estudio de sistemas dinámicos que describen la evolución de sistemas biológicos, para esto se toma el modelo matemático tradicional de Lotka-Volterra (1) - (2) conformado por dos Ecuaciones Diferenciales Ordinarias de primer-orden, los objetivos de esta actividad son los siguientes:

1. Desarrollar un algoritmo computacional para solucionar el sistema con el método de Euler.
2. Desarrollar un algoritmo computacional para solucionar el sistema con el método de Heun.
3. Verificar la convergencia de las soluciones obtenidas reduciendo lo más posible el tamaño del paso de integración y comparando gráficamente los resultados obtenidos para cada nuevo conjunto de soluciones.
4. Construir el diagrama de bloques necesario para solucionar el sistema con diferentes funciones ode (ode15s, ode23s, ode23t, ode23tb, ode45, ode23 y ode113).
5. Ilustrar las soluciones y trayectoria del sistema de Lotka-Volterra.
6. Desarrollar un algoritmo computacional para calcular los puntos de equilibrio, la matriz jacobiana y los valores propios para analizar la estabilidad asintótica local del sistema.
7. Ilustrar el campo direccional del sistema.

Para los parámetros utilizar los valores alpha = 0.5, beta = 0.01, delta = 0.005 y gamma = 0.2, para las condiciones iniciales x(0) = 75 y y(0) = 75 en un intervalo del tiempo de t∈[0,100].

Materiales de Enseñanza

Se utiliza MATLAB para desarrollar algoritmos que permitan solucionar sistemas de Ecuaciones Diferenciales Ordinarias de primer orden que describen la evolución de sistemas biológicos. Se ilustran resultados correspondientes a las trayectorias en el plano de fase y las soluciones en el tiempo, se ilustra también el concepto de estabilidad de los puntos de equilibrio del sistema, adicionalmente, se guardan las figuras en formato PDF para su posterior visualización.

Archivo con la descripción completa de la práctica: Manual (Acrobat (PDF) 857kB Sep11 23)

Archivo Live Script con el código principal: lotka_volterra.mlx (MATLAB Live Script 700kB Dec1 23)

Archivo de Simulink con el modelo construido mediante diagramas de bloques: lotka_volterra_sim.slx ( 35kB Sep11 23) or lotka_volterra_sim.mdl ( 69kB Sep11 23)

Método de Euler

El método de Euler es una herramienta numérica que se utiliza para aproximar la solución phi(t,x₀) de un sistema dinámico no lineal definido por la ecuación de estado dx/dt = f(x), con una condición inicial x(0) = x₀. Al aplicar el método de Euler se calculan valores aproximados de un estado futuro de la solución con base en un valor inicial en intervalos equidistantes en el tiempo. Entonces, sea

dx/dt = f(x),

un sistema dinámico descrito por EDOs de primer orden, con el método de Euler se aproxima su solución al considerar que

dx/dt = f(x) ≈ Δx/Δt,

y si Δt define intervalos equidistantes, entonces se tiene lo siguiente

Δx/Δt = (x_i+1-x_i)/Δt ≈ f(x),

donde el primer valor de la solución es la condición inicial, por lo tanto, en la primera iteración se tiene que x_i=x₀, entonces el valor futuro de la solución que se desea aproximar está dado por

x_i+1 ≈ x_i + f(x_i)Δt, (3)

y esta ecuación se conoce como el método de Euler, la cual es una ecuación de diferencias finitas que se resuelve de manera iterativa en intervalos definidos por Δt en un dominio del tiempo dado por t∈[t₀,t_f], dondet₀=0 y t_f es el valor final en el tiempo en que se desea aproximar la solución phi(t,x₀) del sistema dinámico dx/dt = f(x). Cabe destacar que entre más pequeño sea el valor de Δt, el error resultante obtenido al resolver la Ecuación (3) será menor, sin embargo, esto implica un poder de cómputo mayor que no todos los ordenadores pueden soportar, por lo tanto, el valor de Δt se debe escoger con base en la información del sistema dinámico cuya solución se desea aproximar.

Método de Euler mejorado (Método de Heun)

El método de Euler mejorado, también conocido como el método de Heun, ofrece una mayor precisión al utilizar un promedio ponderado de las pendientes en dos puntos adyacentes, esta técnica se ha utilizado ampliamente en la resolución numérica de EDOs y ha demostrado ser una herramienta valiosa en la simulación y en el modelizado en una variedad de campos y disciplinas científicas, reduciendo significativamente el costo computacional en la aproximación de soluciones de sistemas dinámicos complejos. El proceso necesario para aplicar el método de Heun y obtener una mejor aproximación de x_i+1 consta de tres etapas descritas a continuación:

1. Primero, se aplica el método de Euler tradicional, sin embargo, este paso es solo una estimación provisional de x_i+1, denominada x_n, por lo tanto,

x_n = x_i + f(x_i)Δt,

donde x_n ∈Rⁿ, en este procedimiento es necesario recordar el valor obtenido en la primera pendiente de la estimación inicial, es decir, f(x_i).

2. Segundo, se toma el resultado obtenido para x_n  y se calcula la segunda pendiente definida por f(x_n ), esto implica que se debe evaluar el valor x_n  en el campo vectorial del sistema dinámico no lineal definido por la ecuación de estado dx/dt = f(x), la cual se utilizará para mejorar la aproximación inicial proporcionada por el método de Euler tradicional.
3. Finalmente, se utilizan los valores de calculados para ambas pendientes, f(x_i) y f(x_n), y se formula la siguiente ecuación para obtener una mejor estimación del valor futuro de la solución phi(t,x₀), es decir, x_i+1 :

x_i+1 = x_i + Δt(f(x_i) + f(x_n))/2. (4)

El método de Heun es una mejora con respecto al método de Euler tradicional en términos de precisión, la cual aumenta cuadráticamente en relación con la reducción del tamaño del paso Δt, sin embargo, sigue siendo un método de primer orden en términos de convergencia; por lo tanto, a medida que aumenta la complejidad de las ecuaciones o se busca una mayor precisión, existen otros métodos numéricos más avanzados, como el método de Runge-Kutta de cuarto orden que puede ser más apropiado.

Notas para los Educadores usando la Actividad

Es importante que los estudiantes tengan un conocimiento básico de sistemas dinámicos, ecuaciones diferenciales, métodos numéricos y programación (construcción de funciones, arreglos, ciclos y diagramas de bloques).

Evaluación

Para esta actividad, es indispensable que todos los alumnos obtengan los resultados indicados en cada objetivo, debido a que es la base para continuar con sistemas más complejos en el curso.

Bibliografía

[1] A. Garfinkel, J. Shevtsov, and Y. Guo, Modeling life: the mathematics of biological systems. New York, New York, USA: Springer, 2017, doi: https://doi.org/10.1007/978-3-319-59731-7.
[2] K. Bryan. (2021, March) SIMIODE EXPO 2021 - (p-r3) - Differential Equations: A Toolbox for Modeling the World. https://www.simiode.org/resources/8307.
[3] H. K. Khalil, Nonlinear systems, 3rd ed. Hoboken, New Jersey, USA: Prentice-Hall, 2002.
[4] H. Motulsky, Intuitive biostatistics: A nonmathematical guide to statistical thinking. Oxford, New York, USA: Oxford University Press, 2018.
[5] P. A. Valle, L. N. Coria, C. Plata & Y. Salazar, "CAR-T cell therapy for the treatment of ALL: eradication conditions and in silico experimentation," Hemato, vol. 2, issue 3, pp. 441-462, Jul 2021. https://doi.org/10.3390/hemato2030028

Esta actividad fue creada como parte del Taller con MATLAB Septiembre 2023.