Computational Cell Biology: Spatiotemporal Simulation of Cellular Events
Slepchenko et al. (2002): Computational Cell Biology: Spatiotemporal Simulation of Cellular Events

This is an introduction to Computational Cell Biology focusing on the system the authors developed, which is called "Virtual Cell". It also mentions several other programs, in particular StochSim and MCell. To illustrate their ideas, the authors provide examples respective to "RNA trafficking" and "neuronal calcium dynamics".

The paper first mentions a couple of pieces of technology that have contributed to the progress of Cell Biology in general in the past twenty years:

Confocal and two-photon excited fluorescence microscopies permit investigators to study the structure and dynamics of living cells with submicrometer three-dimensional (3D) spatial resolution and with time resolutions as fast as milliseconds. These quantitative microscopies can be combined with fluorescent indicators and fluorescent protein constructs to enable the study of the spatiotemporal behavior of individual molecules in cells. Patch clamp electrophysiological recording can be used to study ion currents through single-channel proteins or across the entire cell membrane. All these techniques can be further combined with methods to impart specific perturbations to cells such as photorelease of caged compounds to deliver controlled doses of second messengers or laser tweezer manipulations to determine the response of cells to mechanical stresses.

With all these advances, scientists have gained the following data:

Massive structural biology efforts have produced extensive databases of 3D protein structures. High-throughput molecular biology and molecular genetics technologies have led to descriptions of the full genomes of several organisms, including, of course, the human genome. More recently, highthroughput proteomics technologies promise to catalog, for a given state of a given cell, the dynamic levels of and interactions between all proteins and their posttranslational modifications.

To "link all the molecular-level data to the cellular processes that can be probed with the microscope", computational approaches are needed.

Regarding the mathematical knowledge required to implement these approaches, the authors write:

The concentrations of reacting molecular species as a function of time in a well-mixed reactor can be obtained by solving ordinary differential equations (ODEs) that specify the rate of change of each species as a function of the concentrations of the molecules in the system. If membrane transport and electrical potential are to be included in the model, the rate expressions can become more complex but can still be formulated in terms of a system of ODEs. However, when diffusion of molecules within the complex geometry of a cell is also considered, the resultant “reaction/diffusion” system requires the solution of partial differential equations (PDEs) that describe variations in concentration over space and time.
The finite volume method, developed originally for problems in heat transfer, is especially well-suited to simulations in cell biological systems. It is closely related to finite difference methods but allows for good control of boundary conditions and surface profile assumptions while preserving the conservative nature of the equations. Most importantly, the finite volume formalism accommodates the heterogeneous spatial organization of cellular compartments. [...] Within such elements, the rate of change of the concentration of a given molecular species is simply the sum of fluxes entering the volume element from its adjacent neighbors plus the rate of production of the given species via reactions. [...] Linear solvers based on Krylov space approximations, such as the conjugate gradient method, in conjunction with a preconditioner (an operator that approximates the inverse of the matrix but can be applied at a low computational cost), become powerful and robust. There are commercial packages that implement a range of Krylov space methods, as well as many of the well-known preconditioners (e.g., PCGPAK, Scientific Computing Associates, New Haven, Connecticut).

When can we use deterministic models and when do we have to use stochastic models instead? The authors write:

If the number of molecules involved in a process is relatively small, the fluctuations can become important. In this case, the continuous description is no longer sufficient and stochastic effects have to be included in a model. Single-channel ionic currents are one such example. [...] Stochastic fluctuations of macromolecules are crucial for understanding the dynamics of vesicles and granules driven by competing molecular motors. In the case of a relatively small number of participating particles, a system that would be described deterministically by reaction-diffusion PDEs requires fully stochastic treatment. In this approach, diffusion is described as Brownian random walks of individual particles, and chemical kinetics is simulated as stochastic reaction events. Numerical stochastic simulations in this case are based on pseudo-random-number generation. They are often called Monte Carlo simulations (the term, originally introduced by Ulam and von Neumann in the days of the Manhattan Project) since throwing a dice is actually a way to generate a random number.

They also provide an example of a stochastic model:

As an example, in the Hodgkin-Huxley model, the membrane voltage is treated as a continuous deterministic variable described through a set of differential equations, whereas the single channel behavior is random. A natural way to introduce stochasticity in the model is to replace open probabilities by the actual numbers of open channels. In fact, Hodgkin and Huxley introduced variables in their model to represent the proportion of open gates for various ions. The number of open channels is random and is governed by a corresponding Markov kinetic model that explicitly incorporates the internal workings of the ion channels. Mathematically, the membrane potential is now described by a stochastic differential equation with a discrete random process.

The authors further mention two papers on stochastic methods from Gillespie, which he deems "especially relevant" for Computational Cell Biology (Gillespie (1977): Exact stochastic simulation of coupled chemical reactions ; Gillespie (2001): Approximate accelerated stochastic simulation of chemically reacting systems ). Regarding the pros and cons of the algorithm described in these papers, the authors write:

The extraordinary efficiency of the Gillespie stochastic kinetics algorithm is achieved by restricting the decision process to selecting which reaction will occur and adjusting the time step accordingly. Focusing exclusively on the reaction avoids consideration of the properties of individual reactive species as discrete entities, which minimizes processing time when the number of reacting species is large. However, processing time increases in proportion to the number of different reactions. Furthermore, the Gillespie approach does not easily accommodate the existence of multiple states of different substrates, which may affect their reactivities, and since individual reactive species are not identified as discrete elements, their states, positions, and velocities within the reaction volume cannot be followed over time.
This type of approach has been utilized in the Virtual Cell to combine the deterministic description of a continuously distributed species (RNA) with the stochastic treatment of discrete particles (RNA granules)[.]

What follows is a review of programs used in Computational Neuroscience. The authors mention the programs NEURON and GENESIS, the two of which "use cable theory to treat the dynamics of electrical signals in the complex geometries of neurons", which "solves the equation for membrane potential in a series of connected segments with the overall topology of the neuron". Further, he mentions the model description language NMODL which has been added to NEURON and the interface KINETIKIT which makes GENESIS work with chemical kinetics.

The authors also write about software that is supposed "to build complex biochemical reaction pathways and numerically simulate the time course of the individual molecular species within them", such as GEPASI, Jarnac/Scamp, DBSolve, Berekeley Madonna, ECELL, BioSpice and JSIM.

Then, the authors introduce StochSim:

In this program individual molecules or molecular complexes are represented as discrete software objects or intracellular automata. The time step is set to accommodate the most rapid reaction in the system. [...] When a reaction occurs the system is updated according to the stoichiometry of the reaction. Molecules that exist in more than one state are encoded as “multi-state molecules” using a series of binary flags to represent different states of the molecule such as conformation, ligand binding, or covalent modification. The flags can modify the reactivity of the molecule, and reactions can modify the flags associated with a multi-state molecule.

Compared to the Gillespie algorithm, StochSim is supposed to be faster "in systems where molecules can exist in multiple states".

Next, they write about MCell, "a general Monte Carlo simulator of cellular microphysiology":

MCell utilizes Monte Carlo randomwalk and chemical reaction algorithms using pseudo-randomnumber generation. One of MCell’s convenient features is checkpointing, which involves stopping and restarting a simulation as many times as desired. [...] To speed up simulations, MCell is optimized by using 3D spatial partitioning that makes computing speed virtually independent of microdomain geometric complexity. Running parallel computations, another way to speed up Monte Carlo simulations, is also being pursued in MCell.

The paper mentions "microphysiology of synaptic transmission, [...] statistical chemistry, diffusion theory, single-channel simulation and data analysis, noise analysis, and Markov processes" as possible applications of MCell.

Finally comes the main part of the publication, which is about the authors' own program, Virtual Cell:

Simulations of both nonspatial (i.e., ODEs) and spatial (PDEs) models can be performed. For nonspatial models, compartments are assigned appropriate volume fractions relative to their parents in the model topology and surface-to-volume ratios for the proper treatment of membrane fluxes. In spatial models, the segmented regions within a 1D, 2D, or 3D image are connected to the corresponding compartments in the topology. The geometry is prepared for a model in a separate Geometry workspace and can come from a segmented experimental image or can be defined analytically.
The Virtual Cell software displays spatial and nonspatial simulation solutions for the variables over time. The spatial data viewer displays a single plane section of a 3D data set and can sample the solution along an arbitrary curve (piecewise linear or Bezier spline) or at a set of points. Membranes are displayed as curves superimposed on the volume mesh, and membrane variables are displayed along these curves. The nonspatial data viewer plots any number of variables over time on the same plot.

The authors summarize two studies conducted using Virtual Cell, which are about a "model of Calcium dynamics in a neuronal cell" and "stochastic models for RNA trafficking", and conclude:

The Virtual Cell program has several important advantages for stochastic modeling in eukaryotic cells. First, realistic image-based cell geometries are used to define intracellular reaction volumes, which constrain the stochastic behavior of intracellular reactants in unexpected ways. Second, definitions of reactive species can include multiple states described as either discrete parameters or continuous variables, which provide extraordinary contextual richness and behavioral versatility. Third, dynamic transformation and translocation of multiple individual reactive species can be tracked over time, facilitating integration of spatially heterogeneous stochastic models with simultaneous deterministic reaction/diffusion models. A major future challenge for the Virtual Cell will be to integrate dynamic shape changes in the reaction volume within the powerful and flexible stochastic modeling platform already developed. If this can be accomplished, the holy grail of stochastic modeling of cell motility may be attainable using the Virtual Cell.

In the last chapter of the publication, the authors address future challenges for Computational Cell Biology. Among other things, they write:

To improve stability, accuracy, and overall efficiency of numerical simulations, the issues of reaction stiffness in the PDEs, more accurate representation of irregular boundaries, and choice of effective linear solvers need to be addressed. [...] [A]dditional features are being developed, including modeling membrane potential, stochastic processes, lateral diffusion in membranes, and one-dimensional structures such as microtubules and microfilaments. [...] Also needed are computational tools to treat cell structural dynamics to enable the construction of models of such processes as cell migration or mitosis.