Home | About | Journals | Submit | Contact Us | Français |

**|**BMC Bioinformatics**|**v.11; 2010**|**PMC2921409

Formats

Article sections

- Abstract
- Background
- Results
- Conclusions
- Availability and requirements
- Authors' contributions
- Supplementary Material
- References

Authors

Related links

BMC Bioinformatics. 2010; 11: 404.

Published online 2010 July 30. doi: 10.1186/1471-2105-11-404

PMCID: PMC2921409

Reviewed by Joshua Colvin,^{1} Michael I Monine,^{2} Ryan N Gutenkunst,^{2} William S Hlavacek,^{2,}^{3} Daniel D Von Hoff,^{1} and Richard G Posner^{}^{1,}^{4}

Joshua Colvin: gro.negt@nivlocj; Michael I Monine: moc.liamg@eninomim; Ryan N Gutenkunst: vog.lnal@gnayr; William S Hlavacek: vog.lnal@hsiw; Daniel D Von Hoff: gro.negt@hvd; Richard G Posner: gro.negt@rensopr

Received 2010 June 21; Accepted 2010 July 30.

Copyright ©2010 Colvin et al; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article has been cited by other articles in PMC.

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.

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.

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.

A great deal of knowledge about signal transduction, which is mediated largely by the interactions of proteins, has been built up over the years [1,2], in part because molecular changes that affect signal-transduction systems play a role in many diseases, such as cancer [3]. Signal-transduction systems are exceedingly complex, and as a result, our ability to manipulate the behaviors of these systems (e.g., through therapies based on molecularly targeted drugs [4]) is limited. To extend our understanding beyond that reachable through intuition alone, researchers are attempting to develop predictive mathematical models for signal-transduction systems [5-7]. These systems are difficult to model for a variety of reasons [8], and new modeling approaches, such as rule-based modeling [9,10], are needed. Rule-based modeling is notable because it provides a solution to the problem of combinatorial complexity [11].

Central to rule-based modeling is the concept of reaction rules, which will be elaborated in detail below. One can think of a rule as a generalized reaction or a coarse-grained description of a molecular interaction and its consequences. The power of a rule is that it can concisely capture the site-specific details and consequences of an interaction of a structured molecule, such as a protein [9,12]. Functionally, a rule defines 1) the necessary and sufficient properties of reactants in a class of reactions arising from a molecular interaction, 2) the products of a reaction in this class given any particular set of reactants, and 3) a rate law that governs all reactions defined by the rule.

A set of rules for the protein interactions in a signal-transduction system typically implies a much larger set of reactions. The reason is that processes such as multivalent binding and multisite phosphorylation can yield a combinatorial number of chemically distinct protein complexes and phosphoforms [11,13]. The sheer size of such networks, as well as the cost of simulating large-scale chemical reaction networks using conventional methods, has posed a formidable barrier to the development and analysis of models for signal-transduction systems that account for site-specific details of protein interactions (in terms of rules).

To address the problem of simulating models defined by rules that imply large-scale (bio)chemical reaction networks, Colvin et al. [14] generalized the stochastic simulation method of STOCHSIM[15,16] and implemented this method in software called DYNSTOC, which is freely available [17]. The STOCHSIM/DYNSTOC simulation method involves the use of rules to determine if randomly selected components of molecules undergo a reaction. The method has a computational cost that is independent of the size of the reaction network implied by a set of rules. In this respect, the method is similar to methods reported by Danos et al. [18], Yang et al. [19], and Yang and Hlavacek [20]. These methods have been called "network-free" methods [21], in part because the efficiency of these methods is independent of the size of the reaction network implied by rules.

A key feature of DYNSTOC is its ability to simulate models encoded in the BioNetGen language (BNGL), a general-purpose model-specification language that has been described in detail and that can be interpreted by a number of software tools [21]. Thus, DYNSTOC can be used to simulate a wide array of rule sets that imply large-scale chemical reaction networks. Unfortunately, the network-free simulation method implemented in DYNSTOC is a null-event stochastic simulation method, and as a result, simulation of models with fast reactions is inefficient, because the size of the fixed time step used in the simulation procedure is determined by the rate of the fastest reaction.

Here, we present a software tool, RuleMonkey, which implements a network-free stochastic simulation method for determining the system-level dynamics of molecular interactions represented by rules. In this method, unlike in the null-event method of STOCHSIM/DYNSTOC, the size of the time step is variable, and each time step results in a change of the state of the system being simulated. Like DYNSTOC, RuleMonkey is capable of simulating models defined using BNGL.

The algorithm presented below is intended for simulation of rule-based models that can be specified using BNGL [21], particularly models for the system-level dynamics of protein-protein interactions in signal-transduction systems. In such a model (for examples, see [22-24]), molecules and molecular complexes are represented as graphs, and molecular interactions are represented as graph-rewriting rules. It will be useful to briefly review the graphical formalism underlying BNGL [21,25-27].

Molecules (e.g., proteins) are taken to be the building blocks of chemical species (e.g., protein complexes) and to be composed of reactive components, such as amino acid residues subject to post-translational modifications, linear motifs (e.g., proline-rich sequences of the form PxxP), and protein interaction domains [e.g., Src homology 3 (SH3) domains]. Thus, in a rule-based model, molecules and their components comprise a set of chemical species *S *= {*S*_{1}, *S*_{2}, ...}, and each species *S _{i }*is represented by a graph

Rules are used to implicitly represent reactions arising from molecular interactions. A rule *g*_{LHS }→ *g*_{RHS }is composed of 1) sets of left-hand side (LHS) and right-hand side (RHS) pattern graphs, *g*_{LHS }and *g*_{RHS}, which can be viewed as subgraphs of the graphs used to represent chemical species or species graphs, which are denoted *G *= {*G*_{1}, *G*_{2}, ...}; 2) a mapping of the vertices of LHS pattern graphs to the vertices of RHS pattern graphs, which defines a graph transformation; 3) a rate law; and 4) application conditions.

Application conditions are optional but almost all rules encoded in BNGL are endowed with an application condition that restricts molecularity [21,27]. A particular chemical species is a reactant in a rule-defined reaction if its corresponding species graph is matched by a LHS pattern graph (i.e., the species graph contains a subgraph that is isomorphic to the LHS pattern graph) and, in addition, it and any potential co-reactant satisfy the application conditions associated with the rule (if any). The set of molecular components (vertices) directly affected by a transformation will be called a reaction center. Because many species graphs can potentially be matched by a single LHS pattern graph, a rule generally defines not a single reaction, but an entire class of reactions, all of which involve a common transformation. This transformation affects the same reaction center but can take place in multiple molecular contexts, depending on the LHS pattern graph(s) of the rule and the application conditions (if any) of the rule.

Rules are typically associated with application conditions that impose constraints on the molecularity of reactions and the number of reaction products [19,27]. In BNGL, separation of two LHS pattern graphs by a "+" symbol indicates a molecularity of 2, and a rule with a single LHS pattern graph indicates a molecularity of 1. Similarly, separation of two RHS pattern graphs by a "+" symbol indicates two reaction products, and a single RHS pattern graph indicates a single reaction product. Application conditions can be introduced for a variety of other reasons. For example, Monine et al. [28] recently modeled steric effects on multivalent ligand-receptor binding using rules with application conditions that place constraints on the allowable geometric configurations of ligand-induced receptor aggregates. Here, we will only consider application conditions related to molecularity and number of reaction products.

Above, we have implicitly assumed that rules define unidirectional reactions. In BNGL, a short hand notation (< - >) can be used to indicate that the pattern graphs of a rule can be reversed to define reverse reactions [21]. This notation is convenient, as it allows one to avoid redundancy in a model specification. However, the distinction between LHS and RHS graphs is obfuscated. In the discussion that follows, for clarity, we will continue to assume that rules define unidirectional reactions only, which is equivalent to assuming that only the - > notation within BNGL is available. From this point of view, two rules are required to define reversible transformations. For example, A - > B and B - > A must be used instead of simply A <-> B where A and B are pattern graphs. In short, we will refer to those graphs in rules that are used to identify reactants as LHS pattern graphs and those graphs that are used to define transformations as RHS pattern graphs.

The advantage of the rule-based modeling formalism is that a modeler can use a rule to concisely and comprehensivley account for the site-specific details and consequences of a molecular interaction, no matter how many distinct reactions might arise from the interaction. However, this expressiveness comes at a cost. Each individual reaction defined by a rule is assigned the same rate law up to a statistical factor, which can vary from reaction to reaction within a reaction class [21,27]. (Statistical factors are discussed further below.) In reality, each reaction in a reaction class may be characterized by a unique rate law. Thus, the rate law associated with a rule provides only a coarse-grained description of the kinetics of the reactions within the rule-defined reaction class. However, because a rule can specify the molecular context that is permissive for reaction, the coarseness of a rule can be adjusted by tuning the contextual elements of the rule. At the finest level, the contextual elements are highly specific and a rule defines a unique reaction. At the coarsest level, the rule is free of contextual constraints, meaning that the rule indicates that a reaction center can undergo a reaction regardless of the molecular context in which that reaction center is found.

Having presented an overview of the graphical formalism underlying BNGL, we can now introduce the assumptions upon which the algorithm implemented in RuleMonkey is based.

We consider a well-mixed isothermal reaction compartment of volume *V *containing a set of *M *molecules *P *= {*P*_{1}, ..., *P*_{M}}, which we take to be proteins. Each molecule *P*_{i }is associated with a set of components *c*_{i }= {*c*_{i1}, *c*_{i2}, ...}. These components, each of which are optionally associated with an internal state, undergo reactions according to a set of *N *reaction rules *R *= {*R*_{1}, ..., *R*_{N}}. Each rule *R*_{i }defines a class of reactions, and each member of this class is characterized by a single rate law up to a statistical factor, as mentioned earlier. Given this rate law, a cumulative rate *r*_{i }can be determined for the class of reactions defined by *R*_{i}, at least in principle. This cumulative rate corresponds to the sum of the rates for the individual reactions in the class of reactions defined by *R*_{i}. By convention, the rate of each reaction will be taken to be given in terms of the number of molecules undergoing the reaction per unit time in *V*. In other words, each *r*_{i }has the same units as 1*/t*.

As a simplification, we will consider only mass-action rate laws. Thus, for example, we will assume that a bimolecular reaction *A *+ *B *→ product(s) is characterized by a rate law of the form *fk _{V }n_{A}n_{B}*, where

As another simplification, we will assume that internal-state dynamics are described by two-state trajectories. Thus, a component *c _{ij }*of molecule

For purposes of discussion, we will focus on rules specifying reactions that result in the association or dissociation of two components or that change the internal state of a component (i.e., reactions that can be modeled as graph transformations involving the addition or removal of an edge in a graph or the change of a vertex attribute in a graph). Note that a component is taken to have at most one binding partner at a time, as usual [21]. Extension of the discussion that follows to account for other types of rules, such as rules for synthesis, degradation and trafficking between compartments, is straightforward. RuleMonkey is capable of handling these types of reactions, which can be specified in accordance with the conventions of BNGL [21]. For a system in which only association, dissociation, and state-change reactions are possible, counts of molecules and components are conserved quantities. Moreover, for such a system, the LHS pattern graph(s) of a rule will serve to identify either one or two reactants for each reaction defined by the rule.

The state of a system governed by a rule set *R *can be taken to be given by the vector **x**(*t*) = (*x*_{1}(*t*), ..., *x _{n}*(

(1)

where the vector field **f**(**x**) is composed of mass-action terms derived from the rate laws associated with the rule set *R*. For a rule set that implies *m *unidirectional reactions, there are *m *mass-action terms.

In principle, Eq. 1 can be derived from *R *[25,29]. However, derivation and numerical integration of Eq. 1 is impractical when a rule set implies a large-scale chemical reaction network (i.e., large *m *and *n*), which is often the case [9,14,18,19]. Reduced forms of Eq. 1, involving transformations of the variables x, can sometimes be found through analysis of *R *[30-34], and on-the-fly stochastic simulation procedures [25,35,36] can sometimes be used to limit enumeration of the reactions implied by *R*, which is expensive. Nevertheless, these methods fail to be practical for many rule sets.

To overcome this problem, kinetic Monte Carlo (KMC) procedures have been developed in which *R *is used directly to generate stochastic trajectories consistent with the chemical master equation that corresponds to Eq. 1 [14,18-20]. These procedures avoid enumerating the reaction network implied by *R*. Thus, the term "network-free" is apt for these simulation procedures. The algorithm presented below is a network-free procedure.

The RuleMonkey algorithm, like other network-free simulation procedures, is particle based, meaning that each and every material component of a system is tracked, regardless of whether the material components are chemically distinguishable. In other words, given a set of rules *R *and a set of instances of chemical species graphs, each potentially present in multiple copies, these graphs and their constituent parts are tracked individually as the graph transformations defined by *R *are applied to them (as part of a simulation procedure) and the constituent parts move through the space of possible chemical species. Let us use to denote the set of instances of chemical species in a system of volume *V *at time *t*. These instances of chemical species can, in principle, be partitioned into disjoint sets, the set of *n *populated chemical species:

(2)

where is the *j*th copy of species *S _{i }*and

Instead of relying on a partition of as indicated in Eq. 2 or knowledge of x(*t*) to advance a simulation, a network-free algorithm relies on knowledge of the complete component state of a system, i.e., the state of each component in the system. The state of an individual component *c *is the information required to uniquely identify 1) the component, including its label and the label of the molecule of which it is a part; 2) the internal state of the component (if any); and 3) the bound state of the component, including its binding partner (if any). Let us use ∑* _{c}*(

There is a one-to-one relationship between ∑* _{c}*(

The output of a network-free simulation procedure is determined by a specification of a set of *ω *observables, *O *= {*O*_{1}, ..., *O _{ω}*}. Each observable

(3)

(4)

and *w*(*g _{i}*,

Given the background presented above, we can now present an outline of the simulation procedure implemented in RuleMonkey. This procedure has more or less already been described in earlier work [18-20], although key implementation details, described below, are novel. The procedure can be used to produce stochastic trajectories consistent with the chemical master equation corresponding to a system of ODEs of the form of Eq. 1. A comparison of methods is included below.

We take as given a set of *N *rules *R *= {*R*', *R*"}, where *R*' denotes the set of rules with a single LHS pattern graph and *R*" denotes the set of rules with a pair of LHS pattern graphs. We also take as given an initial time *t*(= 0), a system volume *V *, and a set of initial species instances , from which the component state of the system, ∑* _{c}*(

We will also assume that there are procedures that allow us to identify the potential reactants in rule-defined reactions and the reaction centers in these reactants. Thus, for a rule *R _{u }*in

We can now outline the key steps of the network-free simulation procedure implemented in RuleMonkey, which is similar to Gillespie's method [36].

Specify *V*, *t*, *R*, including the single-site rate constants in rate laws associated with rules, and a seed set of species instances , which determines ∑* _{c}*(

Calculate*τ*, the waiting time to the next reaction event in the system, using

(5)

where *ρ*_{1 } (0,1) is a uniform deviate and Recall that the rule rate *r _{i}*, which is taken to have units of inverse time in Eq. 5, is the cumulative rate of all possible reactions defined by rule

Determine the class of the next reaction event by finding the smallest integer *I *that satisfies

(6)

where *ρ*_{2 } (0,1) is a second uniform deviate. The next reaction event will be in the class defined by rule *R _{I}*, where

Additional uniform deviates are needed in the simulation of a reaction event. If *R _{I }*

Increment time *t *← +*τ *and update the sets of reactants associated with rules. Likewise, update the match numbers associated with observables. The steps outlined above that follow initialization are iterated until a stopping criterion is satisfied.

The simulation method described above has been implemented in software called RuleMonkey. RuleMonkey is available as a stand-alone application, which is freely available via the RuleMonkey web site [40], and as a simulation engine within the GetBonNie environment [41,42]. The input for a simulation is a BNGL-encoded model and simulation specification. In other words, RuleMonkey reads and interprets standard BioNetGen input files, like DYNSTOC [14] and BioNetGen [25,29,43]. The conventions of BNGL are described in detail elsewhere [21]. The novel details of algorithm implementation, which are related to the calculation of rule rates, are described below. Usage considerations are also discussed.

Below, equations are given for exactly calculating the rates associated with five types of reaction rules, which generate reactions of the types illustrated in Fig. Fig.1.1. The five rule types are sufficient to specify a wide range of models. For efficiency, the equations given below are used only in the initialization step of the simulation algorithm. After initialization, as described in detail for rejection-free simulation of multivalent ligand-receptor interactions [20], the rates given by these equations are updated on the fly as reactions occur, which requires that matches of LHS pattern graphs in rules be tracked during a simulation.

We will use *R _{a }*to denote a rule that defines unimolecular state-change reactions, such as that illustrated in Fig. 1(a). We will assume that

Given a set of species instances at time *t*, the instantaneous rate *r _{a}*(

(7)

where *k _{a }*is the single-site rate constant associated with the mass-action rate law of

In the course of a RuleMonkey simulation, rates of the form of *r _{a }*are updated on the fly as reactions occur and are never recalculated

We will use *R _{b }*to denote a rule that defines unimolecular dissociation reactions, such as those illustrated in Fig. 1(b). We will assume that

Given a set of species instances at time *t*, the instantaneous rate *r _{b}*(

(8)

where *k _{b }*is the single-site rate constant associated with the mass-action rate law of

As noted above for rates of the form of *r _{a}*, rates of the form of

We will use *R _{c }*to denote a rule that defines bimolecular association reactions, such as those illustrated in Fig. 1(c). We will assume that

Given a set of species instances at time *t*, the instantaneous rate *r _{c}*(

(9)

where

(10)

and

(11)

In the above equations, *k _{c }*is the single-site volume-dependent rate constant associated with the mass-action rate law of

On-the-fly updates of rates of the form of *r _{c }*are fairly complicated but follow from Eq. 9. Basically, when the system state changes,

We will use *R _{d }*to denote a rule that defines unimolecular monogamous ring-closure reactions, such as that illustrated in Fig. 1(d). We assume that

Given a set of species instances at time *t*, the instantaneous rate *r _{d}*(

(12)

where

(13)

and

(14)

In the above equations, *k _{d }*is the single-site rate constant associated with the mass-action rate law of

We will use *R _{e }*to denote a rule that defines unimolecular non-monogamous ring-closure reactions, such as that illustrated in Fig. 1(e). We assume

Given a set of species instances at time *t*, the instantaneous rate *r _{e}*(

(15)

where

(16)

and

(17)

In the above equations, *k _{e }*is the single-site rate constant associated with the mass-action rate law of

In the Actions block of a BioNetGen input file [21], a RuleMonkey-specific command, simulate_{-}rm, is used to initiate a simulation. This command takes the same arguments as analogous commands that invoke other simulation engines, such as DYNSTOC [14]. The output of a RuleMonkey simulation is a .info file and a .gdat file. The .info file reports metadata about a simulation run, such as the execution time. The .gdat file reports a time course for each observable quantity defined in an input file (a file with a .bngl extension). Observables are sums or weighted sums of population levels of chemical species (Eq. 3), and they are defined according to the conventions of BNGL [21]. The time points at which results are reported in the .gdat file are specified by a user using the simulate_{-}rm command. These results are generated through evaluation of the system state at each of the time points of interest.

We validated RuleMonkey by comparing simulation results against those obtained using BioNetGen [25,29], DYNSTOC [14], and problem-specific codes [19,28,44]. The following models were considered (Table (Table1;1; Figs. Figs.2,2, ,3,3, ,4):4): 1) the multisite phosphorylation model introduced by Colvin et al. [14], testcase1.bngl (see Additional file 1); 2) the TLBR (trivalent ligand-bivalent receptor) model introduced by Yang et al. [19] and considered by Colvin et al. [14] and here in Fig. Fig.2B,2B, testcase2a.bngl (see Additional file 2) and testcase2b.bngl (see Additional file 3); 3) a model introduced here with reaction events occurring on disparate time scales, stiff.bngl (see Additional file 4); 4) a model introduced here for interaction of a pentavalent ligand with a trivalent cell-surface receptor, pltr.bngl (see Additional file 5); 5) the model of Blinov et al. [45] for early events in epidermal growth factor receptor (EGFR) signaling, egfr net.bngl (see Additional file 6); 6) the model of Goldstein et al. [46] and Faeder et al. [47] for early events in IgE receptor (FcϵRI) signaling, `fceri.bngl` (see Additional file 7); and 7) the model of Nag et al. [44] for crosslinking of phosphorylated LAT molecules by intracellular Grb2-Sos1 complexes, `lat.bngl` (see Additional file 8). The BioNetGen input files that specify these models and RuleMonkey simulations of them are available as Additional files 12345678. They are also available at the RuleMonkey web site [40].

Each model is specified as a BioNetGen input file, which can be processed by BioNetGen, DYNSTOC, and RuleMonkey. The input files for these different software tools are the same, except for the simulation commands in the Actions blocks [21]. As noted above, the simulation command for RuleMonkey is simulate_{-}rm. BioNetGen accepts two simulation commands, simulate_{-}ode and simulate_{-}ssa. In calculations with BioNetGen, we used simulate_{-}ssa, which invokes an efficient extension of Gillespie's method, preceded by a generate_{-}network command, which prompts BioNetGen to implement the procedures needed to explicitly derive the reaction network implied by a set of rules. In other words, we only used BioNetGen to perform simulations via the so-called generate-first approach [9,21,48].

The reason we only consider generate-first simulations using BioNetGen is that simulation cost and network-generation cost can be easily separated when this simulation approach is taken, as it consists of two steps. The first step is network generation. The second step is simulation. In the on-the-fly approach, execution switches back and forth between network generation and simulation. It has been observed that stochastic simulation via the on-the-fly approach tends to be faster than stochastic simulation via the generate-first approach [35]. Thus, if we find that generate-first simulation is faster than network-free simulation, we can assume with some confidence that on-the-fly simulation is also faster than network-free simulation.

Typical validation results are shown in Fig. Fig.2.2. In all cases, RuleMonkey was found to produce results consistent with those obtained using the other methods/codes considered.

As can be seen from the results of Table Table1,1, RuleMonkey is more efficient than DYNSTOC, but less efficient than BioNetGen and problem-specific codes (for cases where BioNetGen and problem-specific codes can be applied). BioNetGen is unable to simulate three of the test models via the generate-first approach because the reaction networks implied by the rule sets of these models cannot be generated exhaustively. Thus, of the available general-purpose simulation codes compliant with the conventions of BNGL, RuleMonkey and DYNSTOC are the most widely applicable (because these tools implement general-purpose network-free simulation methods), and RuleMonkey is more efficient than DYNSTOC (Table (Table1).1). It should be noted that we would expect the generate-first approach to simulation [9,21,48], which may involve either stochastic simulation or numerical integration of ODEs (although we consider only stochastic simulation here), to be more efficient than network-free simulation for cases where the generate-first approach is feasible because of the bookkeeping costs of network-free simulation, which requires data structures and update schemes to track the component state of a system. Unfortunately, simulation methods that rely on derivation of the reaction network implied by a rule set, which is computationally expensive, are often not feasible for rule sets that imply large-scale chemical reaction networks. Feasibility depends on effective network size, which is a function of model parameter values. For some parameter values, network generation can be halted arbitrarily or limited as in the on-the-fly approach to simulation [9,21,48]. For other parameter values, as demonstrated by Yang et al. [19], for example, network generation is a limiting factor, and network-free simulation is essential.

The performance of RuleMonkey scales with problem size in a way that is similar to that of DYNSTOC [14] and the problem-specific codes of Yang et al. [19] and Yang and Hlavacek [20]. Results obtained from simulation of the TLBR model of Yang et al. [19] are shown in Fig. Fig.3.3. The TLBR model, a kinetic version of the model of Goldstein and Perelson [49], characterizes the interactions of trivalent ligands with bivalent cell-surface receptors. In Fig. 3(a), the cost of simulation is shown as a function of the number of receptors *N _{R}*. As

Finally, we compared the performances of RuleMonkey and DYNSTOC for problems with reactions occurring on different time scales (Fig. (Fig.4).4). As can be seen, as the rate of the fastest reaction in a system increases, the efficiency of DYNSTOC degrades relative to that of RuleMonkey.

RuleMonkey is a general-purpose simulator for rule-based models encoded in BNGL. RuleMonkey complements BioNetGen [21,25,29], the first software tool to have the capability to simulate rule-based models encoded in BNGL. The rules that comprise BNGL-encoded models can be expanded by BioNetGen into the set of reactions implied by the rules, and BioNetGen can simulate the kinetics of these reactions using conventional methods, such as numerical integration of the corresponding system of ordinary differential equations (ODEs) for mass-action kinetics. Although BioNetGen has been used to study a number of systems [22-24,45,50], the reaction network implied by a set of rules is often so large as to make the simulation approaches implemented in BioNetGen impractical. Network generation is expensive, as is simulation of a large-scale reaction network using a conventional method. Such methods have a computational cost that depends on the size of the network being simulated. In the generate-first and on-the-fly approaches to simulating rule-based models [9,21,25,35,48], which are implemented in BioNetGen, there is both a network generation cost and a network simulation cost, because network generation is a prerequisite for network simulation via these approaches. In contrast, in network-free simulation approaches [14,15,18-20], such as the method presented here and implemented in RuleMonkey, there is only a network simulation cost. However, it should be noted that network-free simulation procedures have a relatively high cost per reaction event, as is evident in Table Table1.1. Thus, one should choose such a simulation procedure only when the cost of an alternative procedure that relies on network-generation is prohibitive.

The cost of network generation is often prohibitive, which has motivated the development of several network-free simulation methods [14,15,18-20]. Like these methods, the method presented here and implemented in RuleMonkey avoids network generation and, consequently, has a computational cost independent of the number of reactions implied by a set of rules (Figs. (Figs.22 and and3).3). We note that Yang and Hlavacek [20] have described a method that is essentially the same as that implemented in RuleMonkey but without providing details about how the method could be implemented for general use (i.e., beyond problem-specific demonstrations of the method) and without providing general-purpose software.

RuleMonkey and DYNSTOC [14] are similar in that they can both simulate BNGL-encoded models while avoiding network generation. However the method implemented in RuleMonkey is a rejection-free method, whereas the method implemented in DYNSTOC is a null-event method. Both types of methods are capable of producing correct simulation results but null-event methods tend to have a higher cost for simulation of systems with fast reactions [51]. For cases that we have examined, RuleMonkey typically performs better than DYNSTOC (Table (Table1)1) and has a decided advantage for systems with fast reactions (Fig. (Fig.4).4). This advantage is perhaps significant because practical problems tend to involve reactions occurring on disparate time scales.

The capability of RuleMonkey to correctly process BNGL-encoded model and simulation specifications (as illustrated in Fig. Fig.2)2) is significant in that this capability allows RuleMonkey to be used to simulate a wide array of models that account for the site-specific details of protein-protein interactions (e.g., multisite phosphorylation) as well as the connectivity of proteins in protein complexes (e.g., cyclic aggregates of proteins can be distinguished from acyclic aggregates). RuleMonkey also distinguishes between intra- and intermolecular reactions, which is important for correct simulation results in many cases [19,28]. We have demonstrated that RuleMonkey can be used to efficiently simulate rule sets that imply large-scale chemical reaction networks (Table (Table1;1; Figs. Figs.22 and and3).3). The capabilities of RuleMonkey complement those of BioNetGen and other available software tools for simulation of rule-based modeling [14,35,52-56]. In the future, we expect RuleMonkey will prove useful for simulating models of signal-transduction systems that are significantly larger and more comprehensive than the rule-based models of such systems that have so far been considered [22-24,45-47].

• **Project name**: RuleMonkey

• **Project home page**: http://public.tgen.org/rulemonkey

• **Operating system(s)**: Platform independent

• **Programming language**: C

• **Other requirements**: None

• **License**: GNU General Public License (GPL), version 3

• **Any restrictions to use by non-academics**: If the terms of the GPL are unacceptable, it may be possible to provide the software under a different license.

JC implemented the algorithm and designed the tool. MIM and RNG tested and validated the tool and evaluated its performance. DVH, RGP and WSH planned the project. RGP and WSH conceived the tool and wrote the manuscript. All authors have read and approved the final version of the manuscript.

**This plain-text file specifies a simulation of the multisite phosphorylation model introduced by Colvin et al**. [14]. This simulation is considered in Fig. Fig.2A2A and Table Table1.1. The file can be processed by RuleMonkey and other BNGL-compatible tools for simulating rule-based models.

Click here for file^{(2.3K, BNGL)}

**This plain-text file specifies a simulation of a model for trivalent ligand-bivalent cell-surface receptor interaction (the so-called TLBR model) introduced by Yang et al**. [19]. This simulation is considered in Fig. Fig.3.3. The file can be processed by RuleMonkey and other BNGL-compatible tools for simulating rule-based models.

Click here for file^{(7.6K, BNGL)}

**This plain-text file specifies a simulation of a model for trivalent ligand-bivalent cell-surface receptor interaction (the so-called TLBR model) introduced by Yang et al**. [19]. This simulation is considered in Fig. Fig.2B2B and Table Table1.1. The difference between the files testcase2a.bngl and testcase2b.bngl is in parameter values. The file can be processed by RuleMonkey and other BNGL-compatible tools for simulating rule-based models.

Click here for file^{(7.7K, BNGL)}

**This plain-text file specifies a simulation of reaction events occurring on disparate time scales**. This simulation is considered in Fig. Fig.44 and Table Table1.1. The file can be processed by RuleMonkey and other BNGL-compatible tools for simulating rule-based models.

Click here for file^{(469 bytes, BNGL)}

**This plain-text file specifies a simulation of a model for pentavalent ligand-trivalent cell-surface receptor interactions**. The model is based on assumptions similar to those upon which the TLBR model [19] and the model of Goldstein and Perelson [49] are based. This simulation is considered in Table Table1.1. The file can be processed by RuleMonkey and other BNGL-compatible tools for simulating rule-based models.

Click here for file^{(7.9K, BNGL)}

**This plain-text file specifies a simulation of the model of Blinov et al**. [45] for early events in epidermal growth factor receptor (EGFR) signaling. This simulation is considered in Table Table1.1. The file can be processed by RuleMonkey and other BNGL-compatible tools for simulating rule-based models.

Click here for file^{(6.1K, BNGL)}

**This plain-text file specifies a simulation of the model of Goldstein et al**. [46] and Faeder et al. [47] for early events in IgE receptor (FcϵRI) signaling. This simulation is considered in Table Table1.1. The file can be processed by RuleMonkey and other BNGL-compatible tools for simulating rule-based models.

Click here for file^{(3.4K, BNGL)}

**This plain-text file specifies a simulation of the model of Nag et al**. [44] for crosslinking of phosphorylated LAT molecules by intracellular Grb2-Sos1 complexes. This simulation is considered in Table Table1.1. The file can be processed by RuleMonkey and other BNGL-compatible tools for simulating rule-based models.

Click here for file^{(1.7K, BNGL)}

We thank J. Yang for helpful discussions and B. Hu for integrating RuleMonkey into GetBonNie. This work was supported by National Institutes of Health grants AI35997, CA109552, NS042262, GM076570; DOE contract DE-AC52-06NA25396; and the Arizona Biomedical Research Commission.

- Hunter T. Signaling--2000 and beyond. Cell. 2000;100:113–127. doi: 10.1016/S0092-8674(00)81688-8. [PubMed] [Cross Ref]
- Scott JD, Pawson T. Cell signaling in space and time: where proteins come together and when they're apart. Science. 2009;326:1220–1224. doi: 10.1126/science.1175668. [PMC free article] [PubMed] [Cross Ref]
- Hanahan D, Weinberg RA. The hallmarks of cancer. Cell. 2000;100:57–70. doi: 10.1016/S0092-8674(00)81683-9. [PubMed] [Cross Ref]
- Hunter T. Treatment for chronic myelogenous leukemia: the long road to imatinib. J Clin Invest. 2007;117:2036–2043. doi: 10.1172/JCI31691. [PMC free article] [PubMed] [Cross Ref]
- Kholodenko BN. Cell-signalling dynamics in time and space. Nat Rev Mol Cell Biol. 2006;7:165–176. doi: 10.1038/nrm1838. [PMC free article] [PubMed] [Cross Ref]
- Aldridge BB, Burke JM, Lauffenburger DA, Sorger PK. Physicochemical modelling of cell signaling pathways. Nat Cell Biol. 2006;8:1195–1203. doi: 10.1038/ncb1497. [PubMed] [Cross Ref]
- Chakraborty AK, Das J. Pairing computation with experimentation: a powerful coupling for understanding T cell signalling. Nat Rev Immunol. 2010;10:59–71. doi: 10.1038/nri2688. [PubMed] [Cross Ref]
- Breitling R, Hoeller D. Current challenges in quantitative modeling of epidermal growth factor signaling. FEBS Lett. 2005;579:6289–6294. doi: 10.1016/j.febslet.2005.10.034. [PubMed] [Cross Ref]
- Hlavacek WS, Faeder JR, Blinov ML, Posner RG, Hucka M, Fontana W. Rules for modeling signal-transduction systems. Sci STKE. 2006;2006:re6. doi: 10.1126/stke.3442006re6. [PubMed] [Cross Ref]
- Hlavacek WS, Faeder JR. The complexity of cell signaling and the need for a new mechanics. Sci Signal. 2009;2:pe46. doi: 10.1126/scisignal.281pe46. [PubMed] [Cross Ref]
- Hlavacek WS, Faeder JR, Blinov ML, Perelson AS, Goldstein B. The complexity of complexes in signal transduction. Biotechnol Bioeng. 2003;84:783–794. doi: 10.1002/bit.10842. [PubMed] [Cross Ref]
- Danos V, Feret J, Fontana W, Harmer R, Krivine J. Rule-based modelling of cellular signalling. Lect Notes Comput Sci. 2007;4703:17–41. full_text.
- Mayer BJ, Blinov ML, Loew LM. Molecular machines or pleiomorphic ensembles: signaling complexes revisited. J Biol. 2009;8:81. doi: 10.1186/jbiol185. [PMC free article] [PubMed] [Cross Ref]
- Colvin J, Monine MI, Faeder JR, Hlavacek WS, Von Hoff DD, Posner RG. Simulation of large-scale rule-based models. Bioinformatics. 2009;25:910–917. doi: 10.1093/bioinformatics/btp066. [PMC free article] [PubMed] [Cross Ref]
- Morton-Firth CJ, Bray D. Predicting temporal fluctuations in an intracellular signalling pathway. J Theor Biol. 1998;192:117–128. doi: 10.1006/jtbi.1997.0651. [PubMed] [Cross Ref]
- Shimizu TS, Bray D. In: Foundations of Systems Biology. Kitano H, editor. Ch 10. Cambridge, MA: MIT Press; 2001. Computational cell biology--the stochastic approach.
- The DYNSTOC web site. http://public.tgen.org/dynstoc/
- Danos V, Feret J, Fontana W, Krivine J. Scalable simulation of cellular signaling networks. Lect Notes Comput Sci. 2007;4807:139–157. full_text.
- Yang J, Monine MI, Faeder JR, Hlavacek WS. Kinetic Monte Carlo method for rule-based modeling of biochemical networks. Phys Rev E. 2008;78:031910. doi: 10.1103/PhysRevE.78.031910. [PMC free article] [PubMed] [Cross Ref]
- Yang J, Hlavacek WS. Rejection-free kinetic Monte Carlo simulation of multivalent biomolecular interactions. http://arxiv.org/abs/0812.4619
- Faeder JR, Blinov ML, Hlavacek WS. Rule-based modeling of biochemical systems with BioNetGen. Methods Mol Biol. 2009;500:113–167. [PubMed]
- Barua D, Faeder JR, Haugh JM. Structure-based kinetic models of modular signaling protein function: focus on Shp2. Biophys J. 2007;92:2290–2300. doi: 10.1529/biophysj.106.093484. [PubMed] [Cross Ref]
- Barua D, Faeder JR, Haugh JM. Computational models of tandem Src homology 2 domain interactions and application to phosphoinositide 3-kinase. J Biol Chem. 2008;283:7338–7345. doi: 10.1074/jbc.M708359200. [PubMed] [Cross Ref]
- Barua D, Faeder JR, Haugh JM. A bipolar clamp mechanism for activation of Jak-family protein tyrosine kinases. PLoS Comput Biol. 2009;5:e1000364. doi: 10.1371/journal.pcbi.1000364. [PMC free article] [PubMed] [Cross Ref]
- Faeder JR, Blinov ML, Goldstein B, Hlavacek WS. Rule-based modeling of biochemical networks. Complexity. 2005;10:22–41. doi: 10.1002/cplx.20074. [Cross Ref]
- Faeder JR, Blinov ML, Hlavacek WS. In: Proceedings of the 2005 ACM Symposium on Applied Computing: 13-17 March 2005; Santa Fe, NM. Liebrock LM, editor. ACM Press; 2005. Graphical rule-based representation of signal-transduction networks; pp. 133–140.
- Blinov ML, Yang J, Faeder JR, Hlavacek WS. Graph theory for rule-based modeling of biochemical networks. Lect Notes Comput Sci. 2006;4230:89–106. full_text.
- Monine MI, Posner RG, Savage PB, Faeder JR, Hlavacek WS. Modeling multivalent ligand-receptor interactions with steric constraints on configurations of cell-surface receptor aggregates. Biophys J. 2010;98:48–56. doi: 10.1016/j.bpj.2009.09.043. [PubMed] [Cross Ref]
- Blinov ML, Faeder JR, Goldstein B, Hlavacek WS. BioNetGen: software for rule-based modeling of signal transduction based on the interactions of molecular domains. Bioinformatics. 2006;20:3289–3291. doi: 10.1093/bioinformatics/bth378. [PubMed] [Cross Ref]
- Koschorreck M, Conzelmann H, Ebert S, Ederer M, Gilles ED. Reduced modeling of signal transduction--a modular approach. BMC Bioinformatics. 2007;8:336. doi: 10.1186/1471-2105-8-336. [PMC free article] [PubMed] [Cross Ref]
- Conzelmann H, Fey D, Gilles ED. Exact model reduction of combinatorial reaction networks. BMC Syst Biol. 2008;2:78. [PMC free article] [PubMed]
- Conzelmann H, Gilles ED. Dynamic pathway modeling of signal transduction networks: a domain-oriented approach. Methods Mol Biol. 2008;484:559–578. full_text. [PubMed]
- Borisov NM, Chistopolsky AS, Faeder JR, Kholodenko BN. Domain-oriented reduction of rule-based network models. IET Syst Biol. 2008;2:342–351. doi: 10.1049/iet-syb:20070081. [PMC free article] [PubMed] [Cross Ref]
- Feret J, Danos V, Krivine J, Harmer R, Fontana W. Internal coarse-graining of molecular systems. Proc Natl Acad Sci USA. 2009;106:6453–6458. doi: 10.1073/pnas.0809908106. [PubMed] [Cross Ref]
- Lok L, Brent R. Automatic generation of cellular reaction networks with Moleculizer 1.0. Nat Biotechnol. 2005;23:131–136. doi: 10.1038/nbt1054. [PubMed] [Cross Ref]
- Gillespie DT. Stochastic simulation of chemical kinetics. Annu Rev Phys Chem. 2007;58:35–55. doi: 10.1146/annurev.physchem.58.032806.104637. [PubMed] [Cross Ref]
- Voter AF. In: Radiation Effects in Solids. Sickafus KE, Kotomin EA, Uberuaga BP, editor. Ch 1. Dordrecht, The Netherlands: Springer; 2007. Introduction to the kinetic Monte Carlo method.
- Blue JL, Beichl I, Sullivan F. Faster Monte Carlo simulations. Phys Rev E. 1995;51:R867-R868. doi: 10.1103/PhysRevE.51.R867. [PubMed] [Cross Ref]
- Gibson MA, Bruck J. Efficient exact stochastic simulation of chemical systems with many species and many channels. J Phys Chem A. 2000;104:1876–1889. doi: 10.1021/jp993732q. [Cross Ref]
- The RuleMonkey web site. http://public.tgen.org/rulemonkey/
- Hu B, Fricke GM, Faeder JR, Posner RG, Hlavacek WS. GetBonNie for building, analyzing and sharing rule-based models. Bioinformatics. 2009;25:1457–1460. doi: 10.1093/bioinformatics/btp173. [PMC free article] [PubMed] [Cross Ref]
- The GetBonNie web site. http://getbonnie.org
- The BioNetGen web site. http://bionetgen.org
- Nag A, Monine MI, Faeder JR, Goldstein B. Aggregation of membrane proteins by cytosolic cross-linkers: theory and simulation of the LAT-Grb2-SOS1 system. Biophys J. 2009;96:2604–2623. doi: 10.1016/j.bpj.2009.01.019. [PubMed] [Cross Ref]
- Blinov ML, Faeder JR, Goldstein B, Hlavacek WS. A network model of early events in epidermal growth factor receptor signaling that accounts for combinatorial complexity. BioSystems. 2006;83:136–151. doi: 10.1016/j.biosystems.2005.06.014. [PubMed] [Cross Ref]
- Goldstein B, Faeder JR, Hlavacek WS, Blinov ML, Redondo A, Wofsy C. Modeling the early signaling events mediated by FcϵRI. Mol Immunol. 2002;38:1213–1219. doi: 10.1016/S0161-5890(02)00066-4. [PubMed] [Cross Ref]
- Faeder JR, Hlavacek WS, Reischl I, Blinov ML, Metzger H, Redondo A, Wofsy C, Goldstein B. Investigation of early events in FcϵRI-mediated signaling using a detailed mathematical model. J Immunol. 2003;170:3769–3781. [PubMed]
- Blinov ML, Faeder JR, Yang J, Goldstein B, Hlavacek WS. 'On-the-fly' or 'generate-first' modeling? Nat Biotechnol. 2005;23:1344–1345. doi: 10.1038/nbt1105-1344. [PubMed] [Cross Ref]
- Goldstein B, Perelson AS. Equilibrium theory for the clustering of bivalent cell surface receptors by trivalent ligands. Application to histamine release from basophils. Biophys J. 1984;45:1109–1123. doi: 10.1016/S0006-3495(84)84259-9. [PubMed] [Cross Ref]
- Faeder JR, Blinov ML, Goldstein B, Hlavacek WS. Combinatorial complexity and dynamical restriction of network flows in signal transduction. Syst Biol. 2005;2:5–15. doi: 10.1049/sb:20045031. [PubMed] [Cross Ref]
- Chatterjee A, Vlachos DG. An overview of spatial microscopic and accelerated kinetic Monte Carlo methods. J Computer-Aided Mater Des. 2007;14:253–308. doi: 10.1007/s10820-006-9042-9. [Cross Ref]
- Meier-Schellersheim M, Xu X, Angermann B, Kunkel EJ, Jin T, Germain RN. Key role of local regulation in chemosensing revealed by a new molecular interaction-based modeling method. PLoS Comput Biol. 2006;2:e82. doi: 10.1371/journal.pcbi.0020082. [PubMed] [Cross Ref]
- Moraru II, Schaff JC, Slepchenko BM, L BM, Morgan F, Lakshminarayana A, Gao F, Li Y, Loew LM. Virtual Cell modelling and simulation software environment. IET Syst Biol. 2008;2:352–362. doi: 10.1049/iet-syb:20080102. [PMC free article] [PubMed] [Cross Ref]
- Mallavarapu A, Thomson M, Ullian B, Gunawardena J. Programming with models: modularity and abstraction provide powerful capabilities for systems biology. J R Soc Interface. 2009;6:257–270. doi: 10.1098/rsif.2008.0205. [PMC free article] [PubMed] [Cross Ref]
- Lis M, Artyomov MN, Devadas S, Chakraborty AK. Efficient stochastic simulation of reaction-diffusion processes via direct compilation. Bioinformatics. 2009;25:2289–2291. doi: 10.1093/bioinformatics/btp387. [PMC free article] [PubMed] [Cross Ref]
- Andrews SS, Addy NJ, Brent R, Arkin AP. Detailed simulations of cell biology with Smoldyn 2.1. PLoS Comput Biol. 2010;6:e1000705. doi: 10.1371/journal.pcbi.1000705. [PMC free article] [PubMed] [Cross Ref]

Articles from BMC Bioinformatics are provided here courtesy of **BioMed Central**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |