SEDML_Tellurium
Tellurium is a Python platform to model, simulate, and analyse biochemical systems [1]. It includes a number of pre-installed Python libraries, plugins and tools for biological modelling and can be used via any Python frontend.
Since the course was already set up in the Google environment, we used Google Colaboratory (Colab), a cloud-based Jupyter notebook service for collaborative Python programming.
Guided by examples in the Tellurium documentation, we built Colab notebooks to reproduce several experiments by Bachmann et al. [2], described in more detail below.
- open Colab notebook (Example)
- install + import tellurium and further required packages (os, phrasedml)
# install Tellurium
!pip3 install tellurium
# import required packages
import os
import tellurium as te
import phrasedml
- mount Google Drive and set working directory (corresponding to folder structure of COMBINE archive)
# mount Google Drive if necessary
from google.colab import drive
drive.mount('/content/drive')
workingDir = _directory where to store the SED-ML file_
os.chdir(workingDir)
- load SBML model from ../models (Bachmann2011.xml) and convert to string (note: models can also be loaded via URLs; however, for the purpose of the COMBINE archive we decided to load it from the SBML file stored locally)
# retrieve SBML model from model folder and set as reference for phraSEDML
sbml = "../model/Bachmann2011.xml"
ant_str = te.sbmlToAntimony(sbml)
sbml_str = te.antimonyToSBML(ant_str)
phrasedml.setReferencedSBML('../model/Bachmann2011.xml', sbml_str)
- create phraSEDML [3] string as template for SED-ML file - this is where the contents of the SED-ML (model +/- modifications, simulations, tasks, plots, reports) are specified:
- model: specify model and any modifications (e.g. changes to parameters)
- simulation: specify time course simulation (start time, end time, intervals)
- task: specify tasks, i.e. running a simulation on a model
- plot: specify plot(s) as output of the simulation
- report: specify parameters to include in model data
# create phraSEDML string as template for SED-ML
phrasedml_str = """
model1 = model "{}"
model2 = model model1 with SOCS3oe=1
model3 = model model1 with CISoe=1
sim1 = simulate uniform(0, 80, 1000)
task1 = run sim1 on model1
task2 = run sim1 on model2
task3 = run sim1 on model3
plot "Figure 5A Nuclear pSTAT5" task1.time vs task1.npSTAT5, task2.npSTAT5, task3.npSTAT5
report "Report Figure 5A" task1.time, task1.npSTAT5, task2.npSTAT5, task3.npSTAT5
""".format('../model/Bachmann2011.xml')
- convert phraSEDML to SED-ML and write to indicated file in ../experiment
# convert to SED-ML and create SED-ML file in working directory
sedml_str = phrasedml.convertString(phrasedml_str)
if sedml_str == None:
raise RuntimeError(phrasedml.getLastError())
print(sedml_str)
with open(os.path.join(workingDir, 'Fig5.sedml'), 'wb') as f:
f.write(sedml_str.encode('utf-8'))
- run the SED-ML file using the referenced model file to check simulation
# run the SED-ML file based on the referenced SBML file in ../models/
workingDir = _directory where to store the SED-ML file_
os.chdir(workingDir)
te.executeSEDML(sedml_str, workingDir)
- validate SED-ML file in SED-ML WebTools (described here) and download amended SED-ML file from SED-ML WebTools
- manually change IDs of plots and reports to include a reference to the figure (e.g. Fig5_npSTAT5, Fig5_Report) - otherwise, when simulating the whole COMBINE archive in Tellurium, the outputs are overwritten by the next figure
Figure 3A and 3C show model simulations for various parameters (pJAK2, pEPOR, tSTAT5, pSTAT5, CIS, SOCS3, CISRNA, SOCS3RNA) in combination with experimental data used for model calibration.
- The model could be used without any modifications (model1).
- A time course simulation from 0-250 min (sim1) was defined and applied to the model, generating task1.
- Plots and reports for each molecule of interest were defined based on the plots presented in the Bachmann paper.
- The plots for pJAK2, pEPOR and tSTAT5 could not be reproduced since these parameters are not included in the model.
Example plot showing time course of pSTAT5
Figure 4 shows model simulations for various parameters (SOCS3, CIS, pJAK2, pEPOR, SHP1, pSTAT5) in different experimental conditions (wildtype, SOCS3 overexpression, CIS overexpression, SHP1 overexpression).
- To simulate the wildtype condition, the model could be used without any modifications (model1).
- To simulate the overexpression of SOCS3 and CIS, two new models were created based on model1, where the parameters SOCS3oe and CISoe were set to 1 to create model2 and model3, respectively. In addition, the initial concentrations for SOCS3 and CIS were roughly set to maximum values for each molecule based on the plots generated.
- The overexpression of SHP1 could not be reproduced since the corresponding parameter (SHP1oe) was not available in the model (nor in any other Bachmann model files).
- A time course simulation from 0-120 min was defined (sim1) and applied to each model, generating task1-3.
- Plots and reports for each molecule of interest were defined based on the plots presented in the Bachmann paper.
- The plots for pJAK2 and pEPOR could not be reproduced since these parameters are not included in the model.
Example plot showing pSTAT5 concentrations under wildtype conditions (task1 - blue line), with SOCS3 overexpression (left, task2 - orange line) or with CIS overexpression (right, task3 - orange line)
Figure 5A shows model simulations for npSTAT5 under different experimental conditions (wildtype, SOCS3 overexpression, CIS overexpression) at Epo levels of 10^−6.78 U/cell.
- To simulate the wildtype condition, the model could be used without any modifications (model1).
- To simulate the overexpression of SOCS3 and CIS, two new models were created based on model1, where the parameters SOCS3oe and CISoe were set to 1 to create model2 and model3, respectively. In addition, the initial concentrations for SOCS3 and CIS were roughly set to maximum values for each molecule based on the plots generated.
- A time course simulation from 0-80 min was defined (sim1) and applied to each model, generating task1-3.
- Plots and reports for npSTAT5 were defined based on the plots presented in the Bachmann paper.
- It was not clear whether and how to modify Epo levels for this experiment. We tried to set
epo_level=1.2E-7as indicated in the supplementary material but got similar results to using the model as is (not shown), which is not identical to the published figure.
Example plot showing npSTAT5 concentrations under wildtype conditions (task1 - blue line), with SOCS3 overexpression (task2 - orange line) or with CIS overexpression (task3 - green line)
Figure 6A (bottom) shows model simulations for npSTAT5 for two exemplary Epo concentrations (I - basal, II - acute) in different experimental conditions (wildtype, SOCS3 knockout, CIS knockout, SOCS3/CIS double knockout). We could not reproduce the experiment at this stage, because we were not sure how to create a knockout simulation or how to simulated changing Epo levels.
- The model could be used without any modifications (model1).
- A time course simulation from 0-250 min (sim1) was defined and applied to the model, generating task1.
- Plots and reports for each molecule of interest were defined based on the plots presented in the Bachmann paper.
- We did not attempt to simulate the "extrinsic noise" described in the Supplementary Information, page 12.
Example plot showing pSTAT5 dynamics
- The plots for pJAK2, pEPOR and tSTAT5 (e.g., Fig3, Fig4) could not be generated since these parameters are missing from the model itself. This is most likely due to the fact that, for instance, pEPOR is composed of several model parameters - p1EpoRpJAK2, p2EpoRpJAK2, p12EpoRpJAK2. Information on how to derive these observables from existing model parameters is included in the supplementary material (e.g., Supplementary Information, page 21). However, given the time frame we did not attempt to include this in the simulation at this stage.
- In contrast to CIS and SOCS3, the overexpression of SHP1 could not be reproduced since the corresponding parameter (SHP1oe) was not included in the model (nor in any other Bachmann model files available online).
- We did not modify the model in order to simulate different level of Epo as described in the paper. While the corresponding figure legends and supplementary material contain some information on Epo concentrations used for these experiments, it was not clear ho to simulate this in the model.
- Many figures (e.g. Figure 4) contain experimental data in addition to model simulations. Although this data is referenced as supplementary material in the online version of the paper, we did not include this in our plots.
- The ranges of the y axes do not correspond to the figures in the paper. This may be due to the fact that the paper shows model outputs based on observation functions (Supplementary Information, page 21), but this information is not included in the model itself.
[1] Choi, K. et al. Tellurium: An extensible python-based modeling environment for systems and synthetic biology. Biosystems 171, 74-79 (2018), https://doi.org/10.1101/054601
[2] Bachmann, J. et al. Division of labor by dual feedback regulators controls JAK2/STAT5 signaling over broad ligand range. Molecular Systems Biology 7, 516 (2011). https://doi.org/10.1038/msb.2011.50
[3] Choi, K.et al. phraSED-ML: A paraphrased, human-readable adaptation of SED-ML. J Bioinform Comput Biol 14, 1650035 (2016), https://doi:10.1142/S0219720016500359