Ordinary Differential Equations: Conversion Reaction

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

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.

In [1]:
%matplotlib inline

from pyabc import (ABCSMC,
                   RV, Distribution,
from pyabc.visualization import plot_kde_2d
import matplotlib.pyplot as plt
import os
import tempfile
import scipy as sp

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


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

In [2]:
theta1_true, theta2_true = sp.exp([-2.5, -2])

and the measurement data

In [3]:
measurement_data = sp.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

In [4]:
measurement_times = sp.arange(len(measurement_data))

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

In [5]:
init = sp.array([1, 0])

Define the ODE model

In [6]:
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,
    return {"X_2": sol[:,1]}

Integration of the ODE model for the true parameter values

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

Let’s visualize the results

In [8]:
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')
In [9]:
def distance(simulation, data):
    return sp.absolute(data["X_2"] - simulation["X_2"]).sum()

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

In [10]:
parameter_prior = Distribution(theta1=RV("uniform", 0, 1),
                               theta2=RV("uniform", 0, 1))
['theta1', 'theta2']
In [11]:
abc = ABCSMC(models=model,
             eps=MedianEpsilon(500, median_multiplier=0.7))
In [12]:
abc.new(db_path, {"X_2": measurement_data});
INFO:History:Start <ABCSMC(id=1, start_time=2018-05-08 15:22:45.838120, end_time=None)>
INFO:Epsilon:initial epsilon is 500
In [13]:
h = abc.run(minimum_epsilon=0.1, max_nr_populations=5)
INFO:ABC:t:0 eps:500
INFO:ABC:t:1 eps:1.393207770332699
INFO:ABC:t:2 eps:0.5522857161704047
INFO:ABC:t:3 eps:0.29564846797006095
INFO:ABC:t:4 eps:0.1617727948665631
INFO:History:Done <ABCSMC(id=1, start_time=2018-05-08 15:22:45.838120, end_time=2018-05-08 15:22:50.357498)>

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

In [14]:
for t in range(h.max_t+1):
    ax = plot_kde_2d(*h.get_distribution(m=0, t=t),
                     "theta1", "theta2",
                xmin=0, xmax=1, numx=300,
                ymin=0, ymax=1, numy=300)
    ax.scatter([theta1_true], [theta2_true],
                label='$\Theta$ true = {:.3f}, {:.3f}'.format(
                    theta1_true, theta2_true))
    ax.set_title("Posterior t={}".format(t))