PEtab extention for FitMultiCell

In this example, we will show how a PEtab problem can be defined in FitMultiCell pipeline. For this, we will use a 1D activator-inhibitor model1. First, we will try to define the problem without the use of PEtab format.

Regular FitMultiCell problem

At first, lets’s import the required packages.

import petab_MS
from fitmulticell.PEtab.base import PetabImporter
from fitmulticell.model import MorpheusModel as morpheus_model
from fitmulticell.model import MorpheusModels as morpheus_models
from fitmulticell.sumstat import SummaryStatistics
from pyabc.sampler import RedisEvalParallelSampler
from pyabc.sampler import MulticoreEvalParallelSampler
from pyabc import QuantileEpsilon
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import math
import pyabc

import matplotlib.pylab as plt
from pathlib import Path
import os
import tempfile

Now, let’s define the xpath for each parameter.

par_map = {'rho_a': './Global/System/Constant[@symbol="rho_a"]',
           'mu_i': './Global/System/Constant[@symbol="mu_i"]',
           'mu_a': './Global/System/Constant[@symbol="mu_a"]',
           }

Here, we define the ground truth value for each parameter which we used to generate the senthitic data.

obs_pars = {'rho_a': 0.01,
            'mu_i': 0.03,
            'mu_a': 0.02,
           }
obs_pars_log = {key: math.log10(val) for key, val in obs_pars.items()}

Next, we will import the senthitic data and give the path to Morpheus model.

condition1_obs = str((Path(os.getcwd())) / 'PEtab_problem_1' / 'Small.csv')

model_file = str((Path(os.getcwd())) / 'PEtab_problem_1' / 'ActivatorInhibitor_1D.xml')

data = pd.read_csv(condition1_obs, sep='\t')
dict_data = {}
for col in data.columns:
    dict_data[col] = data[col].to_numpy()

Then, we will construct Morpheus model and SummaryStatistic object.

sumstat = SummaryStatistics(ignore=["time", "loc"])

model = morpheus_model(
    model_file, par_map=par_map, clean_simulation=True,
    par_scale="lin",
    sumstat=sumstat)

Let’s now run a simulation with the ground truth parameters and plot it to see how it looks.

tryjectory = model.sample(obs_pars)

color=['b','r']

color_index=0
size = 100
for key, val in tryjectory.items():
    if key == "loc": continue
    alpha_index = 1
    if key == "space.x": continue
    for i,time in zip(range (0,int(len(val)),50*size),range(0,4000,50*25)):
        plt.plot(tryjectory["space.x"][i:i+size], val[i:i+size],label=key +" - time= "+str(time),color=color[color_index],alpha=1/alpha_index)
        alpha_index += 1
    color_index+=1
plt.legend()

<matplotlib.legend.Legend at 0x7f509964ff10>

Then, we will further define the limits and the scale for parameters space.

limits = {key: (-3,-1) for
          key, val in obs_pars.items()}

model.par_scale="log10"

prior = pyabc.Distribution(**{key: pyabc.RV("uniform", lb, ub - lb)
                              for key, (lb, ub) in limits.items()})

Now, we will define our objective function:

def eucl_dist(sim, obs):
    total = 0
    for key in sim:
        if key in (
            'loc', "time", "space.x"):
                continue
        x = np.array(sim[key])
        y = np.array(obs[key])
        if x.size != y.size:
            return np.inf

        total += np.sum((x - y) ** 2)
    return total

Now, we are ready to start the fitting.

abc = pyabc.ABCSMC(model, prior, eucl_dist, population_size=100)

ABC.Sampler INFO: Parallelize sampling on 12 processes.

db_path = "sqlite:///" + os.path.join(tempfile.gettempdir(), "ActivatorInhibitor_1D.db")
history = abc.new(db_path, dict_data)
history = abc.run(max_nr_populations=8)

ABC.History INFO: Start <ABCSMC id=2, start_time=2022-01-20 11:45:18>
ABC INFO: Calibration sample t = -1.
ABC INFO: t: 0, eps: 6.92881712e+04.
ABC INFO: Accepted: 100 / 215 = 4.6512e-01, ESS: 1.0000e+02.
ABC INFO: t: 1, eps: 4.14571844e+04.
ABC INFO: Accepted: 100 / 254 = 3.9370e-01, ESS: 9.0927e+01.
ABC INFO: t: 2, eps: 2.72167705e+04.
ABC INFO: Accepted: 100 / 216 = 4.6296e-01, ESS: 9.0174e+01.
ABC INFO: t: 3, eps: 1.67243558e+04.
ABC INFO: Accepted: 100 / 244 = 4.0984e-01, ESS: 8.9091e+01.
ABC INFO: t: 4, eps: 1.10596384e+04.
ABC INFO: Accepted: 100 / 209 = 4.7847e-01, ESS: 7.9877e+01.
ABC INFO: t: 5, eps: 7.73285579e+03.
ABC INFO: Accepted: 100 / 290 = 3.4483e-01, ESS: 9.1969e+01.
ABC INFO: t: 6, eps: 6.57135021e+03.
ABC INFO: Accepted: 100 / 278 = 3.5971e-01, ESS: 8.3232e+01.
ABC INFO: t: 7, eps: 6.24093899e+03.
ABC INFO: Accepted: 100 / 265 = 3.7736e-01, ESS: 7.2012e+01.
ABC INFO: Stop: Maximum number of generations.
ABC.History INFO: Done <ABCSMC id=2, duration=0:05:08.607414, end_time=2022-01-20 11:50:26>

Let’s now plot the epsilon over different iteration.

pyabc.visualization.plot_epsilons(history)

<AxesSubplot:title={'center':'Epsilon values'}, xlabel='Population index', ylabel='Epsilon'>

Also, let’s plot the kernal density function for our parameters of interest.

df, w = history.get_distribution(t=history.max_t)
pyabc.visualization.plot_kde_matrix(df, w, limits=limits, refval=obs_pars_log)

array([[<AxesSubplot:ylabel='mu_a'>, <AxesSubplot:>, <AxesSubplot:>],
       [<AxesSubplot:ylabel='mu_i'>, <AxesSubplot:>, <AxesSubplot:>],
       [<AxesSubplot:xlabel='mu_a', ylabel='rho_a'>,
        <AxesSubplot:xlabel='mu_i'>, <AxesSubplot:xlabel='rho_a'>]],
      dtype=object)

Define the problem using PEtab extention

Next, we want to fit the same model but now we want to use PEtab format to construct our problem.

We will start by import a PEtab problem from the yamel file that we already define according to PEtab standared.

petab_problem_path = str((Path(os.getcwd())) / 'PEtab_problem_1' / 'ActivatorInhibitor_1D.yaml')
petab_problem = petab_MS.Problem.from_yaml(petab_problem_path)

importer = PetabImporter(petab_problem)

Then, we will generate all required component from the PEtab problem, such as prior, measurement, parameters’ scale, etc.

PEtab_prior = importer.create_prior()

/home/emad/Insync/blackhand.3@gmail.com/Google_Drive/Bonn/Github/libpetab-python-MS/petab_MS/parameters.py:375: RuntimeWarning: invalid value encountered in log10
  return np.log10(parameter)

par_map_imported = importer.get_par_map()

obs_pars_imported = petab_problem.get_x_nominal_dict(scaled=True)

PEtab_par_scale = petab_problem.get_optimization_parameter_scales()

dict_data_imported = petab_problem.get_measurement_dict()

After that, we will be able to construct our model form the PEtab importer

PEtab_model = importer.create_model()

Then, we can modify any of the model’s attributes.

PEtab_model.sumstat.ignore = ["time", "loc"]
PEtab_model.sumstat.sum_stat_calculator = None

let’s now run a simulation tryjectory for the imported PEtab problem.

PEtab_tryjectory = PEtab_model.sample(obs_pars_imported)

Let’s now visulise the model tryjectory for both regular and PEtab problem

fig, axes = plt.subplots(1,2, figsize=(16, 6))


color=['b','r']


ax = axes[0]
color_index=0
size = 100
for key, val in PEtab_tryjectory.items():
    if key == "loc": continue
    alpha_index = 1
    if key == "condition1__space.x": continue
    for i,time in zip(range (0,int(len(val)),50*size),range(0,4000,50*25)):
        ax.plot( PEtab_tryjectory["condition1__space.x"][i:i+size], val[i:i+size],label=key +" - time= "+str(time),color=color[color_index],alpha=1/alpha_index)
        alpha_index += 1
    color_index+=1
ax.set_title("PEtab problem")

ax = axes[1]
color_index=0
for key, val in tryjectory.items():
    if key == "loc": continue
    alpha_index = 1
    if key == "space.x": continue
    for i,time in zip(range (0,int(len(val)),50*size),range(0,4000,50*25)):
        ax.plot( tryjectory["space.x"][i:i+size], val[i:i+size],label=key +" - time= "+str(time),color=color[color_index],alpha=1/alpha_index)
        alpha_index += 1
    color_index+=1
ax.legend()
ax.set_title("Regular problem")

Text(0.5, 1.0, 'Regular problem')

Let’s get one final component before starting the fitting, which is the objective function

PEtab_obj_function = importer.get_objective_function()

Now, we can go ahead and start the fitting of the PEtab problem.

abc = pyabc.ABCSMC(PEtab_model, PEtab_prior, PEtab_obj_function, population_size=100)
db_path = "sqlite:///" + os.path.join(tempfile.gettempdir(), "PEtab_ActivatorInhibitor_1D.db")
PEtab_history = abc.new(db_path, dict_data_imported)
PEtab_history = abc.run(max_nr_populations=8)

ABC.Sampler INFO: Parallelize sampling on 12 processes.
ABC.History INFO: Start <ABCSMC id=3, start_time=2022-01-19 12:18:04>
ABC INFO: Calibration sample t = -1.
ABC INFO: t: 0, eps: 1.86342779e+05.
ABC INFO: Accepted: 100 / 209 = 4.7847e-01, ESS: 1.0000e+02.
ABC INFO: t: 1, eps: 4.15795513e+04.
ABC INFO: Accepted: 100 / 243 = 4.1152e-01, ESS: 8.5515e+01.
ABC INFO: t: 2, eps: 2.38253314e+04.
ABC INFO: Accepted: 100 / 257 = 3.8911e-01, ESS: 8.2928e+01.
ABC INFO: t: 3, eps: 1.50475712e+04.
ABC INFO: Accepted: 100 / 263 = 3.8023e-01, ESS: 7.3097e+01.
ABC INFO: t: 4, eps: 1.01178034e+04.
ABC INFO: Accepted: 100 / 267 = 3.7453e-01, ESS: 7.1077e+01.
ABC INFO: t: 5, eps: 7.30028266e+03.
ABC INFO: Accepted: 100 / 274 = 3.6496e-01, ESS: 7.4931e+01.
ABC INFO: t: 6, eps: 6.36639092e+03.
ABC INFO: Accepted: 100 / 350 = 2.8571e-01, ESS: 8.0178e+01.
ABC INFO: t: 7, eps: 6.11468392e+03.
ABC INFO: Accepted: 100 / 344 = 2.9070e-01, ESS: 7.8303e+01.
ABC INFO: Stop: Maximum number of generations.
ABC.History INFO: Done <ABCSMC id=3, duration=0:07:13.584303, end_time=2022-01-19 12:25:17>

Now, let us compare the regular fitting problem to that that was initialized by PEtab. lets start by plotting the epsilon of the two problems.

ax= pyabc.visualization.plot_epsilons(PEtab_history, colors=['C1'])
pyabc.visualization.plot_epsilons(history,ax=ax, colors=['C0'])
ax.legend(["PEtab problem","Regular problem"])

<matplotlib.legend.Legend at 0x7fb39a53db90>

Finally, let’s plot and compare the final distribution of the two problems.

fig, axes = plt.subplots(3,1, figsize=(6, 6))
ax =axes[0]
df, w = history.get_distribution(t=history.max_t)
pyabc.visualization.plot_kde_1d(df, w, 'mu_a',title="mu_a", refval=obs_pars_log, refval_color='C0', ax=ax)
df, w = PEtab_history.get_distribution(t=PEtab_history.max_t)
pyabc.visualization.plot_kde_1d(df, w, 'mu_a', refval_color='C1', ax=ax)
ax.legend(["regular","Ground truth","PEtab"])

ax =axes[1]
df, w = history.get_distribution(t=history.max_t)
pyabc.visualization.plot_kde_1d(df, w, 'mu_i', title="mu_i", refval=obs_pars_log, refval_color='C0', ax=ax)
df, w = PEtab_history.get_distribution(t=PEtab_history.max_t)
pyabc.visualization.plot_kde_1d(df, w, 'mu_i', refval_color='C1', ax=ax)
ax.legend(["regular","Ground truth","PEtab"])

ax =axes[2]
df, w = history.get_distribution(t=history.max_t)
pyabc.visualization.plot_kde_1d(df, w, 'rho_a', title="rho_a", refval=obs_pars_log, refval_color='C0', ax=ax)
df, w = PEtab_history.get_distribution(t=PEtab_history.max_t)
pyabc.visualization.plot_kde_1d(df, w, 'rho_a', refval_color='C1', ax=ax)
ax.legend(["regular","Ground truth","PEtab"])

<matplotlib.legend.Legend at 0x7fb396198210>

References:

1 A. Gierer, H. Meinhardt: A Theory of Biological Pattern Formation. Kybernetik 12: 30-39, 1972.