Simulaciones estocásticas utilizando el algoritmo de Gillespie

2025 Workshop en Modelamiento de Sistemas Biológicos

Autor: Bsc. Frank Britto Bisso | Revisado por: PhD. Christian Cuba Samaniego

Última fecha de modificación: 06 de Julio, 2025

En esta sección vamos a explorar cómo implementar simulaciones estocásticas usando el paquete biocircuits. Para una revisión y más detalles sobre la lógica y las matemáticas detrás del algoritmo, los autores del paquete tienen una excelente explicación que se puede consultar aquí.

Matriz de estequiometría (stoichiometry)

La matriz de estequiometría describe cómo cambia el número de moléculas de cada especie química involucrada en la reacción. Para entender cómo escribir esta matriz, recordemos nuestro modelo de transcripción-traducción (TXTL) de un gen reportero del Technical Note #1. Si $m_{\text{gfp}}$ representa el producto de la transcripción (es decir, el mRNA), y $\text{GFP}$ es el producto de la traducción (es decir, la proteína), las reacciones químicas que describen el sistema TXTL son las siguientes:

\begin{eqnarray} \text{DNA} &\xrightarrow{\theta} \text{DNA} + M_{\text{gfp}} && \quad M_{\text{gfp}} \xrightarrow{\delta} \emptyset && \quad \text{Transcripción y degradación de mRNA}\\ M_{\text{gfp}} &\xrightarrow{\rho} M_{\text{gfp}} + \text{GFP} && \text{GFP} \xrightarrow{\phi} \emptyset && \quad \text{Traducción y degradación de la proteina} \end{eqnarray}

Para formar nuestra matriz de estequiometría, reescribamos las reacciones químicas de modo que formen una columna. Al lado, llevaremos la cuenta si la reacción química se encuentra produciendo una molécula (para lo cual escribiremos un $+1$), o degradando una molécula (para lo cual escribiremos un $-1$). Si una especie química no participa en una reacción, le asignaremos un cero (0).

\begin{align*} && M_{\text{gfp}} && \text{GFP} \\ \text{DNA} &\xrightarrow{\theta} \text{DNA} + M_{\text{gfp}} & +1 && 0 \\ M_{\text{gfp}} & \xrightarrow{\delta} \emptyset & -1 && 0\\ M_{\text{gfp}} &\xrightarrow{\rho} M_{\text{gfp}} + \text{GFP} & 0 && +1\\ \text{GFP} &\xrightarrow{\phi} \emptyset & 0 && -1 \end{align*}


import numpy as np

# Stoichiometry matrix
simple_update = np.array(
    [
        # mgfp, gfp
        [ 1,     0], # DNA -> DNA + mRNA        (producción por transcripción)
        [-1,     0], # mRNA -> 0                (degradación)
        [ 0,     1], # mRNA -> mRNA + protein   (producción por traslación)
        [ 0,    -1], # protein -> 0             (degradación)
    ]
)


import numpy as np

# Stoichiometry matrix
simple_update = np.array(
    [
        # mgfp, gfp
        [ 1,     0], # DNA -> DNA + mRNA        (production via transcription)
        [-1,     0], # mRNA -> 0                (degradation)
        [ 0,     1], # mRNA -> mRNA + protein   (production via translation)
        [ 0,    -1], # protein -> 0             (degradation)
    ]
)

La matriz de probabilidades (propensities)

La matriz de probabilidades, como bien describe su nombre, representa las probabilidades de que cada reacción química ocurra en un determinado instante de tiempo. Dicha probabilidad está dada por los elementos de la parte derecha de las ecuaciones diferenciales que escribimos para modelar las reacciones químicas. Una buena práctica antes de escribir la matriz de probabilidades es analizar la dimensionalidad de las tasas descritas en las reacciones químicas. Esto debido a que cuando usamos el algoritmo de Gillespie del paquete biocircuits, todas las concentraciones se interpretan como número de moléculas. Esto significa que debemos asegurarnos de que los parámetros de nuestro modelo estén expresados en las unidades correctas.

\begin{align*} && M_{\text{gfp}} && \text{GFP} && \text{Tasa}\\ \text{DNA} &\xrightarrow{\theta} \text{DNA} + M_{\text{gfp}} & +1 && 0 && \theta \left [\frac{\mu M}{h} \right ] \\ M_{\text{gfp}} &\xrightarrow{\delta} \emptyset & -1 && 0 && \delta \left [\frac{1}{h} \right ]\\ M_{\text{gfp}} &\xrightarrow{\rho} M_{\text{gfp}} + \text{GFP} & 0 && +1 && \rho \left [\frac{1}{h} \right ]\\ \text{GFP} &\xrightarrow{\phi} \emptyset & 0 && -1 && \phi \left [\frac{1}{h} \right ] \end{align*}

Las tasas de traducción, degradación del mRNA y degradación de la proteína ($\rho$, $\delta$ y $\phi$) ya tienen unidades de 'por tiempo' ($1/h$), por lo que no necesitan convertirse. Sin embargo, la tasa de transcripción $\theta$ está definida en unidades de concentración por tiempo ($\mu M/h$), así que necesitamos convertirla a 'moléculas por unidad de tiempo'. Para hacer esto, debemos multiplicar $\theta$ por un factor de volumen llamado $\Omega$. Su valor depende del volumen efectivo de la célula, y garantiza que la simulación estocástica sea coherente con la escala física del sistema. Por lo tanto, nuestra matriz de probabilidades será la siguiente:

\begin{align*} && M_{\text{gfp}} && \text{GFP} && \text{Tasa} && \text{Probabilidad}\\ \text{DNA} &\xrightarrow{\theta} \text{DNA} + M_{\text{gfp}} & +1 && 0 && \theta \left [\frac{\mu M}{h} \right ] && \theta \times \Omega\\ M_{\text{gfp}} &\xrightarrow{\delta} \emptyset & -1 && 0 && \delta \left [\frac{1}{h} \right ] && \delta \times m_{\text{gfp}}\\ M_{\text{gfp}} &\xrightarrow{\rho} M_{\text{gfp}} + \text{GFP} & 0 && +1 && \rho \left [\frac{1}{h} \right ] && \rho \times m_{\text{gfp}}\\ \text{GFP} &\xrightarrow{\phi} \emptyset & 0 && -1 && \phi \left [\frac{1}{h} \right ] && \phi \times \text{gfp} \end{align*}


# Propensities matrix
def simple_propensity(propensities, population, t, th, rho, d, p, Omega):
  mgfp, gfp = population

  propensities[0] = th * Omega
  propensities[1] = d * mgfp
  propensities[2] = rho * mgfp
  propensities[3] = p * gfp


# Propensities matrix
def simple_propensity(propensities, population, t, th, rho, d, p, Omega):
  mgfp, gfp = population

  propensities[0] = th * Omega
  propensities[1] = d * mgfp
  propensities[2] = rho * mgfp
  propensities[3] = p * gfp

Importante: El orden en el que se escriben las probabilidades es relevante, ya que la matriz de estequiometría debe coincidir con las probabilidades de cada reacción y reflejar correctamente el cambio en el número de moléculas.

Utilizando el algoritmo Gillespie

Como resumen, definimos la matriz de probabilidades para indicarle al algoritmo de Gillespie la probabilidad de que una reacción aumente o disminuya el número de moléculas, lo cual está determinado por la matriz de estequiometría. Para implementar el algoritmo y visualizar los resultados, vamos a usar los paquetes biocircuits y bokeh.


import os, sys, subprocess
if "google.colab" in sys.modules:
  cmd = "pip install --upgrade biocircuits bokeh-catplot watermark blackcellmagic intersect"
  process = subprocess.Popen(cmd.split(), stdout=subprocess.PIPE, stderr=subprocess.PIPE)
  stdout, stderr = process.communicate()
# ------

import biocircuits

import bokeh.io
import bokeh.plotting
import bokeh.palettes
from bokeh.models import Range1d

bokeh.io.output_notebook()


# Kinetic parameters
th = 0.5
rho = 4
d = 0.1
p = 0.05

# Definir el valor de conversión a número de molécules
Omega = 600 # e.g., 600: volumen de una célula E. Coli

# Definir el tiempo de simulación, condiciones iniciales y parámetros del sistema
time_points = np.linspace(0, 80, 801)
population_0 = np.zeros(2, dtype=int) # condiciones iniciales
args1 = (th, rho, d, p, Omega)


# Kinetic parameters
th = 0.5
rho = 4
d = 0.1
p = 0.05

# Define the conversion factor to number of molecules
Omega = 600 # e.g., E. Coli cell volume

# Define simulation time, initial conditions, and system parameters
time_points = np.linspace(0, 80, 801)
population_0 = np.zeros(2, dtype=int) # initial conditions
args1 = (th, rho, d, p, Omega)

La función que estamos usando se llama gillespie_ssa , y recibe como argumentos (no necesariamente en este orden): la matriz de propensiones, la matriz de estequiometría, las condiciones iniciales, el tiempo de simulación, los parámetros cinéticos y dos parámetros adicionales llamados size y n_threads, cuyo producto indica cuántas trayectorias estocásticas se van a simular. Además, como las simulaciones estocásticas suelen tomar algo de tiempo en resolverse, es buena práctica guardar los resultados una vez que terminan.


# Resolviendo utilizando el algoritmo Gillespie
samples = biocircuits.gillespie_ssa(simple_propensity, simple_update, population_0, time_points, size=20, args=args1, n_threads=2, progress_bar=True)

# Guardando los resultados
np.save("TXTL_stochastic.npy", samples)

 100%|██████████| 20/20 [11:35<00:00, 34.77s/it]
 100%|██████████| 20/20 [11:35<00:00, 34.77s/it]


# Solve using the Gillespie algorithm
samples = biocircuits.gillespie_ssa(simple_propensity, simple_update, population_0, time_points, size=20, args=args1, n_threads=2, progress_bar=True)

# Save results
np.save("TXTL_stochastic.npy", samples)

 100%|██████████| 20/20 [11:35<00:00, 34.77s/it]
 100%|██████████| 20/20 [11:35<00:00, 34.77s/it]

Finalmente, visualizamos los resultados:


# Configuración de la visualización con bokeh
fig_size = [320, 250]
plots = []

for i in range(2):
  plots.append(bokeh.plotting.figure(width = fig_size[0], height = fig_size[1],
                                     x_axis_label = 'Tiempo (horas)', y_axis_label = 'Numéro de moléculas'))

for i in range(len(samples)):
  if i == 0:
    plots[0].line(time_points, samples[i,:,0], line_width = 3, color = 'grey', legend_label = 'mRNA')
    plots[1].line(time_points, samples[i,:,1], line_width = 3, color = 'orange', legend_label = 'proteina')
  else:
    plots[0].line(time_points, samples[i,:,0], line_width = 3, color = 'grey', legend_label = 'mRNA')
    plots[1].line(time_points, samples[i,:,1], line_width = 3, color = 'orange', legend_label = 'proteina')

# Visualización
for i in range(2):
  plots[i].legend.location = 'bottom_right'
bokeh.io.show(bokeh.layouts.row(plots))


# Visualization Setup with Bokeh
fig_size = [320, 250]
plots = []

for i in range(2):
  plots.append(bokeh.plotting.figure(width = fig_size[0], height = fig_size[1],
                                     x_axis_label = 'Time (hours)', y_axis_label = 'Number of molecules'))

for i in range(len(samples)):
  if i == 0:
    plots[0].line(time_points, samples[i,:,0], line_width = 3, color = 'grey', legend_label = 'mRNA')
    plots[1].line(time_points, samples[i,:,1], line_width = 3, color = 'orange', legend_label = 'protein')
  else:
    plots[0].line(time_points, samples[i,:,0], line_width = 3, color = 'grey', legend_label = 'mRNA')
    plots[1].line(time_points, samples[i,:,1], line_width = 3, color = 'orange', legend_label = 'protein')

# Visualization
for i in range(2):
  plots[i].legend.location = 'bottom_right'
bokeh.io.show(bokeh.layouts.row(plots))