Tellurium notebooks-An environment for reproducible dynamical modeling in systems biology

Abstract

The considerable difficulty encountered in reproducing the results of published dynamical models limits validation, exploration and reuse of this increasingly large biomedical research resource. To address this problem, we have developed Tellurium Notebook, a software system for model authoring, simulation, and teaching that facilitates building reproducible dynamical models and reusing models by 1) providing a notebook environment which allows models, Python code, and narrative to be intermixed, 2) supporting the COMBINE archive format during model development for capturing model information in an exchangeable format and 3) enabling users to easily simulate and edit public COMBINE-compliant models from public repositories to facilitate studying model dynamics, variants and test cases. Tellurium Notebook, a Python-based Jupyter-like environment, is designed to seamlessly inter-operate with these community standards by automating conversion between COMBINE standards formulations and corresponding in-line, human-readable representations. Thus, Tellurium brings to systems biology the strategy used by other literate notebook systems such as Mathematica. These capabilities allow users to edit every aspect of the standards-compliant models and simulations, run the simulations in-line, and re-export to standard formats. We provide several use cases illustrating the advantages of our approach and how it allows development and reuse of models without requiring technical knowledge of standards. Adoption of Tellurium should accelerate model development, reproducibility and reuse.

Fig 1. A demonstration of Tellurium’s SBO syntax.

This figure shows a model of mitotic exit in budding yeast [38] available via the Biomodels repository entry BIOMD0000000370 [39]. When the SBML for this model is imported, Tellurium automatically extracts SBO identifiers for species, reactions, compartments, etc. and embeds the identifiers in the Antimony code. These identifiers point to specific physical, biological, or mathematical entities recorded in the ontology. For example, the first identifier, SBO:0000290, refers to a compartment in physical space [40]. Other identifiers refer to polypeptide chains (SBO:0000252) and protein complexes (SBO:0000297). These identifiers appear inline in the notebook cell and can be edited by the user. The right panel shows the transient response for this model from 0 to 120 minutes. This example is included with the Tellurium notebook viewer version 2.0.14 and later (File→Open Example Notebook→Mitotic Exit (Vinod)).

Fig 2. A comparison of Tellurium’s human–readable representation of a COMBINE archive shown in a Tellurium notebook (A) and excerpts from the equivalent SBML (B) and SED–ML (C) encodings.

Tellurium’s in-line OMEX format contains human–readable representations of both SBML and SED–ML (A). Here, SBML is represented by Antimony code (with the definition of a single reaction in blue) and PhraSEDML (in red). (B) shows the SBML encoding for a single reaction. The single–line human–readable form of this reaction is highlighted in part (A) for comparison. The components of the Antimony syntax are as follows: R23 is the reaction label, the reactant $UbE, with a dollar sign indicating a boundary species, a => symbol, which indicates an irreversible reaction (reversible reactions can be indicated with ->), the product UbE_star, and the kinetic law comprised of everything following the semicolon. Using the SBML encoding, it is difficult to modify the reaction stoichiometry or kinetic law, whereas this task is easy in Tellurium. Finally, (C) shows the SED–ML encoding corresponding to the human–readable simulation portion of this COMBINE archive. The simulation portion performs the following functions: first, two SBML models are instantiated with different sets of parameters (in the original publicatuion [44], the authors provided one set of parameters for oocyte extract and a different set for intact embryos). Second, a timecourse simulation with an adaptive step size is attained by setting the variable_step_size property of the simulation as well as appropriate tolerance values. Finally, the simulation is run with the two different model instantiations and plotted with representative state variables (active MPF, doubly phosphorylated/inactive MPF, and total cyclin) to show the behavior of the two parameter sets.

Fig 3. A demo of an SBML/SED–ML encoding contained in a COMBINE archive showing two useful features of the encoding: Multiple parameter sets (A) and post–processing (B).

(A) Transient responses of M phase control [44]. This model was published with two parameter sets. One set is based on measurements from Xenopus oocyte extracts (top) whereas the other is based on measuremetns from intact embryos (bottom). (B) Phase portraits of representative state variables in the model. These variables are chosen after [44] and are as follows: total cyclin, doubly–phosphorylated MPF (PPMPF, the predominant inactive form of MPF [44]), active MPF, and time. Each pair of variables is plotted in this matrix. Y–axis variables are indicated in the rows of the plot and x–axis variables are indicated in the columns. The title of each subplot is given in terms of x vs y, e.g., the top left subplot shows total cyclin on the x–axis vs active MPF on the y–axis. Phase portraits can show transients (such as the initial response of total cyclin in the upper left corner in blue, which starts at 100 and decreases to its normal range) as well as limit cycles (exhibited by all three phase portraits in the upper part of (B)). The slope of a given region of the phase portrait is useful for showing the relative rate of change of two quantities. The green and orange curves show regions where one quantity changes rapidly with respect to another. These regions correspond to the rapid rise in active MPF due to positive feedback from MPF to its own self–activation, and the subsequent falloff of total cyclin due to cyclin degradation via a ubiquitin pathway activated by MPF. The plot in part (B) is derived from the parameter set for oocyte extract, corresponding to the top plot in part (A).

Fig 4. Using Tellurium to reproduce model variants in [47] and repackage as a COMBINE archive.

To demonstrate the use of COMBINE archives for encoding model variants, we began with a COMBINE archive describing a single variant of this model without String synthesis or degradation [48], which reproduces Fig 1B of [47] (plots B and C here). We then used Tellurium to add a variant describing String degradation, which reproduces Fig 3 of [47] (plots D through G here). Panel (A) shows the inline OMEX cell with the Antimony code elided (it would belong in the model Model_generated_by_BIOCHAM…end block, where the ellipsis is currently shown). Everything after the end instruction is thus PhraSEDML. Plots B and D show the transient response of the cytoplasmic compartment of the model. Plots C and E show the nuclear compartment (defined as the spatial region around the mitotic spindle). Plot F shows the levels of total String and its phosphorylated state. Plot G shows the level of String mRNA and protein factor X, which degrades String mRNA. Note the y–axis scale on plot G was manually adjusted to show the mRNA dynamics. The subplots in this figure intentionally have different durations, after Calzone et al [47]. The model in [47] was authored using BIOCHAM [53]. Our model reproductions that reproduce these plots are available as a COMBINE archive [54].

Fig 5. Comparison of logarithmic (bottom row of plots) vs. linear (top row of plots) plotting of the transient response in Fig 4B, 4C and 4G, along with an excerpt of the relevant PhraSEDML code for plotting on logarithmic axes.

Fig 4G, which is a plot of String mRNA and hypothetical factor X, which degrades String mRNA, exhibits a large dynamic range. Logarithmic plotting helps visualize the dynamic range of these quantities. This is achieved in PhraSEDML by wrapping the quantities for x and y axes inside a log10 operation.

Fig 6. Testing the shift in regulatory mechanism of mitotic oscillations.

To verify the observation [47] that the number of mitotic divisions in the Drosophila embryo is governed by a shift from negative to positive feedback, we first removed all discrete events and introduced the variable C such that N = 1.95^C. We then compared the limit cycles produced by this eventless model (left) with those produced by a variant with attenuated positive feedback from the regulators Wee1 and String (right). Attenuation was achieved by decreasing the rates of the phosphorylation and dephosphorylation of Wee1 and String. The original model exhibits stable limit cycle oscillations for both early cycles (C), which are putatively dominated by negative feedback, and late cycles (E), which are putatively dominated by positive feedback. The attenuated model only exhibits stable oscillations at early cycles (D), suggesting that positive feedback does indeed play a role in late cycle oscillations (F). Our model reuse and modification study is available as a COMBINE archive that reproduces the figure shown and facilitates further modification and reuse [55].