Kinetic Monte Carlo method for rule-based modeling of biochemical networks

Abstract

We present a kinetic Monte Carlo method for simulating chemical transformations specified by reaction rules, which can be viewed as generators of chemical reactions, or equivalently, definitions of reaction classes. A rule identifies the molecular components involved in a transformation, how these components change, conditions that affect whether a transformation occurs, and a rate law. The computational cost of the method, unlike conventional simulation approaches, is independent of the number of possible reactions, which need not be specified in advance or explicitly generated in a simulation. To demonstrate the method, we apply it to study the kinetics of multivalent ligand-receptor interactions. We expect the method will be useful for studying cellular signaling systems and other physical systems involving aggregation phenomena.

FIG. 1

TLBR model. (a) A ligand with three identical binding sites and a mobile cell-surface receptor with two identical binding sites. The ligand mediates cross-linking of receptors as shown. (b) Rules representing capture of a freely diffusing ligand by a receptor (R₁), ligand-mediated receptor cross-linking (R₂), and ligand-receptor dissociation (R₃). Parameters of the rate laws associated with these rules are single-site rate constants: k₊₁, k₊₂, and k_off, respectively. An empty (filled) circle indicates a free (bound) site, a line connecting circles indicates a bond, and an empty box or wedge indicates a site that may be either free or bound. In BNGL [29], the rules are specified as follows: R₁ is L(r, r, r) + R(l) -> L(r!1, r, r).R(l!1), R₂ is L(r!+, r) + R(l) -> L(r!+, r!1).R(l!1), and R₃ is L(r!1).R(l!1) -> L(r) + R(l), where land rare used to represent binding sites of the receptor ( R) and ligand ( L), respectively.

FIG. 2

Percolation transition between sol and sol-gel regions in the space of c and β. The curve marks the percolation transition boundary according to the equilibrium continuum model of Goldstein and Perelson [28]. Using the rule-based KMC method, we simulated the TLBR model to determine the steady-state value of f_g, the fraction of receptors in the gel phase (i.e., in the largest receptor aggregate), as a function of c and β. At points marked by dots, f_g ≥ 0.05, whereas at points marked by circles, f_g < 0.05. To adjust the values of c and β, we varied k₊₁ and k₊₂ and held other parameters constant at the following values: N_R = 3, 000, N_L = 42, 000, and k_off = 0.01 s⁻¹.

FIG. 3

Efficiency of simulation of the TLBR model. (a) Dependence of CPU time per reaction event for rule-based KMC simulation (solid line) vs. on-the-fly simulation [7, 13, 14] (dashed line). (b) Effective network size as a function of β. The solid and dashed lines indicate the numbers of species populated and reactions fired, respectively, in on-the-fly simulation. Calculations were performed using BioNetGen [6, 29]. (c) Dependence of CPU time per reaction event on N_R for β = 50 (solid line), β = 0.1 (dashed line), and β = 50 without connectivity checks (dotted line). For β = 50, the fraction of KMC steps that result in null events is approximately 0.6 for any value of N_R. The fraction is essentially 0.0 for β = 0.1. Note that the system is above (below) the percolation transition at β = 50 (β = 0.1). (d) Importance of null events. The solid and dashed lines are calculated using auxiliary non-local component state information to minimize the cost of null events for β = 50 and β = 0.1, respectively. The line broken in a dash-dot pattern and the dotted line are calculated using a problem-specific rejection-free procedure for β = 50 and β = 0.1, respectively. Additional simulation parameters: (a) and (b) N_R = 300, N_L = 4, 200, and c = 0.84; (c) and (d) N_L = 14N_R and c = 0.84. The value of k_off was held fixed at 0.01 s⁻¹ in all simulations. All reported results are based on simulation for 3,000 s after equilibration.

FIG. 4

Kinetics of the TLBR model. (a) Fraction of receptors in aggregates with 1, 3, 5, 7, or 9 receptors or in the largest aggregate as a function of time in the sol-gel coexistence phase (N_L = 50, 000 and c = 2.7). (b) Mean aggregate size as a function of time for same conditions as (a) (solid line) and at a lower ligand concentration (dashed line) that gives the same mean size at equilibrium (N_L = 2, 000 and c = 0.11). Additional simulation parameters: N_R = 3, 000, β = 16.8, and k_off = 0.01 s⁻¹. Results are averaged over 40 simulation runs. Mean aggregate size is determined by $\langle S\rangle ={\sum }_{i=2}^{N_R}i{n}_i/{\sum }_{i=2}^{N_R}{n}_i$, where n_i is the number of aggregates containing i receptors. Parameter values were chosen arbitrarily for the purpose of demonstrating the rule-based KMC method, but they are expected to be somewhat reasonable for the case of a population of ligands, each with three 2, 4-dinitrophenol (DNP) hapten groups, interacting with a population of monoclonal cell-surface anti-DNP IgE antibodies, each with two antigen-combining sites [24, 26].

FIG. 5

Generation of the reaction network implied by the rules of Fig. 1. Starting from two speed species (free ligand and free receptor), successive rounds of rule application generate new chemical species and reactions. In the process of network generation, species are represented by graphs and rule application is comprised of graph rewriting operations [15]. The two seed species and the four species generated in the first two rounds of rule application are illustrated using the conventions of Fig. 1. White bars indicate the number of species in the partially generated network at each step in the process of network generation. Black bars indicate the number of reactions. Indicated at top is the total CPU time required to perform each of the first four rounds of rule application using BioNetGen [6, 29] running on a desktop workstation. CPU time is not reported for the fifth round of rule application, which was performed over the course of several days.