Solución Numérica de EDOs

Modelando la dinámica de reacciones bioquímicas

A grandes rasgos, la investigación en biología inicia con la recolección de datos experimentales, para luego utilizar herramientas matemáticas para describir cómo se relacionan distintas variables — por ejemplo, mediante una regresión lineal. Lograr que un modelo se ajuste bien a los datos es un gran paso, ya que puede revelar aspectos clave de los mecanismos regulatorios que determinan el comportamiento del sistema. Hay mucho trabajo interesante en esta área, y el campo de la identificación de sistemas en particular ofrece herramientas útiles para descubrir estas relaciones.

En biología sintética también usamos modelos matemáticos, pero muchas veces con otro objetivo. En lugar de usarlos solo para describir relaciones entre variables, los usamos para diseñar circuitos biológicos de manera dirigida (forward engineering). Es decir, los modelos nos ayudan a crear sistemas que se comporten de una forma específica. En este contexto, los modelos predictivos resultan especialmente útiles, ya que nos permiten simular el comportamiento de un circuito (in silico) y reducir la cantidad de iteraciones experimentales necesarias para hacerlo funcionar. No esperamos que nuestros modelos representen todos los detalles de la biología. Más bien, buscamos que capturen las características clave que son relevantes para el problema de diseño que estamos abordando.

Ecuaciones cinéticas (o de velocidad de reacciones químicas)

Un marco matemático muy útil para modelar sistemas biológicos es el de las ecuaciones cinéticas, las cuales describen cómo cambian las concentraciones de las especies químicas a lo largo del tiempo. Este enfoque se basa en conceptos de la mecánica estadística, donde las interacciones moleculares se tratan como eventos estocásticos (aleatorios) que obedecen las leyes de la termodinámica. Para una explicación más detallada, recomendamos la sección introductoria del libro de Del Vecchio y Murray.

Dentro de este marco, usamos Ecuaciones Diferenciales Ordinarias (EDOs) para modelar un conjunto de reacciones químicas que representan las interacciones entre los componentes de un determinado circuito biológico. Como ejemplo, vamos a modelar la dinámica de un gen reportero. Supongamos que tenemos un gen regulado por un promotor constitutivo, que puede transcribirse para generar la molécula de ARNm llamada $m_{\text{gfp}}$, y que puede traducirse en la proteína llamada ${\text{gfp}}$.

Escribiendo las reacciones químicas

Los eventos bioquímicos involucrados en este sistema son la transcripción y la traducción. Vamos a asumir que la transcripción ocurre a una tasa $\theta$, y la traducción a una tasa $\rho$. Además, tanto el ARNm como la proteína están sujetos a degradación: el ARNm se degrada a una tasa $\delta$, y la proteína a una tasa $\phi$. Con estas suposiciones, podemos escribir las reacciones químicas correspondientes de la siguiente manera:

\begin{align} \text{DNA} &\xrightarrow{\theta} \text{DNA} + m_{\text{gfp}} && m_{\text{gfp}} \xrightarrow{\delta} \emptyset && \text{Transcripcion y degradación del mRNA}\\ m_{\text{gfp}} &\xrightarrow{\rho} m_{\text{gfp}} + \text{gfp} && \text{gfp} \xrightarrow{\phi} \emptyset && \text{Traducción y degradación de la proteina} \end{align}

Escribiendo las EDOs

Utilizando la ley de acción de masas, podemos modelar las reacciones químicas con las siguientes EDOs:

\begin{eqnarray} \frac{d}{dt}m_{\text{gfp}} &=& \theta - \delta m_{\text{gfp}} \\ \frac{d}{dt} \text{gfp} &=& \rho m_{\text{gfp}} - \phi \text{gfp} \end{eqnarray}

Resolviendo las EDOs: métodos simbólicos y numéricos

En este sistema, $m_{\text{gfp}}$ y $\text{gfp}$ son las variables de estado, lo que significa que representan cantidades internas cuyos valores cambian con el tiempo de acuerdo con la dinámica del sistema. El sistema es lineal porque cada ecuación es una combinación lineal de estas variables: aparecen solo multiplicadas por constantes, no están elevadas a potencias distintas de uno, y no hay productos entre variables ni funciones no lineales. Entonces, como el sistema es lo suficientemente simple, podemos usar una librería de Python llamada sympy para encontrar las soluciones de estas EDOs.

import sympy as sp

# Paso #1: definir los símbolos
t = sp.symbols('t')
theta, delta, rho, phi = sp.symbols('theta delta rho phi', positive=True, real=True)
m = sp.Function('m_gfp')
gfp = sp.Function('gfp')

# Paso #2: escribir la EDO que describe la dinámica del ARNm
ode1 = sp.Eq(m(t).diff(t), theta - delta * m(t))
sol_m = sp.dsolve(ode1, m(t), ics={m(0): 0})
m_sol = sol_m.rhs

# Paso #3: escribir la EDO que describe la dinámica de la proteína
ode2 = sp.Eq(gfp(t).diff(t), rho * m_sol - phi * gfp(t))
sol_gfp = sp.dsolve(ode2, gfp(t), ics={gfp(0): 0})
gfp_sol = sol_gfp.rhs

# Imprimir resultados
print("Solución de m_gfp(t):")
print(sp.latex(m_sol))
print("\nSolución de gfp(t):")
print(sp.latex(gfp_sol))

  Solución de m_gfp(t):
  \frac{\theta}{\delta} - \frac{\theta e^{- \delta t}}{\delta}

  Solución de gfp(t):
  - \frac{\rho \theta e^{- \phi t}}{\delta \phi - \phi^{2}} + \frac{\rho \theta e^{- \delta t}}{\delta \left(\delta - \phi\right)} + \frac{\rho \theta}{\delta \phi}

import sympy as sp

# Step #1: defining the symbols
t = sp.symbols('t')
theta, delta, rho, phi = sp.symbols('theta delta rho phi', positive=True, real=True)
m = sp.Function('m_gfp')
gfp = sp.Function('gfp')

# Step #2: writing the ODE describing mRNA dynamics, considering mRNA(0) = 0 as the initial condition
ode1 = sp.Eq(m(t).diff(t), theta - delta * m(t))
sol_m = sp.dsolve(ode1, m(t), ics={m(0): 0})  # Initial condition: mRNA(0) = 0
m_sol = sol_m.rhs

# Now, we use the previous solution for the ODE describing protein dynamics, for the same initial condition
ode2 = sp.Eq(gfp(t).diff(t), rho * m_sol - phi * gfp(t))
sol_gfp = sp.dsolve(ode2, gfp(t), ics={gfp(0): 0})  # Initial condition: gfp(0) = 0
gfp_sol = sol_gfp.rhs

# Print results in Latex
print("m_gfp(t) solution in LaTeX:")
print(sp.latex(m_sol))

print("\ngfp(t) solution in LaTeX:")
print(sp.latex(gfp_sol))

  m_gfp(t) solution in LaTeX:
  \frac{\theta}{\delta} - \frac{\theta e^{- \delta t}}{\delta}

  gfp(t) solution in LaTeX:
  - \frac{\rho \theta e^{- \phi t}}{\delta \phi - \phi^{2}} + \frac{\rho \theta e^{- \delta t}}{\delta \left(\delta - \phi\right)} + \frac{\rho \theta}{\delta \phi}

Utilicemos LaTeX para mostrar las soluciones de forma más entendible:

\begin{eqnarray} m_{\text{gfp}}(t) &=& \frac{\theta}{\delta} - \frac{\theta e^{- \delta t}}{\delta} \\ \text{gfp}(t) &=& - \frac{\rho \theta e^{- \phi t}}{\delta \phi - \phi^{2}} + \frac{\rho \theta e^{- \delta t}}{\delta \left(\delta - \phi\right)} + \frac{\rho \theta}{\delta \phi} \end{eqnarray}

Ahora vamos a visualizar la dinámica de ambas variables de estado. A lo largo de este workshop, usaremos la librería bokeh en lugar de matplotlib (u otras similares), simplemente porque permite generar visualizaciones más interactivas.

import numpy as np
import bokeh.io
import bokeh.plotting
import bokeh.palettes
from bokeh.models import Range1d

bokeh.io.output_notebook()

# Parámetros cinéticos
th = 1
rho = 2
d = 0.1
p = 0.05

# Tiempo de simulación
t = np.linspace(0,80,801)

# Escribir la solución de las EDOs
mgfp = (th/d) - (th * np.exp(-d*t))/d
gfp = -(rho*th*np.exp(-p*t))/(d*p - p**2) + (rho*th*np.exp(-d*t))/(d*(d - p)) + rho*th/(d*p)

# Configurar la visualización con bokeh
plot = bokeh.plotting.figure(width = 320, height = 225,
                             x_axis_label = 'Tiempo (minutos)', y_axis_label = 'Nivel de expresión')

# Graficar las soluciones
plot.line(t, mgfp, line_width = 3, color = 'grey', legend_label = 'mRNA')
plot.line(t, gfp, line_width = 3, color = 'orange', legend_label = 'proteína')

# Visualización
plot.legend.location = 'top_left'
plot.legend.background_fill_alpha = 0.0
bokeh.io.show(plot)

import numpy as np
import bokeh.io
import bokeh.plotting
import bokeh.palettes
from bokeh.models import Range1d

bokeh.io.output_notebook()

# Kinetic parameters
th = 1
rho = 2
d = 0.1
p = 0.05

# Simulation time
t = np.linspace(0,80,801)

# Writing the solution to the ODEs
mgfp = (th/d) - (th * np.exp(-d*t))/d
gfp = -(rho*th*np.exp(-p*t))/(d*p - p**2) + (rho*th*np.exp(-d*t))/(d*(d - p)) + rho*th/(d*p)

# Setting the bokeh visualization
plot = bokeh.plotting.figure(width = 320, height = 225,
                             x_axis_label = 'Time (minutes)', y_axis_label = 'Expression level')

# Plotting the solutions
plot.line(t, mgfp, line_width = 3, color = 'grey', legend_label = 'mRNA')
plot.line(t, gfp, line_width = 3, color = 'orange', legend_label = 'protein')

# Visualization
plot.legend.location = 'top_left'
plot.legend.background_fill_alpha = 0.0
bokeh.io.show(plot)

Sin embargo, en la práctica es bastante raro encontrarse con un sistema de EDOs para el cual podamos escribir una solución simbólica. En biología, los sistemas suelen ser no lineales, lo que hace que resolver las ecuaciones a mano sea casi imposible, incluso usando software especializados que trabajan con álgebra simbólica, como MatLab o Mathematica. Lo más común es aplicar métodos numéricos. En nuestro caso, vamos a usar la función integrate.odeint del paquete scipy. Es rápida, fácil de usar, y como los sistemas con los que estamos trabajando son pequeños, no necesitamos nada más avanzado.


import scipy.integrate

# Necesitamos definir el sistema de EDOs que queremos resolver
def TXTL(X, t, th, rho, d, p):

  mgfp, gfp = X

  return(
      [
          th - d*mgfp,
          rho*mgfp - p*gfp
      ]
  )

# Parámetros cinéticos
th = 1
rho = 2
d = 0.1
p = 0.05

# Tiempo de simulación y condiciones iniciales
t = np.linspace(0,80,801)
x0 = np.zeros(2)

# Resolviendo numéricamente las EDOs
args1 = (th, rho, d, p)
sol = scipy.integrate.odeint(TXTL, x0, t, args = args1)


import scipy.integrate

# We need to define the system of ODEs that we want to solve
def TXTL(X, t, th, rho, d, p):

  mgfp, gfp = X

  return(
      [
          th - d*mgfp,
          rho*mgfp - p*gfp
      ]
  )

# Kinetic parameters
th = 1
rho = 2
d = 0.1
p = 0.05

# Simulation time and initial conditions
t = np.linspace(0,80,801)
x0 = np.zeros(2)

# Numerically solving the ODEs
args1 = (th, rho, d, p)
sol = scipy.integrate.odeint(TXTL, x0, t, args = args1)

La solución de nuestras EDOs ahora está almacenada en la variable sol, que es un vector de tamaño $(n,m)$. Aquí, $n$ corresponde a la resolución de la simulación, que fue definida por el array de tiempo. En nuestro ejemplo, el array de tiempo es $t$ y tiene una resolución (o número de pasos) de 801. Por otro lado, $m$ es el número de variables de estado en nuestra función TXTL. Para acceder a la solución, debemos tener en cuenta el orden en que nuestra función retorna las variables de estado. Por ejemplo, sol[:,0] contiene la solución para $m_{\text{gfp}}$, mientras que sol[:,1] contiene la solución para $\text{gfp}$.

Visualizando los resultados


# Configurar la visualización con bokeh
plot = bokeh.plotting.figure(width = 320, height = 225,
                             x_axis_label = 'Tiempo (minutos)', y_axis_label = 'Nivel de expresión')

# Graficar las soluciones
plot.line(t, sol[:,0], line_width = 3, color = 'grey', legend_label = 'mRNA')
plot.line(t, sol[:,1], line_width = 3, color = 'orange', legend_label = 'proteína')

# Visualización
plot.legend.location = 'top_left'
plot.legend.background_fill_alpha = 0.0
bokeh.io.show(plot)


# Setting the bokeh visualization
plot = bokeh.plotting.figure(width = 320, height = 225,
                             x_axis_label = 'Time (minutes)', y_axis_label = 'Expression level')

# Plotting the solutions
plot.line(t, sol[:,0], line_width = 3, color = 'grey', legend_label = 'mRNA')
plot.line(t, sol[:,1], line_width = 3, color = 'orange', legend_label = 'protein')

# Visualization
plot.legend.location = 'top_left'
plot.legend.background_fill_alpha = 0.0
bokeh.io.show(plot)

Solución numérica de EDOs utilizando SciPy