Signaling code¶
In the following we will implement and analyse the simple signaling motives from the lecture - phosphorylation cycle - two-component system
In addition we will have a look at a real world example. A MAP kinase model.
We will use tellurium
to encode the model, i.e., the ODEs and for
simulation. For documentation see
https://tellurium.readthedocs.io/en/latest/notebooks.html#modeling-case-studies
%matplotlib inline
1. Building ODE models¶
1.1 Modeling in a nutshell¶
Basic steps: - write down the players (species Si) in your system - write down the processes vi which influence your species and the stoichiometry (number of species consumed/produced via the process) - write down mathematical equations for the processes - equations define the rates of the processes, defining how fast the process changes the species
Processes - require stoichiometry and rate law (kinetics)
System of ordinary differential equations (ODEs)
The result is a list of mathematical equations (differential equations) * can be solved with numerical methods * results in time courses, i.e., the time-dependent
\[\frac{dS1}{dt} = -v1 = -k1 \cdot S1\]
\[\frac{dS2}{dt} = v1 = k1 \cdot S1\]
Main types of equations * Mass-Action kinetics * Michaelis-Menten Kinetics (irreversible & reversible) * Hill equations (cooperativity) * (Non-)Competitive Inhibition
1.2 Model definition¶
import tellurium as te
import pandas as pd
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
<ipython-input-2-8d1c52be6743> in <module>
----> 1 import tellurium as te
2 import pandas as pd
ModuleNotFoundError: No module named 'tellurium'
# model definition
model = te.loada(f"""
model example1()
// compartments and species
species S1, S2;
// processes
v1: S1 -> S2; k1*S1;
// initial conditions
S1 = 10.0;
// parameters
k1 = 0.1;
end
""")
print(te.sbmlToAntimony(model.getSBML()))
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-3-b51a36b51f47> in <module>
1 # model definition
----> 2 model = te.loada(f"""
3 model example1()
4 // compartments and species
5 species S1, S2;
NameError: name 'te' is not defined
1.3 simulation¶
model.reset()
s = model.simulate(start=0, end=100, steps=100)
model.plot(s)
s = pd.DataFrame(s, columns=s.colnames)
display(s)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-4-748f20177d7d> in <module>
----> 1 model.reset()
2 s = model.simulate(start=0, end=100, steps=100)
3 model.plot(s)
4 s = pd.DataFrame(s, columns=s.colnames)
5 display(s)
NameError: name 'model' is not defined
2. Phosphorylation cycle¶
One of the most important signaling motives are protein phosphorylation cycles, i.e., a post-translational modification of a protein in which an amino acid residue is phosphorylated by a protein kinase, and dephosphorylated by a protein phosphatase.
Simple models follow mass-action kinetics. For example the following simple phosphorylation cycle where the kinase activity represents the signal \(S\), and the activity of the phosphatase is assumed to be constant (and included in the rate constant \(k_2\))
The system exhibits mass conservation \(R_p + R = R^T\), where \(R^T\) denotes the amount of total protein.
2.1 Build model phosphorylation cycle¶
import tellurium as te
import pandas as pd
from IPython.display import display, HTML
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
<ipython-input-5-68a3a525987d> in <module>
----> 1 import tellurium as te
2 import pandas as pd
3 from IPython.display import display, HTML
ModuleNotFoundError: No module named 'tellurium'
model_pcycle = te.loada(f"""
model pcycle()
// compartments and species
species R, Rp;
// reactions
v1: R -> Rp; k1*S*R;
v2: Rp -> R; k2*Rp;
// initial conditions
R = 10.0; Rp = 0.0;
// parameters
k1 = 0.1; k2 = 0.1;
S = 1.0;
end
""")
# print(te.sbmlToAntimony(model_pcycle.getSBML()))
s = model_pcycle.simulate(start=0, end=50, steps=200)
model_pcycle.plot(s)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-6-0e52bacfcfb1> in <module>
----> 1 model_pcycle = te.loada(f"""
2 model pcycle()
3 // compartments and species
4 species R, Rp;
5
NameError: name 'te' is not defined
2.2 Check mass balance¶
- to check the mass balance we add an assignment rule to the model
- to be able to access the variables in the results we
Rt
to the selections
model_pcycle = te.loada(f"""
model pcycle()
// compartments and species
species R, Rp;
// reactions
v1: R -> Rp; k1*S*R;
v2: Rp -> R; k2*Rp;
// initial conditions
R = 10.0; Rp = 0.0;
// parameters
k1 = 0.1; k2 = 0.1;
S = 1.0;
// rules
Rt := R + Rp
end
""")
# print(te.sbmlToAntimony(model_pcycle.getSBML()))
selections = ["time", "R", "Rp", "S", "Rt"]
s = model_pcycle.simulate(start=0, end=50, steps=200, selections=selections)
model_pcycle.plot(s)
s = pd.DataFrame(s, columns=s.colnames)
display(s)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-7-a4374d7c48b1> in <module>
----> 1 model_pcycle = te.loada(f"""
2 model pcycle()
3 // compartments and species
4 species R, Rp;
5
NameError: name 'te' is not defined
2.3 Parameter scan¶
A typical task in a simulation study is to scan parameters and check the influence of parameters on the model.
In the following we will change the parameter of the kinase k1
import numpy as np
from matplotlib import pyplot as plt
model.resetToOrigin()
k1_vec = np.linspace(0, 1, num=21)
# k1_vec = np.logspace(-10, 1, num=21)
# run parameter scan
results = []
for k1 in k1_vec:
# print(k1)
model_pcycle.reset()
model_pcycle["k1"] = k1
s = model_pcycle.simulate(start=0, end=50, steps=200)
s = pd.DataFrame(s, columns=s.colnames)
results.append(s)
# create plot
f, ax = plt.subplots()
kwargs = {"alpha": 0.8}
for k, s in enumerate(results):
ax.plot(s.time, s.R, color="blue", **kwargs)
ax.plot(s.time, s.Rp, color="darkorange", **kwargs)
ax.plot(s.time, s.Rt, color="red", **kwargs)
ax.set_xlabel("time")
ax.set_ylabel("amount")
plt.show()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-9-bc393d3f761e> in <module>
----> 1 model.resetToOrigin()
2 k1_vec = np.linspace(0, 1, num=21)
3 # k1_vec = np.logspace(-10, 1, num=21)
4
5 # run parameter scan
NameError: name 'model' is not defined
Important things - numerical simulations (with floating point values and certain accuracy) - steady state concentration depend on the \(k_1\) parameter - the higher the phosphorylation rate the faster the steady state is reached - mass balance is conserved for changing of parameters
Often one is interested in changes in steady state (or other system properties with parameters). - for instance how do the steady state values of R and Rp change with \(k_1\)
# collect steady state values
R_ss = np.zeros_like(k1_vec)
Rp_ss = np.zeros_like(k1_vec)
Rt_ss = np.zeros_like(k1_vec)
for k, s in enumerate(results):
R_ss[k] = s.R.values[-1]
Rp_ss[k] = s.Rp.values[-1]
Rt_ss[k] = s.Rt.values[-1]
f, ax = plt.subplots()
kwargs = {"marker": "s"}
ax.plot(k1_vec, R_ss, color="blue", **kwargs)
ax.plot(k1_vec, Rp_ss, color="darkorange", **kwargs)
ax.plot(k1_vec, Rt_ss, color="red", **kwargs)
ax.set_xlabel("k1")
ax.set_ylabel("amount")
# ax.set_xscale("log")
plt.show()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-10-2eec67600bbe> in <module>
1 # collect steady state values
----> 2 R_ss = np.zeros_like(k1_vec)
3 Rp_ss = np.zeros_like(k1_vec)
4 Rt_ss = np.zeros_like(k1_vec)
5
NameError: name 'k1_vec' is not defined
2.4 Steady state (dependency on signal)¶
Steady state of the system is given by
Note that the dependence on the kinase activity (signal) is hyperbolic, whereas the dependence on total protein is linear.
In the following we compare the analytical solution against the numerical solution (vs actual steady state simulation).
# make sure the model is reset
model_pcycle.resetToOrigin()
# print model (to check parameter values)
print(te.sbmlToAntimony(model_pcycle.getSBML()))
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-11-e479a66a6a9f> in <module>
1 # make sure the model is reset
----> 2 model_pcycle.resetToOrigin()
3 # print model (to check parameter values)
4 print(te.sbmlToAntimony(model_pcycle.getSBML()))
NameError: name 'model_pcycle' is not defined
S_vec = np.linspace(0, 10, num=21)
# S_vec = np.logspace(-3, 3, num=100)
R_ss = np.zeros_like(S_vec)
Rp_ss = np.zeros_like(S_vec)
Rt_ss = np.zeros_like(S_vec)
# run parameter scan
results = []
for k, S in enumerate(S_vec):
model_pcycle.reset()
model_pcycle["S"] = S
s = model_pcycle.simulate(start=0, end=50, steps=100)
s = pd.DataFrame(s, columns=s.colnames)
# store results
results.append(s)
# collect steady state values
R_ss[k] = s.R.values[-1]
Rp_ss[k] = s.Rp.values[-1]
Rt_ss[k] = s.Rt.values[-1]
f, ax = plt.subplots()
kwargs = {"marker": "s", "alpha": 0.8}
ax.plot(S_vec, R_ss, color="blue", label="R_ss", **kwargs)
ax.plot(S_vec, Rp_ss, color="darkorange", label="Rp_ss", **kwargs)
ax.plot(S_vec, Rt_ss, color="red", label="Rt_ss", **kwargs)
# analytical solution
Rp_ssf = model_pcycle.Rt * (1-S_vec/(S_vec + model_pcycle.k2/model_pcycle.k1))
ax.plot(S_vec, Rp_ssf, "-k", label="Rp_ss exact", linewidth="2")
ax.legend()
ax.set_xlabel("S")
ax.set_ylabel("amount")
# ax.set_xscale("log")
plt.show()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-12-4f1ffbc467d7> in <module>
9 results = []
10 for k, S in enumerate(S_vec):
---> 11 model_pcycle.reset()
12 model_pcycle["S"] = S
13 s = model_pcycle.simulate(start=0, end=50, steps=100)
NameError: name 'model_pcycle' is not defined
2.5 Steady state sensitivity (signal)¶
Steady state of the system is given by
Dependency on signal
S_vec = np.linspace(0, 10, num=21)
# S_vec = np.logspace(-3, 3, num=100)
dR_dRt = np.zeros_like(S_vec)
dRp_dRt = np.zeros_like(S_vec)
ln_dR_dRt = np.zeros_like(S_vec)
ln_dRp_dRt = np.zeros_like(S_vec)
# change
delta = 0.1 # increase parameter 1%
# run parameter scan
results = []
for k, S in enumerate(S_vec):
# complete reset
model_pcycle.resetToOrigin()
model_pcycle["S"] = S
s = model_pcycle.simulate(start=0, end=50, steps=100)
s = pd.DataFrame(s, columns=s.colnames)
# collect steady state values
R_ss = s.R.values[-1]
Rp_ss = s.Rp.values[-1]
# calculate sensitivity (right sided)
model_pcycle.resetToOrigin()
# small parameter change delta S (only one-sided sensitivity)
model_pcycle["S"] = S * (1 + delta)
s = model_pcycle.simulate(start=0, end=50, steps=100)
s = pd.DataFrame(s, columns=s.colnames)
# collect steady state values
R_ss_delta = s.R.values[-1]
Rp_ss_delta = s.Rp.values[-1]
# dR/dRt ~ (R(S+delta S)-R(Rt))/(delta*S)
dR_dRt[k] = (R_ss_delta-R_ss)/(S*delta)
dRp_dRt[k] = (Rp_ss_delta-Rp_ss)/(S*delta)
# log sensitivities
ln_dR_dRt[k] = S/R_ss * dR_dRt[k]
ln_dRp_dRt[k] = S/Rp_ss * dRp_dRt[k]
# plot the sensitivities
f, (ax1, ax2) = plt.subplots(ncols=2, nrows=1, figsize=(10, 5))
f.subplots_adjust(wspace=0.3)
kwargs = {"marker": "s", "alpha": 0.8}
ax1.set_ylabel("dR(p)/dS")
ax1.plot(S_vec, dR_dRt, color="blue", label="dR/dS", **kwargs)
ax1.plot(S_vec, dRp_dRt, color="darkorange", label="dRp/dS", **kwargs)
ax2.set_ylabel("dlnR(p)/dlnS")
ax2.plot(S_vec, ln_dR_dRt, color="blue", label="dlnR/dlnS", **kwargs)
ax2.plot(S_vec, ln_dRp_dRt, color="darkorange", label="dRp/dS", **kwargs)
# analytical solutions
Rp_ssf = model_pcycle.k2/model_pcycle.k1 * model_pcycle.Rt/((S_vec + (model_pcycle.k2/model_pcycle.k1))**2)
ax1.plot(S_vec, Rp_ssf, "-k", label="dRp_dRt exact", linewidth="2")
Rp_ssf = model_pcycle.k2/model_pcycle.k1/(S_vec + (model_pcycle.k2/model_pcycle.k1))
ax2.plot(S_vec, Rp_ssf, "-k", label="dRp_dRt exact", linewidth="2")
for ax in (ax1, ax2):
ax.legend()
ax.set_xlabel("S")
# ax.set_xscale("log")
plt.show()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-13-b660ac94f742> in <module>
14 for k, S in enumerate(S_vec):
15 # complete reset
---> 16 model_pcycle.resetToOrigin()
17 model_pcycle["S"] = S
18 s = model_pcycle.simulate(start=0, end=50, steps=100)
NameError: name 'model_pcycle' is not defined
3. Two-component system¶
Two-component signaling systems typically consist of - (membrane-bound) histidine kinase (HK) that senses a specific environmental stimulus (typically homodimeric transmembrane proteins containing a histidine phosphotransfer domain and an ATP binding domain)
- mass action kinetics
- \(H\): histidine kinase
- \(R\): response regulator
mass conservation: \(H + H_p = H^T\) and \(R + R_p = R^T\)
model_twocomp = te.loada(f"""
model pcycle()
v1: H -> Hp; k1*S*H;
v2: Hp + R -> H + Rp; k2*R*Hp;
v3: Rp -> R; k3 * Rp;
// species
species H, Hp, R, Rp;
// initial values
H = 10.0; Hp = 0.0;
R = 10.0; Rp = 0.0;
// parameters
S = 1.0;
k1 = 0.1; k2 = 0.1; k3 = 0.1;
// rules
Ht := H + Hp
Rt := R + Rp
end
""")
model_twocomp.selections = ["time", "H", "Hp", "R", "Rp"]
model_twocomp.simulate(start=0, end=100, steps=200)
model_twocomp.plot()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-14-d8866ae53e1f> in <module>
----> 1 model_twocomp = te.loada(f"""
2 model pcycle()
3 v1: H -> Hp; k1*S*H;
4 v2: Hp + R -> H + Rp; k2*R*Hp;
5 v3: Rp -> R; k3 * Rp;
NameError: name 'te' is not defined
Steady state solution can be calculated, but lengthy quadratic equation.
# S_vec = np.linspace(0, 10, num=21)
S_vec = np.logspace(-3, 3, num=100)
H_ss = np.zeros_like(S_vec)
Hp_ss = np.zeros_like(S_vec)
R_ss = np.zeros_like(S_vec)
Rp_ss = np.zeros_like(S_vec)
# run parameter scan
results = []
for k, S in enumerate(S_vec):
model_twocomp.reset()
model_twocomp["S"] = S
s = model_twocomp.simulate(start=0, end=500, steps=100)
s = pd.DataFrame(s, columns=s.colnames)
# store results
results.append(s)
# collect steady state values
R_ss[k] = s.R.values[-1]
Rp_ss[k] = s.Rp.values[-1]
H_ss[k] = s.H.values[-1]
Hp_ss[k] = s.Hp.values[-1]
f, ax = plt.subplots()
kwargs = {"marker": "s", "alpha": 0.8}
ax.plot(S_vec, H_ss, color="blue", label="H_ss", **kwargs)
ax.plot(S_vec, Hp_ss, color="darkorange", label="Hp_ss", **kwargs)
ax.plot(S_vec, R_ss, color="darkgreen", label="R_ss", **kwargs)
ax.plot(S_vec, Rp_ss, color="red", label="Rp_ss", **kwargs)
ax.legend()
ax.set_xlabel("S")
ax.set_ylabel("amount")
ax.set_xscale("log")
plt.show()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-15-c00e474149a9> in <module>
10 results = []
11 for k, S in enumerate(S_vec):
---> 12 model_twocomp.reset()
13 model_twocomp["S"] = S
14 s = model_twocomp.simulate(start=0, end=500, steps=100)
NameError: name 'model_twocomp' is not defined
3.2 Perfect adaption¶
model_tcperfect = te.loada(f"""
model pcycle()
v1: H -> Hp; k1*S*H;
v2: Hp + R -> H + Rp; k2*R*Hp;
v3: Rp -> R; k3 * Rp * H;
// initial values
H = 10.0; Hp = 0.0;
R = 10.0; Rp = 0.0;
// parameters
S = 1.0;
k1 = 0.1; k2 = 0.1; k3 = 0.1;
// rules
Ht := H + Hp
Rt := R + Rp
end
""")
model_tcperfect.selections = ["time", "H", "Hp", "R", "Rp"]
model_tcperfect.simulate(start=0, end=50, steps=200)
model_tcperfect.plot()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-16-8ed5631fdc32> in <module>
----> 1 model_tcperfect = te.loada(f"""
2 model pcycle()
3 v1: H -> Hp; k1*S*H;
4 v2: Hp + R -> H + Rp; k2*R*Hp;
5 v3: Rp -> R; k3 * Rp * H;
NameError: name 'te' is not defined
She steady state solution for the response regulator is
The resulting expression is independent of the expression of the
proteins :math:R
and :math:H
. This is often termed perfect
adaption or integral feedback.
S_vec = np.linspace(0, 10, num=100)
# S_vec = np.logspace(-3, 3, num=100)
H_ss = np.zeros_like(S_vec)
Hp_ss = np.zeros_like(S_vec)
R_ss = np.zeros_like(S_vec)
Rp_ss = np.zeros_like(S_vec)
# run parameter scan
results = []
for k, S in enumerate(S_vec):
model_tcperfect.reset()
model_tcperfect["Hp"] = 20
model_tcperfect["Rp"] = 150
model_tcperfect["S"] = S
s = model_tcperfect.simulate(start=0, end=500, steps=100)
s = pd.DataFrame(s, columns=s.colnames)
# collect steady state values
R_ss[k] = s.R.values[-1]
Rp_ss[k] = s.Rp.values[-1]
H_ss[k] = s.H.values[-1]
Hp_ss[k] = s.Hp.values[-1]
f, ax = plt.subplots()
kwargs = {"marker": "s", "alpha": 0.8}
ax.plot(S_vec, H_ss, color="blue", label="H_ss", **kwargs)
ax.plot(S_vec, Hp_ss, color="darkorange", label="Hp_ss", **kwargs)
ax.plot(S_vec, R_ss, color="darkgreen", label="R_ss", **kwargs)
ax.plot(S_vec, Rp_ss, color="red", label="Rp_ss", **kwargs)
# analytical solution
Rp_ssf = model_tcperfect.k1/model_tcperfect.k3 * S_vec
ax.plot(S_vec, Rp_ssf, "-k", label="dRp_dRt exact", linewidth="2")
ax.legend()
ax.set_xlabel("S")
ax.set_ylabel("amount")
# ax.set_xscale("log")
plt.show()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-17-266b2cc18fe9> in <module>
10 results = []
11 for k, S in enumerate(S_vec):
---> 12 model_tcperfect.reset()
13 model_tcperfect["Hp"] = 20
14 model_tcperfect["Rp"] = 150
NameError: name 'model_tcperfect' is not defined
4. MAPK (real-world example)¶
Kholodenko2000 - Ultrasensitivity and negative feedback bring oscillations in MAPK cascade¶
Abstract:
Functional organization of signal transduction into protein phosphorylation cascades, such as the mitogen-activated protein kinase (MAPK) cascades, greatly enhances the sensitivity of cellular targets to external stimuli. The sensitivity increases multiplicatively with the number of cascade levels, so that a tiny change in a stimulus results in a large change in the response, the phenomenon referred to as ultrasensitivity. In a variety of cell types, the MAPK cascades are imbedded in long feedback loops, positive or negative, depending on whether the terminal kinase stimulates or inhibits the activation of the initial level. Here we demonstrate that a negative feedback loop combined with intrinsic ultrasensitivity of the MAPK cascade can bring about sustained oscillations in MAPK phosphorylation. Based on recent kinetic data on the MAPK cascades, we predict that the period of oscillations can range from minutes to hours. The phosphorylation level can vary between the base level and almost 100% of the total protein. The oscillations of the phosphorylation cascades and slow protein diffusion in the cytoplasm can lead to intracellular waves of phospho-proteins.
https://www.ebi.ac.uk/biomodels-main/BIOMD0000000010
import tellurium as te
# Load model from biomodels (may not work with https).
r = te.loadSBMLModel("https://www.ebi.ac.uk/biomodels-main/download?mid=BIOMD0000000010")
result = r.simulate(start=0, end=3000, steps=3000)
r.plot(result)
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
<ipython-input-18-9a85773b2df0> in <module>
----> 1 import tellurium as te
2
3 # Load model from biomodels (may not work with https).
4 r = te.loadSBMLModel("https://www.ebi.ac.uk/biomodels-main/download?mid=BIOMD0000000010")
5 result = r.simulate(start=0, end=3000, steps=3000)
ModuleNotFoundError: No module named 'tellurium'
print(te.sbmlToAntimony(r.getSBML()))
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-19-4cfabc9f32e0> in <module>
----> 1 print(te.sbmlToAntimony(r.getSBML()))
NameError: name 'te' is not defined