RuleMonkey: software for stochastic simulation of rule-based models

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.

Programmatic modeling for biological systems

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

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.

Analysis of Fn14-NF-κB signaling response dynamics using a mechanistic model

Fn14 is a transmembrane receptor protein belonging to the tumor necrosis factor receptor (TNFR) superfamily. Many experimental reports have shown that crosslinking of the receptor by its extracellular ligand TWEAK induces prolonged activation of transcription factor NF-κB. This behavior is distinct from TNF-α receptor, which is a more well-characterized member of the TNFR family. TNF-α receptor, despite sharing many similar molecular interactions with Fn14, only transiently activates NF-κB in response to TNF-α stimulation. Here, we investigate molecular mechanisms that enable Fn14 to display such distinctive behavior. In particular, we focus on two specific features of the Fn14 pathway that potentially give rise to a positive feedback regulation and differentiate it from the TNF-α receptor signaling. By developing a mechanistic model, we analyze how these features may determine the dynamics of an Fn14-NF-κB response. Our analysis reveals that stimulation of Fn14 by TWEAK may generate highly non-linear dynamics, including stable limit cycles and bistable responses. The type of response depends both on the strength and duration of a TWEAK signal. Our predictions and analyses also show that the molecular interactions underlying the positive feedback explain the prolonged activation of NF-κB under certain parameter regimes. In light of the model predictions, we propose possible deregulations of Fn14 leading to its overexpression in solid tumors and tissue injuries.

Modeling cell line-specific recruitment of signaling proteins to the insulin-like growth factor 1 receptor

Receptor tyrosine kinases (RTKs) typically contain multiple autophosphorylation sites in their cytoplasmic domains. Once activated, these autophosphorylation sites can recruit downstream signaling proteins containing Src homology 2 (SH2) and phosphotyrosine-binding (PTB) domains, which recognize phosphotyrosine-containing short linear motifs (SLiMs). These domains and SLiMs have polyspecific or promiscuous binding activities. Thus, multiple signaling proteins may compete for binding to a common SLiM and vice versa. To investigate the effects of competition on RTK signaling, we used a rule-based modeling approach to develop and analyze models for ligand-induced recruitment of SH2/PTB domain-containing proteins to autophosphorylation sites in the insulin-like growth factor 1 (IGF1) receptor (IGF1R). Models were parameterized using published datasets reporting protein copy numbers and site-specific binding affinities. Simulations were facilitated by a novel application of model restructuration, to reduce redundancy in rule-derived equations. We compare predictions obtained via numerical simulation of the model to those obtained through simple prediction methods, such as through an analytical approximation, or ranking by copy number and/or K_D value, and find that the simple methods are unable to recapitulate the predictions of numerical simulations. We created 45 cell line-specific models that demonstrate how early events in IGF1R signaling depend on the protein abundance profile of a cell. Simulations, facilitated by model restructuration, identified pairs of IGF1R binding partners that are recruited in anti-correlated and correlated fashions, despite no inclusion of cooperativity in our models. This work shows that the outcome of competition depends on the physicochemical parameters that characterize pairwise interactions, as well as network properties, including network connectivity and the relative abundances of competitors.

Cells rely on networks of interacting biomolecules to sense and respond to environmental perturbations and signals. However, it is unclear how information is processed to generate appropriate and specific responses to signals, especially given that these networks tend to share many components. For example, receptors that detect distinct ligands and regulate distinct cellular activities commonly interact with overlapping sets of downstream signaling proteins. Here, to investigate the downstream signaling of a well-studied receptor tyrosine kinase (RTK), the insulin-like growth factor 1 (IGF1) receptor (IGF1R), we formulated and analyzed 45 cell line-specific mathematical models, which account for recruitment of 18 different binding partners to six sites of receptor autophosphorylation in IGF1R. The models were parameterized using available protein copy number and site-specific affinity measurements, and restructured to allow for network generation. We find that recruitment is influenced by the protein abundance profile of a cell, with different patterns of recruitment in different cell lines. Furthermore, in a given cell line, we find that pairs of IGF1R binding partners may be recruited in a correlated or anti-correlated fashion. We demonstrate that the simulations of the model have greater predictive power than protein copy number and/or binding affinity data, and that even a simple analytical model cannot reproduce the predicted recruitment ranking obtained via simulations. These findings represent testable predictions and indicate that the outputs of IGF1R signaling depend on cell line-specific properties in addition to the properties that are intrinsic to the biomolecules involved.

MCell-R: A Particle-Resolution Network-Free Spatial Modeling Framework

Spatial heterogeneity can have dramatic effects on the biochemical networks that drive cell regulation and decision-making. For this reason, a number of methods have been developed to model spatial heterogeneity and incorporated into widely used modeling platforms. Unfortunately, the standard approaches for specifying and simulating chemical reaction networks become untenable when dealing with multistate, multicomponent systems that are characterized by combinatorial complexity. To address this issue, we developed MCell-R, a framework that extends the particle-based spatial Monte Carlo simulator, MCell, with the rule-based model specification and simulation capabilities provided by BioNetGen and NFsim. The BioNetGen syntax enables the specification of biomolecules as structured objects whose components can have different internal states that represent such features as covalent modification and conformation and which can bind components of other molecules to form molecular complexes. The network-free simulation algorithm used by NFsim enables efficient simulation of rule-based models even when the size of the network implied by the biochemical rules is too large to enumerate explicitly, which frequently occurs in detailed models of biochemical signaling. The result is a framework that can efficiently simulate systems characterized by combinatorial complexity at the level of spatially resolved individual molecules over biologically relevant time and length scales.

Rule-based modelling provides an extendable framework for comparing candidate mechanisms underpinning clathrin polymerisation

Polymerisation of clathrin is a key process that underlies clathrin-mediated endocytosis. Clathrin-coated vesicles are responsible for cell internalization of external substances required for normal homeostasis and life –sustaining activity. There are several hypotheses describing formation of closed clathrin structures. According to one of the proposed mechanisms cage formation may start from a flat lattice buildup on the cellular membrane, which is later transformed into a curved structure. Creation of the curved surface requires rearrangement of the lattice, induced by additional molecular mechanisms. Different potential mechanisms require a modeling framework that can be easily modified to compare between them. We created an extendable rule-based model that describes polymerisation of clathrin molecules and various scenarios of cage formation. Using Global Sensitivity Analysis (GSA) we obtained parameter sets describing clathrin pentagon closure and the emergence/production and closure of large-size clathrin cages/vesicles. We were able to demonstrate that the model can reproduce budding of the clathrin cage from an initial flat array.

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.

Quantitative computational models of molecular self-assembly in systems biology

Molecular self-assembly is the dominant form of chemical reaction in living systems, yet efforts at systems biology modeling are only beginning to appreciate the need for and challenges to accurate quantitative modeling of self-assembly. Self-assembly reactions are essential to nearly every important process in cell and molecular biology and handling them is thus a necessary step in building comprehensive models of complex cellular systems. They present exceptional challenges, however, to standard methods for simulating complex systems. While the general systems biology world is just beginning to deal with these challenges, there is an extensive literature dealing with them for more specialized self-assembly modeling. This review will examine the challenges of self-assembly modeling, nascent efforts to deal with these challenges in the systems modeling community, and some of the solutions offered in prior work on self-assembly specifically. The review concludes with some consideration of the likely role of self-assembly in the future of complex biological system models more generally.

eSTGt: a programming and simulation environment for population dynamics

Background

We have previously presented a formal language for describing population dynamics based on environment-dependent Stochastic Tree Grammars (eSTG). The language captures in broad terms the effect of the changing environment while abstracting away details on interaction among individuals. An eSTG program consists of a set of stochastic tree grammar transition rules that are context-free. Transition rule probabilities and rates, however, can depend on global parameters such as population size, generation count and elapsed time. In addition, each individual may have an internal state, which can change during transitions.

Results

This paper presents eSTGt (eSTG tool), an eSTG programming and simulation environment. When executing a program, the tool generates the corresponding lineage trees as well as the internal states values, which can then be analyzed either through the tool’s GUI or using MATLAB’s command-line environment.

Conclusions

The presented tool allows researchers to use existing biological knowledge in order to model the dynamics of a developmental process and analyze its behavior throughout the historical events. Simulated lineage trees can be used to validate various hypotheses in silico and to predict the behavior of dynamical systems under various conditions. Written under MATLAB environment, the tool also enables to easily integrate the output data within the user’s downstream analysis.

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.

SpringSaLaD: A Spatial, Particle-Based Biochemical Simulation Platform with Excluded Volume

We introduce Springs, Sites, and Langevin Dynamics (SpringSaLaD), a comprehensive software platform for spatial, stochastic, particle-based modeling of biochemical systems. SpringSaLaD models biomolecules in a coarse-grained manner as a group of linked spherical sites with excluded volume. This mesoscopic approach bridges the gap between highly detailed molecular dynamics simulations and the various methods used to study network kinetics and diffusion at the cellular level. SpringSaLaD is a standalone tool that supports model building, simulation, visualization, and data analysis, all through a user-friendly graphical user interface that should make it more accessible than tools built into more comprehensive molecular dynamics infrastructures. Importantly, for bimolecular reactions we derive an exact expression relating the macroscopic on-rate to the various microscopic parameters with the inclusion of excluded volume; this makes SpringSaLaD more accurate than other tools, which rely on approximate relationships between these parameters.

Analytical reduction of combinatorial complexity arising from multiple protein modification sites

Combinatorial complexity is a major obstacle to ordinary differential equation (ODE) modelling of biochemical networks. For example, a protein with 10 sites that can each be unphosphorylated, phosphorylated or bound to adaptor protein requires 3(10) ODEs. This problem is often dealt with by making ad hoc assumptions which have unclear validity and disallow modelling of site-specific dynamics. Such site-specific dynamics, however, are important in many biological systems. We show here that for a common biological situation where adaptors bind modified sites, binding is slow relative to modification/demodification, and binding to one modified site hinders binding to other sites, for a protein with n modification sites and m adaptor proteins the number of ODEs needed to simulate the site-specific dynamics of biologically relevant, lumped bound adaptor states is independent of the number of modification sites and equal to m + 1, giving a significant reduction in system size. These considerations can be relaxed considerably while retaining reasonably accurate descriptions of the true system dynamics. We apply the theory to model, using only 11 ODEs, the dynamics of ligand-induced phosphorylation of nine tyrosines on epidermal growth factor receptor (EGFR) and primary recruitment of six signalling proteins (Grb2, PI3K, PLCγ1, SHP2, RasA1 and Shc1). The model quantitatively accounts for experimentally determined site-specific phosphorylation and dephosphorylation rates, differential affinities of binding proteins for the phosphorylated sites and binding protein expression levels. Analysis suggests that local concentration of site-specific phosphatases such as SHP2 in membrane subdomains by a factor of approximately 10(7) is critical for effective site-specific regulation. We further show how our framework can be extended with minimal effort to consider binding cooperativity between Grb2 and c-Cbl, which is important for receptor trafficking. Our theory has potentially broad application to reduce combinatorial complexity and allow practical simulation of a variety ODE models relevant to systems biology and pharmacology applications to allow exploration of key aspects of complexity that control signal flux.

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.

Pleomorphic ensembles: formation of large clusters composed of weakly interacting multivalent molecules

Molecular interactions of importance to cell biology are subject to sol-gel transitions: large clusters of weakly interacting multivalent molecules (gel phase) are produced at a critical concentration of monomers. Examples include cell-cell and cell-matrix adhesions, nucleoprotein bodies, and cell signaling platforms. We use the term pleomorphic ensembles (PEs) to describe these clusters, because they have dynamic compositions and sizes and have rapid turnover of their molecular constituents; this plasticity can be highly responsive to cellular signals. The classical polymer physical chemistry theory developed by Flory and Stockmayer provides a brilliant framework for treating multivalent interactions for simple idealized systems. But the complexity and variability of PEs challenges existing modeling approaches. Here we describe and validate a computational algorithm that extends the Flory-Stockmayer formalism to overcome the limitations of analytic theories. We divide the problem by deterministically calculating the fraction of bound sites for each type of binding site, followed by the stochastic assignment of the bonds to a finite number of molecules. The method allows for high valency within many different kinds of interacting molecules and site types, permits simulation of steady-state distributions, as well as assembly kinetics, and can treat cooperative binding within one of the interacting molecules. We then apply our method to the analysis of interactions in the nephrin-Nck-N-Wasp signaling system, demonstrating how multivalent layered scaffolds produce PEs at low monomer concentrations despite weak binding interactions. We show how the experimental data for this system are most consistent with synergistic cooperative interactions between Nck and N-Wasp.

From proteomes to complexomes in the era of systems biology

Protein complexes carry out almost the entire signaling and functional processes in the cell. The protein complex complement of a cell, and its network of complex-complex interactions, is referred to here as the complexome. Computational methods to predict protein complexes from proteomics data, resulting in network representations of complexomes, have recently being developed. In addition, key advances have been made toward understanding the network and structural organization of complexomes. We review these bioinformatics advances, and their discovery-potential, as well as the merits of integrating proteomics data with emerging methods in systems biology to study protein complex signaling. It is envisioned that improved integration of proteomics and systems biology, incorporating the dynamics of protein complexes in space and time, may lead to more predictive models of cell signaling networks for effective modulation.

Phosphorylation site dynamics of early T-cell receptor signaling

In adaptive immune responses, T-cell receptor (TCR) signaling impacts multiple cellular processes and results in T-cell differentiation, proliferation, and cytokine production. Although individual protein-protein interactions and phosphorylation events have been studied extensively, we lack a systems-level understanding of how these components cooperate to control signaling dynamics, especially during the crucial first seconds of stimulation. Here, we used quantitative proteomics to characterize reshaping of the T-cell phosphoproteome in response to TCR/CD28 co-stimulation, and found that diverse dynamic patterns emerge within seconds. We detected phosphorylation dynamics as early as 5 s and observed widespread regulation of key TCR signaling proteins by 30 s. Development of a computational model pointed to the presence of novel regulatory mechanisms controlling phosphorylation of sites with central roles in TCR signaling. The model was used to generate predictions suggesting unexpected roles for the phosphatase PTPN6 (SHP-1) and shortcut recruitment of the actin regulator WAS. Predictions were validated experimentally. This integration of proteomics and modeling illustrates a novel, generalizable framework for solidifying quantitative understanding of a signaling network and for elucidating missing links.

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).

The Simmune Modeler visual interface for creating signaling networks based on bi-molecular interactions

MOTIVATION

Biochemical modeling efforts now frequently take advantage of the possibility to automatically create reaction networks based on the specification of pairwise molecular interactions. Even though a variety of tools exist to visualize the resulting networks, defining the rules for the molecular interactions typically requires writing scripts, which impacts the non-specialist accessibility of those approaches. We introduce the Simmune Modeler that allows users to specify molecular complexes and their interactions as well as the reaction-induced modifications of the molecules through a flexible visual interface. It can take into account the positions of the components of trans-membrane complexes relative to the embedding membranes as well as symmetry aspects affecting the reactions of multimeric molecular structures. Models created with this tool can be simulated using the Simmune Simulator or be exported as SBML code or as files describing the reaction networks as systems of ODEs for import into Matlab.

AVAILABILITY

The Simmune Modeler and the associated simulators as well as extensive additional documentation and tutorials are freely available for Linux, Mac and Windows: http://go.usa.gov/QeH (Note shortened case-sensitive URL!).

Chemical kinetic mechanistic models to investigate cancer biology and impact cancer medicine

Traditional experimental biology has provided a mechanistic understanding of cancer in which the malignancy develops through the acquisition of mutations that disrupt cellular processes. Several drugs developed to target such mutations have now demonstrated clinical value. These advances are unequivocal testaments to the value of traditional cellular and molecular biology. However, several features of cancer may limit the pace of progress that can be made with established experimental approaches alone. The mutated genes (and resultant mutant proteins) function within large biochemical networks. Biochemical networks typically have a large number of component molecules and are characterized by a large number of quantitative properties. Responses to a stimulus or perturbation are typically nonlinear and can display qualitative changes that depend upon the specific values of variable system properties. Features such as these can complicate the interpretation of experimental data and the formulation of logical hypotheses that drive further research. Mathematical models based upon the molecular reactions that define these networks combined with computational studies have the potential to deal with these obstacles and to enable currently available information to be more completely utilized. Many of the pressing problems in cancer biology and cancer medicine may benefit from a mathematical treatment. As work in this area advances, one can envision a future where such models may meaningfully contribute to the clinical management of cancer patients.

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.

The Process-Interaction-Model: a common representation of rule-based and logical models allows studying signal transduction on different levels of detail

Background

Signaling systems typically involve large, structured molecules each consisting of a large number of subunits called molecule domains. In modeling such systems these domains can be considered as the main players. In order to handle the resulting combinatorial complexity, rule-based modeling has been established as the tool of choice. In contrast to the detailed quantitative rule-based modeling, qualitative modeling approaches like logical modeling rely solely on the network structure and are particularly useful for analyzing structural and functional properties of signaling systems.

Results

We introduce the Process-Interaction-Model (PIM) concept. It defines a common representation (or basis) of rule-based models and site-specific logical models, and, furthermore, includes methods to derive models of both types from a given PIM. A PIM is based on directed graphs with nodes representing processes like post-translational modifications or binding processes and edges representing the interactions among processes. The applicability of the concept has been demonstrated by applying it to a model describing EGF insulin crosstalk. A prototypic implementation of the PIM concept has been integrated in the modeling software ProMoT.

Conclusions

The PIM concept provides a common basis for two modeling formalisms tailored to the study of signaling systems: a quantitative (rule-based) and a qualitative (logical) modeling formalism. Every PIM is a compact specification of a rule-based model and facilitates the systematic set-up of a rule-based model, while at the same time facilitating the automatic generation of a site-specific logical model. Consequently, modifications can be made on the underlying basis and then be propagated into the different model specifications – ensuring consistency of all models, regardless of the modeling formalism. This facilitates the analysis of a system on different levels of detail as it guarantees the application of established simulation and analysis methods to consistent descriptions (rule-based and logical) of a particular signaling system.

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 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.

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.

Toward a comprehensive language for biological systems

Rule-based modeling has become a powerful approach for modeling intracellular networks, which are characterized by rich molecular diversity. Truly comprehensive models of cell behavior, however, must address spatial complexity at both the intracellular level and at the level of interacting populations of cells, and will require richer modeling languages and tools. A recent paper in BMC Systems Biology represents a signifcant step toward the development of a unified modeling language and software platform for the development of multi-level, multiscale biological models.

See research article: http://www.biomedcentral.com/1752-0509/5/166

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.

RuleBender: a visual interface for rule-based modeling

SUMMARY

Rule-based modeling (RBM) is a powerful and increasingly popular approach to modeling intracellular biochemistry. Current interfaces for RBM are predominantly text-based and command-line driven. Better visual tools are needed to make RBM accessible to a broad range of users, to make specification of models less error prone and to improve workflows. We present RULEBENDER, an open-source visual interface that facilitates interactive debugging, simulation and analysis of RBMs.

AVAILABILITY

RULEBENDER is freely available for Mac, Windows and Linux at http://rulebender.org.

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.