Simulation of large-scale rule-based models

Abstract

Motivation: Interactions of molecules, such as signaling proteins, with multiple binding sites and/or multiple sites of post-translational covalent modification can be modeled using reaction rules. Rules comprehensively, but implicitly, define the individual chemical species and reactions that molecular interactions can potentially generate. Although rules can be automatically processed to define a biochemical reaction network, the network implied by a set of rules is often too large to generate completely or to simulate using conventional procedures. To address this problem, we present DYNSTOC, a general-purpose tool for simulating rule-based models.

Results: DYNSTOC implements a null-event algorithm for simulating chemical reactions in a homogenous reaction compartment. The simulation method does not require that a reaction network be specified explicitly in advance, but rather takes advantage of the availability of the reaction rules in a rule-based specification of a network to determine if a randomly selected set of molecular components participates in a reaction during a time step. DYNSTOC reads reaction rules written in the BioNetGen language which is useful for modeling protein-protein interactions involved in signal transduction. The method of DYNSTOC is closely related to that of StochSim. DYNSTOC differs from StochSim by allowing for model specification in terms of BNGL, which extends the range of protein complexes that can be considered in a model. DYNSTOC enables the simulation of rule-based models that cannot be simulated by conventional methods. We demonstrate the ability of DYNSTOC to simulate models accounting for multisite phosphorylation and multivalent binding processes that are characterized by large numbers of reactions.

Availability: DYNSTOC is free for non-commercial use. The C source code, supporting documentation and example input files are available at http://public.tgen.org/dynstoc/.

Supplementary information: Supplementary data are available at Bioinformatics online.

Fig. 1.

The steady-state fraction of receptors in ligand-induced receptor aggregates containing three or more receptors, S₃, correlates with the secretory response of RBL cells at different doses of the trivalent ligand of Posner et al. (2002). Points represent measurements of secretion reported in Figure 4 of Posner et al. (2002). The line is obtained by using DYNSTOC to simulate the TLBR model [Equations (10a-c)] with the parameter values indicated below. There are six parameters in this model: ligand dose; N_R, the number of receptors per cell; V, the volume of extracellular fluid surrounding a single cell; k₊₁, the single-site rate constant for ligand–receptor binding when ligand is freely diffusing; k₊₂, the single-site rate constant for ligand–receptor binding when ligand is tethered to the cell surface; and k_off, the single-site rate constant for ligand–receptor dissociation. We set V=10⁻⁹ L, which corresponds to a cell density of 10⁶ cells/ml; we set k_off=0.01 s⁻¹, which is consistent with assays of binding of monovalent hapten to anti-DNP IgE (Erickson et al., ; Xu et al., 1998); and we set N_R=300 000/cell, which is consistent with assays of the number of FcɛRI on RBL cells (Erickson et al., 1987). We systematically varied k₊₁ and k₊₂ in a grid search to find values for which S₃ is maximal at the same ligand dose that yields maximal secretory response. The line shown in this figure is calculated using k₊₁=3.6 × 10⁵ M⁻¹ s⁻¹ and k₊₂N_R=9 × 10⁻⁴ s⁻¹. To speed calculations, we considered only 1% of the volume surrounding a cell in simulations.

Fig. 2.

Validation of DYNSTOC simulation results. (A) Time courses of tyrosine phosphorylation calculated by using numerical integration to solve Equations (8a-c) and (9) for n=6 (dotted lines) and by using DYNSTOC to simulate the corresponding set of rules given by Equations (6a-c) (solid lines). Note the transformation of time t (s). The initial data are [A]=[A_T]=0.12 μM and [U_ij]=1.2 μM, [P_ij]=0, and [AP_ij]=0 for i=1,…, 6 and j=1, 2. Additional parameters have the following values: V=1.4 × 10⁻¹² L (cell volume); k₊₁=k₊₂=7.6 × 10⁶ M⁻¹s⁻¹, k₊₃=k₊₄=1.7 × 10⁶ M⁻¹s⁻¹, k₊₅=k₊₆=6.7 × 10⁶ M⁻¹s⁻¹; k_−i=0.3 s⁻¹ for i=1,…, 6; φ₁=φ₂=0.01 s⁻¹, φ₃=φ₄=0.1 s⁻¹, φ₅=φ₆=0.8 s⁻¹; and δ_I=0.01 s⁻¹ for i=1,…, 6. (B) Time courses of ligand-induced receptor aggregation after the addition of ligand (6 nM) according to Equations (10a-c). Time courses are calculated by using a problem-specific implementation of the YMFH method (dotted lines) and DYNSTOC (solid lines). Parameter values are the same as for Figure 1. Time courses are shown for aggregates containing 2, 3 and 4 receptors. To speed calculations, we considered only 1% of the relevant volume in DYNSTOC simulations; we considered a larger volume in simulations using the YMFH method to reduce noise. The model/simulation-specification files processed by DYNSTOC to produce the results shown here (testcase1.bngl and testcase2.bngl) are available at the DYNSTOC web site.

Fig. 3.

Efficiency of DYNSTOC. We compare DYNSTOC (solid lines) and a problem-specific implementation of the YMFH method (dotted lines); these methods are used to simulate Equations (10a-c). (A) Scaling of computational cost with system size, where size is measured by N_R, the number of cell-surface receptors. Additional parameter values: V=10⁻¹² L; k₊₁=1.8 × 10⁵ M⁻¹s⁻¹; k₊₂N_R^ref=0.03 s⁻¹, where N_R^ref=300; k_off=0.01 s⁻¹; and N_L=4200. (B) Scaling of computational cost with dimensionless parameter β=N_Rk₊₂/k_off, which controls the (equilibrium) extent of ligand-induced receptor aggregation. The value of β was adjusted by varying k₊₂ while holding N_R=300, and k_off=0.01 s⁻¹ fixed. Additional parameter values: V=10⁻¹² L, N_L=4200 and k₊₁=1.8 × 10⁷ M⁻¹s⁻¹. In each panel, the y-axis indicates the total CPU time per reaction event required to simulate the kinetics of the TLBR model from time t=0s to 1000 s with all ligand initially free. Parameter values used in simulations are the same as those used by Yang et al. (2008). See the file tlbr.bngl at the DYNSTOC web site.