PyBioNetFit and the Biological Property Specification Language

Reusable rule-based cell cycle model explains compartment-resolved dynamics of 16 observables in RPE-1 cells

The mammalian cell cycle is regulated by a well-studied but complex biochemical reaction system. Computational models provide a particularly systematic and systemic description of the mechanisms governing mammalian cell cycle control. By combining both state-of-the-art multiplexed experimental methods and powerful computational tools, this work aims at improving on these models along four dimensions: model structure, validation data, validation methodology and model reusability. We developed a comprehensive model structure of the full cell cycle that qualitatively explains the behaviour of human retinal pigment epithelial-1 cells. To estimate the model parameters, time courses of eight cell cycle regulators in two compartments were reconstructed from single cell snapshot measurements. After optimisation with a parallel global optimisation metaheuristic we obtained excellent agreements between simulations and measurements. The PEtab specification of the optimisation problem facilitates reuse of model, data and/or optimisation results. Future perturbation experiments will improve parameter identifiability and allow for testing model predictive power. Such a predictive model may aid in drug discovery for cell cycle-related disorders.

Nonself RNA rewires IFN-β signaling: A mathematical model of the innate immune response

Type I interferons (IFNs) are key coordinators of the innate immune response to viral infection, which, through activation of the transcriptional regulators STAT1 and STAT2 (STAT1/2) in bystander cells, induce the expression of IFN-stimulated genes (ISGs). Here, we showed that in cells transfected with poly(I:C), an analog of viral RNA, the transcriptional activity of STAT1/2 was terminated because of depletion of the interferon-β (IFN-β) receptor, IFNAR. Activation of RNase L and PKR, products of two ISGs, not only hindered the replenishment of IFNAR but also suppressed negative regulators of IRF3 and NF-κB, consequently promoting IFNB transcription. We incorporated these findings into a mathematical model of innate immunity. By coupling signaling through the IRF3-NF-κB and STAT1/2 pathways with the activities of RNase L and PKR, the model explains how poly(I:C) switches the transcriptional program from being STAT1/2 induced to being IRF3 and NF-κB induced, which converts IFN-β-responding cells to IFN-β-secreting cells.

Epigenetic targeting of autophagy for cancer: DNA and RNA methylation

Autophagy, a crucial cellular mechanism responsible for degradation and recycling of intracellular components, is modulated by an intricate network of molecular signals. Its paradoxical involvement in oncogenesis, acting as both a tumor suppressor and promoter, has been underscored in recent studies. Central to this regulatory network are the epigenetic modifications of DNA and RNA methylation, notably the presence of N6-methyldeoxyadenosine (6mA) in genomic DNA and N6-methyladenosine (m6A) in eukaryotic mRNA. The 6mA modification in genomic DNA adds an extra dimension of epigenetic regulation, potentially impacting the transcriptional dynamics of genes linked to autophagy and, especially, cancer. Conversely, m6A modification, governed by methyltransferases and demethylases, influences mRNA stability, processing, and translation, affecting genes central to autophagic pathways. As we delve deeper into the complexities of autophagy regulation, the importance of these methylation modifications grows more evident. The interplay of 6mA, m6A, and autophagy points to a layered regulatory mechanism, illuminating cellular reactions to a range of conditions. This review delves into the nexus between DNA 6mA and RNA m6A methylation and their influence on autophagy in cancer contexts. By closely examining these epigenetic markers, we underscore their promise as therapeutic avenues, suggesting novel approaches for cancer intervention through autophagy modulation.

BioNetGMMFit: estimating parameters of a BioNetGen model from time-stamped snapshots of single cells

Mechanistic models are commonly employed to describe signaling and gene regulatory kinetics in single cells and cell populations. Recent advances in single-cell technologies have produced multidimensional datasets where snapshots of copy numbers (or abundances) of a large number of proteins and mRNA are measured across time in single cells. The availability of such datasets presents an attractive scenario where mechanistic models are validated against experiments, and estimated model parameters enable quantitative predictions of signaling or gene regulatory kinetics. To empower the systems biology community to easily estimate parameters accurately from multidimensional single-cell data, we have merged a widely used rule-based modeling software package BioNetGen, which provides a user-friendly way to code for mechanistic models describing biochemical reactions, and the recently introduced CyGMM, that uses cell-to-cell differences to improve parameter estimation for such networks, into a single software package: BioNetGMMFit. BioNetGMMFit provides parameter estimates of the model, supplied by the user in the BioNetGen markup language (BNGL), which yield the best fit for the observed single-cell, time-stamped data of cellular components. Furthermore, for more precise estimates, our software generates confidence intervals around each model parameter. BioNetGMMFit is capable of fitting datasets of increasing cell population sizes for any mechanistic model specified in the BioNetGen markup language. By streamlining the process of developing mechanistic models for large single-cell datasets, BioNetGMMFit provides an easily-accessible modeling framework designed for scale and the broader biochemical signaling community.

A plausible identifiable model of the canonical NF-κB signaling pathway

An overwhelming majority of mathematical models of regulatory pathways, including the intensively studied NF-κB pathway, remains non-identifiable, meaning that their parameters may not be determined by existing data. The existing NF-κB models that are capable of reproducing experimental data contain non-identifiable parameters, whereas simplified models with a smaller number of parameters exhibit dynamics that differs from that observed in experiments. Here, we reduced an existing model of the canonical NF-κB pathway by decreasing the number of equations from 15 to 6. The reduced model retains two negative feedback loops mediated by IκBα and A20, and in response to both tonic and pulsatile TNF stimulation exhibits dynamics that closely follow that of the original model. We carried out the sensitivity-based linear analysis and Monte Carlo-based analysis to demonstrate that the resulting model is both structurally and practically identifiable given measurements of 5 model variables from a simple TNF stimulation protocol. The reduced model is capable of reproducing different types of responses that are characteristic to regulatory motifs controlled by negative feedback loops: nearly-perfect adaptation as well as damped and sustained oscillations. It can serve as a building block of more comprehensive models of the immune response and cancer, where NF-κB plays a decisive role. Our approach, although may not be automatically generalized, suggests that models of other regulatory pathways can be transformed to identifiable, while retaining their dynamical features.

Verifiable biology

The formalization of biological systems using computational modelling approaches as an alternative to mathematical-based methods has recently received much interest because computational models provide a deeper mechanistic understanding of biological systems. In particular, formal verification, complementary approach to standard computational techniques such as simulation, is used to validate the system correctness and obtain critical information about system behaviour. In this study, we survey the most frequently used computational modelling approaches and formal verification techniques for computational biology. We compare a number of verification tools and software suites used to analyse biological systems and biochemical networks, and to verify a wide range of biological properties. For users who have no expertise in formal verification, we present a novel methodology that allows them to easily apply formal verification techniques to analyse their biological or biochemical system of interest.

Control of cell state transitions

Understanding cell state transitions and purposefully controlling them is a longstanding challenge in biology. Here we present cell state transition assessment and regulation (cSTAR), an approach for mapping cell states, modelling transitions between them and predicting targeted interventions to convert cell fate decisions. cSTAR uses omics data as input, classifies cell states, and develops a workflow that transforms the input data into mechanistic models that identify a core signalling network, which controls cell fate transitions by influencing whole-cell networks. By integrating signalling and phenotypic data, cSTAR models how cells manoeuvre in Waddington's landscape1 and make decisions about which cell fate to adopt. Notably, cSTAR devises interventions to control the movement of cells in Waddington's landscape. Testing cSTAR in a cellular model of differentiation and proliferation shows a high correlation between quantitative predictions and experimental data. Applying cSTAR to different types of perturbation and omics datasets, including single-cell data, demonstrates its flexibility and scalability and provides new biological insights. The ability of cSTAR to identify targeted perturbations that interconvert cell fates will enable designer approaches for manipulating cellular development pathways and mechanistically underpinned therapeutic interventions.

Combining hypothesis- and data-driven neuroscience modeling in FAIR workflows

Modeling in neuroscience occurs at the intersection of different points of view and approaches. Typically, hypothesis-driven modeling brings a question into focus so that a model is constructed to investigate a specific hypothesis about how the system works or why certain phenomena are observed. Data-driven modeling, on the other hand, follows a more unbiased approach, with model construction informed by the computationally intensive use of data. At the same time, researchers employ models at different biological scales and at different levels of abstraction. Combining these models while validating them against experimental data increases understanding of the multiscale brain. However, a lack of interoperability, transparency, and reusability of both models and the workflows used to construct them creates barriers for the integration of models representing different biological scales and built using different modeling philosophies. We argue that the same imperatives that drive resources and policy for data - such as the FAIR (Findable, Accessible, Interoperable, Reusable) principles - also support the integration of different modeling approaches. The FAIR principles require that data be shared in formats that are Findable, Accessible, Interoperable, and Reusable. Applying these principles to models and modeling workflows, as well as the data used to constrain and validate them, would allow researchers to find, reuse, question, validate, and extend published models, regardless of whether they are implemented phenomenologically or mechanistically, as a few equations or as a multiscale, hierarchical system. To illustrate these ideas, we use a classical synaptic plasticity model, the Bienenstock-Cooper-Munro rule, as an example due to its long history, different levels of abstraction, and implementation at many scales.

Implementation of a practical Markov chain Monte Carlo sampling algorithm in PyBioNetFit

SUMMARY

Bayesian inference in biological modeling commonly relies on Markov chain Monte Carlo (MCMC) sampling of a multidimensional and non-Gaussian posterior distribution that is not analytically tractable. Here, we present the implementation of a practical MCMC method in the open-source software package PyBioNetFit (PyBNF), which is designed to support parameterization of mathematical models for biological systems. The new MCMC method, am, incorporates an adaptive move proposal distribution. For warm starts, sampling can be initiated at a specified location in parameter space and with a multivariate Gaussian proposal distribution defined initially by a specified covariance matrix. Multiple chains can be generated in parallel using a computer cluster. We demonstrate that am can be used to successfully solve real-world Bayesian inference problems, including forecasting of new Coronavirus Disease 2019 case detection with Bayesian quantification of forecast uncertainty.

AVAILABILITY AND IMPLEMENTATION

PyBNF version 1.1.9, the first stable release with am, is available at PyPI and can be installed using the pip package-management system on platforms that have a working installation of Python 3. PyBNF relies on libRoadRunner and BioNetGen for simulations (e.g. numerical integration of ordinary differential equations defined in SBML or BNGL files) and Dask.Distributed for task scheduling on Linux computer clusters. The Python source code can be freely downloaded/cloned from GitHub and used and modified under terms of the BSD-3 license (https://github.com/lanl/pybnf). Online documentation covering installation/usage is available (https://pybnf.readthedocs.io/en/latest/). A tutorial video is available on YouTube (https://www.youtube.com/watch?v=2aRqpqFOiS4&t=63s).

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Bayesian Inference of State-Level COVID-19 Basic Reproduction Numbers across the United States

Although many persons in the United States have acquired immunity to COVID-19, either through vaccination or infection with SARS-CoV-2, COVID-19 will pose an ongoing threat to non-immune persons so long as disease transmission continues. We can estimate when sustained disease transmission will end in a population by calculating the population-specific basic reproduction number ℛ0, the expected number of secondary cases generated by an infected person in the absence of any interventions. The value of ℛ0 relates to a herd immunity threshold (HIT), which is given by 1-1/ℛ0. When the immune fraction of a population exceeds this threshold, sustained disease transmission becomes exponentially unlikely (barring mutations allowing SARS-CoV-2 to escape immunity). Here, we report state-level ℛ0 estimates obtained using Bayesian inference. Maximum a posteriori estimates range from 7.1 for New Jersey to 2.3 for Wyoming, indicating that disease transmission varies considerably across states and that reaching herd immunity will be more difficult in some states than others. ℛ0 estimates were obtained from compartmental models via the next-generation matrix approach after each model was parameterized using regional daily confirmed case reports of COVID-19 from 21 January 2020 to 21 June 2020. Our ℛ0 estimates characterize the infectiousness of ancestral strains, but they can be used to determine HITs for a distinct, currently dominant circulating strain, such as SARS-CoV-2 variant Delta (lineage B.1.617.2), if the relative infectiousness of the strain can be ascertained. On the basis of Delta-adjusted HITs, vaccination data, and seroprevalence survey data, we found that no state had achieved herd immunity as of 20 September 2021.

Efficient gradient-based parameter estimation for dynamic models using qualitative data

MOTIVATION

Unknown parameters of dynamical models are commonly estimated from experimental data. However, while various efficient optimization and uncertainty analysis methods have been proposed for quantitative data, methods for qualitative data are rare and suffer from bad scaling and convergence.

RESULTS

Here, we propose an efficient and reliable framework for estimating the parameters of ordinary differential equation models from qualitative data. In this framework, we derive a semi-analytical algorithm for gradient calculation of the optimal scaling method developed for qualitative data. This enables the use of efficient gradient-based optimization algorithms. We demonstrate that the use of gradient information improves performance of optimization and uncertainty quantification on several application examples. On average, we achieve a speedup of more than one order of magnitude compared to gradient-free optimization. In addition, in some examples, the gradient-based approach yields substantially improved objective function values and quality of the fits. Accordingly, the proposed framework substantially improves the parameterization of models from qualitative data.

AVAILABILITY AND IMPLEMENTATION

The proposed approach is implemented in the open-source Python Parameter EStimation TOolbox (pyPESTO). pyPESTO is available at https://github.com/ICB-DCM/pyPESTO. All application examples and code to reproduce this study are available at https://doi.org/10.5281/zenodo.4507613.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Bayesian Inference of State-Level COVID-19 Basic Reproduction Numbers across the United States

Although many persons in the United States have acquired immunity to COVID-19, either through vaccination or infection with SARS-CoV-2, COVID-19 will pose an ongoing threat to non-immune persons so long as disease transmission continues. We can estimate when sustained disease transmission will end in a population by calculating the population-specific basic reproduction number ℛ ₀ , the expected number of secondary cases generated by an infected person in the absence of any interventions. The value of ℛ ₀ relates to a herd immunity threshold (HIT), which is given by 1 - 1/ℛ ₀ . When the immune fraction of a population exceeds this threshold, sustained disease transmission becomes exponentially unlikely (barring mutations allowing SARS-CoV-2 to escape immunity). Here, we report state-level ℛ ₀ estimates obtained using Bayesian inference. Maximum a posteriori estimates range from 7.1 for New Jersey to 2.3 for Wyoming, indicating that disease transmission varies considerably across states and that reaching herd immunity will be more difficult in some states than others. ℛ ₀ estimates were obtained from compartmental models via the next-generation matrix approach after each model was parameterized using regional daily confirmed case reports of COVID-19 from 21-January-2020 to 21-June-2020. Our ℛ ₀ estimates characterize infectiousness of ancestral strains, but they can be used to determine HITs for a distinct, currently dominant circulating strain, such as SARS-CoV-2 variant Delta (lineage B.1.617.2), if the relative infectiousness of the strain can be ascertained. On the basis of Delta-adjusted HITs, vaccination data, and seroprevalence survey data, we find that no state has achieved herd immunity as of 20-September-2021.

Significance Statement

COVID-19 will continue to threaten non-immune persons in the presence of ongoing disease transmission. We can estimate when sustained disease transmission will end by calculating the population-specific basic reproduction number ℛ ₀ , which relates to a herd immunity threshold (HIT), given by 1 - 1/ℛ ₀ . When the immune fraction of a population exceeds this threshold, sustained disease transmission becomes exponentially unlikely. Here, we report state-level ℛ ₀ estimates indicating that disease transmission varies considerably across states. Our ℛ ₀ estimates can also be used to determine HITs for the Delta variant of COVID-19. On the basis of Delta-adjusted HITs, vaccination data, and serological survey results, we find that no state has yet achieved herd immunity.

Programmatic modeling for biological systems

Computational modeling has become an established technique to encode mathematical representations of cellular processes and gain mechanistic insights that drive testable predictions. These models are often constructed using graphical user interfaces or domain-specific languages, with community standards used for interchange. Models undergo steady state or dynamic analysis, which can include simulation and calibration within a single application, or transfer across various tools. Here, we describe a novel programmatic modeling paradigm, whereby modeling is augmented with software engineering best practices. We focus on Python - a popular programming language with a large scientific package ecosystem. Models can be encoded as programs, adding benefits such as modularity, testing, and automated documentation generators, while still being extensible and exportable to standardized formats for use with external tools if desired. Programmatic modeling is a key technology to enable collaborative model development and enhance dissemination, transparency, and reproducibility.

Relating simulation studies by provenance-Developing a family of Wnt signaling models

For many biological systems, a variety of simulation models exist. A new simulation model is rarely developed from scratch, but rather revises and extends an existing one. A key challenge, however, is to decide which model might be an appropriate starting point for a particular problem and why. To answer this question, we need to identify entities and activities that contributed to the development of a simulation model. Therefore, we exploit the provenance data model, PROV-DM, of the World Wide Web Consortium and, building on previous work, continue developing a PROV ontology for simulation studies. Based on a case study of 19 Wnt/β-catenin signaling models, we identify crucial entities and activities as well as useful metadata to both capture the provenance information from individual simulation studies and relate these forming a family of models. The approach is implemented in WebProv, a web application for inserting and querying provenance information. Our specialization of PROV-DM contains the entities Research Question, Assumption, Requirement, Qualitative Model, Simulation Model, Simulation Experiment, Simulation Data, and Wet-lab Data as well as activities referring to building, calibrating, validating, and analyzing a simulation model. We show that most Wnt simulation models are connected to other Wnt models by using (parts of) these models. However, the overlap, especially regarding the Wet-lab Data used for calibration or validation of the models is small. Making these aspects of developing a model explicit and queryable is an important step for assessing and reusing simulation models more effectively. Exposing this information helps to integrate a new simulation model within a family of existing ones and may lead to the development of more robust and valid simulation models. We hope that our approach becomes part of a standardization effort and that modelers adopt the benefits of provenance when considering or creating simulation models.

We revise a provenance ontology for simulation studies of cellular biochemical models. Provenance information is useful for understanding the creation of a simulation model because it not only contains information about the entities and activities that have led to a simulation model but also their relations, all of which can be visualized. It provides additional structure by explicitly recording research questions, assumptions, and requirements and relating them along with data, qualitative models, simulation models, and simulation experiments through a small set of predefined but extensible activities.

We have applied our concept to a family of 19 Wnt signaling models and implemented a web-based tool (WebProv) to store the provenance information from these studies. The resulting provenance graph visualizes the story line of simulation studies and demonstrates the creation and calibration of simulation models, the successive attempts of validation and extension, and shows, beyond an individual simulation study, how the Wnt models are related. Thereby, the steps and sources that contributed to a simulation model are made explicit.

Our approach complements other approaches aimed at facilitating the reuse and assessment of simulation products in systems biology such as model repositories as well as annotation and documentation guidelines.

Modeling the Nonlinear Dynamics of Intracellular Signaling Networks

This protocol illustrates a pipeline for modeling the nonlinear behavior of intracellular signaling pathways. At fixed spatial points, nonlinear signaling dynamics are described by ordinary differential equations (ODEs). At constant parameters, these ODEs may have multiple attractors, such as multiple steady states or limit cycles. Standard optimization procedures fine-tune the parameters for the system trajectories localized within the basin of attraction of only one attractor, usually a stable steady state. The suggested protocol samples the parameter space and captures the overall dynamic behavior by analyzing the number and stability of steady states and the shapes of the assembly of nullclines, which are determined as projections of quasi-steady-state trajectories into different 2D spaces of system variables. Our pipeline allows identifying main qualitative features of the model behavior, perform bifurcation analysis, and determine the borders separating the different dynamical regimes within the assembly of 2D parametric planes. Partial differential equation (PDE) systems describing the nonlinear spatiotemporal behavior are derived by coupling fixed point dynamics with species diffusion.

Efficient gradient-based parameter estimation for dynamic models using qualitative data

Motivation

Unknown parameters of dynamical models are commonly estimated from experimental data. However, while various efficient optimization and uncertainty analysis methods have been proposed for quantitative data, methods for qualitative data are rare and suffer from bad scaling and convergence.

Results

Here, we propose an efficient and reliable framework for estimating the parameters of ordinary differential equation models from qualitative data. In this framework, we derive a semi-analytical algorithm for gradient calculation of the optimal scaling method developed for qualitative data. This enables the use of efficient gradient-based optimization algorithms. We demonstrate that the use of gradient information improves performance of optimization and uncertainty quantification on several application examples. On average, we achieve a speedup of more than one order of magnitude compared to gradient-free optimization. In addition, in some examples, the gradient-based approach yields substantially improved objective function values and quality of the fits. Accordingly, the proposed framework substantially improves the parameterization of models from qualitative data.

Availability and implementation

The proposed approach is implemented in the open-source Python Parameter EStimation TOolbox (pyPESTO). pyPESTO is available at https://github.com/ICB-DCM/pyPESTO. All application examples and code to reproduce this study are available at https://doi.org/10.5281/zenodo.4507613.

Supplementary information

Supplementary data are available at Bioinformatics online.

Publishing reproducible dynamic kinetic models

Publishing repeatable and reproducible computational models is a crucial aspect of the scientific method in computational biology and one that is often forgotten in the rush to publish. The pressures of academic life and the lack of any reward system at institutions, granting agencies and journals means that publishing reproducible science is often either non-existent or, at best, presented in the form of an incomplete description. In the article, we will focus on repeatability and reproducibility in the systems biology field where a great many published models cannot be reproduced and in many cases even repeated. This review describes the current landscape of software tooling, model repositories, model standards and best practices for publishing repeatable and reproducible kinetic models. The review also discusses possible future remedies including working more closely with journals to help reviewers and editors ensure that published kinetic models are at minimum, repeatable. Contact: hsauro@uw.edu.

Atlas: automatic modeling of regulation of bacterial gene expression and metabolism using rule-based languages

MOTIVATION

Cells are complex systems composed of hundreds of genes whose products interact to produce elaborated behaviors. To control such behaviors, cells rely on transcription factors to regulate gene expression, and gene regulatory networks (GRNs) are employed to describe and understand such behavior. However, GRNs are static models, and dynamic models are difficult to obtain due to their size, complexity, stochastic dynamics and interactions with other cell processes.

RESULTS

We developed Atlas, a Python software that converts genome graphs and gene regulatory, interaction and metabolic networks into dynamic models. The software employs these biological networks to write rule-based models for the PySB framework. The underlying method is a divide-and-conquer strategy to obtain sub-models and combine them later into an ensemble model. To exemplify the utility of Atlas, we used networks of varying size and complexity of Escherichia coli and evaluated in silico modifications, such as gene knockouts and the insertion of promoters and terminators. Moreover, the methodology could be applied to the dynamic modeling of natural and synthetic networks of any bacteria.

AVAILABILITY AND IMPLEMENTATION

Code, models and tutorials are available online (https://github.com/networkbiolab/atlas).

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Epidemiological modeling in StochSS Live!

Summary

We present StochSS Live!, a web-based service for modeling, simulation and analysis of a wide range of mathematical, biological and biochemical systems. Using an epidemiological model of COVID-19, we demonstrate the power of StochSS Live! to enable researchers to quickly develop a deterministic or a discrete stochastic model, infer its parameters and analyze the results.

Availability and implementation

StochSS Live! is freely available at https://live.stochss.org/

Supplementary information

Supplementary data are available at Bioinformatics online.

PEtab-Interoperable specification of parameter estimation problems in systems biology

Reproducibility and reusability of the results of data-based modeling studies are essential. Yet, there has been-so far-no broadly supported format for the specification of parameter estimation problems in systems biology. Here, we introduce PEtab, a format which facilitates the specification of parameter estimation problems using Systems Biology Markup Language (SBML) models and a set of tab-separated value files describing the observation model and experimental data as well as parameters to be estimated. We already implemented PEtab support into eight well-established model simulation and parameter estimation toolboxes with hundreds of users in total. We provide a Python library for validation and modification of a PEtab problem and currently 20 example parameter estimation problems based on recent studies.

Quantifying the roles of space and stochasticity in computer simulations for cell biology and cellular biochemistry

Most of the fascinating phenomena studied in cell biology emerge from interactions among highly organized multimolecular structures embedded into complex and frequently dynamic cellular morphologies. For the exploration of such systems, computer simulation has proved to be an invaluable tool, and many researchers in this field have developed sophisticated computational models for application to specific cell biological questions. However, it is often difficult to reconcile conflicting computational results that use different approaches to describe the same phenomenon. To address this issue systematically, we have defined a series of computational test cases ranging from very simple to moderately complex, varying key features of dimensionality, reaction type, reaction speed, crowding, and cell size. We then quantified how explicit spatial and/or stochastic implementations alter outcomes, even when all methods use the same reaction network, rates, and concentrations. For simple cases, we generally find minor differences in solutions of the same problem. However, we observe increasing discordance as the effects of localization, dimensionality reduction, and irreversible enzymatic reactions are combined. We discuss the strengths and limitations of commonly used computational approaches for exploring cell biological questions and provide a framework for decision making by researchers developing new models. As computational power and speed continue to increase at a remarkable rate, the dream of a fully comprehensive computational model of a living cell may be drawing closer to reality, but our analysis demonstrates that it will be crucial to evaluate the accuracy of such models critically and systematically.

MASSpy: Building, simulating, and visualizing dynamic biological models in Python using mass action kinetics

Mathematical models of metabolic networks utilize simulation to study system-level mechanisms and functions. Various approaches have been used to model the steady state behavior of metabolic networks using genome-scale reconstructions, but formulating dynamic models from such reconstructions continues to be a key challenge. Here, we present the Mass Action Stoichiometric Simulation Python (MASSpy) package, an open-source computational framework for dynamic modeling of metabolism. MASSpy utilizes mass action kinetics and detailed chemical mechanisms to build dynamic models of complex biological processes. MASSpy adds dynamic modeling tools to the COnstraint-Based Reconstruction and Analysis Python (COBRApy) package to provide an unified framework for constraint-based and kinetic modeling of metabolic networks. MASSpy supports high-performance dynamic simulation through its implementation of libRoadRunner: the Systems Biology Markup Language (SBML) simulation engine. Three examples are provided to demonstrate how to use MASSpy: (1) a validation of the MASSpy modeling tool through dynamic simulation of detailed mechanisms of enzyme regulation; (2) a feature demonstration using a workflow for generating ensemble of kinetic models using Monte Carlo sampling to approximate missing numerical values of parameters and to quantify biological uncertainty, and (3) a case study in which MASSpy is utilized to overcome issues that arise when integrating experimental data with the computation of functional states of detailed biological mechanisms. MASSpy represents a powerful tool to address challenges that arise in dynamic modeling of metabolic networks, both at small and large scales.

Co-regulation of the transcription controlling ATF2 phosphoswitch by JNK and p38

Transcription factor phosphorylation at specific sites often activates gene expression, but how environmental cues quantitatively control transcription is not well-understood. Activating protein 1 transcription factors are phosphorylated by mitogen-activated protein kinases (MAPK) in their transactivation domains (TAD) at so-called phosphoswitches, which are a hallmark in response to growth factors, cytokines or stress. We show that the ATF2 TAD is controlled by functionally distinct signaling pathways (JNK and p38) through structurally different MAPK binding sites. Moreover, JNK mediated phosphorylation at an evolutionarily more recent site diminishes p38 binding and made the phosphoswitch differently sensitive to JNK and p38 in vertebrates. Structures of MAPK-TAD complexes and mechanistic modeling of ATF2 TAD phosphorylation in cells suggest that kinase binding motifs and phosphorylation sites line up to maximize MAPK based co-regulation. This study shows how the activity of an ancient transcription controlling phosphoswitch became dependent on the relative flux of upstream signals.

Best Practices for Making Reproducible Biochemical Models

Like many scientific disciplines, dynamical biochemical modeling is hindered by irreproducible results. This limits the utility of biochemical models by making them difficult to understand, trust, or reuse. We comprehensively list the best practices that biochemical modelers should follow to build reproducible biochemical model artifacts-all data, model descriptions, and custom software used by the model-that can be understood and reused. The best practices provide advice for all steps of a typical biochemical modeling workflow in which a modeler collects data; constructs, trains, simulates, and validates the model; uses the predictions of a model to advance knowledge; and publicly shares the model artifacts. The best practices emphasize the benefits obtained by using standard tools and formats and provides guidance to modelers who do not or cannot use standards in some stages of their modeling workflow. Adoption of these best practices will enhance the ability of researchers to reproduce, understand, and reuse biochemical models.

Parameterization of mechanistic models from qualitative data using an efficient optimal scaling approach

Quantitative dynamical models facilitate the understanding of biological processes and the prediction of their dynamics. These models usually comprise unknown parameters, which have to be inferred from experimental data. For quantitative experimental data, there are several methods and software tools available. However, for qualitative data the available approaches are limited and computationally demanding. Here, we consider the optimal scaling method which has been developed in statistics for categorical data and has been applied to dynamical systems. This approach turns qualitative variables into quantitative ones, accounting for constraints on their relation. We derive a reduced formulation for the optimization problem defining the optimal scaling. The reduced formulation possesses the same optimal points as the established formulation but requires less degrees of freedom. Parameter estimation for dynamical models of cellular pathways revealed that the reduced formulation improves the robustness and convergence of optimizers. This resulted in substantially reduced computation times. We implemented the proposed approach in the open-source Python Parameter EStimation TOolbox (pyPESTO) to facilitate reuse and extension. The proposed approach enables efficient parameterization of quantitative dynamical models using qualitative data.

SBML Level 3: an extensible format for the exchange and reuse of biological models

Systems biology has experienced dramatic growth in the number, size, and complexity of computational models. To reproduce simulation results and reuse models, researchers must exchange unambiguous model descriptions. We review the latest edition of the Systems Biology Markup Language (SBML), a format designed for this purpose. A community of modelers and software authors developed SBML Level 3 over the past decade. Its modular form consists of a core suited to representing reaction-based models and packages that extend the core with features suited to other model types including constraint-based models, reaction-diffusion models, logical network models, and rule-based models. The format leverages two decades of SBML and a rich software ecosystem that transformed how systems biologists build and interact with models. More recently, the rise of multiscale models of whole cells and organs, and new data sources such as single-cell measurements and live imaging, has precipitated new ways of integrating data with models. We provide our perspectives on the challenges presented by these developments and how SBML Level 3 provides the foundation needed to support this evolution.

A Probabilistic Approach to Explore Signal Execution Mechanisms With Limited Experimental Data

Mathematical models of biochemical reaction networks are central to the study of dynamic cellular processes and hypothesis generation that informs experimentation and validation. Unfortunately, model parameters are often not available and sparse experimental data leads to challenges in model calibration and parameter estimation. This can in turn lead to unreliable mechanistic interpretations of experimental data and the generation of poorly conceived hypotheses for experimental validation. To address this challenge, we evaluate whether a Bayesian-inspired probability-based approach, that relies on expected values for quantities of interest calculated from available information regarding the reaction network topology and parameters can be used to qualitatively explore hypothetical biochemical network execution mechanisms in the context of limited available data. We test our approach on a model of extrinsic apoptosis execution to identify preferred signal execution modes across varying conditions. Apoptosis signal processing can take place either through a mitochondria independent (Type I) mode or a mitochondria dependent (Type II) mode. We first show that in silico knockouts, represented by model subnetworks, successfully identify the most likely execution mode for specific concentrations of key molecular regulators. We then show that changes in molecular regulator concentrations alter the overall reaction flux through the network by shifting the primary route of signal flow between the direct caspase and mitochondrial pathways. Our work thus demonstrates that probabilistic approaches can be used to explore the qualitative dynamic behavior of model biochemical systems even with missing or sparse data.

Bayesian inference using qualitative observations of underlying continuous variables

MOTIVATION

Recent work has demonstrated the feasibility of using non-numerical, qualitative data to parameterize mathematical models. However, uncertainty quantification (UQ) of such parameterized models has remained challenging because of a lack of a statistical interpretation of the objective functions used in optimization.

RESULTS

We formulated likelihood functions suitable for performing Bayesian UQ using qualitative observations of underlying continuous variables or a combination of qualitative and quantitative data. To demonstrate the resulting UQ capabilities, we analyzed a published model for immunoglobulin E (IgE) receptor signaling using synthetic qualitative and quantitative datasets. Remarkably, estimates of parameter values derived from the qualitative data were nearly as consistent with the assumed ground-truth parameter values as estimates derived from the lower throughput quantitative data. These results provide further motivation for leveraging qualitative data in biological modeling.

AVAILABILITY AND IMPLEMENTATION

The likelihood functions presented here are implemented in a new release of PyBioNetFit, an open-source application for analyzing Systems Biology Markup Language- and BioNetGen Language-formatted models, available online at www.github.com/lanl/PyBNF.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Multisite EGFR phosphorylation is regulated by adaptor protein abundances and dimer lifetimes

Differential epidermal growth factor receptor (EGFR) phosphorylation is thought to couple receptor activation to distinct signaling pathways. However, the molecular mechanisms responsible for biased signaling are unresolved due to a lack of insight into the phosphorylation patterns of full-length EGFR. We extended a single-molecule pull-down technique previously used to study protein-protein interactions to allow for robust measurement of receptor phosphorylation. We found that EGFR is predominantly phosphorylated at multiple sites, yet phosphorylation at specific tyrosines is variable and only a subset of receptors share phosphorylation at the same site, even with saturating ligand concentrations. We found distinct populations of receptors as soon as 1 min after ligand stimulation, indicating early diversification of function. To understand this heterogeneity, we developed a mathematical model. The model predicted that variations in phosphorylation are dependent on the abundances of signaling partners, while phosphorylation levels are dependent on dimer lifetimes. The predictions were confirmed in studies of cell lines with different expression levels of signaling partners, and in experiments comparing low- and high-affinity ligands and oncogenic EGFR mutants. These results reveal how ligand-regulated receptor dimerization dynamics and adaptor protein concentrations play critical roles in EGFR signaling.

Parameter Estimation and Uncertainty Quantification for Systems Biology Models

Mathematical models can provide quantitative insights into immunoreceptor signaling, and other biological processes, but require parameterization and uncertainty quantification before reliable predictions become possible. We review currently available methods and software tools to address these problems. We consider gradient-based and gradient-free methods for point estimation of parameter values, and methods of profile likelihood, bootstrapping, and Bayesian inference for uncertainty quantification. We consider recent and potential future applications of these methods to systems-level modeling of immune-related phenomena.

Pleione: A tool for statistical and multi-objective calibration of Rule-based models

Mathematical models based on Ordinary Differential Equations (ODEs) are frequently used to describe and simulate biological systems. Nevertheless, such models are often difficult to understand. Unlike ODE models, Rule-Based Models (RBMs) utilise formal language to describe reactions as a cumulative number of statements that are easier to understand and correct. They are also gaining popularity because of their conciseness and simulation flexibility. However, RBMs generally lack tools to perform further analysis that requires simulation. This situation arises because exact and approximate simulations are computationally intensive. Translating RBMs into ODEs is commonly used to reduce simulation time, but this technique may be prohibitive due to combinatorial explosion. Here, we present the software called Pleione to calibrate RBMs. Parameter calibration is essential given the incomplete experimental determination of reaction rates and the goal of using models to reproduce experimental data. The software distributes stochastic simulations and calculations and incorporates equivalence tests to determine the fitness of RBMs compared with data. The primary features of Pleione were thoroughly tested on a model of gene regulation in Escherichia coli. Pleione yielded satisfactory results regarding calculation time and error reduction for multiple simulators, models, parameter search strategies, and computing infrastructures.

PyBioNetFit and the Biological Property Specification Language

In systems biology modeling, important steps include model parameterization, uncertainty quantification, and evaluation of agreement with experimental observations. To help modelers perform these steps, we developed the software PyBioNetFit, which in addition supports checking models against known system properties and solving design problems. PyBioNetFit introduces Biological Property Specification Language (BPSL) for the formal declaration of system properties. BPSL allows qualitative data to be used alone or in combination with quantitative data. PyBioNetFit performs parameterization with parallelized metaheuristic optimization algorithms that work directly with existing model definition standards: BioNetGen Language (BNGL) and Systems Biology Markup Language (SBML). We demonstrate PyBioNetFit's capabilities by solving various example problems, including the challenging problem of parameterizing a 153-parameter model of cell cycle control in yeast based on both quantitative and qualitative data. We demonstrate the model checking and design applications of PyBioNetFit and BPSL by analyzing a model of targeted drug interventions in autophagy signaling.