Kinetic Monte Carlo Simulation in Biophysics and Systems Biology

1. Introduction

Generation of large amount of biological data has provided us with some of its unique challenges in modeling of biological systems. Recent advances in imaging techniques, such as confocal, FRET, multi-photon among others, led systems level information for biological processes and accumulation of even more biological data. There is a clear need for quantitative (mathematical and computational) tools that are capable of analyzing such large amount of information. Equally pressing is the need for theoretical models that can elucidate the cellular and molecular mechanisms of complex biological processes. Systems biology modeling (differential equations, stochastic simulations, agent based models) is emerging as a powerful tool that can provide mechanistic insight into workings of complex biological systems as well as allow analysis of large datasets. Such computational models can be validated using synergistic experiments (in an iterative manner) and also possess predictive power. Monte Carlo (MC) simulation is a stochastic modeling approach that can be utilized to carry out in silico studies of biological problems.

Monte Carlo simulation is now an established technique to solve problems in a wide variety of disciplines. It is an essential tool in statistical mechanics and has been frequently used for problems in material science and engineering. Monte Carlo method relies on selecting an initial configuration state and then iterate through a series of steps during which one or more moves are selected and evaluated for acceptance. The number of moves per step is often taken to be the number of available degrees of freedom in that step so as to ensure probing (on the average) of each degree of freedom once per step. However, other types of implementations (such as a rejection-free algorithm) are also possible. Use of Monte Carlo method to solve biological problems has its unique challenges, but if successfully addressed, it can become a powerful tool. Monte Carlo models (for biological systems) can be developed based on statistical mechanics and thus it is expected that it would work well for probing the underlying biophysics of a complex biological process. Of particular interest is a class of Monte Carlo simulations that are based on a set of pre-defined probabilistic rate constants and can be called kinetic Monte Carlo models. Probabilistic rate constant based kinetic Monte Carlo simulations can be developed for a wide variety of biological processes that span multiple length and time scales. We demonstrate that controlled kinetic Monte Carlo experiments are capable of revealing fundamental regulatory mechanisms in a complex biological system.

The organization of the chapter is as follows. We begin with a discussion of the random walk problem in section 2. It elucidates the theoretical basis of Monte Carlo simulations in terms of master equations. It is also an instructive example problem to elucidate some key features, such as the concept of probabilistic simulation parameters (constants), of the kinetic Monte Carlo method. In section 3 Monte Carlo simulation is discussed in the context of physical sciences and materials engineering. This allowed us to discuss how Monte Carlo models can be developed based on statistical mechanics (detailed balance and ergodicity). We also briefly discuss implicit and explicit free energy based implementations of Monte Carlo simulations as applied in physical sciences and materials engineering. In section 4 we describe how Monte Carlo simulations can be used to solve problems in biophysics and biology. A probabilistic rate constant based kinetic Monte Carlo algorithm is discussed in detail. We developed this algorithm to solve the following biological problems: (i) the problem of B cell affinity discrimination and B cell activation (section 5), (ii) systems biology of apoptotic cell death signaling (section 6). We conclude with a discussion on emergence of Monte Carlo simulation as a tool to solve problems of biomedical relevance (section 7). The discussion is focused on our recent efforts to address some key issues in cancer biology and cancer therapy.

2. Monte Carlo simulation of random walk

Random walk is a probabilistic process that is analytically tractable and can also be easily simulated on a computer. Analytical approach is based on master equations that describe the system by equations in probability space. Computer simulation for a random walk is carried out using a set of pre-defined probability constants that are taken from the master equation [1-2]. Random walk based kinetic Monte Carlo simulations can be developed that can describe diffusion and trafficking in biological systems. Below we discuss a simple one-dimensional random walk for which results obtained from Monte Carlo simulations can be readily compared with exact results obtained by analytically solving random walk master equations.

2.1. Master equation description for one-dimensional random walk

At each step, a random walker executing a one-dimensional random walk (on a lattice) jumps to one of the two neighboring lattice sites. One-dimensional random walk is described by a master equation (discrete time and discrete space formulation) of the form: P_n,N+1 = (1/2) P_n-1,N + (1/2) P_n+1,N. P_n,Nis the probability of finding the random walker at site n at time step N [1]. The above master equation describes probability conservation for a symmetric random walk having equal (=1/2) probability of moving left or right. The random walker can reach the site n (at time N+1) from either of the two neighboring sites n-1 or n+1.

Random walk master equation can be solved by defining a characteristic function G = ΣP_n,Ne^ins (the sum is over all possible values of n [-∞,∞]) [1]. By making this change in variable n → s (through the function G), the random walk master equation becomes an algebraic equation G(s,N) = α G(s,N-1), where α = (e^is + e^-is). The equation in G can be solved recursively G(s,N) = α^N G(s,0). For a random walk starting at the origin (n=0 at N=0), G(s,0) = 1, leading to G(s,N) = (e^is + e^-is)^N. The probability function for the random walker P_n,N can then be obtained by expanding G(s,N) in a power series (binomial expansion) and comparing with the starting definition G = ΣP_n,Ne^ins. Various moments of the random walk, such as <n> = 0 and <n²> = N can be obtained directly from G(s,N).

A more general form of the random walk master equation is: P_n,N+1 = p_r P_n-1,N + p_l P_n+1,N + p_wP_n,N; p_r,p_l, and p_w are right jump, left jump, and waiting (at site n) probabilities respectively. This equation is still amenable to exact solution through the characteristic function G(s,N).

2.2. Monte Carlo simulation of a one-dimensional random walk

First, we need to decide on the allowed moves, which are called Monte Carlo moves, in the simulation. For a simple one-dimensional random walk there are two possible Monte Carlo moves: right and left jumps. At each Monte Carlo time step, those two moves are sampled with equal probability and then the chosen move is performed by increasing or decreasing the position, n, by unit length (n → n+1 or n → n-1). We store the information for the position of the random walker in a pre-defined variable (memory location on a computer). The process is then repeated many times by the use of a loop (such as the for loop). Position data (n) is recorded at every time-step or every M (integer) number of MC time steps. Often data is stored in a separate output file location that can be processed later for analyzing simulation data.

Sampling of a probability distribution, such as the equal probability sampling of right / left jumps, can be carried out on a computer using pseudo random number generators. Standard uniform random number generators in a programming platform, such as the rand function in the C library, are often used for this purpose. Uniform random number generators provide a number by sampling a uniform probability density function in the interval (0,1). Equal sampling in case of a one-dimensional symmetric random walk can be achieved by diving the interval (0,1) into two equal parts. At each MC step, the random number is called once, and right / left jump is chosen depending on whether the random number is > or < 0.5 (Head or Tail). Other types of probability distributions can be sampled by utilizing the uniform random number generator.

We can run the random walk simulation many times by generating a different sequence of random numbers in each run. Thus each simulation run corresponds to a distinct random walk trajectory. One can solve the random walk problem by calculating P_n,Nfrom many runs of such random walk simulations; statistical quantities such as <n> or <n²> can also be estimated. A computer program for simulating two-dimensional random walk is provided in the Appendix. Results from simulations should agree with the exact results obtained from the solution of the random walk master equation. For a slightly generalized random walk, having asymmetric right and left jumps and also a waiting probability, the uniform distribution in (0,1) is sampled accordingly (based on P_r, P_l, and P_w). The probability constants P_r, P_l, and P_w can be taken as an input of the simulation. This type of random walk simulations provides a simple example of the probabilistic constant based kinetic Monte Carlo method. Each Monte Carlo step in such a simulation can be assigned a physical time (depending on the system under study); the probability constants P_r, P_l, and P_w would vary with the physical time-scale assigned to a MC step. Periodic boundary conditions can be used to avoid finite size effects.

3. Monte Carlo simulation in physical sciences and materials engineering

Monte Carlo simulation is now an established method in statistical mechanics having wide applicability [2-4]. Monte Carlo moves in such systems typically correspond to changes in configuration space and hence the total energy of the system under consideration. The two key requirements for this class of simulations are ergodicity (of the MC moves) and the detailed balance condition of the acceptance criteria [3].

3.1. Ergodicity and detailed balance

Ergodicity demands that any state in the configuration space of the system can be reached in a finite number of Monte Carlo moves starting from an arbitrary initial state. Such a condition ensures that the configuration space of the system is adequately explored in a given time.

Detailed balance maintains the equilibrium probability distribution for different states by setting a condition for equilibrium. The requirement of detailed balance can be elucidated by considering a two state system. Suppose N₁ is the number of particles in state 1 (or the number of one particle system in state 1 for an ensemble of systems) and N₂ is the same for state 2. S_1→2 represents the probability that state 1 is proposed to move to state 2 (similarly define S_2→1 for 2 → 1) and P_1→2 represents the acceptance probability (similarly P_2→1 for 2 → 1). The condition of detailed balance ensures equilibrium in the following manner:N₁S_1→2P_1→2 = N₂S_2→1P_2→1. Assuming each degree of freedom is sampled once (on average) per Monte Carlo step: S_1→2 = S_2→1. In equilibrium, the particles (or the ensemble) should be distributed according to the Gibbs-Boltzmann distribution so that N_i ~ exp(- E_i/ K_BT). Hence P_1→2 / P_2→1 = exp{-(E₂ – E₁)/K_BT}. In the Metropolis Monte Carlo algorithm, P_2→1 = 1 (assuming E₂> E₁). Thus for any move that results in an increase in free energy the acceptance probability is given by P_1→2= exp{-(E₂ – E₁)/K_BT}.

3.2. Various implementations of Monte Carlo algorithm

1. Explicit free energy simulations: Algorithms with explicit consideration of thermodynamic free energy, such as the Metropolis acceptance criteria, have been widely used to calculate physical properties of a system at equilibrium [2-6]. For systems that are not analytically tractable, such as the Ising model in 3 dimensions and certain disordered systems, Monte Carlo simulation remains a powerful approach. For systems that are driven out-of-equilibrium, Monte Carlo simulation is an essential tool to study the approach-to-equilibrium dynamics in such systems [3].

2. Implicit free energy simulations: Various types of discrete particle based simulations have been employed to study physical properties of systems that remain out-of-equilibrium. Examples of such non-equilibrium systems include surface growth models that show kinetic roughening [7-8]. Non-equilibrium condition is generated in a wide variety of systems such as in molecular beam epitaxy where a surface is grown by vapor deposition. The discrete particle based simulations that are developed to study surface growth and kinetic roughening in non-equilibrium systems can be considered as atomistic (or discrete particle based) kinetic Monte Carlo simulations. These simulations are often based on a set of rules for Monte Carlo moves and any consideration of free energy (if needed) remain implicit. Some of the initial simulations were developed to study universal properties of out-of-equilibrium systems; such universal properties were also studied utilizing coarse-grained field theoretic models [7-10]. However, more complex Monte Carlo models can be developed for realistic simulations of systems in material science, engineering, and nanotechnology [11-15].

4. Monte Carlo simulations in biology

Metropolis Monte Carlo, an algorithm based on explicit consideration of free energy and detailed balance, has been utilized in the past to simulate biological processes such as receptor clustering and cellular signaling. One of the early examples include a model for TNFR1 receptor clustering that is induced by ligand binding [16]. In that study, a Hamiltonian was constructed with a linear chemical energy and binding energy term for each lattice site, as well as an interaction term involving molecules present in adjoining lattice sites. The interaction energy considered in their model was similar to the energy function used in lattice gas models in statistical mechanics [3]. Another study modeled more complex receptor clustering at the T cell synapse [17], in which two types of receptor-ligand pairs cluster at the cell-cell contact but also known to segregate driven by receptor-ligand pair length differences [18-23]. These two types of receptor-ligand pairs, which possess different equilibrium bond lengths, cause membrane bending and tension; thus the Hamiltonian is a sum of elastic energy from the membrane plus the interaction energies of receptors and ligands. In a separate study, the functional role of T cell synapse was explored by a Monte Carlo model of membrane proximal signaling events in T cells [24]. An internal spin state {0 or 1} dependent energy term was considered in the Hamiltonian to model phosphorylation-dephosphorylation type signaling reactions. Examples of Monte Carlo modeling include rolling and adhesion of leukocytes [25] and dynamics of G-protein activation [26].

Another type of Monte Carlo approach, known as Gillespie’s stochastic simulation algorithm (SSA), has been widely used to simulate biochemical signaling networks, especially for spatially homogeneous systems [27-28]. Several attempts have been made to introduce diffusion into Gillespie’s SSA [29]. Some of the initial works that explored stochastic fluctuations in genetic regulatory networks used the SSA for its efficient sampling scheme [28, 30-31]. A similar approach, known as continuous time Monte Carlo, has been employed to study low temperature systems [3, 32]. The advantage of this type of rejection-free algorithms arises from frequent sampling of high probability reaction moves (not wasting computational time by repeatedly rejecting low probability events). When large numbers of signaling reactions are present, such as in a complex cell signaling network, Gillespie’s SSA method might become computationally expensive.

4.1. Study of complex systems in biology: Probabilistic rate constant based kinetic Monte Carlo simulations

Probabilistic rate constant based kinetic Monte Carlo simulations are emerging as a powerful tool to solve biological problems. Though similar implicit free energy simulations have been used in random walk problems as well as in non-equilibrium statistical physics, solving complex biological problems by the Monte Carlo method pose its unique challenges. Finding a suitable parameter mapping scheme, between experimentally measured rate constants and simulation parameters, is an example of such issues that needs to be addressed when simulating biological systems [33]. Below we discuss a kinetic Monte Carlo method that has been developed by us to solve the following problems: (i) antigen affinity discrimination by B lymphocytes, (ii) systems biology of cell death (apoptosis) signaling [34].

In our developed kinetic Monte Carlo method a set of MC moves are carried out with pre-defined probabilistic rate constants. All possible moves are randomly sampled by first randomly sampling individual objects (or agents). Once a move is selected, acceptance / rejection for that move is determined by a pre-defined probability constant. In this method it is possible to satisfy the conditions of detailed balance and ergodicity, thus the model can be placed on a firm basis of statistical mechanics (necessary when thermodynamic considerations are relevant). A parameter mapping scheme was designed that can be used to obtain the probabilistic simulation parameters from known rate constants. Assigning physical time to a MC step is part of this parameter mapping scheme. For lattice simulation, individual nodes are occupied by discrete particles (cells / molecules / agents); such explicit simulation of individual objects (in a temporal model) can be considered agent-based simulations. In our simulations, lattice size and lattice spacing are defined by relevant length scales of the problem. For cell-signaling problems in biology, lattice size will be set by the size of a cell and lattice spacing by size of protein molecules. Our simulation method can be modified, when necessary, to perform off-lattice simulations. At the beginning of the simulation, molecules are either (i) uniformly and randomly distributed, or, (ii) distributed according to their known spatial localization. Even when the molecules are randomly distributed initially, spatial heterogeneity can emerge during the course of simulation.

4.2. Monte Carlo in biophysics and systems biology: Simulation method

1. Random sampling of all possible MC moves and acceptance / rejection criteria: This can be implemented in the following manner. After an object (e.g. molecule) is first randomly sampled, all possible diffusion / reaction moves (degrees of freedom) for that object are randomly sampled to choose one particular move, which is then accepted / rejected based on pre-defined probability constants. Thus the effective probability for a Monte Carlo move is given by P_eff = S P_reaction, where S = 1/N is the sampling probability for N number of allowed Monte Carlo moves and P_reaction= acceptance / rejection probability for the chosen reaction. If we decide to attempt N Monte Carlo moves in a MC step (Δt), then S = 1.

2. Detailed balance and ergodicity: The probabilistic constants (for reaction / diffusion moves) contain the information for free energy changes, associated with a reaction or a diffusion move, in an implicit manner. Let us first consider a bi-molecular reaction of the type A + B → C. We need two probability constants, P_on and P_off, corresponding to the forward and the backward reactions respectively. The ratio P_on / P_off should be equal to exp(-βΔE)to satisfy the detailed balance condition at each spatial point. β = 1/K_BT and ΔE is the free energy change associated with the reaction move (ΔE = E_bound – E_free,E_bound< E_free). Clearly, the reaction probabilities P_onand P_off should be connected to the kinetic reaction rate constants K_onand K_off, the exact relation for which is discussed later in the parameter mapping section. For a class of problems in systems biology, thermodynamic considerations are not necessary and the detailed balance condition to achieve Gibbs-Boltzmann equilibrium distribution need not be satisfied. Rather direct estimation of reaction probabilities (P_reaction) from known rate constants (K_reaction) is more important in such systems biology problems. Ergodicity is satisfied as any state in the configuration space can be reached in a finite number of MC steps starting from an arbitrary initial state of the system.

3. Parameter mapping scheme: Time-scale of the simulation is estimated by considering the time-scale for diffusive motion in the system. We associate P_diffusion = 1 for the fastest diffusing species and determine Δt (MC time step) by matching the diffusion constant (for that species) to the experimentally measured diffusion co-efficient. An approximate estimate for the diffusion constant D ~P_diffusion (Δx)²/ (Δt) can be used, where P_diffusion is the diffusion probability of a given species, (Δx) is the lattice spacing, and, (Δt) is the MC time-step used in the simulation. Reaction probabilities are then estimated using t.

For the A + B → C type reaction, the dissociation reaction probability P_off can be calculated directly from the backward reaction rate constant K_off by using P_off = K_off Δt. Such a mapping would work for a class of reaction rates (expressed in sec^-1 unit) such as the degradation rate or the catalytic conversion rate. Estimation of the association reaction probability P_on (for A + B → C type reaction) is not straightforward as the forward reaction rate K_on has the unit [area sec^-1]. Note, the on rate K_on has two components to it, one of it captures the association rate for two free reactant molecules and a second factor corresponds to the reaction probability [35]. One can use simulations to find a relation between the ratios (P_on / P_off) and (K_on / K_off). Typically, P_off can be varied in a simulation keeping P_on fixed, and the affinity constant K_A (= K_on / K_off) can be obtained from the average steady-state (or equilibrium) value of the reactant molecules. A linear relationship in the form of (P_on / P_off) = α (K_on / K_off) was found, which provides an estimation for the association reaction probability P_on. The constant α captures the effect of diffusion-limited association of free reactant molecules and clearly its value should depend on the spatial dimension of the system. If the value obtained for P_on or P_off turns out to be >1, the MC time-step (Δt) needs to be readjusted accordingly.

1. Simulating complexity in biological systems: Complexities involved in biological processes, such as complex rules of cellular signaling, are handled well by Monte Carlo simulations. In contrast, increasing complexity of a biological system can make developing the correct set of differential equations for that system increasingly difficult. Other examples of such complexity will include formation of multi-molecular complexes (hetero-oligomers) such as the assembly of an apoptosome during apoptotic cell death activation. Our previous studies have shown the effectiveness of kinetic Monte Carlo methods in capturing the complexities of a biological system. In addition, such probabilistic rate constant based methods can be successfully combined with other types of simulation or differential equation based models. We could combine the dynamical equation for membrane shape fluctuations with a kinetic Monte Carlo method describing receptor-ligand binding at the cell-cell contact region [36].We can also couple our probabilistic rate constant based method to explicit free energy based Monte Carlo techniques (such as the Metropolis scheme). Such a hybrid simulation technique was found to be key to developing a Monte Carlo model for receptor-lipid raft formation during B lymphocyte activation.

2. Boundary conditions: Appropriate boundary conditions, depending on the biological problem under consideration, need to be employed in the Monte Carlo simulation code. For example, when simulating an intra-cellular signaling pathway signaling molecules are restricted to diffuse within the cellular volume. Such no flux boundary conditions are implemented by reflective boundary conditions in a lattice simulation. Linear dimensions of the simulation lattice are determined by the size of a cell. Scaled-down versions with smaller lattice sizes are also used in the cases of computationally expensive simulations.

5. B lymphocyte activation and antigen affinity discrimination

Only a few specific lymphocytes (T and B cells), out of billions of such cells, are known to get activated and mount a successful immune response against a particular type of pathogen. How the adaptive immune system carry out such a task remains a problem of considerable interest. B lymphocytes have an ability to adapt to the changing pathogenic load by continuously generating diversity in the variable region of the receptor (somatic hypermutation leading to variations in antigen binding affinity) and selecting the high affinity clones, a feature known as affinity maturation. It is a strategy by which B lymphocytes optimize the immune response given a pathogenic load. Affinity maturation is also key to design of vaccines that rely on B cell mediated antibody production. Earlier studies indicated factors such as competition for antigens could lead to selection of high affinity receptors during the process of affinity maturation. More recent observations, however, focus on antigen affinity discrimination at the level of single cell activation [37-42]. In B cells, affinity discrimination implies increasing level of signaling response as the affinity (K_A=K_on /K_off) for antigen is increased. Such a monotonically increasing response, however, is not obvious as low affinity antigens clearly have the ability to quickly dissociate from a B cell receptor and serially activate a large number of such receptors in a small amount of time. We use Monte Carlo simulations to explore the molecular mechanism of B cell affinity discrimination [36, 43-46].

Mathematical and computational modeling, often in synergy with biological experiments, has been increasingly utilized to solve immunological problems. Previous studies in T cells assumed a series of conformational changes in the T cell receptor upon antigen binding, a concept known as kinetic proofreading, to explain affinity discrimination in T cells [47-49]. Kinetic proofreading was thought to compete against loss in serial triggering (the ability of a single antigen to serially trigger many T cell receptors [50]), as the affinity for antigen is increased. Our initial Monte Carlo studies also utilized such kinetic proofreading requirements in an ad hoc manner. However, explicit simulation of molecular level events, such as B cell receptor oligomerization and lipid raft formation, can allow kinetic proofreading to emerge naturally from simulations. Monte Carlo models seem to be particularly suitable for simulating such molecular level details and biological complexity. Affinity discrimination can be studied in a controlled manner as only the BCR-antigen binding affinity can be varied in our simulations by keeping other parameters fixed.

5.1. Brief description of simulation set up

B cell receptors are placed on a two-dimensional lattice that mimics the B cell surface. Anitgens are also placed on a two-dimensional lattice that represents a lipid bilayer (surrogate for an antigen presenting cell). Diffusion move consists of moving a molecule to one of its four nearest neighbor nodes thus displacing the molecule by unit lattice spacing. Reaction between a B cell receptor and an antigen molecule may take place when they are at the same spatial location (x,y coordinates) on two opposing lattices. The probability of binding / unbinding between a BCR and an antigen should depend on the spatial location of the molecule to take into account the spherical curvature of a B cell. The center of the cell-cell contact region is assumed to adjust itself to the equilibrium bond length of the BCR-antigen pair; thus the probability for the association reaction at that point is maximal. Below we discuss how increasing level of biological complexity can be simulated by our Monte Carlo model:

1. Simulating formation of B cell receptor (BCR) oligomers and oligomerization mediated signaling [46]

   1. We allow two types of MC moves: diffusion and reaction. Once a molecule is sampled randomly, diffusion or reaction is chosen with equal probability. Reaction moves include binding / unbinding between two molecules (such as between a BCR and an antigen), BCR-BCR oligomer formation, phosphorylation by signaling kinases. Simulation procedure is similar to that described in section #.4.2.

   2. Two antigen bound BCRs, when in sufficient spatial proximity, can form dimers and dimerization can eventually lead to formation of higher order oligomers. It is simulated based on the recent experimental observation that BCRs oligomerize due to opening of their Cμ4 domains upon antigen binding [51-52]. We also allow dissociation of such oligomers such as due to antigen dissociation from a dimerized BCR.

   3. In our model, the tyrosine residues in the signaling chains (ITAMs) of BCRs can get phosphorylated only when they are part of an oligomer; Signaling kinase Lyn can then bind and phosphorylate the tyrosine residues. Such a requirement for oligomerization introduces a delay time before activating signal can propagate downstream and thus provides a basis for kinetic proofreading.

   4. Syk (Spleen tyrosine kinase) can only bind to phosphorylated signaling chains of BCR molecules and get phosphorylated. Syk phosphorylation propagates the activation signal downstream.

   5. We are also able to simulate diffusion of various types of molecules (free and antigen bound) in a realistic manner by varying the probability of diffusion move. BCR oligomers, for example, are assumed to be immobile.

   6. Affinity discrimination can be studied by running Monte Carlo experiments where BCR-antigen affinity is varied but all the other parameters remain constant. Simulation results indicate that formation of BCR oligomers is affinity dependent in a manner that can lead to enhanced signaling with increasing affinity.

The complexity of lipid-mediated interactions cannot be captured by the oligomerization model, however, such interactions are instrumental to understanding how signaling kinase Lyn can access antigen bound BCRs (upon antigen stimulation). Lipid rafts are sphingolipid and cholesterol enriched membrane microdomains [53]. Large (> 100 nm) and stable BCR-lipid rafts are known to form upon antigen binding [54-57], whereas in resting B cells lipid rafts are small and transient.

1. Simulating lipid mediated interactions and formation of BCR-lipid rafts

   1. BCR-antigen binding is simulated by the probabilistic rate constant based kinetic Monte Carlo model. We couple the kinetic MC model with an explicit free-energy based model for lipid mediated interactions [58].

   2. In the explicit free-energy based model, we simulate BCR-BCR and BCR-lipid interactions through two energy-based parameters BB and BL. BB captures pairwise BCR-BCR interaction energy such as that is needed to form BCR oligomers. BL denotes the interaction between a BCR molecule and a raft-forming sphingolipid. We can increase the strength of those parameters upon antigen binding and follow the dynamics of BCR-BCR and BCR-lipid cluster formation. Variations in BB and BL parameters, which depend on the state of antigen binding, create a coupling between the implicit free-energy and explicit free-energy Monte Carlo models.

   3. For antigen bound BCRs, origin of large BB lies in the opening of the Cμ4 domains of BCRs upon antigen binding. As for increased BL (upon antigen binding), in addition to energetics and entropy, effects of membrane curvature can be important.

   4. Lipid-lipid interactions, such as one that arises due to hydrophobic mismatch, are captured by an energy-based parameter LL.

   5. In resting B cells, Lyn is already in sphingolipid rich regions due to its dual acyl chains. We can explicitly simulate the spatial localization of Lyn kinase in sphingolipid rich regions through a KL parameter. When there is no antigen, BB and BL are small, and BCRs partition into non-raft regions of the membrane. Thus BCR and Lyn are segregated into distinct spatial regions in resting B cells. This provides a mechanism for inhibition of spontaneous activation in resting B cells.

   6. Antigen binding leads to enhanced BB and BL, generation of BCR-sphingolipid rafts, and co-localization of BCRs with signaling kinase Lyn. The delay in BCR-Lyn association, upon antigen binding, provides a key step in kinetic proofreading mechanism. Thus selective partitioning and recruitment of signaling molecules, in different spatial regions on the B cell membrane, emerge as a basis for kinetic proofreading. Interestingly, the requirement of kinetic proofreading time decreases with increasing affinity thus favoring increasingly higher affinity for antigens.

It becomes evident that Monte Carlo simulations can handle spatial heterogeneity and clustering, a key aspect of biological complexity, in a rigorous manner. In addition, complex regulatory mechanisms of B cell signaling, such as the signal transduction through distinct phosphorylation sites on Syk, can be incorporated in a Monte Carlo simulation. Monte Carlo method is also suitable for capturing stochastic effects that are unavoidable when antigen concentration is low. Such antigen limiting conditions might arise, for example, during the initial phase of an infection. Stochastic effects are also going to be important when BCR concentration is low, such as in pre-plasma cells. Moreover, Monte Carlo models can be used to probe the underlying biophysics of BCR-lipid raft domain formation (such as phase separation, scaling behavior in domain growth, and the nature of criticality).

5.2. Analysis of single cell data

Phosphorylation of BCR signaling chains (pBCR) and Syk (pSyk) are measured in our simulations as readouts for membrane proximal signaling. Each run of Monte Carlo simulations provide us activation data at the level of single cells and thus can be readily compared with data obtained from single cell experiments. Standard statistical measures, such as the average and the second moment (or standard deviation), can be calculated for relevant variables from many runs of the simulation. If simulation data for large number of runs is available, then the full probability distribution for a physical variable can be calculated and compared with experimental histograms obtained from, for example, flow cytometric measurements. In our analysis, probability distributions for pBCR (Figure 1) and pSyk are generated for various affinity values and affinity discrimination is obtained from the separation of probability distributions.

Figure 1.

Histogram plots of the number of BCR molecules with at least one ITAM phosphorylated by Lyn (figure re-produced from our earlier published article in Journal of Theoretical Biology 2012;307: 174-182; reference [46]). Data from one thousand single cell runs were used (at t =, 10⁵ MC step); simulation parameters are listed in Table 2 of [46]. The number of Lyn-phosphorylated BCRs, a measure of early-time membrane proximal signaling, increases with affinity.

We proposed the following metric that can be used to analyze the probability distributions obtained from simulation data and quantitatively characterize affinity discrimination: Δ = overlap area / (m₁ – m₂); the area of overlap between the histograms for two affinity values is divided by m₁ and m₂, the histograms’ mean values for those two affinities. Lower Δ values correspond to better affinity discrimination, with the best discrimination occurring at Δ=0 (no overlap between histograms). When Δ=0, it becomes necessary to compare the mean values of the histograms. From Fig. 1, significant affinity discrimination is obtained even between affinity values K_A = 10⁸ M^-1 and K_A = 10⁹ M^-1 (Δ values are provided in Fig. 6 of reference [46]). The quantitative metric Δ can also be used to analyze single cell experimental data for B cell activation (such as the pSyk data in fig 8 in reference [42]).

5.3. Hybrid simulation methods

Our probabilistic rate constant based kinetic Monte Carlo method has the advantage that it can be coupled to other types of Monte Carlo or differential equation based models. Monte Carlo simulations can be computationally expensive, especially when large-scale intracellular signaling pathways are considered. For many such signaling pathways, presence of large number of molecules make stochastic fluctuations less significant and differential equation based approaches are often sufficient to describe such systems [24, 59]. It is also assumed that differential equation based approaches would be able to handle the biological complexity. Whether differential equation based models are valid descriptions also depends on the biological question one wants to address. In situations where membrane proximal events activate the downstream signaling, such as antigen binding mediated activation of B cell intracellular signaling pathways, spatial heterogeneity and other complexities might demand application of Monte Carlo type approaches for the membrane events but differential equations might be sufficient for describing the intracellular signaling processes. In such cases, a hybrid simulation method, which combines a Monte Carlo simulation with an ODE / PDE based model, can be an efficient modeling approach. We have developed two hybrid simulations in the context of studying B cell activation:

1. Modeling membrane deformation as a result of receptor-ligand binding: A time dependent Landau-Ginzburg type equation, which describes membrane shape fluctuations, was coupled to a probabilistic rate constant based kinetic MC model that simulates receptor-ligand binding dynamics. The most widespread approach for modeling membrane deformation due to receptor-ligand binding is based on modeling the free energy of the membrane as a function of receptor-ligand bond stretching and mechanical restitution forces [17, 20-23].The change in the local intermembrane membrane separation distance can be assumed to evolve (towards the free-energy minimum) according to a time-dependent Landau-Ginzburg formulation. A random noise term in the equation captures the effect of thermal fluctuations (variance of noise ~ K_BT). For the purpose of calculating local concentrations of receptor-ligand complexes, the membrane surface is discretized into square subdomains over which membrane separation remains approximately constant. The spatially averaged concentration of receptor-ligand complexes in each of these subdomains is then calculated from the kinetic Monte Carlo model and entered in the discrete form of the membrane equation [44-45]. In our Monte Carlo model, memebrane separation is updated at the end of each MC time step by solving the membrane equation.

2. Modeling lipid mediated interactions and BCR-lipid raft (microcluster) formation:Here, a kinetic MC model (for receptor-ligand binding) is linked to an explicit free-energy based MC model for lipid mediated interactions [58]. B cell receptors (BCRs) need to be modeled as a two-layer structure: (i) the bivalent Fab domains move in the first layer where it can bind monovalent antigens and (ii) the transmembrane part diffuse in the second layer with lipids and other proteins. Diffusion moves in this second layer is accepted / rejected based on an explicit free-energy based mechanism (such as the Metropolis scheme). Src family kinase Lyn is also simulated in a manner that can interact with lipids. The hybrid Monte Carlo simulation proceeds as follows:

   1. At the beginning of the simulation, all the protein molecules and lipids are distributed uniformly and randomly on their respective 2d lattice grids. Any clustering of proteins and lipids emerges from our simulation and MC simulations are well suited to capture such spatial heterogeneity.

   2. Lipids are placed on a smaller sized lattice grid than that for protein molecules such as BCRs. Lipids, however, are going to be sampled n times more frequently than protein molecules (n is determined by relative diffusion). Diffusion constants for lipids and BCRs (or other proteins) will be used to set the time-step Δt of MC simulations. Lipid diffusion moves are determined by lipid-lipid and lipid-protein free energy of interactions.

   3. Diffusion of antigen molecules is carried out only with the constraint of mutual physical exclusion. Diffusion of BCRs, however, has additional free energy considerations arising from BCR-BCR and BCR-lipid interactions.

   4. Diffusion of Lyn molecules is also governed by free energy changes associated with Lyn’s affinity to raft forming sphingolipids. Lyn prefers to cluster with spingolipids due to increased hydrophobicity caused by its doubly acylated (palmitoylation and myristoylation) form. We assume interdigitation between the two leaflets will allow Lyn to get associated with BCR-lipid raft domains [60].

   5. It is technically challenging to simulate two disparate size objects (lipids and BCRs) on the same 2-dimensional lattice. Two possible ways to implement the simulation:

      1. One BCR and a significantly larger number of lipids will be placed on a relatively large lattice block (approximately, 20 nm x 20 nm). We assume the lone BCR in a lattice block will get a chance to interact with all other lipids in the same block through short-time (<ΔT of MC step) rotational/translational displacements.

      2. Another potential solution is to define BCR and lipids blocks (approximately 10 nm x 10 nm) on a 2-d lattice and allow exchange of blocks (Monte Carlo moves) in a manner that will (i) allow one BCR molecule to interact with several lipid molecules and (ii) satisfy detailed balance. Such an approximation is justified by noting that the relevant length scale of the problem, namely the size of generated raft domains (> 100 nm), is significantly larger than size of the unit blocks (~ 10 nm). A schematic of the simulation is provided on the right. Initial results from simulations of two distinct affinity values (K_A= 10⁷ M^-1 and K_A = 10⁹ M^-1) are shown in Fig 2.

Figure 2.

(A) Schematic of the BCR-lipid raft formation simulation. (B) Time-course of BCR-lipid raft formation is simulated for two affinity values, low affinity (top): K_A = 10⁷ M^-1 (P_on = 10^-3 and P_off=10^-6) and high affinity (bottom): K_A = 10⁹ M^-1 (P_on = 10^-2 and P_off=10^-7). Affinity discrimination has been studied experimentally for similar affinity values: K_A = 9.9 x 10⁶ M^-1 (B1-8-Low antibody) and K_A = 5.2 x 10⁸ M^-1 (B1-8-High antibody) for the antigen hapten NIP [42]. We use the following energy-based parameter values in our simulation: BB = 3.0, BL = 2.0 and LL = 2.0 (in K_BT); BB and BL values are for antigen bound BCRs. Both unbound (blue) and bound BCRs (red) are shown at the cell-cell contact region. Raft forming lipids (not shown) are co-clustered with the B cell receptors. Significantly larger amount of BCRs are clustered in the high affinity case.

The underlying physics of BCR-lipid raft formation is similar to the phase separation process in a two liquid system. BCRs and raft forming lipids can be treated as one type of liquid molecules that phase separates from other lipids and proteins. In the case of two liquids, phase separation is typically driven by a low temperature quench that initially drives the system out-of-equilibrium. Some aspects of the out-of-equilibrium dynamics (towards an equilibrium state) studied in the context of liquid-liquid phase separations [3] can be useful here. However, complexity of the BCR-lipid raft formation makes any simulation and its analysis more challenging. Here, antigen binding drives the system out of its steady state (small unstable rafts). As long as we do not consider any active processes, the problem of BCR-lipid raft formation can be considered as out-of-equilibrium dynamics towards an effective equilibrium state. BCRs are assumed to have significantly increased mutual attractive interactions (BB) upon antigen binding, presumably due to opening of the Cμ4 domains, introducing affinity into the problem of raft formation. As expected, the phase diagram is also affinity dependent elucidating the physical basis for early-time affinity discrimination in B cells.

5.4. Simulating B cell affinity discrimination: Role of spatial heterogeneity

Recent work by us and others have demonstrated that early time (~ 100 seconds) membrane proximal events, BCR oligomerization, BCR-lipid raft formation, and B cell immune synapse formation, hold the key to B cell affinity discrimination problem. BCR-lipid raft structures start forming within tens of seconds after antigen binding and precede the formation of B cell immune synapse (~ 100 seconds). Synaptic structure consists of a central cluster of BCR-antigen complexes surrounded by a peripheral ring of integrin complexes (LFA1 –ICAM1) [61-62]. Spatial clustering and heterogeneity, in the form of lipid rafts or immunological synapses, presumably modulate serial activation (serial triggering) and kinetic proofreading effects and determine the nature of affinity discrimination in B cells. It is apparent that Monte Carlo models can capture the relevant details of spatial clustering and spatial heterogeneity, even when antigen or receptor concentration is low. We discussed how to combine implicit and explicit free energy MC models in a hybrid simulation scheme to simulate BCR-lipid raft formation. Detailed method for simulating B cell immune synapses can be found in our earlier work [36, 43-44]. We showed that, a directed transport (of antigen bound receptors) based mechanism is needed for synaptic pattern formation, as the size difference between BCR-antigen (at least for IgM) and integrin complexes is not significant. Further work needs to done that can link BCR-lipid raft formation to immunological synapse formation in B cells. Whether affinity dependent signaling through rafts, as found in our simulations, can generate an affinity dependent directed transport mechanism remains to be explored.

6. Systems biology of apoptotic cell death signaling

Programmed cell death, apoptosis, is one of the most important cellular processes. Apoptotic death is critical to a wide variety of cellular and physiological phenomena ranging from the normal development of multicellular organisms to maintaining homeostasis [63-64]. Dysregulated apoptosis has been implicated in a large number of diseases including in cancer and degenerative disorders. Targeting the apoptosis pathway is emerging as a new frontier in the therapeutic approaches for those diseases. Large number of signaling species and an intricate network structure generate complexity in the apoptotic cell death signaling pathway and make any computational study of apoptosis a challenging task. Apoptosis pathway has evolved to sense and respond to a wide variety of stimuli through structurally similar signaling molecules having similar pro- or anti- apoptotic functions. The system level regulation of apoptosis signaling is achieved through a loop network structure that combines two distinct pathways known as Type 1 (extrinsic) and Type 2 (intrinsic or mitochondrial) pathways (Figure 3). In addition, several smaller loop network structures exist within the Type 2 pathway that seem capable of generating nonlinear effects in the signaling response. Such loop networks are different from well studied feed forward and feed back loops that are frequently encountered in signaling pathways. Additional complexity in the apoptotic pathway arises due to the formation of multi-molecular complex apoptosome. Translocation of pro- and anti- apoptotic proteins to the mitochondrial membrane, such as that happen for activated Bax dimers, generate spatial heterogeneity in the system.

6.1. Monte Carlo model for the apoptosis signaling pathway

We developed a fairly detailed Monte Carlo model to study the systems biology of apoptotic cell death signaling [65]. We simulate a single cell using a cubic lattice in which all the molecules were placed at different nodes of the lattice. Membrane bound molecules such as the death inducing signaling complex (DISC) are confined to one surface of the cubic box and allowed to diffuse only in two dimensions. Intracellular signaling molecules are allowed to diffuse inside the cubic box when they are in free state; multi-molecular complexes are assumed to be immobile.

We simulate the signaling network shown in Fig. 3. Our model considers apoptotic signaling through two distinct pathways: direct activation of caspase 3 by caspase 8 (type 1) and activation of caspase 3 through mitochondrial cytochrome c release and apoptosome formation (type 2) (Fig. 3). In our model, intracellular apoptosis signalling was triggered by the activation of caspase 8 molecules (at the cell surface) that, in turn, diffused in the cytosol and activated both pathways of apoptosis signaling. In the type 1 pathway, caspase 8 molecules directly catalyze the cleavage reaction of procaspase 3 to generate active form of caspase 3. In the type 2 pathway, caspase 8 binds with Bid and catalyzes its truncation to form tBid that, in turn, binds to Bax to generate Bax2 homodimers. Bcl2, an antiapoptotic molecule, can inhibit both tBid and Bax, thus creates a local loop structure in the type 2 signaling cascade. Formation of the active Bax2 complex leads to cytochrome c release from mitochondria and formation of multi-molecular Cytochorme c-Apaf1-ATP complex apoptosome, which ultimately induces activation of caspase 3. Hence, caspase 3 activation at the end of both the type 1 and the type 2 pathways creates a global loop structure in apoptosis signaling. Values for kinetic rate constants and molecular concentrations are taken from previous measurements (and/or used in previous theoretical studies) [66-69].

Figure 3.

Schematic of a simplified signalling network for the apoptotic cell death pathway.

6.2. Simulation code for apoptosis signaling

We define several generalized functions to carry out diffusion or reaction moves. Those functions are called when a molecule is randomly sampled and a MC move (diffusion or reaction) is attempted. Separate diffusion functions are defined depending on, for example, whether a molecule is diffusing in the intracellular space (within cytosol) or on the mitochondrial surface. Similarly, separate functions are assigned for distinct reaction types such as for binding and unbinding reactions. When a given molecule is known to stay bound with more than one molecule (of the same type or different types), we assign multiple binding sites for that molecule. Ability to bind to multiple partners (at the same time) allows us to simulate oligomeric complexes such as the apoptosome. At each MC time step, N random sampling of molecules are carried out, where N is the total number of molecules in the system; hence, on average one molecule is sampled once within a MC time step. We define separate functions for each molecular species. Once a molecule is randomly chosen, diffusion or reaction move is attempted by calling the function for that molecule. As an example, we provide the function for Caspase 8 molecule in the Appendix.

Simulating the type 2 pathway can be computationally expensive due to the presence of low probability reactions, such as the Bid-Bax reaction or the formation of apoptosome. One potential modification to improve the speed of the simulation is to partition the signaling network into groups of fast and slow reactions and then sample the slow reactions less frequently than the fast reactions [70-71].

6.3. In silico experiments of apoptosis signaling

We wanted to study the mechanism of system level regulations in apoptosis signaling. In silico experiments were designed to address the following biological question: how do the two distinct apoptotic pathways get activated? A related question is why two very distinct time-scales are observed in biological experiments of apoptosis [72]. In order to address these questions we designed three sets of experiments

1. Study the type 1 pathway only (set the rate constant for the type 2 pathway to zero).

2. Study the type 2 pathway only (set the rate constant for the type 1 pathway to zero).

3. Keep both the pathways open for activation but vary the strength of apoptotic stimuli (by varying the level of pro-caspase 8 or caspase 8).

One of the advantages of the Monte Carlo simulation approach is that clean experiments can be carried out in silico, such as selective activation of one of the two pathways of apoptosis (by simply setting the rate constant for other pathway to zero) to study its signaling behavior. In silico experiments revealed some of the fundamental regulatory mechanisms in apoptosis biology [65].

Pure type 1 activation: For the type 1 pathway, population averaged (over many cells) activation data could capture the fast (seconds – minutes) caspase 3 activation observed for all the cells in our simulations. Decreasing strength of stimuli, i.e. decreasing caspase 8 concentrations, leads to slower activation of caspase 3 (Figure 4A). Stochastic fluctuations are starting to be seen for very low (< nanomolar) caspase 8 concentrations.

Pure type 2 activation: In contrast to type 1 behavior, caspase 3 activation in individual cell simulations for type 2 showed slow activation with large cell-to-cell stochastic variability. Such stochastic variability in type 2 apoptosis activation was characterized by all-or-none type activation (of capsase 3) at the level of single cells with large variability in the activation time (minutes – hours) (Figure 4B). Stochastic fluctuations in type 2 activation is caused by (i) Bcl2 inhibition of tBid-Bax reaction and (ii) low probability of apoptosome formation. When caspase 8 concentration is low, stochastic fluctuations are more pronounced due to additional stochastic variability arising from the reaction Bid truncation to tBid. Population average (over many cells) behavior cannot capture the all-or-none type rapid caspase 3 activation observed at the level of single cells. Also note, the information of the strength of apoptotic stimuli is lost at the level of single cell capsase3 activation, but could be captured by the probability distribution of the time-scale of caspase 3 activation.

Figure 4.

Time course of caspase 3 activation, as readout for apoptotic activation, is shown for single cells. Data is shown for simulation of the pure type 1 activation (4A) and the pure type 2 activation (4B). 1 Monte Carlo step = 10^-4 sec. Normalized data is shown for caspase-8 = 1 (left) and caspase-8 = 10 (right). Nano-molar concentration is obtained by multiplying the number of molecules with 1.67. Data is shown for 10 individual single cell runs. Data is similar to our earlier published result (Figure 2) in Biophysical Journal 2008;95: 3559-3562 [65].

Combined type 1 and 2 activation: Under weak apoptotic stimuli, such as for very low concentrations of active caspase 8, stronger binding affinity between caspase 8 and Bid, compared with that between caspase 8 and pro-caspase 3 in the type 1 pathway, leads to selective activation of the type 2 pathway with its characteristic slow activation and large cell-to-cell variability (as seen in pure type 2 activation). Clearly, the all-or-none type behavior in caspase 3 activation, and its rapid completion, shows that the type 2 pathway is designed to amplify an initially weak signaling. However, type 1 activation (accompanied by partial activation of the type 2 pathway) is observed in our simulations, even for low caspase 8 concentrations, if the following conditions are met: [1] active caspase 3 can feedback its activation by directly activating caspase9, (ii) Xiap is inhibited by Smac released in the type 2 pathway, (iii) rate of multi-molecular complex apoptosome formation is very slow. How these cellular factors impact the value of caspase 8 at which the type2 → type 1 transition occurs (depending on the cell type) needs to be explored further. For example, presence of apoptosome inhibiting anti-apoptotic proteins in certain cell types (such as neuroglobin in neural cells [73]) and the cytosolic condition (such as the oxidative condition or the pH level) would affect the rate of apoptosome formation and hence type 2 activation. For strong apoptotic stimuli, large number of available caspase 8 molecules (even after binding with Bid) directly binds and activates pro-caspase 3 to active capsase 3 in a fast deterministic manner. For intermediate strength of apoptotic stimuli, initial slow and gradual activation of capsase 3 through the type 1 pathway can suddenly change to a mode of rapid activation through the type 2 pathway (once the apoptosome has formed). Time-scale of apoptotic activation, as emerges from our single cell simulations, is comparable to that observed in biological experiments. Better estimates can be obtained as quantitative information for the kinetic rate constants and the expression levels of signaling proteins are gathered.

Thus, in silico experiments reveal a unique system level design mechanism for apoptosis signaling; weak apoptotic stimuli activates the type 2 pathway in a slow stochastic manner, whereas under strong stimuli the type 1 pathway gets activated in a fast deterministic manner. We also used a minimal network based approach, to provide further insight into the fundamental design mechanism of apoptosis signaling, and to demonstrate robustness of some of our results [74]. Previously it was thought that DISC (death-inducing-signaling-complex) formation and caspase 8 activation could lead to activation of either the type 1 or the type 2 pathway in a cell type intrinsic manner; cells were labeled either as a type 1 cell or as a type 2 cell depending on their mode of activation [72]. More recent observations, however, indicate that both the pathways can be activated irrespective of cell types [75], as observed in our earlier Monte Carlo studies [65] and clearly elucidated by a minimal network based analysis [74]. Cell type specificity (in type 1 / type 2 choice) arises from differential expression of death receptors in different cell types [75], which presumably modulates the activation level of caspase 8, a key determinant of type2-type1 switch [65,72,74-75]. Whether variations in the concentration of death ligands, even when the death receptor expressions are comparable, can determine the mode of activation (type 1 or type 2), needs to be investigated. Additional complexities (in type 1 / type 2 choice) are generated due to, for example, variations in cellular protein levels such as the Xiap concentration (depending on cell types and cellular conditions) [69]. Lipid mediated interactions can also be important in determining the cell type specificity of type 1 / type 2 choice in apoptosis signaling. These additional complexities are currently being incorporated into our Monte Carlo model. Stochastic type 2 to deterministic type 1 transition has important ramifications for cancer biology and cancer therapy [76] and will be discussed further in the next section.

6.4. Analysis of single cell simulation data for apoptosis signaling

Probability distributions, such as that for activated caspase 3, are obtained from many runs of single cell simulations. Interestingly, probability distributions for type 1 and type 2 signaling show two distinct functional behaviors. For signaling through the type 1 pathway, probability distribution shows a single peak that shifts towards larger caspase 3 value as time increases. The value of caspase 3 corresponding to the peak of the probability distribution indicates the average level of caspase 3 activation. For type 2 signaling, large cell-to-cell variability coupled with rapid all-or-none type activation of caspase 3 leads to a bi-modal distribution [65].

Assuming a perfect bi-modal distribution for caspase 3, the ratio variance/average for capsase 3 activity can be estimated by the metric C3[1-f(e,t), where C3 is the initial number of pro-caspase 3 molecules and f(e,t) is the fraction of cells in which caspase 3 has undergone complete activation within a given time t [77-78]. Another quantitative measure that can be used to assess the single cell variability in apoptotic activation is cell-to-cell variability in time-to-death. Such a quantitative measure is also useful in estimating the relative contributions of inherent stochastic variability and stochastic variability arising from cellular variations in protein levels. Cellular variability in protein levels can arise due to genetic and/or epigenetic variations. Stochastic gene expression can also generate cellular variability in protein levels.

7. Applications of Monte Carlo simulation in systems biology and systems medicine

Systems level approaches can be utilized to study complex biological processes and have been increasingly used in the recent past. Additional impetus to develop such systems level models comes from the recent generation of large amount of biological data and an urgent need for systems level computational models that are capable of analyzing this huge amount of biological data. Monte Carlo models are suitable for handling large number of degrees of freedom in a system and it is becoming an essential tool in systems biology. Monte Carlo simulations could capture the following key aspects of a biological process: (i) stochastic variability, (ii) spatial heterogeneity, and (iii) biological complexity (as demonstrated for the problems discussed in the previous section). For a class of systems biology problems that do not require consideration of free energy, instead a coarse-grained phenomenological approach is more appropriate, suitable kinetic Monte Carlo simulations can be developed for such systems. Examples of such problems would include multi-cell simulations in a tissue, and, some of the diffusion and trafficking problems in biology (such as trafficking of immune cells). Probabilistic rate constants that need to be used in such type of Monte Carlo models can be extracted from biological data. For problems in systems biology, Monte Carlo technique can also be used for rapidly carrying out parametric variation study once the MC model has been calibrated against experimental data. Such parametric variation analysis is also essential in cases where the kinetic rate constants are only approximately known.

7.1. Cancer biology and cancer therapy: Insights from Monte Carlo studies

Modeling diseases, as well as designing therapeutic approaches, is an area where systems level approaches and Monte Carlo simulation will become increasingly important. As a prominent example, here, we discuss the impact of systems level computational modeling in cancer biology and cancer therapy. It is now established that a hallmark of cancer is loss of apoptotic regulation in cancerous cells (along with aberrant growth signaling). Thus it is expected that the system level regulatory mechanisms and cell-to-cell stochastic variability in apoptosis, as emerged from recent studies by us and other labs, will have significant ramifications in cancer biology. Major findings from our Monte Carlo modeling of apoptosis and cancer are summarized below.

Bcl-2 overexpression increases cell-to-cell stochastic variability in apoptotic activation: This is an example of parametric variation study that can be carried out using Monte Carlo simulations. Increase in Bcl-2 concentration leads to increasing inhibition of both tBid and Bax, and leaves only a few number of available reactant molecules for the tBid-Bax reaction; such low reactant concentrations generates stochastic fluctuations and slow activation in the pre-mitochondrial signaling module of the type 2 pathway. In addition, simulation results indicate non-linear effect of Bcl-2 variation due to its simultaneous inhibition of tBid and Bax through the tBid-Bax-Bcl2 loop network. Such slow activation of the death pathway under Bcl-2 inhibition, as well as its non-linear effect, explains apoptosis resistance of cancer cells equipped with higher levels of Bcl-2 proteins [77]. Bcl-2 over-expression has been observed in a variety of cancer cells and is a marker for poor prognosis [79-80]. Increase in cell-to-cell variability in apoptotic activation, due to increased expression of Bcl-2, might allow a normal cell (out of many such cells) that is particularly slow to activate the apoptotic pathway to acquire tumor initiating features. In a similar manner, a cancer cell equipped with higher Bcl-2 level might acquire additional mutations to generate a more malignant genetic subclone.

A low probability Bid-Bax reaction can allow selective killing of cancer cells: Finding the mechanisms for inducing selective apoptotic death of cancer cells, which would leave normal cells unharmed, can be key to designing successful anti-cancer therapy. Monte Carlo models can simulate activation of the apoptotic death pathway, under the action of a chemotherapeutic agent (such as a Bcl-2 inhibitor), for both (i) normal cells and (ii) cancer cells. Such simulations can provide information regarding specificity (selective killing of cancer cells); simulation results can also be used to determine the optimal dose for a specific anti-cancer agent. We have studied the mechanism of BH3 mimetic Bcl-2 inhibitors that are considered as potential anti-cancer agents [81-82]. Monte Carlo simulations indicate a low-probability for the direct Bid-Bax type reaction could allow selective killing of cancer cells, which are presumably equipped with higher levels of Bid and Bax molecules, under the action of a Bcl-2 inhibitor (Figure 5). Such a low probability Bid-Bax reaction is also a mechanism for generating cell-to-cell variability in apoptotic activation and highly relevant for apoptotic death of cancer cells.

Figure 5.

Apoptotic activation, under the action of 3 μm HA14-1 (a BH3 mimetic Bcl-2 protein inhibitor [77-78]) is shown for normal (left) and cancer (right) cells. Time course of caspase 3 activation, as readout for downstream apoptotic activation, is shown for single cells. 1 Monte Carlo step = 10^-4 sec. Nano-molar concentration is obtained by multiplying the number of molecules with 1.67. Data is shown for 10 individual single cell runs (only a fraction of cells show activation within the given simulation time).

Simultaneous inhibition of Bcl-2 and Xiap in cancer therapy: This is an example of applying systems level computational approaches in designing cancer therapy. A key challenge that cancer therapy needs to address is the issue of cell-to-cell stochastic variability in apoptotic death of cancer cells. Cell-to-cell variability is linked to the apoptosis resistance of cancer cells and modulated by key anti-apoptotic proteins such as Bcl-2 (pre-mitochondrial) and Xiap (post-mitochondrial). Thus combinatorial targeting of the apoptotic pathway at different signaling modules can be a novel strategy in cancer therapy, and, systems level approaches can elucidate potential options. We carried out Monte Carlo simulations of apoptotic activation for various combinations of Bcl-2 and Xiap inhibitor concentrations. Such a study can provide us with an estimate of optimal inhibitor concentrations for reducing (i) cell-to-cell stochastic variability (in the type 2 pathway) and (ii) toxicity to healthy cells [78].

7.2. Systems level Monte Carlo modeling: a tool to determine optimal strategies in cancer therapy

Monte Carlo simulation is a tool to determine optimal strategies to perturb the apoptotic network at a system level; controlled Monte Carlo experiments of cancer cell death can be carried out in silico. The following key issues needs be addressed:

1. How to maximize specificity (maximize selectivity for cancer cells and minimize toxicity to healthy cells)

2. How to minimize inherent cell-to-cell variability

3. Compare various potential options

Our recent studies indicate that a stochastic to deterministic transition can be achieved selectively for cancer cells, which are equipped with over-expressed pro-apoptotic proteins, in the type 2 pathway of apoptosis. Alternatively, switching the activation from type 2 to type 1 can also lead to deterministic activation of the apoptotic pathway (in cancer cells), as revealed by our earlier systems level studies of apoptosis signaling. In both cases, inherent cell-to-cell variability would protect normal cells from harmful toxicity. Cell-to-cell stochastic fluctuations in the type 2 pathway can also be exploited to protect cells in degenerative disorders. In cancer therapy, it is also important to compare among various possible options and the optimal strategy might depend on the cancer type under consideration. For example, type 2 to type 1 transition in TRAIL resistant cancer cells can be achieved by hedgehog inhibition [76] or possibly by a Xiap inhibitor. Monte Carlo simulations can be carried out for various possible options to determine the optimal one.

Appendix

Simulation code for 2-dimensional random walk

Caspase-8 function (Apoptosis cell signaling simulation)

Acknowledgments

The author acknowledges Somkanya C. Das (for Bid-Bax work in apoptosis) and Paul Yu-Yang (for BCR-lipid raft simulations in B cell affinity discrimination). Work carried out by former lab members and rotation students is also acknowledged.