Ordinary differential equations: Conversion reaction

This example was kindly contributed by Lukas Sandmeir and Elba Raimundez. It can be downloaded here: Ordinary Differential Equations.

Note: Before you use pyABC to parametrize your ODE, please be aware of potential errors introduced by inadequately representing the data generation process, see also the “Measurement noise assessment” notebook. For deterministic models, there are often more efficient alternatives to ABC, check out for example our tool pyPESTO.

This example provides a model for the interconversion of two species (\(X_1\) and \(X_2\)) following first-order mass action kinetics with the parameters \(\Theta_1\) and \(\Theta_2\) respectively:

\[X_1 \rightarrow X_2, \quad\text{rate} = \Theta_1 \cdot [X_1]\]

\[X_2 \rightarrow X_1, \quad\text{rate} = \Theta_2 \cdot [X_2]\]

Measurement of \([X_2]\) is provided as \(Y = [X_2]\).

We will show how to estimate \(\Theta_1\) and \(\Theta_2\) using pyABC.

# install if not done yet
!pip install pyabc --quiet

%matplotlib inline

import os
import tempfile

import matplotlib.pyplot as plt
import numpy as np
import scipy as sp

from pyabc import ABCSMC, RV, Distribution, LocalTransition, MedianEpsilon
from pyabc.visualization import plot_data_callback, plot_kde_2d

db_path = "sqlite:///" + os.path.join(tempfile.gettempdir(), "test.db")

Data

We use an artificial data set which consists of a vector of time points \(t\) and a measurement vector \(Y\). This data was created using the parameter values which are assigned to \(\Theta_{\text{true}}\) and by adding normaly distributed measurement noise with variance \(\sigma^2 = 0.015^2\).

ODE model

\[\begin{split}\begin{align*} \frac{dX_1}{dt} &= -\Theta_1 \cdot X_1 + \Theta_2 \cdot X_2\\ \frac{dX_2}{dt} &= \Theta_1 \cdot X_1 - \Theta_2 \cdot X_2 \end{align*}\end{split}\]

Define the true parameters

theta1_true, theta2_true = np.exp([-2.5, -2])

and the measurement data

measurement_data = np.array(
    [
        0.0244,
        0.0842,
        0.1208,
        0.1724,
        0.2315,
        0.2634,
        0.2831,
        0.3084,
        0.3079,
        0.3097,
        0.3324,
    ]
)

as well as the time points at whith to evaluate

measurement_times = np.arange(len(measurement_data))

and the initial conditions for \(X_1\) and \(X_2\)

init = np.array([1, 0])

Define the ODE model

def f(y, t0, theta1, theta2):
    x1, x2 = y
    dx1 = -theta1 * x1 + theta2 * x2
    dx2 = theta1 * x1 - theta2 * x2
    return dx1, dx2


def model(pars):
    sol = sp.integrate.odeint(
        f, init, measurement_times, args=(pars["theta1"], pars["theta2"])
    )
    return {"X_2": sol[:, 1]}

Integration of the ODE model for the true parameter values

true_trajectory = model({"theta1": theta1_true, "theta2": theta2_true})["X_2"]

Let’s visualize the results

plt.plot(true_trajectory, color="C0", label='Simulation')
plt.scatter(measurement_times, measurement_data, color="C1", label='Data')
plt.xlabel('Time $t$')
plt.ylabel('Measurement $Y$')
plt.title('Conversion reaction: True parameters fit')
plt.legend()
plt.show()

def distance(simulation, data):
    return np.absolute(data["X_2"] - simulation["X_2"]).sum()

Define the prior for \(\Theta_1\) and \(\Theta_2\)

parameter_prior = Distribution(
    theta1=RV("uniform", 0, 1), theta2=RV("uniform", 0, 1)
)
parameter_prior.get_parameter_names()

['theta1', 'theta2']

abc = ABCSMC(
    models=model,
    parameter_priors=parameter_prior,
    distance_function=distance,
    population_size=50,
    transitions=LocalTransition(k_fraction=0.3),
    eps=MedianEpsilon(500, median_multiplier=0.7),
)

ABC.Sampler INFO: Parallelize sampling on 4 processes.

abc.new(db_path, {"X_2": measurement_data});

ABC.History INFO: Start <ABCSMC id=4, start_time=2021-11-18 17:40:36>

h = abc.run(minimum_epsilon=0.1, max_nr_populations=5)

ABC INFO: t: 0, eps: 5.00000000e+02.
ABC INFO: Accepted: 50 / 53 = 9.4340e-01, ESS: 5.0000e+01.
ABC INFO: t: 1, eps: 1.89305859e+00.
ABC INFO: Accepted: 50 / 164 = 3.0488e-01, ESS: 4.4195e+01.
ABC INFO: t: 2, eps: 9.64400074e-01.
ABC INFO: Accepted: 50 / 159 = 3.1447e-01, ESS: 3.7093e+01.
ABC INFO: t: 3, eps: 4.96905744e-01.
ABC INFO: Accepted: 50 / 329 = 1.5198e-01, ESS: 3.5372e+01.
ABC INFO: t: 4, eps: 2.94398923e-01.
ABC INFO: Accepted: 50 / 448 = 1.1161e-01, ESS: 4.4592e+01.
ABC INFO: Stop: Maximum number of generations.
ABC.History INFO: Done <ABCSMC id=4, duration=0:00:03.032233, end_time=2021-11-18 17:40:39>

Visualization of the probability density functions for \(\Theta_1\) and \(\Theta_2\)

fig = plt.figure(figsize=(10, 8))
for t in range(h.max_t + 1):
    ax = fig.add_subplot(3, int(np.ceil(h.max_t / 3)), t + 1)

    ax = plot_kde_2d(
        *h.get_distribution(m=0, t=t),
        "theta1",
        "theta2",
        xmin=0,
        xmax=1,
        numx=200,
        ymin=0,
        ymax=1,
        numy=200,
        ax=ax,
    )
    ax.scatter(
        [theta1_true],
        [theta2_true],
        color="C1",
        label=r'$\Theta$ true = {:.3f}, {:.3f}'.format(
            theta1_true, theta2_true
        ),
    )
    ax.set_title(f"Posterior t={t}")

    ax.legend()
fig.tight_layout()

We can also plot the simulated trajectories:

_, ax = plt.subplots()


def plot_data(sum_stat, weight, ax, **kwargs):
    """Plot a single trajectory"""
    ax.plot(measurement_times, sum_stat['X_2'], color='grey', alpha=0.1)


def plot_mean(sum_stats, weights, ax, **kwargs):
    """Plot mean over all samples"""
    weights = np.array(weights)
    weights /= weights.sum()
    data = np.array([sum_stat['X_2'] for sum_stat in sum_stats])
    mean = (data * weights.reshape((-1, 1))).sum(axis=0)
    ax.plot(measurement_times, mean, color='C2', label='Sample mean')


ax = plot_data_callback(h, plot_data, plot_mean, ax=ax)

plt.plot(true_trajectory, color="C0", label='Simulation')
plt.scatter(measurement_times, measurement_data, color="C1", label='Data')
plt.xlabel('Time $t$')
plt.ylabel('Measurement $Y$')
plt.title('Conversion reaction: Simulated data fit')
plt.legend()
plt.show()