RoadRunner in pyPESTO

After going through this notebook, you will be able to…

  • … create a pyPESTO problem using RoadRunner as a simulator directly from a PEtab problem.

  • … perform a parameter estimation using pyPESTO with RoadRunner as a simulator, setting advanced simulator features.

[1]:
# install pyPESTO with roadrunner support
# %pip install pypesto[roadrunner,petab] --quiet
[2]:
# import
import matplotlib as mpl
import numpy as np
import petab
import pypesto.objective
import pypesto.optimize as optimize
import pypesto.objective.roadrunner as pypesto_rr
import pypesto.sample as sample
import pypesto.visualize as visualize
import pypesto.profile as profile
from IPython.display import Markdown, display
from pprint import pprint

mpl.rcParams["figure.dpi"] = 100
mpl.rcParams["font.size"] = 18

np.random.seed(1912)


# name of the model that will also be the name of the python module
model_name = "boehm_JProteomeRes2014"

# output directory
model_output_dir = "tmp/" + model_name

Creating pyPESTO problem from PEtab

The PEtab file format stores all the necessary information to define a parameter estimation problem. This includes the model, the experimental data, the parameters to estimate, and the experimental conditions. Using the pypesto_rr.PetabImporterRR class, we can create a pyPESTO problem directly from a PEtab problem.

[3]:
petab_yaml = f"./{model_name}/{model_name}.yaml"

petab_problem = petab.Problem.from_yaml(petab_yaml)
importer = pypesto_rr.PetabImporterRR(petab_problem)
problem = importer.create_problem()

We now have a pyPESTO problem that we can use to perform parameter estimation. We can get some information on the RoadRunnerObjective and access the RoadRunner model.

[4]:
pprint(problem.objective.get_config())
{'roadrunner_instance': '<roadrunner.RoadRunner() { \n'
                        "'this' : 0x5600e823b4c0\n"
                        "'modelLoaded' : true\n"
                        "'modelName' : FullModel\n"
                        "'libSBMLVersion' : LibSBML Version: 5.19.5\n"
                        "'jacobianStepSize' : 1e-05\n"
                        "'steadyStateThreshold' : 0.01\n"
                        "'fluxThreshold' : 1e-09\n"
                        "'conservedMoietyAnalysis' : false\n"
                        "'simulateOptions' : \n"
                        '< roadrunner.SimulateOptions() \n'
                        '{ \n'
                        "'this' : 0x5600e8236b28, \n"
                        "'reset' : 0,\n"
                        "'structuredResult' : 0,\n"
                        "'copyResult' : 1,\n"
                        "'steps' : 50,\n"
                        "'start' : 0,\n"
                        "'duration' : 5\n"
                        "'output_file' : \n"
                        '}>, \n'
                        "'integrator' : \n"
                        '< roadrunner.Integrator() >\n'
                        '  name: cvode\n'
                        '  settings:\n'
                        '      relative_tolerance: 1e-06\n'
                        '      absolute_tolerance: 1e-12\n'
                        '                   stiff: true\n'
                        '       maximum_bdf_order: 5\n'
                        '     maximum_adams_order: 12\n'
                        '       maximum_num_steps: 20000\n'
                        '       maximum_time_step: 0\n'
                        '       minimum_time_step: 0\n'
                        '       initial_time_step: 0\n'
                        '          multiple_steps: false\n'
                        '      variable_step_size: false\n'
                        '         max_output_rows: 100000\n'
                        '\n'
                        '}>',
 'solver_options': "SolverOptions({'integrator': 'cvode', "
                   "'relative_tolerance': 1e-06, 'absolute_tolerance': 1e-12, "
                   "'maximum_num_steps': 20000})",
 'type': 'RoadRunnerObjective',
 'x_names': ['Epo_degradation_BaF3',
             'k_exp_hetero',
             'k_exp_homo',
             'k_imp_hetero',
             'k_imp_homo',
             'k_phos',
             'sd_pSTAT5A_rel',
             'sd_pSTAT5B_rel',
             'sd_rSTAT5A_rel']}
[5]:
# direct simulation of the model using roadrunner
sim_res = problem.objective.roadrunner_instance.simulate(
    times=[0, 2.5, 5, 10, 20, 50]
)
pprint(sim_res)
problem.objective.roadrunner_instance.plot();
    time, [STAT5A], [STAT5B],  [pApB],      [pApA],      [pBpB], [nucpApA], [nucpApB], [nucpBpB]
 [[    0,  143.867,  63.7332,       0,           0,           0,         0,         0,         0],
  [  2.5,  53.4606,  28.8659,  17.304, 5.43139e-05, 1.58348e-05,   112.991,   1.44768,   26.5969],
  [    5,  34.0638,  19.6392, 21.0102,  2.0613e-05, 6.85177e-06,   136.157,   3.93134,   33.9425],
  [   10,  21.7736,  12.8932, 22.6399, 7.35936e-06, 2.58051e-06,   149.923,   9.56434,   39.0847],
  [   20,  16.3397,  8.95152, 21.5682, 3.16462e-06, 9.49793e-07,   154.368,    20.913,   41.2089],
  [   50,  18.3666,  7.66901, 15.5584, 1.78015e-06, 3.10368e-07,   146.417,   49.2072,   38.4054]]

../_images/example_roadrunner_7_1.png

For more details on interacting with the roadrunner instance, we refer to the documentation of libroadrunner. However, we point out that while theoretical possible, we strongly advice against changing the model in that manner.

[6]:
ret = problem.objective(
    petab_problem.get_x_nominal(fixed=False,scaled=True),
    mode="mode_fun",
    return_dict=True,
)
pprint(ret)
{'fval': 151.50489751558354,
 'llh': -151.50489751558354,
 'simulation_results': {'model1_data1':     time, pSTAT5A_rel, pSTAT5B_rel, rSTAT5A_rel
 [[    0,           0,          -0,     21.2893],
  [  2.5,      75.155,     40.1655,     34.8294],
  [    5,     85.2164,     54.5033,      39.009],
  [   10,     90.6695,     66.2879,     42.2331],
  [   15,     92.1423,      70.943,     43.5006],
  [   20,     92.5019,      72.959,     44.0944],
  [   30,     91.9629,     73.5974,     44.4535],
  [   40,     90.6253,      72.083,      44.302],
  [   50,     88.7858,     69.4358,     43.8852],
  [   60,     86.5229,     66.0721,      43.299],
  [   80,     80.8019,     58.0459,     41.8054],
  [  100,     73.5773,     49.2225,     40.0835],
  [  120,     65.1055,     40.4541,     38.3234],
  [  160,     46.4623,     25.2907,     35.2469],
  [  200,     29.5146,     14.6585,     33.1842],
  [  240,     17.1454,     8.09325,     32.0669]]
}}

Optimization

To optimize a problem using a RoadRunner objective, we can set additional solver options for the ODE solver.

[7]:
optimizer = optimize.ScipyOptimizer()

solver_options = pypesto_rr.SolverOptions(
    relative_tolerance = 1e-6,
    absolute_tolerance = 1e-12,
    maximum_num_steps = 10000
)
engine = pypesto.engine.SingleCoreEngine()
problem.objective.solver_options = solver_options
[8]:
result = optimize.minimize(
    problem=problem,
    optimizer=optimizer,
    n_starts=5,  # usually a value >= 100 should be used
    engine=engine
)
display(Markdown(result.summary()))

Optimization Result

  • number of starts: 5

  • best value: 159.08150326977574, id=2

  • worst value: 2028.1213376090186, id=3

  • number of non-finite values: 0

  • execution time summary:

    • Mean execution time: 0.401s

    • Maximum execution time: 0.847s, id=2

    • Minimum execution time: 0.135s, id=3

  • summary of optimizer messages:

    Count

    Message

    3

    CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH

    1

    STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT

    1

    ABNORMAL_TERMINATION_IN_LNSRCH

  • best value found (approximately) 1 time(s)

  • number of plateaus found: 1

A summary of the best run:

Optimizer Result

  • optimizer used: <ScipyOptimizer method=L-BFGS-B options={‘disp’: False, ‘maxfun’: 1000}>

  • message: STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT

  • number of evaluations: 1010

  • time taken to optimize: 0.847s

  • startpoint: [ 0.17032373 2.40671846 -0.42272524 -0.68191772 -4.36152045 -0.13272644 -0.41145952 3.23408298 -3.50203955]

  • endpoint: [-1.43879756 2.02948835 -0.74039897 -1.82283304 -1.89156185 3.99816474 0.85597109 0.94022365 -3.50203955]

  • final objective value: 159.08150326977574

  • final gradient value: [-1.4679813 0.01610658 0.16264607 0.11554278 -0.56849672 0.37223345 0.46586024 -0.12937278 0. ]

Disclaimer: Currently there are two main things not yet fully supported with roadrunner objectives. One is parallelization of the optimization using MultiProcessEngine. The other is explicit gradients of the objective function. While the former will be added in a near future version, we will show a workaround for the latter, until it is implemented.

Visualization Methods

In order to visualize the optimization, there are a few things possible. For a more extensive explanation we refer to the “getting started” notebook.

[9]:
visualize.waterfall(result);
../_images/example_roadrunner_15_0.png
[10]:
visualize.parameters(result);
../_images/example_roadrunner_16_0.png
[11]:
visualize.parameters_correlation_matrix(result);
../_images/example_roadrunner_17_0.png

Sensitivities via finite differences

Some solvers need a way to calculate the sensitivities, which currently RoadRunner Objectives do not suport. For this scenario, we can use the FiniteDifferences objective in pypesto, which wraps a given objective into one, that computes sensitivities via finite differences.

[12]:
# no support for sensitivities
try:
    ret = problem.objective(
        petab_problem.get_x_nominal(fixed=False,scaled=True),
        mode="mode_fun",
        return_dict=True,
        sensi_orders=(1,),
    )
    pprint(ret)
except ValueError as e:
    pprint(e)
ValueError('This Objective cannot be called with sensi_orders= (1,) and mode=mode_fun.')
[13]:
objective_fd = pypesto.objective.FD(problem.objective)
# support through finite differences
try:
    ret = objective_fd(
        petab_problem.get_x_nominal(fixed=False,scaled=True),
        mode="mode_fun",
        return_dict=True,
        sensi_orders=(1,),
    )
    pprint(ret)
except ValueError as e:
    pprint(e)
{'grad': array([-3.24436940e-02, -3.59825663e-02,  6.61385968e-02,  5.43096235e-03,
        2.54184527e+01, -2.87776913e-02,  1.06348068e-02,  2.40250699e-02,
        1.91774348e-02])}