Model development

This section walk you through the creation of a mechanistic model from text.

How to use

pasmopy.Text2Model is a useful class to build an ordinary differential equation (ODE) model from a text file describing biochemical systems.

The text file you need to prepare can be divided into three parts:

  1. Reaction layer

  2. Observable layer (Prefix: @obs)

  3. Simulation layer (Prefix: @sim)

Note

To write a comment, simply put the hash mark # before your desired comment.

# This is a comment.

Reaction layer

description | parameters | initial conditions

In the reaction layer, you need to describe biochemical reactions. Each reaction described in the line number i will be converted into ith rate equation. To specify parameters or initial conditions, you can put those information after |.

  • If you don’t specify parameters/initial_conditions, they are initialized to 1 and 0, respectively, and the parameter values will be estimated from experimental data.

  • If you want to set a parameter value to 1 and don’t want to estimate, you can add const prefix:

# The Hill coefficient is fixed to 1.
TF transcribes mRNA | const n=1

  • You can impose parameter constraints by specifying line number in the parameter section.

1# Nucleocytoplasmic Shuttling of DUSP
2DUSPc translocates from the cytoplasm to the nucleus <--> DUSPn
3pDUSPc translocates from the cytoplasm to the nucleus <--> pDUSPn |2|

In the example above, you can assume that import and export rates were identical for DUSP (line 2) and pDUSP (line 3).

  • If the amount of a model species should be held fixed (never consumed) during simulation, you can add fixed prefix:

# [Ligand] will be held fixed to 10.0 during simulation
Ligand binds Receptor <--> LR | kf = 1e-6, kr = 1e-1 | fixed Ligand = 10.0

  • To describe more complex rate equations, you can use @rxn prefix:

@rxn Reactant --> Product: define rate equation here

Please also refer to the following example: cfos_model

Note

The available rules can be found at https://biomass-core.readthedocs.io/en/latest/api/reaction_rules.html.

You can also supply your own terminology in a reaction rule via:

from pasmopy import Text2Model

# Supply "releases" to the reaction rule: "dissociate"
mm_kinetics = Text2Model("michaelis_menten.txt")
mm_kinetics.register_word({"dissociate": ["releases"]})
# Now you can use "releases" in your text, e.g., 'ES releases E and P'
mm_kinetics.convert()

Observable layer (Prefix: @obs)

In the observable layer, you need to specify biomass.observables, which can correlate model simulations and experimental measurements. You can create an observable by using model parameters (p) and species (u). For example, if the total amount of SOS bound to EGFR should be the sum of RGS (EGFR-Grb2-SOS) and RShGS (EGFR-Shc-Grb2-SOS) complexes in your model, then you can write as follows:

@obs Total_SOS_bound_to_EGFR: u[RGS] + u[RShGS]

Simulation layer (Prefix: @sim)

In the simulation layer, you can set simulation conditions, e.g, the simulation time span, the initial concentration of model species, etc.

Example:

@sim tspan: [0, 120]
@sim unperturbed: init[EGF] = 0
@sim condition EGF20nM: init[EGF] = 680
@sim condition EGF2nM: init[EGF] = 68

  • tspan:

    Two element vector [t0, tf] specifying the initial and final times.

  • unperturbed (optional):

    Description of the untreated condition to find the steady state.

  • condition (optional):

    Experimental conditions. Use p and init to modify model parameters and initial conditions, respectively.

Examples

Michaelis-Menten enzyme kinetics

This example shows you how to build a simple Michaelis-Menten two-step enzyme catalysis model with Pasmopy.

E + S ⇄ ES → E + P

An enzyme, E, binding to a substrate, S, to form a complex, ES, which in turn releases a product, P, regenerating the original enzyme.

  1. Prepare a text file describing biochemical reactions (e.g., michaelis_menten.txt)

     1E + S <--> ES | kf=0.003, kr=0.001 | E=100, S=50
 2ES --> E + P | kf=0.002
 3
 4@obs Substrate: u[S]
 5@obs E_free: u[E]
 6@obs E_total: u[E] + u[ES]
 7@obs Product: u[P]
 8@obs Complex: u[ES]
 9
10@sim tspan: [0, 100]

  2. Convert the text into an executable model

    $ python

    >>> from pasmopy import Text2Model
>>> description = Text2Model("michaelis_menten.txt")
>>> description.convert()
Model information
-----------------
2 reactions
4 species
3 parameters

  3. Run simulation

    >>> from pasmopy import create_model, run_simulation
>>> model = create_model("michaelis_menten")
>>> run_simulation(model)
    _images/michaelis_menten_sim.png

EGF signaling

Below is an example of Pasmopy in action to illustrate EGF signalling pathway.

Reference:

Kholodenko, B. N., Demin, O. V, Moehren, G. & Hoek, J. B. Quantification of short term signaling by the epidermal growth factor receptor. J. Biol. Chem. 274, 30169–30181 (1999). https://doi.org/10.1074/jbc.274.42.30169

  1. Prepare a text describing EGF signaling in hepatocytes (Kholodenko1999.txt)

    You can learn how to build the model via Text2Model by comparing the desccription below and the network scheme in Figure 1 in the original paper. The description in the line n denotes the n-th reaction in the scheme.

     1EGF binds EGFR <--> Ra | kf=0.003, kr=0.06 | EGFR=100
 2Ra dimerizes <--> R2 | kf=0.01, kr=0.1
 3R2 is phosphorylated <--> RP | kf=1, kr=0.01
 4RP is dephosphorylated --> R2 | V=450, K=50
 5RP binds PLCg <--> RPL | kf=0.06, kr=0.2 | PLCg=105
 6RPL is phosphorylated <--> RPLP | kf=1, kr=0.05
 7RPLP is dissociated into RP and PLCgP | kf=0.3, kr=0.006
 8PLCgP is dephosphorylated --> PLCg | V=1, K=100
 9RP binds Grb2 <--> RG | kf=0.003, kr=0.05 | Grb2=85
10RG binds SOS <--> RGS | kf=0.01, kr=0.06 | SOS=34
11RGS is dissociated into RP and GS | kf=0.03, kr=4.5e-3
12GS is dissociated into Grb2 and SOS | kf=1.5e-3, kr=1e-4
13RP binds Shc <--> RSh | kf=0.09, kr=0.6 | Shc=150
14RSh is phosphorylated <--> RShP | kf=6, kr=0.06
15RShP is dissociated into ShP and RP | kf=0.3, kr=9e-4
16ShP is dephosphorylated --> Shc | V=1.7, K=340
17RShP binds Grb2 <--> RShG | kf=0.003, kr=0.1
18RShG is dissociated into RP and ShG | kf=0.3, kr=9e-4
19RShG binds SOS <--> RShGS | kf=0.01, kr=2.14e-2
20RShGS is dissociated into ShGS and RP | kf=0.12, kr=2.4e-4
21ShP binds Grb2 <--> ShG | kf=0.003, kr=0.1
22ShG binds SOS <--> ShGS | kf=0.03, kr=0.064
23ShGS is dissociated into ShP and GS | kf=0.1, kr=0.021
24RShP binds GS <--> RShGS | kf=0.009, kr=4.29e-2
25PLCgP is translocated to cytoskeletal or membrane structures <--> PLCgP_I | kf=1, kr=0.03
26
27# observable layer
28@obs Total_phosphorylated_Shc: u[RShP] + u[RShG] + u[RShGS] + u[ShP] + u[ShG] + u[ShGS]
29@obs Total_Grb2_coprecipitated_with_Shc: u[RShG] + u[ShG] + u[RShGS] + u[ShGS]
30@obs Total_phosphorylated_Shc_bound_to_EGFR: u[RShP] + u[RShG] + u[RShGS]
31@obs Total_Grb2_bound_to_EGFR: u[RG] + u[RGS] + u[RShG] + u[RShGS]
32@obs Total_SOS_bound_to_EGFR: u[RGS] + u[RShGS]
33@obs ShGS_complex: u[ShGS]
34@obs Total_phosphorylated_PLCg: u[RPLP] + u[PLCgP]
35
36# simulation layer
37@sim tspan: [0, 120]
38@sim condition EGF20nM: init[EGF] = 680
39@sim condition EGF2nM: init[EGF] = 68
40@sim condition Absence_PLCgP_transloc: init[EGF] = 680; p[kf25] = 0; p[kr25] = 0

  2. Convert the text into an executable model

    $ python

    To display thermodynamic restrictions, set show_restrictions to True.

    >>> from pasmopy import Text2Model
>>> description = Text2Model("Kholodenko_JBC_1999.txt")
>>> description.convert(show_restrictions=True)
Model information
-----------------
25 reactions
23 species
50 parameters

Thermodynamic restrictions
--------------------------
{9, 12, 10, 11}
{15, 18, 21, 17}
{18, 22, 20, 19}
{17, 24, 12, 19}
{23, 24, 20, 15}
{23, 12, 22, 21}

    The output of Thermodynamic restrictions shows the cyclic pathways in the biochemical reaction network. These detailed balance relations require the product of the equilibrium constants along a cycle to be equal to 1, since at equilibrium the net flux through any cycle vanishes.

  3. Run simulation

    >>> from pasmopy import create_model, run_simulation
>>> model = create_model("Kholodenko_JBC_1999")
>>> run_simulation(model)

  4. Plot simulation results

    %matplotlib inline
import os
import matplotlib.pyplot as plt
import numpy as np


def plot_simulation_results(res):

   plt.figure(figsize=(9, 9))
   plt.rcParams['font.family'] = 'Arial'
   plt.rcParams['font.size'] = 12
   plt.rcParams['axes.linewidth'] = 1
   plt.rcParams['lines.linewidth'] = 2

   plt.subplots_adjust(wspace=0.5, hspace=0.4)

   plt.subplot(2, 2, 1)  # ----------------------------------------------------
   for obs_name, color in zip(
      ['Total_phosphorylated_Shc', 'Total_Grb2_coprecipitated_with_Shc'],
      ['g', 'm'],
   ):
      obs_idx = model.observables.index(obs_name)
      for j, condition in enumerate(['EGF20nM', 'EGF2nM']):
            plt.plot(
               model.problem.t,
               res[obs_idx, j],
               color=color,
               alpha=0.5 if condition == 'EGF2nM' else None,
            )
   plt.xlim(0, 120)
   plt.xticks([30*i for i in range(5)])
   plt.ylim(0, 150)
   plt.xlabel("TIME (s)")
   plt.ylabel("Protein concentrations (nM)")

   plt.subplot(2, 2, 2)  # ----------------------------------------------------
   for obs_name, color in zip(
      ['Total_phosphorylated_Shc_bound_to_EGFR', 'Total_Grb2_bound_to_EGFR'],
      ['g', 'm'],
   ):
      obs_idx = model.observables.index(obs_name)
      for j, condition in enumerate(['EGF20nM', 'EGF2nM']):
            plt.plot(
               model.problem.t,
               res[obs_idx, j],
               color=color,
               alpha=0.5 if condition == 'EGF2nM' else None,
            )
   plt.xlim(0, 120)
   plt.xticks([30*i for i in range(5)])
   plt.ylim(0, 25)
   plt.xlabel("TIME (s)")
   plt.ylabel("Protein concentrations (nM)")

   ax1=plt.subplot(2, 2, 3)  # ------------------------------------------------
   ax2 = ax1.twinx()
   for j, condition in enumerate(['EGF20nM', 'EGF2nM']):
      ax1.plot(
            model.problem.t,
            res[model.observables.index('Total_SOS_bound_to_EGFR'), j],
            color='g',
            alpha=0.5 if condition == 'EGF2nM' else None,
      )
      ax2.plot(
            model.problem.t,
            res[model.observables.index('ShGS_complex'), j],
            color='m',
            alpha=0.5 if condition == 'EGF2nM' else None,
      )
   ax1.set_xlim(0, 120)
   ax1.set_xticks([30*i for i in range(5)])
   ax1.set_xlabel("TIME (s)")
   ax1.set_ylim(0, 8)
   ax2.set_ylim(0, 30)
   ax1.set_ylabel("SOS bound to EGFR (nM)")
   ax2.set_ylabel("Concentration of Sh-G-S (nM)")

   ax1=plt.subplot(2, 2, 4)  # ------------------------------------------------
   ax2 = ax1.twinx()
   obs_idx = model.observables.index('Total_phosphorylated_PLCg')
   ax1.plot(
      model.problem.t,
      res[obs_idx, model.problem.conditions.index('EGF20nM')],
      'g',
   )
   ax1.plot(
      model.problem.t,
      res[obs_idx, model.problem.conditions.index('EGF2nM')],
      'g',
      alpha=0.5,
   )
   ax2.plot(
      model.problem.t,
      res[obs_idx, model.problem.conditions.index('Absence_PLCgP_transloc')],
      'g--',
   )
   ax1.set_xlim(0, 120)
   ax1.set_xticks([30*i for i in range(5)])
   ax1.set_ylim(0, 15)
   ax1.set_yticks([5*i for i in range(4)])
   ax1.set_xlabel("TIME (s)")
   ax1.set_ylabel("Total Phosphorylated PLCγ (nM)")
   ax2.set_ylim(0, 105)
   ax2.set_yticks([30*i for i in range(4)])

   plt.show()


if __name__ == '__main__':
   res = np.load(os.path.join(model.path, "simulation_data", "simulations_original.npy"))
   plot_simulation_results(res)
    _images/EGF_signaling_sim.png

c-Fos expression dynamics

Please refer to https://biomass-core.readthedocs.io/en/latest/tutorial/cfos.html.