2025 Workshop en Modelamiento de Sistemas Biológicos
Los sistemas multiestables pueden terminar en uno de los posibles estados, dependiendo de dónde comiencen, lo que se conoce como las condiciones iniciales del sistema. En circuitos biológicos, esto depende de las concentraciones iniciales de las especies químicas involucradas. Pero ¿cómo podemos saber a qué estado llegará el sistema cuando pase suficiente tiempo? Para responder a eso, necesitamos analizar la estabilidad de cada estado. Un estado se considera estable si el sistema tiende a regresar a él después de una pequeña perturbación. En cambio, un estado es inestable si incluso una perturbación leve hace que el sistema se aleje de él y se dirija hacia otro estado.
Una forma de identificar los estados a los que puede converger un sistema multiestable es a través de un análisis de nullclines. Las nullclines son las curvas donde la tasa de cambio de una especie química en particular es cero. Si nos enfocamos en dos especies del circuito, podemos representar sus dinámicas en un plano de fases, donde cada nullclina corresponde a una de las especies. Los puntos donde estas curvas se cruzan se llaman puntos fijos, y representan concentraciones en estado estacionario que definen los estados del sistema.
Para ilustrar estos conceptos, analicemos el toggle switch. De la Lecture Note #2, recordamos las ecuaciones que describen este circuito.
Para ilustrar estos conceptos, analicemos el toggle switch. De la Lecture Note #2, recordamos las ecuaciones que describen este circuito.
Aunque derivar las ecuaciones de las nullclines para este sistema es relativamente sencillo, escribamos un script que utilice sympy para obtener dichas ecuaciones.
import sympy as sym
# Definir los símbolos
x, y = sym.symbols('x y', real=True, positive=True)
alpha, K, m, delta = sym.symbols('alpha K m delta', real=True, positive=True)
# Definir las EDOs
dxdt = alpha * (K**m / (K**m + y**m)) - delta * x
dydt = alpha * (K**m / (K**m + x**m)) - delta * y
# Nullclines: igualar dx/dt = 0 y dy/dt = 0
x_nullcline = sym.Eq(dxdt, 0)
y_nullcline = sym.Eq(dydt, 0)
# Resolver para x y para y
x_nullcline_expr = sym.solve(x_nullcline, x)[0]
y_nullcline_expr = sym.solve(y_nullcline, y)[0]
# Display de nullclines como ecuaciones en formato LaTeX
x_latex = sym.latex(sym.Eq(x, x_nullcline_expr))
y_latex = sym.latex(sym.Eq(y, y_nullcline_expr))
print("Nullclina de x:", x_latex)
print("Nullclina de y:", y_latex)
Nullclina de x: x = \frac{K^{m} \alpha}{\delta \left(K^{m} + y^{m}\right)}
Nullclina de y: y = \frac{K^{m} \alpha}{\delta \left(K^{m} + x^{m}\right)}
import sympy as sym
# Define symbols
x, y = sym.symbols('x y', real=True, positive=True)
alpha, K, m, delta = sym.symbols('alpha K m delta', real=True, positive=True)
# Define the ODEs
dxdt = alpha * (K**m / (K**m + y**m)) - delta * x
dydt = alpha * (K**m / (K**m + x**m)) - delta * y
# Nullclines: set dx/dt = 0 and dy/dt = 0
x_nullcline = sym.Eq(dxdt, 0)
y_nullcline = sym.Eq(dydt, 0)
# Solve for x and y
x_nullcline_expr = sym.solve(x_nullcline, x)[0]
y_nullcline_expr = sym.solve(y_nullcline, y)[0]
# Export nullclines as LaTeX strings
x_latex = sym.latex(sym.Eq(x, x_nullcline_expr))
y_latex = sym.latex(sym.Eq(y, y_nullcline_expr))
print("x-nullcline:", x_latex)
print("y-nullcline:", y_latex)
x-nullcline: x = \frac{K^{m} \alpha}{\delta \left(K^{m} + y^{m}\right)}
y-nullcline: y = \frac{K^{m} \alpha}{\delta \left(K^{m} + x^{m}\right)}
Podemos mostrar estas ecuaciones en formato LaTeX para una mejor legibilidad.
Ahora, vamos a graficar las nullclinas en el plano de fases y a encontrar los puntos fijos utilizando la función intersection.
# Definimos las ecuaciones de nullclines
def null_X(Y_x, a, d, K, m):
X_x = (a/d) * K**m / (K**m + Y_x**m)
return X_x
def null_Y(X_y, a, d, K, m):
Y_y = (a/d) * K**m / (K**m + X_y**m)
return Y_y
# Define the nullclines equations
def null_X(Y_x, a, d, K, m):
X_x = (a/d) * K**m / (K**m + Y_x**m)
return X_x
def null_Y(X_y, a, d, K, m):
Y_y = (a/d) * K**m / (K**m + X_y**m)
return Y_y
import os, sys, subprocess
if "google.colab" in sys.modules:
cmd = "pip install --upgrade bokeh-catplot intersect"
process = subprocess.Popen(cmd.split(), stdout=subprocess.PIPE, stderr=subprocess.PIPE)
stdout, stderr = process.communicate()
# ------
# Librerias para el analisis matematico
import numpy as np
from intersect import intersection
# Librerias para visualizar nuestros resultados
import bokeh.io
import bokeh.plotting
import bokeh.palettes
from bokeh.models import LinearColorMapper, ColorBar
from bokeh.models import Range1d
from bokeh.io import export_svgs
bokeh.io.output_notebook()
import os, sys, subprocess
if "google.colab" in sys.modules:
cmd = "pip install --upgrade bokeh-catplot intersect"
process = subprocess.Popen(cmd.split(), stdout=subprocess.PIPE, stderr=subprocess.PIPE)
stdout, stderr = process.communicate()
# ------
# Libraries for mathematical analysis
import numpy as np
from intersect import intersection
# Libraries to visualize our results
import bokeh.io
import bokeh.plotting
import bokeh.palettes
from bokeh.models import LinearColorMapper, ColorBar
from bokeh.models import Range1d
from bokeh.io import export_svgs
bokeh.io.output_notebook()
# Definir la visualización con Bokeh
p = bokeh.plotting.figure(width = 325, height = 275,
x_range = Range1d(0, 3.1),
y_range = Range1d(0, 3.1),
x_axis_label = 'X',
y_axis_label = 'Y')
# Parámetros cinéticos
a = 2.2
K = 1
d = 1
m = 3
# Definimos las variables de estado
X_y = np.linspace(0, 3, 300)
Y_x = np.linspace(0, 3, 300)
# Calculamos las nullclines
X_x = null_X(X_y, a, d, K, m)
Y_y = null_Y(Y_x, a, d, K, m)
# Encontramos los puntos fijos mediante la función de intersección
fx_points = intersection(X_x, Y_x, X_y, Y_y)
X_i = fx_points[0] # coordenadas en el eje x
Y_i = fx_points[1] # coordenadas en el eje y
# Graficamos las nullclines y encontramos los puntos fijos
p.line(X_y, Y_y, line_width = 3, color = 'black')
p.line(X_x, Y_x, line_width = 3, color = 'black')
p.scatter(X_i, Y_i, size = 10, fill_color = 'black', color = 'black')
# Visualizando
bokeh.io.show(p)
# Define the bokeh visualization
p = bokeh.plotting.figure(width = 325, height = 275,
x_range = Range1d(0, 3.1),
y_range = Range1d(0, 3.1),
x_axis_label = 'X',
y_axis_label = 'Y')
# Kinetic parameters
a = 2.2
K = 1
d = 1
m = 3
# Definimos las variables de estados
X_y = np.linspace(0,3,300)
Y_x = np.linspace(0,3,300)
# Computing the nullclines
X_x = null_X(X_y, a, d, K, m)
Y_y = null_Y(Y_x, a, d, K, m)
# Finding the fixed points with the intersection
fx_points = intersection(X_x, Y_x, X_y, Y_y)
X_i = fx_points[0] # x-axis coordinates
Y_i = fx_points[1] # y-axis coordinates
# Plotting the nullclines and highlight the fixed points
p.line(X_y, Y_y, line_width = 3, color = 'black')
p.line(X_x, Y_x, line_width = 3, color = 'black')
p.scatter(X_i, Y_i, size = 10, fill_color = 'black', color = 'black')
# Visualization
bokeh.io.show(p)
La idea detrás del análisis de estabilidad local es simplificar el comportamiento de un sistema no lineal cerca de un punto fijo mediante su linearización. Las EDOs no lineales pueden ser difíciles de analizar directamente, por lo que usamos una expansión en series de Taylor alrededor de un punto fijo para aproximar el sistema con uno lineal.
Supongamos que tenemos una EDO que describe la dinámica de una especie química $Y$:
La expansión de Taylor de $f(y)$ alrededor de un punto fijo $y_0$ será:
Como $y_0$ es un punto fijo, sabemos por definición que la derivada en ese punto es cero. Por lo tanto,
Así, la expresión se simplifica a:
Esta aproximación lineal nos permite estudiar el comportamiento local del sistema cerca de $y_0$. Para sistemas con más de una variable, la derivada $\frac{df}{dy}\bigg|_{y = y_0}$ se convierte en una matriz Jacobiana, denotada como $\nabla f(y_0)$ y evaluada en el punto fijo. Una vez más, en lugar de calcular esta matriz a mano, vamos a computarla simbólicamente usando el paquete sympy.
def symbolic_jacobian(v_str, f_list):
varss = sym.symbols(v_str)
f = sym.sympify(f_list)
J = sym.zeros(len(f),len(varss))
for i, fi in enumerate(f):
for j, s in enumerate(varss):
J[i,j] = sym.diff(fi, s)
return J
v_str = 'y x'
f_list = [f'(a * K**m / (K**m + y**m) - d*x)',
f'(a * K**m / (K**m + x**m) - d*y)',
]
J = symbolic_jacobian(v_str, f_list)
J
def symbolic_jacobian(v_str, f_list):
varss = sym.symbols(v_str)
f = sym.sympify(f_list)
J = sym.zeros(len(f),len(varss))
for i, fi in enumerate(f):
for j, s in enumerate(varss):
J[i,j] = sym.diff(fi, s)
return J
v_str = 'y x'
f_list = [f'(a * K**m / (K**m + y**m) - d*x)',
f'(a * K**m / (K**m + x**m) - d*y)',
]
J = symbolic_jacobian(v_str, f_list)
J
Ahora, vamos a reescribir nuestra matriz Jacobiana como una función para que podamos calcular su valor numéricamente.
def jacobian(y0, x0, a, d, K, m):
# Calcular la matriz Jacobiana
jmat = np.zeros([2,2])
jmat[0,0] = - (K**m * a * m * y0**(m-1)) / (K**m + y0**m)**2
jmat[0,1] = -d
jmat[1,0] = - (K**m * a * m * x0**(m-1)) / (K**m + x0**m)**2
jmat[1,1] = -d
return jmat
def jacobian(y0, x0, a, d, K, m):
# Calculate the jacobian matrix
jmat = np.zeros([2,2])
jmat[0,0] = - (K**m * a * m * y0**(m-1)) / (K**m + y0**m)**2
jmat[0,1] = -d
jmat[1,0] = - (K**m * a * m * x0**(m-1)) / (K**m + x0**m)**2
jmat[1,1] = -d
return jmat
En este último paso de nuestro análisis de estabilidad local, determinaremos la estabilidad observando la parte real de los valores propios (o eigenvalores) de la matriz Jacobiana. Matemáticamente, esos valores propios son las raíces del polinomio característico. Para calcularlo, si $J$ es nuestra matriz Jacobiana, necesitamos calcular el determinante de $\lambda I - J$, donde $\lambda$ representa los valores propios, e $I$ es la matriz identidad de dimensión apropiada.
# Definir el símbolo del eigenvalores: λ
λ = sym.symbols('λ')
# Computando el polinomio caracteristico
characteristic_poly = sym.Matrix(J - λ * sym.eye(J.shape[0])).charpoly(λ)
characteristic_poly
# Define the symbolic variable for the λ
λ = sym.symbols('λ')
# Compute the characteristic polynomial using SymPy's charpoly function
characteristic_poly = sym.Matrix(J - λ * sym.eye(J.shape[0])).charpoly(λ)
characteristic_poly
Tener una expresión explícita para el polinomio característico podría ayudarnos a entender mejor la relación entre las constantes de reacción que permiten un comportamiento deseado. Pero dejemos eso para un trabajo futuro! Terminemos nuestra nota aprendiendo cómo calcular numéricamente los valores propios.
# Parámetros cinéticos
a = 2.2
K = 1
d = 1
m = 3
# Definimos las variables de estado
X_y = np.linspace(0, 3, 300)
Y_x = np.linspace(0, 3, 300)
# Paso N°1: calcular los puntos fijos
X_x = null_X(X_y, a, d, K, m)
Y_y = null_Y(Y_x, a, d, K, m)
fx_points = intersection(X_x, Y_x, X_y, Y_y)
X_i = fx_points[0] # coordenadas en el eje x
Y_i = fx_points[1] # coordenadas en el eje y
# Calcular la Jacobiana
jmat = jacobian(Y_i[0], X_i[0], a, d, K, m)
# Kinetic parameters
a = 2.2
K = 1
d = 1
m = 3
# Definimos las variables de estados
X_y = np.linspace(0,3,300)
Y_x = np.linspace(0,3,300)
# Step N°1: compute the fixed points
X_x = null_X(X_y, a, d, K, m)
Y_y = null_Y(Y_x, a, d, K, m)
fx_points = intersection(X_x, Y_x, X_y, Y_y)
X_i = fx_points[0] # x-axis coordinates
Y_i = fx_points[1] # y-axis coordinates
# Compute the Jacobian
jmat = jacobian(Y_i[0], X_i[0], a, d, K, m)
Para calcular los valores propios, usamos la función linalg.eigvals del paquete NumPy.
# Obteniendo los eigenvalores
eign = np.linalg.eigvals(jmat)
# Mostrando los componentes reales y complejos de los eigenvalores
print("Real part:")
print(np.real(eign))
print("\nImaginary part:")
print(np.imag(eign))
Real part:
[-1.15652215e-05 -1.24135953e+00]
Imaginary part:
[0. 0.]
# Get the eigenvalues
eign = np.linalg.eigvals(jmat)
# Displaying the real components of eigenvalue
print("Real part:")
print(np.real(eign))
print("\nImaginary part:")
print(np.imag(eign))
Real part:
[-1.15652215e-05 -1.24135953e+00]
Imaginary part:
[0. 0.]
Finalmente, necesitamos evaluar la parte real de nuestros valores propios.
• Si la parte real de todos los valores propios es negativa ($\mathcal{Re} (\lambda) < 0$), entonces el punto fijo donde evaluamos la Jacobiana es estable.
• Si la parte real de al menos un valor propio es positiva ($\mathcal{Re} (\lambda) > 0$), entonces el punto fijo donde evaluamos la Jacobiana es inestable.