Simulation of large-scale rule-based models

A query-response causal analysis of reaction events in biochemical reaction networks

The stochastic kinetics of biochemical reaction networks is described by a chemical master equation (CME) and the underlying laws of mass action. Assuming network-free simulations of the rule-based models of biochemical reaction networks (BRNs), this paper departs from the usual analysis of network dynamics as the time-dependent distributions of chemical species counts, and instead considers statistically evaluating the sequences of reaction events generated from the stochastic simulations. The reaction event-time series can be used for reaction clustering, identifying rare events, and recognizing the periods of increased or steady-state activity. However, the main aim of this paper is to device an effective method for identifying causally and anti-causally related sub-sequences of reaction events using their empirical probabilities. This allows discovering some of the causal dynamics of BRNs as well as uncovering their short-term deterministic behaviors. In particular, it is proposed that the reaction sub-sequences that are conditionally nearly certain or nearly uncertain can be considered as being causally related. Moreover, since the time-ordering of reaction events is locally irrelevant, the reaction sub-sequences can be transformed into the reaction sets or multi-sets. The distance metrics can be then used to define the equivalences among the reaction events. The proposed method for identifying the causally related reaction sub-sequences has been implemented as a computationally efficient query-response mechanism. The method was evaluated for five models of genetic networks in seven defined numerical experiments. The models were simulated in BioNetGen using the open-source network-free simulator NFsim. This simulator had to be modified first to allow recording the traces of reaction events, and it is available in the Github repository, ploskot/nfsim_1.20. The generated event time-series were analyzed with Python and Matlab scripts. The whole process of data generation, analysis and visualization has been nearly fully automated using shell scripts. This demonstrates the opportunities for substantially increasing the research productivity by creating automated data generation and processing pipelines.

Catalyst: Fast and flexible modeling of reaction networks

We introduce Catalyst.jl, a flexible and feature-filled Julia library for modeling and high-performance simulation of chemical reaction networks (CRNs). Catalyst supports simulating stochastic chemical kinetics (jump process), chemical Langevin equation (stochastic differential equation), and reaction rate equation (ordinary differential equation) representations for CRNs. Through comprehensive benchmarks, we demonstrate that Catalyst simulation runtimes are often one to two orders of magnitude faster than other popular tools. More broadly, Catalyst acts as both a domain-specific language and an intermediate representation for symbolically encoding CRN models as Julia-native objects. This enables a pipeline of symbolically specifying, analyzing, and modifying CRNs; converting Catalyst models to symbolic representations of concrete mathematical models; and generating compiled code for numerical solvers. Leveraging ModelingToolkit.jl and Symbolics.jl, Catalyst models can be analyzed, simplified, and compiled into optimized representations for use in numerical solvers. Finally, we demonstrate Catalyst's broad extensibility and composability by highlighting how it can compose with a variety of Julia libraries, and how existing open-source biological modeling projects have extended its intermediate representation.

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.

Applications of personalised signalling network models in precision oncology

As our ability to provide in-depth, patient-specific characterisation of the molecular alterations within tumours rapidly improves, it is becoming apparent that new approaches will be required to leverage the power of this data and derive the full benefit for each individual patient. Systems biology approaches are beginning to emerge within this field as a potential method of incorporating large volumes of network level data and distilling a coherent, clinically-relevant prediction of drug response. However, the initial promise of this developing field is yet to be realised. Here we argue that in order to develop these precise models of individual drug response and tailor treatment accordingly, we will need to develop mathematical models capable of capturing both the dynamic nature of drug-response signalling networks and key patient-specific information such as mutation status or expression profiles. We also review the modelling approaches commonly utilised within this field, and outline recent examples of their use in furthering the application of systems biology for a precision medicine approach to cancer treatment.

Strategy-Driven Exploration for Rule-Based Models of Biochemical Systems with PORGY

This chapter presents PORGY-an interactive visual environment for rule-based modelling of biochemical systems. We model molecules and molecule interactions as port graphs and port graph rewrite rules, respectively. We use rewriting strategies to control which rules to apply, and where and when to apply them. Our main contributions to rule-based modelling of biochemical systems lie in the strategy language and the associated visual and interactive features offered by PORGY. These features facilitate an exploratory approach to test different ways of applying the rules while recording the model evolution, and tracking and plotting parameters. We illustrate PORGY's features with a study of the role of a scaffold protein in RAF/MEK/ERK signalling.

Generalizing Gillespie's Direct Method to Enable Network-Free Simulations

Gillespie's direct method for stochastic simulation of chemical kinetics is a staple of computational systems biology research. However, the algorithm requires explicit enumeration of all reactions and all chemical species that may arise in the system. In many cases, this is not feasible due to the combinatorial explosion of reactions and species in biological networks. Rule-based modeling frameworks provide a way to exactly represent networks containing such combinatorial complexity, and generalizations of Gillespie's direct method have been developed as simulation engines for rule-based modeling languages. Here, we provide both a high-level description of the algorithms underlying the simulation engines, termed network-free simulation algorithms, and how they have been applied in systems biology research. We also define a generic rule-based modeling framework and describe a number of technical details required for adapting Gillespie's direct method for network-free simulation. Finally, we briefly discuss potential avenues for advancing network-free simulation and the role they continue to play in modeling dynamical systems in biology.

An Overview of Network-Based and -Free Approaches for Stochastic Simulation of Biochemical Systems

Stochastic simulation has been widely used to model the dynamics of biochemical reaction networks. Several algorithms have been proposed that are exact solutions of the chemical master equation, following the work of Gillespie. These stochastic simulation approaches can be broadly classified into two categories: network-based and -free simulation. The network-based approach requires that the full network of reactions be established at the start, while the network-free approach is based on reaction rules that encode classes of reactions, and by applying rule transformations, it generates reaction events as they are needed without ever having to derive the entire network. In this study, we compare the efficiency and limitations of several available implementations of these two approaches. The results allow for an informed selection of the implementation and methodology for specific biochemical modeling applications.

Connecting the dots across time: reconstruction of single-cell signalling trajectories using time-stamped data

Single-cell responses are shaped by the geometry of signalling kinetic trajectories carved in a multidimensional space spanned by signalling protein abundances. It is, however, challenging to assay a large number (more than 3) of signalling species in live-cell imaging, which makes it difficult to probe single-cell signalling kinetic trajectories in large dimensions. Flow and mass cytometry techniques can measure a large number (4 to more than 40) of signalling species but are unable to track single cells. Thus, cytometry experiments provide detailed time-stamped snapshots of single-cell signalling kinetics. Is it possible to use the time-stamped cytometry data to reconstruct single-cell signalling trajectories? Borrowing concepts of conserved and slow variables from non-equilibrium statistical physics we develop an approach to reconstruct signalling trajectories using snapshot data by creating new variables that remain invariant or vary slowly during the signalling kinetics. We apply this approach to reconstruct trajectories using snapshot data obtained from in silico simulations, live-cell imaging measurements, and, synthetic flow cytometry datasets. The application of invariants and slow variables to reconstruct trajectories provides a radically different way to track objects using snapshot data. The approach is likely to have implications for solving matching problems in a wide range of disciplines.

Symbolic computation of differential equivalences

Ordinary differential equations (ODEs) are widespread in many natural sciences including chemistry, ecology, and systems biology, and in disciplines such as control theory and electrical engineering. Building on the celebrated molecules-as-processes paradigm, they have become increasingly popular in computer science, with high-level languages and formal methods such as Petri nets, process algebra, and rule-based systems that are interpreted as ODEs. We consider the problem of comparing and minimizing ODEs automatically. Influenced by traditional approaches in the theory of programming, we propose differential equivalence relations. We study them for a basic intermediate language, for which we have decidability results, that can be targeted by a class of high-level specifications. An ODE implicitly represents an uncountable state space, hence reasoning techniques cannot be borrowed from established domains such as probabilistic programs with finite-state Markov chain semantics. We provide novel symbolic procedures to check an equivalence and compute the largest one via partition refinement algorithms that use satisfiability modulo theories. We illustrate the generality of our framework by showing that differential equivalences include (i) well-known notions for the minimization of continuous-time Markov chains (lumpability), (ii)~bisimulations for chemical reaction networks recently proposed by Cardelli et al., and (iii) behavioral relations for process algebra with ODE semantics. With a prototype implementation we are able to detect equivalences in biochemical models from the literature that cannot be reduced using competing automatic techniques.

Modeling for (physical) biologists: an introduction to the rule-based approach

Models that capture the chemical kinetics of cellular regulatory networks can be specified in terms of rules for biomolecular interactions. A rule defines a generalized reaction, meaning a reaction that permits multiple reactants, each capable of participating in a characteristic transformation and each possessing certain, specified properties, which may be local, such as the state of a particular site or domain of a protein. In other words, a rule defines a transformation and the properties that reactants must possess to participate in the transformation. A rule also provides a rate law. A rule-based approach to modeling enables consideration of mechanistic details at the level of functional sites of biomolecules and provides a facile and visual means for constructing computational models, which can be analyzed to study how system-level behaviors emerge from component interactions.

Modeling biomolecular site dynamics in immunoreceptor signaling systems

The immune system plays a central role in human health. The activities of immune cells, whether defending an organism from disease or triggering a pathological condition such as autoimmunity, are driven by the molecular machinery of cellular signaling systems. Decades of experimentation have elucidated many of the biomolecules and interactions involved in immune signaling and regulation, and recently developed technologies have led to new types of quantitative, systems-level data. To integrate such information and develop nontrivial insights into the immune system, computational modeling is needed, and it is essential for modeling methods to keep pace with experimental advances. In this chapter, we focus on the dynamic, site-specific, and context-dependent nature of interactions in immunoreceptor signaling (i.e., the biomolecular site dynamics of immunoreceptor signaling), the challenges associated with capturing these details in computational models, and how these challenges have been met through use of rule-based modeling approaches.

Multi-state modeling of biomolecules

Multi-state modeling of biomolecules refers to a series of techniques used to represent and compute the behavior of biological molecules or complexes that can adopt a large number of possible functional states. Biological signaling systems often rely on complexes of biological macromolecules that can undergo several functionally significant modifications that are mutually compatible. Thus, they can exist in a very large number of functionally different states. Modeling such multi-state systems poses two problems: the problem of how to describe and specify a multi-state system (the “specification problem”) and the problem of how to use a computer to simulate the progress of the system over time (the “computation problem”). To address the specification problem, modelers have in recent years moved away from explicit specification of all possible states and towards rule-based formalisms that allow for implicit model specification, including the κ-calculus [1], BioNetGen [2]–[5], the Allosteric Network Compiler [6], and others [7], [8]. To tackle the computation problem, they have turned to particle-based methods that have in many cases proved more computationally efficient than population-based methods based on ordinary differential equations, partial differential equations, or the Gillespie stochastic simulation algorithm[9], [10]. Given current computing technology, particle-based methods are sometimes the only possible option. Particle-based simulators fall into two further categories: nonspatial simulators, such as StochSim [11], DYNSTOC [12], RuleMonkey [9], [13], and the Network-Free Stochastic Simulator (NFSim) [14], and spatial simulators, including Meredys [15], SRSim [16], [17], and MCell [18]–[20]. Modelers can thus choose from a variety of tools, the best choice depending on the particular problem. Development of faster and more powerful methods is ongoing, promising the ability to simulate ever more complex signaling processes in the future.

Analysis of Wnt signaling β-catenin spatial dynamics in HEK293T cells

Background

Wnt/β-catenin signaling is involved in different stages of mammalian development and implicated in various cancers (e.g. colorectal cancer). Recent experimental and computational studies have revealed characteristics of the pathway, however a cell-specific spatial perspective is lacking. In this study, a novel 3D confocal quantitation protocol is developed to acquire spatial (two cellular compartments: nucleus and cytosol-membrane) and temporal quantitative data on target protein (e.g. β-catenin) concentrations in Human Epithelial Kidney cells (HEK293T) during perturbation (with either cycloheximide or Wnt3A). Computational models of the Wnt pathway are constructed and interrogated based on this data.

Results

A single compartment Wnt pathway model is compared with a simple β-catenin two compartment model to investigate Wnt3A signaling in HEK293T cells. When protein synthesis is inhibited, β-catenin decreases at the same rate in both cellular compartments, suggesting diffusional transport is fast compared to β-catenin degradation in the cytosol. With Wnt3A stimulation, the total amount of β-catenin rises throughout the cell, however the increase is initially (~first hour) faster in the nuclear compartment. While both models were able to reproduce the whole cell changes in β-catenin, only the compartment model reproduced the Wnt3A induced changes in β-catenin distribution and it was also the best fit for the data obtained when active transport was included alongside passive diffusion transport.

Conclusions

This integrated 3D quantitation imaging protocol and computational modeling approach allowed cell-specific compartment models of the signaling pathways to be constructed and analyzed. The Wnt models constructed in this study are the first for HEK293T and have suggested potential roles of inter-compartment transport to the dynamics of signaling.

Exact hybrid particle/population simulation of rule-based models of biochemical systems

Detailed modeling and simulation of biochemical systems is complicated by the problem of combinatorial complexity, an explosion in the number of species and reactions due to myriad protein-protein interactions and post-translational modifications. Rule-based modeling overcomes this problem by representing molecules as structured objects and encoding their interactions as pattern-based rules. This greatly simplifies the process of model specification, avoiding the tedious and error prone task of manually enumerating all species and reactions that can potentially exist in a system. From a simulation perspective, rule-based models can be expanded algorithmically into fully-enumerated reaction networks and simulated using a variety of network-based simulation methods, such as ordinary differential equations or Gillespie's algorithm, provided that the network is not exceedingly large. Alternatively, rule-based models can be simulated directly using particle-based kinetic Monte Carlo methods. This “network-free” approach produces exact stochastic trajectories with a computational cost that is independent of network size. However, memory and run time costs increase with the number of particles, limiting the size of system that can be feasibly simulated. Here, we present a hybrid particle/population simulation method that combines the best attributes of both the network-based and network-free approaches. The method takes as input a rule-based model and a user-specified subset of species to treat as population variables rather than as particles. The model is then transformed by a process of “partial network expansion” into a dynamically equivalent form that can be simulated using a population-adapted network-free simulator. The transformation method has been implemented within the open-source rule-based modeling platform BioNetGen, and resulting hybrid models can be simulated using the particle-based simulator NFsim. Performance tests show that significant memory savings can be achieved using the new approach and a monetary cost analysis provides a practical measure of its utility.

Rule-based modeling is a modeling paradigm that addresses the problem of combinatorial complexity in biochemical systems. The key idea is to specify only those components of a biological macromolecule that are directly involved in a biochemical transformation. Until recently, this “pattern-based” approach greatly simplified the process of model building but did nothing to improve the performance of model simulation. This changed with the introduction of “network-free” simulation methods, which operate directly on the compressed rule set of a rule-based model rather than on a fully-enumerated set of reactions and species. However, these methods represent every molecule in a system as a particle, limiting their use to systems containing less than a few million molecules. Here, we describe an extension to the network-free approach that treats rare, complex species as particles and plentiful, simple species as population variables, while retaining the exact dynamics of the model system. By making more efficient use of computational resources for species that do not require the level of detail of a particle representation, this hybrid particle/population approach can simulate systems much larger than is possible using network-free methods and is an important step towards realizing the practical simulation of detailed, mechanistic models of whole cells.

Rule-based modeling: a computational approach for studying biomolecular site dynamics in cell signaling systems

Rule-based modeling was developed to address the limitations of traditional approaches for modeling chemical kinetics in cell signaling systems. These systems consist of multiple interacting biomolecules (e.g., proteins), which themselves consist of multiple parts (e.g., domains, linear motifs, and sites of phosphorylation). Consequently, biomolecules that mediate information processing generally have the potential to interact in multiple ways, with the number of possible complexes and posttranslational modification states tending to grow exponentially with the number of binary interactions considered. As a result, only large reaction networks capture all possible consequences of the molecular interactions that occur in a cell signaling system, which is problematic because traditional modeling approaches for chemical kinetics (e.g., ordinary differential equations) require explicit network specification. This problem is circumvented through representation of interactions in terms of local rules. With this approach, network specification is implicit and model specification is concise. Concise representation results in a coarse graining of chemical kinetics, which is introduced because all reactions implied by a rule inherit the rate law associated with that rule. Coarse graining can be appropriate if interactions are modular, and the coarseness of a model can be adjusted as needed. Rules can be specified using specialized model-specification languages, and recently developed tools designed for specification of rule-based models allow one to leverage powerful software engineering capabilities. A rule-based model comprises a set of rules, which can be processed by general-purpose simulation and analysis tools to achieve different objectives (e.g., to perform either a deterministic or stochastic simulation).

A predictive mathematical model of the DNA damage G2 checkpoint

A predictive mathematical model of the transition from the G2 phase in the cell cycle to mitosis (M) was constructed from the known interactions of the proteins that are thought to play significant roles in the G2 to M transition as well as the DNA damage- induced G2 checkpoint. The model simulates the accumulation of active cyclin B1/Cdk1 (MPF) complexes in the nucleus to activate mitosis, the inhibition of this process by DNA damage, and transport of component proteins between cytoplasm and nucleus. Interactions in the model are based on activities of individual phospho-epitopes and binding sites of proteins involved in G2/M. Because tracking phosphoforms leads to combinatorial explosion, we employ a rule-based approach using the BioNetGen software. The model was used to determine the effects of depletion or over-expression of selected proteins involved in the regulation of the G2 to M transition in the presence and absence of DNA damage. Depletion of Plk1 delayed mitotic entry and recovery from the DNA damage-induced G2 arrest and over-expression of MPF attenuated the DNA damage-induced G2 delay. The model recapitulates the G2 delay observed in the biological response to varying levels of a DNA damage signal. The model produced the novel prediction that depletion of pkMyt1 results in an abnormal biological state in which G2 cells with DNA damage accumulate inactive nuclear MPF. Such a detailed model may prove useful for predicting DNA damage G2 checkpoint function in cancer and, therefore, sensitivity to cancer therapy.

Leveraging modeling approaches: reaction networks and rules

We have witnessed an explosive growth in research involving mathematical models and computer simulations of intracellular molecular interactions, ranging from metabolic pathways to signaling and gene regulatory networks. Many software tools have been developed to aid in the study of such biological systems, some of which have a wealth of features for model building and visualization, and powerful capabilities for simulation and data analysis. Novel high-resolution and/or high-throughput experimental techniques have led to an abundance of qualitative and quantitative data related to the spatiotemporal distribution of molecules and complexes, their interactions kinetics, and functional modifications. Based on this information, computational biology researchers are attempting to build larger and more detailed models. However, this has proved to be a major challenge. Traditionally, modeling tools require the explicit specification of all molecular species and interactions in a model, which can quickly become a major limitation in the case of complex networks - the number of ways biomolecules can combine to form multimolecular complexes can be combinatorially large. Recently, a new breed of software tools has been created to address the problems faced when building models marked by combinatorial complexity. These have a different approach for model specification, using reaction rules and species patterns. Here we compare the traditional modeling approach with the new rule-based methods. We make a case for combining the capabilities of conventional simulation software with the unique features and flexibility of a rule-based approach in a single software platform for building models of molecular interaction networks.

Specification, annotation, visualization and simulation of a large rule-based model for ERBB receptor signaling

BACKGROUND

Mathematical/computational models are needed to understand cell signaling networks, which are complex. Signaling proteins contain multiple functional components and multiple sites of post-translational modification. The multiplicity of components and sites of modification ensures that interactions among signaling proteins have the potential to generate myriad protein complexes and post-translational modification states. As a result, the number of chemical species that can be populated in a cell signaling network, and hence the number of equations in an ordinary differential equation model required to capture the dynamics of these species, is prohibitively large. To overcome this problem, the rule-based modeling approach has been developed for representing interactions within signaling networks efficiently and compactly through coarse-graining of the chemical kinetics of molecular interactions.

RESULTS

Here, we provide a demonstration that the rule-based modeling approach can be used to specify and simulate a large model for ERBB receptor signaling that accounts for site-specific details of protein-protein interactions. The model is considered large because it corresponds to a reaction network containing more reactions than can be practically enumerated. The model encompasses activation of ERK and Akt, and it can be simulated using a network-free simulator, such as NFsim, to generate time courses of phosphorylation for 55 individual serine, threonine, and tyrosine residues. The model is annotated and visualized in the form of an extended contact map.

CONCLUSIONS

With the development of software that implements novel computational methods for calculating the dynamics of large-scale rule-based representations of cellular signaling networks, it is now possible to build and analyze models that include a significant fraction of the protein interactions that comprise a signaling network, with incorporation of the site-specific details of the interactions. Modeling at this level of detail is important for understanding cellular signaling.

A quick guide to biomolecular network studies: construction, analysis, applications, and resources

Over the past decade, a rapid increase in network data including signaling, transcription regulation, metabolic reaction, protein-protein interaction and genetic interaction has been observed. Many biology issues have been investigated by analyzing these diverse networks, providing new insights into biology. Networks also play an important role in disease studies including disease gene screening and clinical diagnosis. Large amounts of databases and software have been developed to facilitate the storage, exchange, integration, and analysis of network data and network analysis is becoming a routine procedure for biologists to infer biological information. In this review, several main aspects of network studies are discussed, including network construction, analysis, application, and resources.

A computational model for early events in B cell antigen receptor signaling: analysis of the roles of Lyn and Fyn

BCR signaling regulates the activities and fates of B cells. BCR signaling encompasses two feedback loops emanating from Lyn and Fyn, which are Src family protein tyrosine kinases (SFKs). Positive feedback arises from SFK-mediated trans phosphorylation of BCR and receptor-bound Lyn and Fyn, which increases the kinase activities of Lyn and Fyn. Negative feedback arises from SFK-mediated cis phosphorylation of the transmembrane adapter protein PAG1, which recruits the cytosolic protein tyrosine kinase Csk to the plasma membrane, where it acts to decrease the kinase activities of Lyn and Fyn. To study the effects of the positive and negative feedback loops on the dynamical stability of BCR signaling and the relative contributions of Lyn and Fyn to BCR signaling, we consider in this study a rule-based model for early events in BCR signaling that encompasses membrane-proximal interactions of six proteins, as follows: BCR, Lyn, Fyn, Csk, PAG1, and Syk, a cytosolic protein tyrosine kinase that is activated as a result of SFK-mediated phosphorylation of BCR. The model is consistent with known effects of Lyn and Fyn deletions. We find that BCR signaling can generate a single pulse or oscillations of Syk activation depending on the strength of Ag signal and the relative levels of Lyn and Fyn. We also show that bistability can arise in Lyn- or Csk-deficient cells.

RuleBender: integrated modeling, simulation and visualization for rule-based intracellular biochemistry

BACKGROUND

Rule-based modeling (RBM) is a powerful and increasingly popular approach to modeling cell signaling networks. However, novel visual tools are needed in order to make RBM accessible to a broad range of users, to make specification of models less error prone, and to improve workflows.

RESULTS

We introduce RuleBender, a novel visualization system for the integrated visualization, modeling and simulation of rule-based intracellular biochemistry. We present the user requirements, visual paradigms, algorithms and design decisions behind RuleBender, with emphasis on visual global/local model exploration and integrated execution of simulations. The support of RBM creation, debugging, and interactive visualization expedites the RBM learning process and reduces model construction time; while built-in model simulation and results with multiple linked views streamline the execution and analysis of newly created models and generated networks.

CONCLUSION

RuleBender has been adopted as both an educational and a research tool and is available as a free open source tool at http://www.rulebender.org. A development cycle that includes close interaction with expert users allows RuleBender to better serve the needs of the systems biology community.

Rule-based modeling of signal transduction: a primer

Biological cells accomplish their physiological functions using interconnected networks of genes, proteins, and other biomolecules. Most interactions in biological signaling networks, such as bimolecular association or covalent modification, can be modeled in a physically realistic manner using elementary reaction kinetics. However, the size and combinatorial complexity of such reaction networks have hindered such a mechanistic approach, leading many to conclude that it is premature and to adopt alternative statistical or phenomenological approaches. The recent development of rule-based modeling languages, such as BioNetGen (BNG) and Kappa, enables the precise and succinct encoding of large reaction networks. Coupled with complementary advances in simulation methods, these languages circumvent the combinatorial barrier and allow mechanistic modeling on a much larger scale than previously possible. These languages are also intuitive to the biologist and accessible to the novice modeler. In this chapter, we provide a self-contained tutorial on modeling signal transduction networks using the BNG Language and related software tools. We review the basic syntax of the language and show how biochemical knowledge can be articulated using reaction rules, which can be used to capture a broad range of biochemical and biophysical phenomena in a concise and modular way. A model of ligand-activated receptor dimerization is examined, with a detailed treatment of each step of the modeling process. Sections discussing modeling theory, implicit and explicit model assumptions, and model parameterization are included, with special focus on retaining biophysical realism and avoiding common pitfalls. We also discuss the more advanced case of compartmental modeling using the compartmental extension to BioNetGen. In addition, we provide a comprehensive set of example reaction rules that cover the various aspects of signal transduction, from signaling at the membrane to gene regulation. The reader can modify these reaction rules to model their own systems of interest.

Modeling formalisms in Systems Biology

Systems Biology has taken advantage of computational tools and high-throughput experimental data to model several biological processes. These include signaling, gene regulatory, and metabolic networks. However, most of these models are specific to each kind of network. Their interconnection demands a whole-cell modeling framework for a complete understanding of cellular systems. We describe the features required by an integrated framework for modeling, analyzing and simulating biological processes, and review several modeling formalisms that have been used in Systems Biology including Boolean networks, Bayesian networks, Petri nets, process algebras, constraint-based models, differential equations, rule-based models, interacting state machines, cellular automata, and agent-based models. We compare the features provided by different formalisms, and discuss recent approaches in the integration of these formalisms, as well as possible directions for the future.

Guidelines for visualizing and annotating rule-based models

Rule-based modeling provides a means to represent cell signaling systems in a way that captures site-specific details of molecular interactions. For rule-based models to be more widely understood and (re)used, conventions for model visualization and annotation are needed. We have developed the concepts of an extended contact map and a model guide for illustrating and annotating rule-based models. An extended contact map represents the scope of a model by providing an illustration of each molecule, molecular component, direct physical interaction, post-translational modification, and enzyme-substrate relationship considered in a model. A map can also illustrate allosteric effects, structural relationships among molecular components, and compartmental locations of molecules. A model guide associates elements of a contact map with annotation and elements of an underlying model, which may be fully or partially specified. A guide can also serve to document the biological knowledge upon which a model is based. We provide examples of a map and guide for a published rule-based model that characterizes early events in IgE receptor (FcεRI) signaling. We also provide examples of how to visualize a variety of processes that are common in cell signaling systems but not considered in the example model, such as ubiquitination. An extended contact map and an associated guide can document knowledge of a cell signaling system in a form that is visual as well as executable. As a tool for model annotation, a map and guide can communicate the content of a model clearly and with precision, even for large models.

The efficiency of reactant site sampling in network-free simulation of rule-based models for biochemical systems

Rule-based models, which are typically formulated to represent cell signaling systems, can now be simulated via various network-free simulation methods. In a network-free method, reaction rates are calculated for rules that characterize molecular interactions, and these rule rates, which each correspond to the cumulative rate of all reactions implied by a rule, are used to perform a stochastic simulation of reaction kinetics. Network-free methods, which can be viewed as generalizations of Gillespie's method, are so named because these methods do not require that a list of individual reactions implied by a set of rules be explicitly generated, which is a requirement of other methods for simulating rule-based models. This requirement is impractical for rule sets that imply large reaction networks (i.e. long lists of individual reactions), as reaction network generation is expensive. Here, we compare the network-free simulation methods implemented in RuleMonkey and NFsim, general-purpose software tools for simulating rule-based models encoded in the BioNetGen language. The method implemented in NFsim uses rejection sampling to correct overestimates of rule rates, which introduces null events (i.e. time steps that do not change the state of the system being simulated). The method implemented in RuleMonkey uses iterative updates to track rule rates exactly, which avoids null events. To ensure a fair comparison of the two methods, we developed implementations of the rejection and rejection-free methods specific to a particular class of kinetic models for multivalent ligand-receptor interactions. These implementations were written with the intention of making them as much alike as possible, minimizing the contribution of irrelevant coding differences to efficiency differences. Simulation results show that performance of the rejection method is equal to or better than that of the rejection-free method over wide parameter ranges. However, when parameter values are such that ligand-induced aggregation of receptors yields a large connected receptor cluster, the rejection-free method is more efficient.

Hierarchical graphs for rule-based modeling of biochemical systems

Background

In rule-based modeling, graphs are used to represent molecules: a colored vertex represents a component of a molecule, a vertex attribute represents the internal state of a component, and an edge represents a bond between components. Components of a molecule share the same color. Furthermore, graph-rewriting rules are used to represent molecular interactions. A rule that specifies addition (removal) of an edge represents a class of association (dissociation) reactions, and a rule that specifies a change of a vertex attribute represents a class of reactions that affect the internal state of a molecular component. A set of rules comprises an executable model that can be used to determine, through various means, the system-level dynamics of molecular interactions in a biochemical system.

Results

For purposes of model annotation, we propose the use of hierarchical graphs to represent structural relationships among components and subcomponents of molecules. We illustrate how hierarchical graphs can be used to naturally document the structural organization of the functional components and subcomponents of two proteins: the protein tyrosine kinase Lck and the T cell receptor (TCR) complex. We also show that computational methods developed for regular graphs can be applied to hierarchical graphs. In particular, we describe a generalization of Nauty, a graph isomorphism and canonical labeling algorithm. The generalized version of the Nauty procedure, which we call HNauty, can be used to assign canonical labels to hierarchical graphs or more generally to graphs with multiple edge types. The difference between the Nauty and HNauty procedures is minor, but for completeness, we provide an explanation of the entire HNauty algorithm.

Conclusions

Hierarchical graphs provide more intuitive formal representations of proteins and other structured molecules with multiple functional components than do the regular graphs of current languages for specifying rule-based models, such as the BioNetGen language (BNGL). Thus, the proposed use of hierarchical graphs should promote clarity and better understanding of rule-based models.

Rule-based modelling and simulation of biochemical systems with molecular finite automata

The authors propose a theoretical formalism, molecular finite automata (MFA), to describe individual proteins as rule-based computing machines. The MFA formalism provides a framework for modelling individual protein behaviours and systems-level dynamics via construction of programmable and executable machines. Models specified within this formalism explicitly represent the context-sensitive dynamics of individual proteins driven by external inputs and represent protein-protein interactions as synchronised machine reconfigurations. Both deterministic and stochastic simulations can be applied to quantitatively compute the dynamics of MFA models. They apply the MFA formalism to model and simulate a simple example of a signal-transduction system that involves an MAP kinase cascade and a scaffold protein.

Efficient modeling, simulation and coarse-graining of biological complexity with NFsim

Managing the overwhelming numbers of molecular states and interactions is a fundamental obstacle to building predictive models of biological systems. Here we introduce the Network-Free Stochastic Simulator (NFsim), a general-purpose modeling platform that overcomes the combinatorial nature of molecular interactions. Unlike standard simulators that represent molecular species as variables in equations, NFsim uses a biologically intuitive representation: objects with binding and modification sites acted on by reaction rules. During simulations, rules operate directly on molecular objects to produce exact stochastic results with performance that scales independently of the reaction network size. Reaction rates can be defined as arbitrary functions of molecular states to provide powerful coarse-graining capabilities, for example to merge Boolean and kinetic representations of biological networks. NFsim enables researchers to simulate many biological systems that were previously inaccessible to general-purpose software, as we illustrate with models of immune system signaling, microbial signaling, cytoskeletal assembly and oscillating gene expression.

Scalable rule-based modelling of allosteric proteins and biochemical networks

Much of the complexity of biochemical networks comes from the information-processing abilities of allosteric proteins, be they receptors, ion-channels, signalling molecules or transcription factors. An allosteric protein can be uniquely regulated by each combination of input molecules that it binds. This “regulatory complexity” causes a combinatorial increase in the number of parameters required to fit experimental data as the number of protein interactions increases. It therefore challenges the creation, updating, and re-use of biochemical models. Here, we propose a rule-based modelling framework that exploits the intrinsic modularity of protein structure to address regulatory complexity. Rather than treating proteins as “black boxes”, we model their hierarchical structure and, as conformational changes, internal dynamics. By modelling the regulation of allosteric proteins through these conformational changes, we often decrease the number of parameters required to fit data, and so reduce over-fitting and improve the predictive power of a model. Our method is thermodynamically grounded, imposes detailed balance, and also includes molecular cross-talk and the background activity of enzymes. We use our Allosteric Network Compiler to examine how allostery can facilitate macromolecular assembly and how competitive ligands can change the observed cooperativity of an allosteric protein. We also develop a parsimonious model of G protein-coupled receptors that explains functional selectivity and can predict the rank order of potency of agonists acting through a receptor. Our methodology should provide a basis for scalable, modular and executable modelling of biochemical networks in systems and synthetic biology.

The complexity of biochemical networks challenges our ability to create quantitative and predictive models of cellular responses to extracellular changes. In these networks, the regulation of allosteric receptors and proteins by multiple drugs or endogenous ligands introduces “regulatory complexity” because a large number of parameters is required to describe such interactions. Protein interactions also give rise to “combinatorial complexity” by generating large numbers of protein complexes and covalent modification states. To address these twin problems, we propose a modelling framework that combines a modular description of protein structure and function with a rule-based description of protein interactions. We define the input-output function of an allosteric protein through its thermodynamic properties and structural components. We show that our “biomolecule-centric” methodology, in contrast to ad hoc approaches that emphasize the regulatory logic of interactions, can reduce the number of parameters required to model experimental observations. We also demonstrate how the application of our framework gives insights into the assembly of macromolecular complexes and increases the predictive power of a standard model of G protein-coupled receptors. These benefits are possible in many systems, given the ubiquity of allostery in biochemical networks. Our research delineates a fundamental relationship between allostery, modularity, and complexity in biochemical networks.

RuleMonkey: software for stochastic simulation of rule-based models

Background

The system-level dynamics of many molecular interactions, particularly protein-protein interactions, can be conveniently represented using reaction rules, which can be specified using model-specification languages, such as the BioNetGen language (BNGL). A set of rules implicitly defines a (bio)chemical reaction network. The reaction network implied by a set of rules is often very large, and as a result, generation of the network implied by rules tends to be computationally expensive. Moreover, the cost of many commonly used methods for simulating network dynamics is a function of network size. Together these factors have limited application of the rule-based modeling approach. Recently, several methods for simulating rule-based models have been developed that avoid the expensive step of network generation. The cost of these "network-free" simulation methods is independent of the number of reactions implied by rules. Software implementing such methods is now needed for the simulation and analysis of rule-based models of biochemical systems.

Results

Here, we present a software tool called RuleMonkey, which implements a network-free method for simulation of rule-based models that is similar to Gillespie's method. The method is suitable for rule-based models that can be encoded in BNGL, including models with rules that have global application conditions, such as rules for intramolecular association reactions. In addition, the method is rejection free, unlike other network-free methods that introduce null events, i.e., steps in the simulation procedure that do not change the state of the reaction system being simulated. We verify that RuleMonkey produces correct simulation results, and we compare its performance against DYNSTOC, another BNGL-compliant tool for network-free simulation of rule-based models. We also compare RuleMonkey against problem-specific codes implementing network-free simulation methods.

Conclusions

RuleMonkey enables the simulation of rule-based models for which the underlying reaction networks are large. It is typically faster than DYNSTOC for benchmark problems that we have examined. RuleMonkey is freely available as a stand-alone application http://public.tgen.org/rulemonkey. It is also available as a simulation engine within GetBonNie, a web-based environment for building, analyzing and sharing rule-based models.

Modelling the response of FOXO transcription factors to multiple post-translational modifications made by ageing-related signalling pathways

FOXO transcription factors are an important, conserved family of regulators of cellular processes including metabolism, cell-cycle progression, apoptosis and stress resistance. They are required for the efficacy of several of the genetic interventions that modulate lifespan. FOXO activity is regulated by multiple post-translational modifications (PTMs) that affect its subcellular localization, half-life, DNA binding and transcriptional activity. Here, we show how a mathematical modelling approach can be used to simulate the effects, singly and in combination, of these PTMs. Our model is implemented using the Systems Biology Markup Language (SBML), generated by an ancillary program and simulated in a stochastic framework. The use of the ancillary program to generate the SBML is necessary because the possibility that many regulatory PTMs may be added, each independently of the others, means that a large number of chemically distinct forms of the FOXO molecule must be taken into account, and the program is used to generate them. Although the model does not yet include detailed representations of events upstream and downstream of FOXO, we show how it can qualitatively, and in some cases quantitatively, reproduce the known effects of certain treatments that induce various single and multiple PTMs, and allows for a complex spatiotemporal interplay of effects due to the activation of multiple PTM-inducing treatments. Thus, it provides an important framework to integrate current knowledge about the behaviour of FOXO. The approach should be generally applicable to other proteins experiencing multiple regulations.

Modeling multivalent ligand-receptor interactions with steric constraints on configurations of cell-surface receptor aggregates

We use flow cytometry to characterize equilibrium binding of a fluorophore-labeled trivalent model antigen to bivalent IgE-FcepsilonRI complexes on RBL cells. We find that flow cytometric measurements are consistent with an equilibrium model for ligand-receptor binding in which binding sites are assumed to be equivalent and ligand-induced receptor aggregates are assumed to be acyclic. However, this model predicts extensive receptor aggregation at antigen concentrations that yield strong cellular secretory responses, which is inconsistent with the expectation that large receptor aggregates should inhibit such responses. To investigate possible explanations for this discrepancy, we evaluate four rule-based models for interaction of a trivalent ligand with a bivalent cell-surface receptor that relax simplifying assumptions of the equilibrium model. These models are simulated using a rule-based kinetic Monte Carlo approach to investigate the kinetics of ligand-induced receptor aggregation and to study how the kinetics and equilibria of ligand-receptor interaction are affected by steric constraints on receptor aggregate configurations and by the formation of cyclic receptor aggregates. The results suggest that formation of linear chains of cyclic receptor dimers may be important for generating secretory signals. Steric effects that limit receptor aggregation and transient formation of small receptor aggregates may also be important.

The complexity of cell signaling and the need for a new mechanics

Cell signaling systems respond to multiple inputs, such as ligands of cell-surface receptors; and produce multiple outputs, such as changes in gene expression and cellular activities, including motility, proliferation, and death. This "macroscopic" input-output behavior is generated by a web of molecular interactions that can be viewed as taking place at a lower, "microscopic" level. These interactions prominently involve posttranslational modification of proteins and the nucleation of protein complexes. Behaviors at both the micro- and macroscopic levels are complex and must be probed systematically and characterized quantitatively as a prelude to the development of a predictive understanding of a cell signaling system. We must also have a theoretical framework or a mechanics within which we can determine how macroscopic behaviors emerge from known microscopic behaviors or change with manipulations of microscopic behaviors. To connect behaviors at both levels, we suggest that a new mechanics is now required. Newly available data support the idea that this mechanics should enable one to track the site-specific details of molecular interactions in a model, such as the phosphorylation status of individual amino acid residues within a protein.

GetBonNie for building, analyzing and sharing rule-based models

SUMMARY

GetBonNie is a web-based application for building, analyzing and sharing rule-based models encoded in the BioNetGen language (BNGL). Tools accessible within the GetBonNie environment include (i) an applet for drawing graphs that correspond to BNGL code; (ii) a network-generation engine for translating a set of rules into a chemical reaction network; (iii) simulation engines that implement generate-first, on-the-fly and network-free methods for simulating rule-based models; and (iv) a database for sharing models, parameter values, annotations, simulation tasks and results.

AVAILABILITY

GetBonNie is free at (http://getbonnie.org).