Rule-based modeling of biochemical systems with BioNetGen

Leveraging public AI tools to explore systems biology resources in mathematical modeling

Predictive mathematical modeling is an essential part of systems biology and is interconnected with information management. Systems biology information is often stored in specialized formats to facilitate data storage and analysis. These formats are not designed for easy human readability and thus require specialized software to visualize and interpret results. Therefore, comprehending modeling and underlying networks and pathways is contingent on mastering systems biology tools, which is particularly challenging for users with no or little background in data science or system biology. To address this challenge, we investigated the usage of public Artificial Intelligence (AI) tools in exploring systems biology resources in mathematical modeling. We tested public AI's understanding of mathematics in models, related systems biology data, and the complexity of model structures. Our approach can enhance the accessibility of systems biology for non-system biologists and help them understand systems biology without a deep learning curve.

A spatial model of autophosphorylation of Ca2+/calmodulin-dependent protein kinase II (CaMKII) predicts that the lifetime of phospho-CaMKII after induction of synaptic plasticity is greatly prolonged by CaM-trapping

Long-term potentiation (LTP) is a biochemical process that underlies learning in excitatory glutamatergic synapses in the Central Nervous System (CNS). The critical early driver of LTP is autophosphorylation of the abundant postsynaptic enzyme, Ca2+/calmodulin-dependent protein kinase II (CaMKII). Autophosphorylation is initiated by Ca2+ flowing through NMDA receptors activated by strong synaptic activity. Its lifetime is ultimately determined by the balance of the rates of autophosphorylation and of dephosphorylation by protein phosphatase 1 (PP1). Here we have modeled the autophosphorylation and dephosphorylation of CaMKII during synaptic activity in a spine synapse using MCell4, an open source computer program for creating particle-based stochastic, and spatially realistic models of cellular microchemistry. The model integrates four earlier detailed models of separate aspects of regulation of spine Ca2+ and CaMKII activity, each of which incorporate experimentally measured biochemical parameters and have been validated against experimental data. We validate the composite model by showing that it accurately predicts previous experimental measurements of effects of NMDA receptor activation, including high sensitivity of induction of LTP to phosphatase activity in vivo, and persistence of autophosphorylation for a period of minutes after the end of synaptic stimulation. We then use the model to probe aspects of the mechanism of regulation of autophosphorylation of CaMKII that are difficult to measure in vivo. We examine the effects of "CaM-trapping," a process in which the affinity for Ca2+/CaM increases several hundred-fold after autophosphorylation. We find that CaM-trapping does not increase the proportion of autophosphorylated subunits in holoenzymes after a complex stimulus, as previously hypothesized. Instead, CaM-trapping may dramatically prolong the lifetime of autophosphorylated CaMKII through steric hindrance of dephosphorylation by protein phosphatase 1. The results provide motivation for experimental measurement of the extent of suppression of dephosphorylation of CaMKII by bound Ca2+/CaM. The composite MCell4 model of biochemical effects of complex stimuli in synaptic spines is a powerful new tool for realistic, detailed dissection of mechanisms of synaptic plasticity.

Parallelization of particle-based reaction-diffusion simulations using MPI

Particle-based reaction-diffusion models offer a high-resolution alternative to the continuum reaction-diffusion approach, capturing the discrete and volume-excluding nature of molecules undergoing stochastic dynamics. These methods are thus uniquely capable of simulating explicit self-assembly of particles into higher-order structures like filaments, spherical cages, or heterogeneous macromolecular complexes, which are ubiquitous across living systems and in materials design. The disadvantage of these high-resolution methods is their increased computational cost. Here we present a parallel implementation of the particle-based NERDSS software using the Message Passing Interface (MPI) and spatial domain decomposition, achieving close to linear scaling for up to 96 processors in the largest simulation systems. The scalability of parallel NERDSS is evaluated for bimolecular reactions in 3D and 2D, for self-assembly of trimeric and hexameric complexes, and for protein lattice assembly from 3D to 2D, with all parallel test cases producing accurate solutions. We demonstrate how parallel efficiency depends on the system size, the reaction network, and the limiting timescales of the system, showing optimal scaling only for smaller assemblies with slower timescales. The formation of very large assemblies represents a challenge in evaluating reaction updates across processors, and here we restrict assembly sizes to below the spatial decomposition size. We provide the parallel NERDSS code open source, with detailed documentation for developers and extension to other particle-based reaction-diffusion software.

Expressive rule-based modeling and fast simulation for dynamic compartments

Compartmentalization is vital for cell biological processes. The field of rule-based stochastic simulation has acknowledged this, and many tools and methods have capabilities for compartmentalization. However, mostly, this is limited to a static compartmental hierarchy and does not integrate compartmental changes. Integrating compartmental dynamics is challenging for the design of the modeling language and the simulation engine. The language should support a concise yet flexible modeling of compartmental dynamics. Our work is based on ML-Rules, a rule-based language for multi-level cell biological modeling that supports a wide variety of compartmental dynamics, whose syntax we slightly adapt. To develop an efficient simulation engine for compartmental dynamics, we combine specific data structures and new and existing algorithms and implement them in the Rust programming language. We evaluate the concept and implementation using two case studies from existing cell-biological models. The execution of these models outperforms previous simulations of ML-Rules by two orders of magnitude. Finally, we present a prototype of a WebAssembly-based implementation to allow for a low barrier of entry when exploring the language and associated models without the need for local installation.

Spleen Tyrosine Kinase (SYK) negatively regulates ITAM-mediated human NK cell signaling and CD19-CAR NK cell efficacy

NK cells express activating receptors that signal through ITAM-bearing adapter proteins. The phosphorylation of each ITAM creates binding sites for SYK and ZAP70 protein tyrosine kinases to propagate downstream signaling including the induction ofCa 2 + influx. While all immature and mature human NK cells co-express SYK and ZAP70, clonally driven memory or adaptive NK cells can methylate SYK genes and signaling is mediated exclusively using ZAP70. Here, we examined the role of SYK and ZAP70 in a clonal human NK cell line KHYG1 by CRISPR-based deletion using a combination of experiments and mechanistic computational modeling. Elimination of SYK resulted in more robustCa + + influx after cross-linking of the CD16 and NKp30 receptors and enhanced phosphorylation of downstream proteins, whereas ZAP70 deletion diminished these responses. By contrast, ZAP70 depletion increased proliferation of the NK cells. As immature T cells express both SYK and ZAP70 but mature T cells often express only ZAP70, we transduced the human Jurkat cell line with SYK and found that expression of SYK increased proliferation but diminished TCR-inducedCa 2 + flux and activation. We performed transcriptional analysis of the matched sets of variant Jurkat and KHYG1 cells and observed profound alterations caused by SYK expression. As depletion of SYK in NK cells increased their activation, primary human NK cells were transduced with a CD19-targeting CAR and were CRISPR edited to ablate SYK or ZAP70. Deletion of SYK resulted in more robust cytotoxic activity and cytokine production, providing a new therapeutic strategy of NK cell engineering for cancer immunotherapy.

Identifying Bayesian optimal experiments for uncertain biochemical pathway models

Pharmacodynamic (PD) models are mathematical models of cellular reaction networks that include drug mechanisms of action. These models are useful for studying predictive therapeutic outcomes of novel drug therapies in silico. However, PD models are known to possess significant uncertainty with respect to constituent parameter data, leading to uncertainty in the model predictions. Furthermore, experimental data to calibrate these models is often limited or unavailable for novel pathways. In this study, we present a Bayesian optimal experimental design approach for improving PD model prediction accuracy. We then apply our method using simulated experimental data to account for uncertainty in hypothetical laboratory measurements. This leads to a probabilistic prediction of drug performance and a quantitative measure of which prospective laboratory experiment will optimally reduce prediction uncertainty in the PD model. The methods proposed here provide a way forward for uncertainty quantification and guided experimental design for models of novel biological pathways.

bnglViz: online visualization of rule-based models

MOTIVATION

Rule-based modeling is a powerful method to describe and simulate interactions among multi-site molecules and multi-molecular species, accounting for the internal connectivity of molecules in chemical species. This modeling technique is implemented in BioNetGen software that is used by various tools and software frameworks, such as BioNetGen stand-alone software, NFSim simulation engine, Virtual Cell simulation and modeling framework, SmolDyn and PySB software tools. These tools exchange models using BioNetGen scripting language (BNGL). Until now, there was no online visualization of such rule-based models. Modelers and researchers reading the manuscripts describing rule-based models had to learn BNGL scripting or master one of these tools to understand the models.

RESULTS

Here, we introduce bnglViz, an online platform for visualizing BNGL files as graphical cartoons, empowering researchers to grasp the nuances of rule-based models swiftly and efficiently, and making the exploration of complex biological systems more accessible than ever before. The produced visualizations can be used as supplemental figures in publications or as a way to annotate BNGL models on web repositories.

AVAILABILITY AND IMPLEMENTATION

Available at https://bnglviz.github.io/.

An architecture for collaboration in systems biology at the age of the Metaverse

As the current state of the Metaverse is largely driven by corporate interests, which may not align with scientific goals and values, academia should play a more active role in its development. Here, we present the challenges and solutions for building a Metaverse that supports systems biology research and collaboration. Our solution consists of two components: Kosmogora, a server ensuring biological data access, traceability, and integrity in the context of a highly collaborative environment such as a metaverse; and ECellDive, a virtual reality application to explore, interact, and build upon the data managed by Kosmogora. We illustrate the synergy between the two components by visualizing a metabolic network and its flux balance analysis. We also argue that the Metaverse of systems biology will foster closer communication and cooperation between experimentalists and modelers in the field.

A detailed molecular network map and model of the NLRP3 inflammasome

The NLRP3 inflammasome is a key regulator of inflammation that responds to a broad range of stimuli. The exact mechanism of activation has not been determined, but there is a consensus on cellular potassium efflux as a major common denominator. Once NLRP3 is activated, it forms high-order complexes together with NEK7 that trigger aggregation of ASC into specks. Typically, there is only one speck per cell, consistent with the proposal that specks form - or end up at - the centrosome. ASC polymerisation in turn triggers caspase-1 activation, leading to maturation and release of IL-1β and pyroptosis, i.e., highly inflammatory cell death. Several gain-of-function mutations in the NLRP3 inflammasome have been suggested to induce spontaneous activation of NLRP3 and hence contribute to development and disease severity in numerous autoinflammatory and autoimmune diseases. Consequently, the NLRP3 inflammasome is of significant clinical interest, and recent attention has drastically improved our insight in the range of involved triggers and mechanisms of signal transduction. However, despite recent progress in knowledge, a clear and comprehensive overview of how these mechanisms interplay to shape the system level function is missing from the literature. Here, we provide such an overview as a resource to researchers working in or entering the field, as well as a computational model that allows for evaluating and explaining the function of the NLRP3 inflammasome system from the current molecular knowledge. We present a detailed reconstruction of the molecular network surrounding the NLRP3 inflammasome, which account for each specific reaction and the known regulatory constraints on each event as well as the mechanisms of drug action and impact of genetics when known. Furthermore, an executable model from this network reconstruction is generated with the aim to be used to explain NLRP3 activation from priming and activation to the maturation and release of IL-1β and IL-18. Finally, we test this detailed mechanistic model against data on the effect of different modes of inhibition of NLRP3 assembly. While the exact mechanisms of NLRP3 activation remains elusive, the literature indicates that the different stimuli converge on a single activation mechanism that is additionally controlled by distinct (positive or negative) priming and licensing events through covalent modifications of the NLRP3 molecule. Taken together, we present a compilation of the literature knowledge on the molecular mechanisms on NLRP3 activation, a detailed mechanistic model of NLRP3 activation, and explore the convergence of diverse NLRP3 activation stimuli into a single input mechanism.

MolClustPy: a Python package to characterize multivalent biomolecular clusters

SUMMARY

Low-affinity interactions among multivalent biomolecules may lead to the formation of molecular complexes that undergo phase transitions to become supply-limited large clusters. In stochastic simulations, such clusters display a wide range of sizes and compositions. We have developed a Python package, MolClustPy, which performs multiple stochastic simulation runs using NFsim (Network-Free stochastic simulator); MolClustPy characterizes and visualizes the distribution of cluster sizes, molecular composition, and bonds across molecular clusters. The statistical analysis offered by MolClustPy is readily applicable to other stochastic simulation software, such as SpringSaLaD and ReaDDy.

AVAILABILITY AND IMPLEMENTATION

The software is implemented in Python. A detailed Jupyter notebook is provided to enable convenient running. Code, user guide, and examples are freely available at https://molclustpy.github.io/.

Combining Multikinase Tyrosine Kinase Inhibitors Targeting the Vascular Endothelial Growth Factor and Cluster of Differentiation 47 Signaling Pathways Is Predicted to Increase the Efficacy of Antiangiogenic Combination Therapies

Angiogenesis is a critical step in tumor growth, development, and invasion. Nascent tumor cells secrete vascular endothelial growth factor (VEGF) that significantly remodels the tumor microenvironment through interaction with multiple receptors on vascular endothelial cells, including type 2 VEGF receptor (VEGFR2). The complex pathways initiated by VEGF binding to VEGFR2 lead to enhanced proliferation, survival, and motility of vascular endothelial cells and formation of a new vascular network, enabling tumor growth. Antiangiogenic therapies that inhibit VEGF signaling pathways were among the first drugs that targeted stroma rather than tumor cells. Despite improvements in progression-free survival and higher response rates relative to chemotherapy in some types of solid tumors, the impact on overall survival (OS) has been limited, with the majority of tumors eventually relapsing due to resistance or activation of alternate angiogenic pathways. Here, we developed a molecularly detailed computational model of endothelial cell signaling and angiogenesis-driven tumor growth to investigate combination therapies targeting different nodes of the endothelial VEGF/VEGFR2 signaling pathway. Simulations predicted a strong threshold-like behavior in extracellular signal-regulated kinases 1/2 (ERK1/2) activation relative to phosphorylated VEGFR2 levels, as continuous inhibition of at least 95% of receptors was necessary to abrogate phosphorylated ERK1/2 (pERK1/2). Combinations with mitogen-activated protein kinase/ERK kinase (MEK) and spingosine-1-phosphate inhibitors were found to be effective in overcoming the ERK1/2 activation threshold and abolishing activation of the pathway. Modeling results also identified a mechanism of resistance whereby tumor cells could reduce pERK1/2 sensitivity to inhibitors of VEGFR2 by upregulation of Raf, MEK, and sphingosine kinase 1 (SphK1), thus highlighting the need for deeper investigation of the dynamics of the crosstalk between VEGFR2 and SphK1 pathways. Inhibition of VEGFR2 phosphorylation was found to be more effective at blocking protein kinase B, also known as AKT, activation; however, to effectively abolish AKT activation, simulations identified Axl autophosphorylation or the Src kinase domain as potent targets. Simulations also supported activating cluster of differentiation 47 (CD47) on endothelial cells as an effective combination partner with tyrosine kinase inhibitors to inhibit angiogenesis signaling and tumor growth. Virtual patient simulations supported the effectiveness of CD47 agonism in combination with inhibitors of VEGFR2 and SphK1 pathways. Overall, the rule-based system model developed here provides new insights, generates novel hypothesis, and makes predictions regarding combinations that may enhance the OS with currently approved antiangiogenic therapies.

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.

Computational speed-up of large-scale, single-cell model simulations via a fully integrated SBML-based format

Summary

Large-scale and whole-cell modeling has multiple challenges, including scalable model building and module communication bottlenecks (e.g. between metabolism, gene expression, signaling, etc.). We previously developed an open-source, scalable format for a large-scale mechanistic model of proliferation and death signaling dynamics, but communication bottlenecks between gene expression and protein biochemistry modules remained. Here, we developed two solutions to communication bottlenecks that speed-up simulation by ∼4-fold for hybrid stochastic-deterministic simulations and by over 100-fold for fully deterministic simulations. Fully deterministic speed-up facilitates model initialization, parameter estimation and sensitivity analysis tasks.

Availability and implementation

Source code is freely available at https://github.com/birtwistlelab/SPARCED/releases/tag/v1.3.0 implemented in python, and supported on Linux, Windows and MacOS (via Docker).

MolClustPy: A Python Package to Characterize Multivalent Biomolecular Clusters

S ummary

Low-affinity interactions among multivalent biomolecules may lead to the formation of molecular complexes that undergo phase transitions to become extra-large clusters. Characterizing the physical properties of these clusters is important in recent biophysical research. Due to weak interactions such clusters are highly stochastic, demonstrating a wide range of sizes and compositions. We have developed a Python package to perform multiple stochastic simulation runs using NFsim (Network-Free stochastic simulator), characterize and visualize the distribution of cluster sizes, molecular composition, and bonds across molecular clusters and individual molecules of different types.

A vailability and implementation

The software is implemented in Python. A detailed Jupyter notebook is provided to enable convenient running. Code, user guide and examples are freely available at https://molclustpy.github.io/.

C ontact

achattaraj007@gmail.com , blinov@uchc.edu.

S upplementary information

Available at https://molclustpy.github.io/.

McSNAC: A software to approximate first-order signaling networks from mass cytometry data

Background

Mass cytometry (CyTOF) gives unprecedented opportunity to simultaneously measure up to 40 proteins in single cells, with a theoretical potential to reach 100 proteins. This high-dimensional single-cell information can be very useful in dissecting mechanisms of cellular activity. In particular, measuring abundances of signaling proteins like phospho-proteins can provide detailed information on the dynamics of single-cell signaling processes. However, computational analysis is required to reconstruct such networks with a mechanistic model.

Methods

We propose our Mass cytometry Signaling Network Analysis Code (McSNAC), a new software capable of reconstructing signaling networks and estimating their kinetic parameters from CyTOF data. McSNAC approximates signaling networks as a network of first-order reactions between proteins. This assumption often breaks down as signaling reactions can involve binding and unbinding, enzymatic reactions, and other nonlinear constructions. Furthermore, McSNAC may be limited to approximating indirect interactions between protein species, as cytometry experiments are only able to assay a small fraction of protein species involved in signaling.

Results

We carry out a series of in silico experiments here to show (1) McSNAC is capable of accurately estimating the ground-truth model in a scalable manner when given data originating from a first-order system; (2) McSNAC is capable of qualitatively predicting outcomes to perturbations of species abundances in simple second-order reaction models and in a complex in silico nonlinear signaling network in which some proteins are unmeasured.

Conclusions

These findings demonstrate that McSNAC can be a valuable screening tool for generating models of signaling networks from time-stamped CyTOF data.

Targeted Quantification of Protein Phosphorylation and Its Contributions towards Mathematical Modeling of Signaling Pathways

Post-translational modifications (PTMs) are key regulatory mechanisms that can control protein function. Of these, phosphorylation is the most common and widely studied. Because of its importance in regulating cell signaling, precise and accurate measurements of protein phosphorylation across wide dynamic ranges are crucial to understanding how signaling pathways function. Although immunological assays are commonly used to detect phosphoproteins, their lack of sensitivity, specificity, and selectivity often make them unreliable for quantitative measurements of complex biological samples. Recent advances in Mass Spectrometry (MS)-based targeted proteomics have made it a more useful approach than immunoassays for studying the dynamics of protein phosphorylation. Selected reaction monitoring (SRM)-also known as multiple reaction monitoring (MRM)-and parallel reaction monitoring (PRM) can quantify relative and absolute abundances of protein phosphorylation in multiplexed fashions targeting specific pathways. In addition, the refinement of these tools by enrichment and fractionation strategies has improved measurement of phosphorylation of low-abundance proteins. The quantitative data generated are particularly useful for building and parameterizing mathematical models of complex phospho-signaling pathways. Potentially, these models can provide a framework for linking analytical measurements of clinical samples to better diagnosis and treatment of disease.

Excel-SBOL Converter: Creating SBOL from Excel Templates and Vice Versa

Standards support synthetic biology research by enabling the exchange of component information. However, using formal representations, such as the Synthetic Biology Open Language (SBOL), typically requires either a thorough understanding of these standards or a suite of tools developed in concurrence with the ontologies. Since these tools may be a barrier for use by many practitioners, the Excel-SBOL Converter was developed to facilitate the use of SBOL and integration into existing workflows. The converter consists of two Python libraries: one that converts Excel templates to SBOL and another that converts SBOL to an Excel workbook. Both libraries can be used either directly or via a SynBioHub plugin.

Computational analysis of 4-1BB-induced NFκB signaling suggests improvements to CAR cell design

BACKGROUND

Chimeric antigen receptor (CAR)-expressing cells are a powerful modality of adoptive cell therapy against cancer. The potency of signaling events initiated upon antigen binding depends on the costimulatory domain within the structure of the CAR. One such costimulatory domain is 4-1BB, which affects cellular response via the NFκB pathway. However, the quantitative aspects of 4-1BB-induced NFκB signaling are not fully understood.

METHODS

We developed an ordinary differential equation-based mathematical model representing canonical NFκB signaling activated by CD19scFv-4-1BB. After a global sensitivity analysis on model parameters, we ran Monte Carlo simulations of cell population-wide variability in NFκB signaling and quantified the mutual information between the extracellular signal and different levels of the NFκB signal transduction pathway.

RESULTS

In response to a wide range of antigen concentrations, the magnitude of the transient peak in NFκB nuclear concentration varies significantly, while the timing of this peak is relatively consistent. Global sensitivity analysis showed that the model is robust to variations in parameters, and thus, its quantitative predictions would remain applicable to a broad range of parameter values. The model predicts that overexpressing NEMO and disabling IKKβ deactivation can increase the mutual information between antigen levels and NFκB activation.

CONCLUSIONS

Our modeling predictions provide actionable insights to guide CAR development. Particularly, we propose specific manipulations to the NFκB signal transduction pathway that can fine-tune the response of CD19scFv-4-1BB cells to the antigen concentrations they are likely to encounter. Video Abstract.

Vivarium: an interface and engine for integrative multiscale modeling in computational biology

MOTIVATION

This article introduces Vivarium-software born of the idea that it should be as easy as possible for computational biologists to define any imaginable mechanistic model, combine it with existing models and execute them together as an integrated multiscale model. Integrative multiscale modeling confronts the complexity of biology by combining heterogeneous datasets and diverse modeling strategies into unified representations. These integrated models are then run to simulate how the hypothesized mechanisms operate as a whole. But building such models has been a labor-intensive process that requires many contributors, and they are still primarily developed on a case-by-case basis with each project starting anew. New software tools that streamline the integrative modeling effort and facilitate collaboration are therefore essential for future computational biologists.

RESULTS

Vivarium is a software tool for building integrative multiscale models. It provides an interface that makes individual models into modules that can be wired together in large composite models, parallelized across multiple CPUs and run with Vivarium's discrete-event simulation engine. Vivarium's utility is demonstrated by building composite models that combine several modeling frameworks: agent-based models, ordinary differential equations, stochastic reaction systems, constraint-based models, solid-body physics and spatial diffusion. This demonstrates just the beginning of what is possible-Vivarium will be able to support future efforts that integrate many more types of models and at many more biological scales.

AVAILABILITY AND IMPLEMENTATION

The specific models, simulation pipelines and notebooks developed for this article are all available at the vivarium-notebooks repository: https://github.com/vivarium-collective/vivarium-notebooks. Vivarium-core is available at https://github.com/vivarium-collective/vivarium-core, and has been released on Python Package Index. The Vivarium Collective (https://vivarium-collective.github.io) is a repository of freely available Vivarium processes and composites, including the processes used in Section 3. Supplementary Materials provide with an extensive methodology section, with several code listings that demonstrate the basic interfaces.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

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.

Relationship Between Dimensionality and Convergence of Optimization Algorithms: A Comparison Between Data-Driven Normalization and Scaling Factor-Based Methods Using PEPSSBI

Ordinary differential equation models are used to represent intracellular signaling pathways in silico, aiding and guiding biological experiments to elucidate intracellular regulation. To construct such quantitative and predictive models of intracellular interactions, unknown parameters need to be estimated. Most of biological data are expressed in relative or arbitrary units, raising the question of how to compare model simulations with data. It has recently been shown that for models with large number of unknown parameters, fitting algorithms using a data-driven normalization of the simulations (DNS) performs best in terms of the convergence time and parameter identifiability. DNS approach compares model simulations and corresponding data both normalized by the same normalization procedure, without requiring additional parameters to be estimated, as necessary for widely used scaling factor-based methods. However, currently there is no parameter estimation software that directly supports DNS. In this chapter, we show how to apply DNS to dynamic models of systems and synthetic biology using PEPSSBI (Parameter Estimation Pipeline for Systems and Synthetic Biology). PEPSSBI is the first software that supports DNS, through algorithmically supported data normalization and objective function construction. PEPSSBI also supports model import using SBML and repeated parameter estimation runs executed in parallel either on a personal computer or a multi-CPU cluster.

RunBioSimulations: an extensible web application that simulates a wide range of computational modeling frameworks, algorithms, and formats

Comprehensive, predictive computational models have significant potential for science, bioengineering, and medicine. One promising way to achieve more predictive models is to combine submodels of multiple subsystems. To capture the multiple scales of biology, these submodels will likely require multiple modeling frameworks and simulation algorithms. Several community resources are already available for working with many of these frameworks and algorithms. However, the variety and sheer number of these resources make it challenging to find and use appropriate tools for each model, especially for novice modelers and experimentalists. To make these resources easier to use, we developed RunBioSimulations (https://run.biosimulations.org), a single web application for executing a broad range of models. RunBioSimulations leverages community resources, including BioSimulators, a new open registry of simulation tools. These resources currently enable RunBioSimulations to execute nine frameworks and 44 algorithms, and they make RunBioSimulations extensible to additional frameworks and algorithms. RunBioSimulations also provides features for sharing simulations and interactively visualizing their results. We anticipate that RunBioSimulations will foster reproducibility, stimulate collaboration, and ultimately facilitate the creation of more predictive models.

BiPSim: a flexible and generic stochastic simulator for polymerization processes

Detailed whole-cell modeling requires an integration of heterogeneous cell processes having different modeling formalisms, for which whole-cell simulation could remain tractable. Here, we introduce BiPSim, an open-source stochastic simulator of template-based polymerization processes, such as replication, transcription and translation. BiPSim combines an efficient abstract representation of reactions and a constant-time implementation of the Gillespie’s Stochastic Simulation Algorithm (SSA) with respect to reactions, which makes it highly efficient to simulate large-scale polymerization processes stochastically. Moreover, multi-level descriptions of polymerization processes can be handled simultaneously, allowing the user to tune a trade-off between simulation speed and model granularity. We evaluated the performance of BiPSim by simulating genome-wide gene expression in bacteria for multiple levels of granularity. Finally, since no cell-type specific information is hard-coded in the simulator, models can easily be adapted to other organismal species. We expect that BiPSim should open new perspectives for the genome-wide simulation of stochastic phenomena in biology.

A systems-biology model of the tumor necrosis factor (TNF) interactions with TNF receptor 1 and 2

MOTIVATION

Clustering enables TNF receptors to stimulate intracellular signaling. The differential soluble ligand-induced clustering behavior of TNF receptor 1 (TNFR1) and TNFR2 was modeled. A structured, rule-based model implemented ligand-independent pre-ligand binding assembly domain (PLAD)-mediated homotypic low affinity interactions of unliganded and liganded TNF receptors.

RESULTS

Soluble TNF initiates TNFR1 signaling but not TNFR2 signaling despite receptor binding unless it is secondarily oligomerized. We consider high affinity binding of TNF to signaling-incompetent pre-assembled dimeric TNFR1 and TNFR2 molecules and secondary clustering of liganded dimers to signaling competent ligand-receptor clusters. Published receptor numbers, affinities and measured different activities of clustered receptors validated model simulations for a large range of receptor and ligand concentrations. Different PLAD-PLAD affinities and different activities of receptor clusters explain the observed differences in the TNF receptor stimulating activities of soluble TNF.

AVAILABILITY AND IMPLEMENTATION

All scripts and data are in manuscript and supplement at Bioinformatics online.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Model-based optimization of combination protocols for irradiation-insensitive cancers

Alternations in the p53 regulatory network may render cancer cells resistant to the radiation-induced apoptosis. In this theoretical study we search for the best protocols combining targeted therapy with radiation to treat cancers with wild-type p53, but having downregulated expression of PTEN or overexpression of Wip1 resulting in resistance to radiation monotherapy. Instead of using the maximum tolerated dose paradigm, we exploit stochastic computational model of the p53 regulatory network to calculate apoptotic fractions for both normal and cancer cells. We consider combination protocols, with irradiations repeated every 12, 18, 24, or 36 h to find that timing between Mdm2 inhibitor delivery and irradiation significantly influences the apoptotic cell fractions. We assume that uptake of the inhibitor is higher by cancer than by normal cells and that cancer cells receive higher irradiation doses from intersecting beams. These two assumptions were found necessary for the existence of protocols inducing massive apoptosis in cancer cells without killing large fraction of normal cells neighboring tumor. The best found protocols have irradiations repeated every 24 or 36 h with two inhibitor doses per irradiation cycle, and allow to induce apoptosis in more than 95% of cancer cells, killing less than 10% of normal cells.

Interactions between calmodulin and neurogranin govern the dynamics of CaMKII as a leaky integrator

Calmodulin-dependent kinase II (CaMKII) has long been known to play an important role in learning and memory as well as long term potentiation (LTP). More recently it has been suggested that it might be involved in the time averaging of synaptic signals, which can then lead to the high precision of information stored at a single synapse. However, the role of the scaffolding molecule, neurogranin (Ng), in governing the dynamics of CaMKII is not yet fully understood. In this work, we adopt a rule-based modeling approach through the Monte Carlo method to study the effect of Ca2+ signals on the dynamics of CaMKII phosphorylation in the postsynaptic density (PSD). Calcium surges are observed in synaptic spines during an EPSP and back-propagating action potential due to the opening of NMDA receptors and voltage dependent calcium channels. Using agent-based models, we computationally investigate the dynamics of phosphorylation of CaMKII monomers and dodecameric holoenzymes. The scaffolding molecule, Ng, when present in significant concentration, limits the availability of free calmodulin (CaM), the protein which activates CaMKII in the presence of calcium. We show that Ng plays an important modulatory role in CaMKII phosphorylation following a surge of high calcium concentration. We find a non-intuitive dependence of this effect on CaM concentration that results from the different affinities of CaM for CaMKII depending on the number of calcium ions bound to the former. It has been shown previously that in the absence of phosphatase, CaMKII monomers integrate over Ca2+ signals of certain frequencies through autophosphorylation (Pepke et al, Plos Comp. Bio., 2010). We also study the effect of multiple calcium spikes on CaMKII holoenzyme autophosphorylation, and show that in the presence of phosphatase, CaMKII behaves as a leaky integrator of calcium signals, a result that has been recently observed in vivo. Our models predict that the parameters of this leaky integrator are finely tuned through the interactions of Ng, CaM, CaMKII, and PP1, providing a mechanism to precisely control the sensitivity of synapses to calcium signals. Author Summary not valid for PLOS ONE submissions.

Understanding the Dynamics of T-Cell Activation in Health and Disease Through the Lens of Computational Modeling

T cells in the immune system are activated by binding to foreign peptides (from an external pathogen) or mutant peptide (derived from endogenous proteins) displayed on the surface of a diseased cell. This triggers a series of intracellular signaling pathways, which ultimately dictate the response of the T cell. The insights from computational models have greatly improved our understanding of the mechanisms that control T-cell activation. In this review, we focus on the use of ordinary differential equation–based mechanistic models to study T-cell activation. We highlight several examples that demonstrate the models’ utility in answering specific questions related to T-cell activation signaling, from antigen discrimination to the feedback mechanisms that initiate transcription factor activation. In addition, we describe other modeling approaches that can be combined with mechanistic models to bridge time scales and better understand how intracellular signaling events, which occur on the order of seconds to minutes, influence phenotypic responses of T-cell activation, which occur on the order of hours to days. Overall, through concrete examples, we emphasize how computational modeling can be used to enable the rational design and optimization of immunotherapies.

A domain-level DNA strand displacement reaction enumerator allowing arbitrary non-pseudoknotted secondary structures

Information technologies enable programmers and engineers to design and synthesize systems of startling complexity that nonetheless behave as intended. This mastery of complexity is made possible by a hierarchy of formal abstractions that span from high-level programming languages down to low-level implementation specifications, with rigorous connections between the levels. DNA nanotechnology presents us with a new molecular information technology whose potential has not yet been fully unlocked in this way. Developing an effective hierarchy of abstractions may be critical for increasing the complexity of programmable DNA systems. Here, we build on prior practice to provide a new formalization of 'domain-level' representations of DNA strand displacement systems that has a natural connection to nucleic acid biophysics while still being suitable for formal analysis. Enumeration of unimolecular and bimolecular reactions provides a semantics for programmable molecular interactions, with kinetics given by an approximate biophysical model. Reaction condensation provides a tractable simplification of the detailed reactions that respects overall kinetic properties. The applicability and accuracy of the model is evaluated across a wide range of engineered DNA strand displacement systems. Thus, our work can serve as an interface between lower-level DNA models that operate at the nucleotide sequence level, and high-level chemical reaction network models that operate at the level of interactions between abstract species.

NERDSS: A Nonequilibrium Simulator for Multibody Self-Assembly at the Cellular Scale

Currently, a significant barrier to building predictive models of cellular self-assembly processes is that molecular models cannot capture minutes-long dynamics that couple distinct components with active processes, whereas reaction-diffusion models cannot capture structures of molecular assembly. Here, we introduce the nonequilibrium reaction-diffusion self-assembly simulator (NERDSS), which addresses this spatiotemporal resolution gap. NERDSS integrates efficient reaction-diffusion algorithms into generalized software that operates on user-defined molecules through diffusion, binding and orientation, unbinding, chemical transformations, and spatial localization. By connecting the fast processes of binding with the slow timescales of large-scale assembly, NERDSS integrates molecular resolution with reversible formation of ordered, multisubunit complexes. NERDSS encodes models using rule-based formatting languages to facilitate model portability, usability, and reproducibility. Applying NERDSS to steps in clathrin-mediated endocytosis, we design multicomponent systems that can form lattices in solution or on the membrane, and we predict how stochastic but localized dephosphorylation of membrane lipids can drive lattice disassembly. The NERDSS simulations reveal the spatial constraints on lattice growth and the role of membrane localization and cooperativity in nucleating assembly. By modeling viral lattice assembly and recapitulating oscillations in protein expression levels for a circadian clock model, we illustrate the adaptability of NERDSS. NERDSS simulates user-defined assembly models that were previously inaccessible to existing software tools, with broad applications to predicting self-assembly in vivo and designing high-yield assemblies in vitro.

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.

Setting the basis of best practices and standards for curation and annotation of logical models in biology-highlights of the [BC]2 2019 CoLoMoTo/SysMod Workshop

The fast accumulation of biological data calls for their integration, analysis and exploitation through more systematic approaches. The generation of novel, relevant hypotheses from this enormous quantity of data remains challenging. Logical models have long been used to answer a variety of questions regarding the dynamical behaviours of regulatory networks. As the number of published logical models increases, there is a pressing need for systematic model annotation, referencing and curation in community-supported and standardised formats. This article summarises the key topics and future directions of a meeting entitled 'Annotation and curation of computational models in biology', organised as part of the 2019 [BC]2 conference. The purpose of the meeting was to develop and drive forward a plan towards the standardised annotation of logical models, review and connect various ongoing projects of experts from different communities involved in the modelling and annotation of molecular biological entities, interactions, pathways and models. This article defines a roadmap towards the annotation and curation of logical models, including milestones for best practices and minimum standard requirements.

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.

Interactive Multiresolution Visualization of Cellular Network Processes

Visualization plays a central role in the analysis of biochemical network models to identify patterns that arise from reaction dynamics and perform model exploratory analysis. To facilitate these analyses, we developed PyViPR, a visualization tool that generates static and dynamic representations of biochemical network processes within a Python-based environment. PyViPR embeds network visualizations within Jupyter notebooks, thus enabling integration with modeling, simulation, and analysis workflows. To present the capabilities of PyViPR, we explore execution mechanisms of extrinsic apoptosis in HeLa cells. We show that community-detection algorithms identify groups of molecular species that capture key biological functions and ease exploration of the apoptosis network. We then show how different kinetic parameter sets that fit the experimental data equally well exhibit significantly different signal-execution dynamics as the system progresses toward mitochondrial outer-membrane permeabilization. Therefore, PyViPR aids the conceptual understanding of dynamic network processes and accelerates hypothesis generation for further testing and validation.

Potential based, spatial simulation of dynamically nested particles

BACKGROUND

To study cell biological phenomena which depend on diffusion, active transport processes, or the locations of species, modeling and simulation studies need to take space into account. To describe the system as a collection of discrete objects moving and interacting in continuous space, various particle-based reaction diffusion simulators for cell-biological system have been developed. So far the focus has been on particles as solid spheres or points. However, spatial dynamics might happen at different organizational levels, such as proteins, vesicles or cells with interrelated dynamics which requires spatial approaches that take this multi-levelness of cell biological systems into account.

RESULTS

Based on the perception of particles forming hollow spheres, ML-Force contributes to the family of particle-based simulation approaches: in addition to excluded volumes and forces, it also supports compartmental dynamics and relating dynamics between different organizational levels explicitly. Thereby, compartmental dynamics, e.g., particles entering and leaving other particles, and bimolecular reactions are modeled using pair-wise potentials (forces) and the Langevin equation. In addition, forces that act independently of other particles can be applied to direct the movement of particles. Attributes and the possibility to define arbitrary functions on particles, their attributes and content, to determine the results and kinetics of reactions add to the expressiveness of ML-Force. Its implementation comprises a rudimentary rule-based embedded domain-specific modeling language for specifying models and a simulator for executing models continuously. Applications inspired by cell biological models from literature, such as vesicle transport or yeast growth, show the value of the realized features. They facilitate capturing more complex spatial dynamics, such as the fission of compartments or the directed movement of particles, and enable the integration of non-spatial intra-compartmental dynamics as stochastic events.

CONCLUSIONS

By handling all dynamics based on potentials (forces) and the Langevin equation, compartmental dynamics, such as dynamic nesting, fusion and fission of compartmental structures are handled continuously and are seamlessly integrated with traditional particle-based reaction-diffusion dynamics within the cell. Thereby, attributes and arbitrary functions allow to flexibly describe diverse spatial phenomena, and relate dynamics across organizational levels. Also they prove crucial in modeling intra-cellular or intra-compartmental dynamics in a non-spatial manner, and, thus, to abstract from spatial dynamics, on demand which increases the range of multi-compartmental processes that can be captured.

ModelBricks-modules for reproducible modeling improving model annotation and provenance

Most computational models in biology are built and intended for "single-use"; the lack of appropriate annotation creates models where the assumptions are unknown, and model elements are not uniquely identified. Simply recreating a simulation result from a publication can be daunting; expanding models to new and more complex situations is a herculean task. As a result, new models are almost always created anew, repeating literature searches for kinetic parameters, initial conditions and modeling specifics. It is akin to building a brick house starting with a pile of clay. Here we discuss a concept for building annotated, reusable models, by starting with small well-annotated modules we call ModelBricks. Curated ModelBricks, accessible through an open database, could be used to construct new models that will inherit ModelBricks annotations and thus be easier to understand and reuse. Key features of ModelBricks include reliance on a commonly used standard language (SBML), rule-based specification describing species as a collection of uniquely identifiable molecules, association with model specific numerical parameters, and more common annotations. Physical bricks can vary substantively; likewise, to be useful the structure of ModelBricks must be highly flexible-it should encapsulate mechanisms from single reactions to multiple reactions in a complex process. Ultimately, a modeler would be able to construct large models by using multiple ModelBricks, preserving annotations and provenance of model elements, resulting in a highly annotated model. We envision the library of ModelBricks to rapidly grow from community contributions. Persistent citable references will incentivize model creators to contribute new ModelBricks.

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.

Mechanistic Models of Cellular Signaling, Cytokine Crosstalk, and Cell-Cell Communication in Immunology

The cells of the immune system respond to a great variety of different signals that frequently reach them simultaneously. Computational models of signaling pathways and cellular behavior can help us explore the biochemical mechanisms at play during such responses, in particular when those models aim at incorporating molecular details of intracellular reaction networks. Such detailed models can encompass hypotheses about the interactions among molecular binding domains and how these interactions are modulated by, for instance, post-translational modifications, or steric constraints in multi-molecular complexes. In this way, the models become formal representations of mechanistic immunological hypotheses that can be tested through quantitative simulations. Due to the large number of parameters (molecular abundances, association-, dissociation-, and enzymatic transformation rates) the goal of simulating the models can, however, in many cases no longer be the fitting of particular parameter values. Rather, the simulations perform sweeps through parameter space to test whether a model can account for certain experimentally observed features when allowing the parameter values to vary within experimentally determined or physiologically reasonable ranges. We illustrate how this approach can be used to explore possible mechanisms of immunological pathway crosstalk. Probing the input-output behavior of mechanistic pathway models through systematic simulated variations of receptor stimuli will soon allow us to derive cell population behavior from single-cell models, thereby bridging a scale gap that currently still is frequently addressed through heuristic phenomenological multi-scale models.

Rule-Based Modeling Using Wildcards in the Smoldyn Simulator

Many biological molecules exist in multiple variants, such as proteins with different posttranslational modifications, DNAs with different sequences, and phospholipids with different chain lengths. Representing these variants as distinct species, as most biochemical simulators do, leads to the problem that the number of species, and chemical reactions that interconvert them, typically increase combinatorially with the number of ways that the molecules can vary. This can be alleviated by "rule-based modeling methods," in which software generates the chemical reaction network from relatively simple "rules." This chapter presents a new approach to rule-based modeling. It is based on wildcards that match to species names, much as wildcards can match to file names in computer operating systems. It is much simpler to use than the formal rule-based modeling approaches developed previously but can lead to unintended consequences if not used carefully. This chapter demonstrates rule-based modeling with wildcards through examples for signaling systems, protein complexation, polymerization, nucleic acid sequence copying and mutation, the "SMILES" chemical notation, and others. The method is implemented in Smoldyn, a spatial and stochastic biochemical simulator, for both generate-first and on-the-fly expansion, meaning whether the reaction network is generated before or during the simulation.

Exploring the Extracellular Regulation of the Tumor Angiogenic Interaction Network Using a Systems Biology Model

Tumor angiogenesis is regulated by pro- and anti-angiogenic factors. Anti-angiogenic agents target the interconnected network of angiogenic factors to inhibit neovascularization, which subsequently impedes tumor growth. Due to the complexity of this network, optimizing anti-angiogenic cancer treatments requires detailed knowledge at a systems level. In this study, we constructed a tumor tissue-based model to better understand how the angiogenic network is regulated by opposing mediators at the extracellular level. We consider the network comprised of two pro-angiogenic factors: vascular endothelial growth factor (VEGF) and basic fibroblast growth factor (FGF2), and two anti-angiogenic factors: thrombospondin-1 (TSP1) and platelet factor 4 (PF4). The model's prediction of angiogenic factors' distribution in tumor tissue reveals the localization of different factors and indicates the angiogenic state of the tumor. We explored how the distributions are affected by the secretion of the pro- and anti-angiogenic factors, illustrating how the angiogenic network is regulated in the extracellular space. Interestingly, we identified a counterintuitive result that the secretion of the anti-angiogenic factor PF4 can enhance pro-angiogenic signaling by elevating the levels of the interstitial and surface-level pro-angiogenic species. This counterintuitive situation is pertinent to the clinical setting, such as the release of anti-angiogenic factors in platelet activation or the administration of exogenous PF4 for anti-angiogenic therapy. Our study provides mechanistic insights into this counterintuitive result and highlights the role of heparan sulfate proteoglycans in regulating the interactions between angiogenic factors. This work complements previous studies aimed at understanding the formation of angiogenic complexes in tumor tissue and helps in the development of anti-cancer strategies targeting angiogenesis.

Scaling methods for accelerating kinetic Monte Carlo simulations of chemical reaction networks

Various kinetic Monte Carlo algorithms become inefficient when some of the population sizes in a system are large, which gives rise to a large number of reaction events per unit time. Here, we present a new acceleration algorithm based on adaptive and heterogeneous scaling of reaction rates and stoichiometric coefficients. The algorithm is conceptually related to the commonly used idea of accelerating a stochastic simulation by considering a subvolume λΩ (0 < λ < 1) within a system of interest, which reduces the number of reaction events per unit time occurring in a simulation by a factor 1/λ at the cost of greater error in unbiased estimates of first moments and biased overestimates of second moments. Our new approach offers two unique benefits. First, scaling is adaptive and heterogeneous, which eliminates the pitfall of overaggressive scaling. Second, there is no need for an a priori classification of populations as discrete or continuous (as in a hybrid method), which is problematic when discreteness of a chemical species changes during a simulation. The method requires specification of only a single algorithmic parameter, N_c, a global critical population size above which populations are effectively scaled down to increase simulation efficiency. The method, which we term partial scaling, is implemented in the open-source BioNetGen software package. We demonstrate that partial scaling can significantly accelerate simulations without significant loss of accuracy for several published models of biological systems. These models characterize activation of the mitogen-activated protein kinase ERK, prion protein aggregation, and T-cell receptor signaling.

Aanat Knockdown and Melatonin Supplementation in Embryo Development: Involvement of Mitochondrial Function and DNA Methylation

Aims: In addition to pineal gland, many cells, tissues, and organs also synthesize melatonin (N-acetyl-5-methoxytryptamine). Embryos are a group of special cells and whether they can synthesize melatonin is still an open question. However, melatonin application promoted embryo development in many species in in vitro condition. The purpose of this study was to investigate whether embryos can synthesize melatonin; if it is so, what are the impacts of the endogenously produced melatonin on embryo development and the associated molecular mechanisms. These have never been reported previously. Results: Melatonin synthesis was observed at different stages of embryonic development. Aanat (aralkylamine N-acetyltransferase), a rate-limiting enzyme for melatonin production, was found to mostly localize in the mitochondria. Aanat knockdown significantly impeded embryonic development, and melatonin supplementation rescued it. The potential mechanisms might be that melatonin preserved mitochondrial intact and its function, thus providing sufficient adenosine 5'-triphosphate for the embryo development. In addition, melatonin scavenged intracellular reactive oxygen species (ROS) and reduced the DNA mutation induced by oxidative stress. In the molecular level, Aanat knockdown reduced tet methylcytosine dioxygenase 2 (Tet2) expression and DNA demethylation in blastocyst and melatonin supplementation rescued these processes. Innovation: This is the first report to show that embryos synthesize melatonin, and its synthetic enzyme Aanat was located in the mitochondria of embryos. An effect of melatonin is to maintain Tet2 expression and normal methylation status, and thereby promote embryonic development. Conclusion: Embryos can produce melatonin that reduces ROS production, preserves mitochondrial function, and maintains Tet2 expression and the normal DNA methylation.

Using RuleBuilder to Graphically Define and Visualize BioNetGen-Language Patterns and Reaction Rules

RuleBuilder is a tool for drawing graphs that can be represented by the BioNetGen language (BNGL), which is used to formulate mathematical, rule-based models of biochemical systems. BNGL provides an intuitive plain text, or string, representation of such systems, which is based on a graphical formalism. Reactions are defined in terms of graph-rewriting rules that specify the necessary intrinsic properties of the reactants, a transformation, and a rate law. Rules also contain contextual constraints that restrict application of the rule. In some cases, the specification of contextual constraints can be verbose, making a rule difficult to read. RuleBuilder is designed to ease the task of reading and writing individual reaction rules or other BNGL patterns required for model formulation. The software assists in the reading of existing models by converting BNGL strings of interest into a graph-based representation composed of nodes and edges. RuleBuilder also enables the user to construct de novo a visual representation of BNGL strings using drawing tools available in its interface. As objects are added to the drawing canvas, the corresponding BNGL string is generated on the fly, and objects are similarly drawn on the fly as BNGL strings are entered into the application. RuleBuilder thus facilitates construction and interpretation of rule-based models.

A nonlinear mathematical model of cell-mediated immune response for tumor phenotypic heterogeneity

Human cancers display intra-tumor heterogeneity in many phenotypic features, such as expression of cell surface receptors, growth, and angiogenic, proliferative, and immunogenic factors, which represent obstacles to a successful immune response. In this paper, we propose a nonlinear mathematical model of cancer immunosurveillance that takes into account some of these features based on cell-mediated immune responses. The model describes phenomena that are seen in vivo, such as tumor dormancy, robustness, immunoselection over tumor heterogeneity (also called "cancer immunoediting") and strong sensitivity to initial conditions in the composition of tumor microenvironment. The results framework has as common element the tumor as an attractor for abnormal cells. Bifurcation analysis give us as tumor attractors fixed-points, limit cycles and chaotic attractors, the latter emerging from period-doubling cascade displaying Feigenbaum's universality. Finally, we simulated both elimination and escape tumor scenarios by means of a stochastic version of the model according to the Doob-Gillespie algorithm.

Limits to the rate of information transmission through the MAPK pathway

Two important signalling pathways of NF-κB and ERK transmit merely 1 bit of information about the level of extracellular stimulation. It is thus unclear how such systems can coordinate complex cell responses to external cues. We analyse information transmission in the MAPK/ERK pathway that converts both constant and pulsatile EGF stimulation into pulses of ERK activity. Based on an experimentally verified computational model, we demonstrate that, when input consists of sequences of EGF pulses, transmitted information increases nearly linearly with time. Thus, pulse-interval transcoding allows more information to be relayed than the amplitude–amplitude transcoding considered previously for the ERK and NF-κB pathways. Moreover, the information channel capacity C, or simply bitrate, is not limited by the bandwidth B = 1/τ, where τ ≈ 1 h is the relaxation time. Specifically, when the input is provided in the form of sequences of short binary EGF pulses separated by intervals that are multiples of τ/n (but not shorter than τ), then for n = 2, C ≈ 1.39 bit h⁻¹; and for n = 4, C ≈ 1.86 bit h⁻¹. The capability to respond to random sequences of EGF pulses enables cells to propagate spontaneous ERK activity waves across tissue.

The impact of tumor receptor heterogeneity on the response to anti-angiogenic cancer treatment

Multiple promoters and inhibitors mediate angiogenesis, the formation of new blood vessels, and these factors represent potential targets for impeding vessel growth in tumors. Vascular endothelial growth factor (VEGF) is a potent angiogenic factor targeted in anti-angiogenic cancer therapies. In addition, thrombospondin-1 (TSP1) is a major endogenous inhibitor of angiogenesis, and TSP1 mimetics are being developed as an alternative type of anti-angiogenic agent. The combination of bevacizumab, an anti-VEGF agent, and ABT-510, a TSP1 mimetic, has been tested in clinical trials to treat advanced solid tumors. However, the patients' responses are highly variable and show disappointing outcomes. To obtain mechanistic insight into the effects of this combination anti-angiogenic therapy, we have constructed a novel whole-body systems biology model including the VEGF and TSP1 reaction networks. Using this molecular-detailed model, we investigated how the combination anti-angiogenic therapy changes the amounts of pro-angiogenic and anti-angiogenic complexes in cancer patients. We particularly focus on answering the question of how the effect of the combination therapy is influenced by tumor receptor expression, one aspect of patient-to-patient variability. Overall, this model complements the clinical administration of combination anti-angiogenic therapy, highlights the role of tumor receptor variability in the heterogeneous responses to anti-angiogenic therapy, and identifies the tumor receptor profiles that correlate with a high likelihood of a positive response to the combination therapy. Our model provides novel understanding of the VEGF-TSP1 balance in cancer patients at the systems-level and could be further used to optimize combination anti-angiogenic therapy.

Harmonizing semantic annotations for computational models in biology

Life science researchers use computational models to articulate and test hypotheses about the behavior of biological systems. Semantic annotation is a critical component for enhancing the interoperability and reusability of such models as well as for the integration of the data needed for model parameterization and validation. Encoded as machine-readable links to knowledge resource terms, semantic annotations describe the computational or biological meaning of what models and data represent. These annotations help researchers find and repurpose models, accelerate model composition and enable knowledge integration across model repositories and experimental data stores. However, realizing the potential benefits of semantic annotation requires the development of model annotation standards that adhere to a community-based annotation protocol. Without such standards, tool developers must account for a variety of annotation formats and approaches, a situation that can become prohibitively cumbersome and which can defeat the purpose of linking model elements to controlled knowledge resource terms. Currently, no consensus protocol for semantic annotation exists among the larger biological modeling community. Here, we report on the landscape of current annotation practices among the COmputational Modeling in BIology NEtwork community and provide a set of recommendations for building a consensus approach to semantic annotation.