In Silico Biology 5, 0016 (2004); ©2004, Bioinformation Systems e.V.  
Dagstuhl Seminar "Integrative Bioinformatics"

A life-like virtual cell membrane using discrete automata

Gordon Broderick, Melania Ru'aini, Eugene Chan and Michael J. Ellison*

Project CyberCell™, Institute for Biomolecular Design,
3-67, Medical Sciences Bldg, University of Alberta,
Edmonton, Alberta, Canada, T6G 2H7

* Corresponding author

Edited by E. Wingender; received September 17, 2004; revised and accepted November 17, 2004; published November 27, 2004


A framework is presented that captures the discrete and probabilistic nature of molecular transport and reaction kinetics found in a living cell as well as formally representing the spatial distribution of these phenomena. This particle or agent-based approach is computationally robust and complements established methods. Namely it provides a higher level of spatial resolution than formulations based on ordinary differential equations (ODE) while offering significant advantages in computational efficiency over molecular dynamics (MD). Using this framework, a model cell membrane has been constructed with discrete particle agents that respond to local component interactions that resemble flocking or herding behavioural cues in animals. Results from simulation experiments are presented where this model cell exhibits many of the characteristic behaviours associated with its biological counterpart such as lateral diffusion, response to osmotic pressure gradients, membrane growth and cell division. Lateral diffusion rates and estimates for the membrane modulus of elasticity derived from these simple experiments fall well within a biologically relevant range of values. More importantly, these estimates were obtained by applying a simple qualitative tuning of the model membrane. Membrane growth was simulated by injecting precursor molecules into the proto-cell at different rates and produced a variety of morphologies ranging from a single large cell to a cluster of cells. The computational scalability of this methodology has been tested and results from benchmarking experiments indicate that real-time simulation of a complete bacterial cell will be possible within 10 years.

Keywords: mathematical model, agents, flocking, lipid bilayer, cell membrane, discrete automata, stochastic simulation, virtual cell


The recent emergence of integrative or systems biology on the scientific landscape fills an important gap in the drive towards formulating approaches to understand cells as highly intricate and integrated machines. A logical extension of this initiative is the development of comprehensive or whole-cell models. The construction of in silico models for the simulation of cellular metabolism has traditionally drawn on formalisms and concepts used in electrical and chemical engineering. Metabolic reaction pathways continue to be analyzed using circuit models and large sets of coupled ordinary or partial differential equations. Drawing on a rich background in the classical mathematics, these models rely heavily on the assumption that the concentrations of metabolites vary smoothly and continuously throughout the reaction space. On this basis the cell is generally represented as one or several ideally mixed reactors where deterministic reaction kinetics are applied [1, 2]. These regions can be quite sophisticated in shape but the basic assumptions of a continuum still apply within each compartment [2, 3]. Although these models and their underlying assumptions allow for rapid computation, they are nonetheless quite limited in their ability to describe discrete spatially localized reactions and transport phenomena. Furthermore, irregular and discontinuous spatial distributions of molecules and reaction processes are known to exist in vivo; examples include DNA replication, transcription and translation processes. Such methods are not mathematically robust to conditions found at this level of spatial resolution and will in some cases fail to converge to a solution altogether.

Particle-based models preserve the irregular and discontinuous distribution of matter and can be applied at different levels of resolution. Molecular dynamics (MD) for instance, is a fine grain particle-based approach that is suited to atomic-scale resolution. MD is often used to compute the three-dimensional conformation dynamics of complex molecules. Weak and strong atomic interactions are balanced for individual atoms and new conformations are proposed until the overall system energy is minimized. Modelling at this level of detail is a computationally intense optimization problem. At present MD simulations rarely exceed 10,000 to 20,000 particles in size and can require several days of compute time to simulate less than 5 ns of biomolecular activity [4]. Coarse-grain adaptations of MD, such as dissipative particle dynamics (DPD) [5], can be used to extend the practical temporal scale of these simulations somewhat. Unfortunately the population sizes and timescales which result are still only a fraction of those required to model basic biology of even a single cell.

In parallel efforts, computer scientists have long been interested in applying simple interaction rules to populations of discrete particles as a means to simulate the behaviour of complex dynamical systems with properties that do not adhere to first principles. Simple rules that describe the local interactions between individual particles or automata have been shown to produce highly intricate emergent behaviour at the level of the overall system. Discrete cellular automata positioned on a regular lattice have been used to study problems such as non-equilibrium phase transition [6] and two-phase flow in porous media [7] among others. These classical automata are fixed in space and assume a number of discrete states computed by way of a transition function based on the states of the immediate neighbours. In recent work, 2-dimensional automata using simple state transition rules have been applied to examine diffusion phenomena in systems crowded with slow-moving or immobile obstacles [8]. As discrete automata states are often the norm, studying large state spaces can be an issue even for simple automata. For example, the state space defined by the number of chemical species in a typical biological system can lead to levels of combinatorial complexity which are computationally prohibitive.

Estimates of the molecular assay of even a simple organism such as E. coli suggest that a whole-cell simulation conducted at the level of individual molecules is a sizable undertaking. As an example of scale, a typical E. coli cell is estimated to contain in excess of 60 million biomolecules [9]. This estimate excludes water molecules. It is precisely because of this scale and the computational expense of classical particle-based simulations that whole cell simulations have strongly favoured equation-based methods. Whole-cell simulation at mesoscale resolution has so far received little attention. Nonetheless, hybrid schemes have been used to simulate multi-cellular assemblies by using geometric surfaces to effectively reduce the molecular population [10]. Ideally a single scalable methodology would marry the strengths of particle-based methods with those of equation-based methods to provide both the speed and resolution to tackle the large-scale problem associated with a whole-cell model.

Here we describe such a scalable framework directed at constructing a detailed model of a single cell where each biomolecule is explicitly represented. We have used probabilistic methods describe particle dynamics that realistically mimic actual molecular interactions while requiring only a fraction of the computational power. This framework is applied in the construction of a cell membrane prototype where the fluid movement of the idealized membrane components is maintained. Schools of fish and flocks of birds exhibit highly coordinated patterns of movement based solely on simple local interactions with their nearest neighbours as acquired over a very limited range of perception [11, 12]. Molecules also interact over limited ranges with their immediate neighbours. Of course, this interaction is governed by basic physical laws rather than animal intelligence. Nonetheless, we investigate the possibility of formulating standard Lennard-Jones potentials as a set of behavioural cues or interaction potentials. Combined with molecule-level particle resolution and simplified concise state descriptions, the present model demonstrates how this type of formulation inspired from population biology might be used to provide efficient representations of large scale membranes. Simulation experiments are presented that demonstrate how a proto-cell is bounded by such a membrane, can grow, divide and respond to osmotic pressure gradients.


Automata-based formulation

In the present approach we have designed a framework for simulating nonlinear chemical kinetics which is in many ways conceptually rooted in automata theory. In representing individual molecules we maintain the simplistic state behaviour so crucial to the automata's computational efficiency. Furthermore, in an attempt to minimize memory requirements and accelerate calculations the solvent is not represented explicitly in the model. The resulting lattice space is populated by a large proportion of cells that represent void space. Although sharing the simplicity of stationary automata, our model molecules assume levels and patterns of cooperativity more akin to agent-based systems. In this work, our chemical automata or molecular agents are free to navigate through the simulation space using simple behavioural cues or component interactions. These component interactions are computed on the basis of local environmental conditions and combined into conditional probability distributions from which specific navigational course corrections and chemical state transitions are selected. The state of any given molecular agent is defined by a current position, a velocity vector and a chemical species or state. Relative diffusion rates or velocities for different chemical species are translated into a molecular agent having a greater or lesser probability of not moving at all on any given iteration.

At each iteration, molecular agents examine the local space bounded by a cut-off radius of interaction. Should other molecules be present within a critical contact radius their identity is ascertained and a stoichiometric truth table is examined for chemical affinities and the probability of a candidate reaction occurring. In essence, many of the seminal concepts in stochastic kinetics proposed by Gillespie [13] and applied on a global basis are in this work being applied in a discrete and very local manner. Validation of this framework was first carried out by simulating simple solution chemistry and comparing results to those obtained from classical reaction rate equations. Simulation of the dynamics of basic enzymatic inhibition reaction networks confirmed this framework of discrete molecular agents will indeed express the expected Michaelis-Menten dynamics. Furthermore, as the size of molecular populations is increased the averaging effect becomes clear and the system behaviour approaches that expected from a deterministic set of rate equations for simple reaction networks.

Design of a cell proto-membrane

To overcome some of the computational requirements, a single discrete agent is used to represent a pair of opposing phospholipid molecules which constitute a membrane bilayer element. For similar reasons, solvent is not explicitly modeled in this work. The overall geometry of this model membrane is an emergent property of the interaction between these simple automata or so-called molecular agents. The motion of each membrane particle is driven by three main types of interaction with its neighbouring particles: (1) attraction, (2) dispersion and (3) alignment. This last interaction potential is required because the agents being used in this work are spherical and devoid of polarity. As a result, attraction and repulsion interactions alone are not sufficient to promote the formation of a monolayer. Instead, given a sufficient population of neighbours, each agent computes a small local plane and attempts to move within this local surface tile. This local plane serves as a frame of reference for the interaction potentials. As described in Equation (1), an overall interaction potential is computed for each particle agent using a weighted sum of in-plane and out-of-plane components. Cohesion of the membrane tile results from a balance between attraction and repulsion applied to the particles within this two-dimensional plane. These interactions are modeled after the Lennard-Jones potential and their respective magnitudes are computed based on the projected in-plane distance as shown in Equations (2) and (3). A separate component of interaction, one acting normal to the local plane, is now required to align the molecular agents as a monolayer. The magnitude of this alignment interaction is proportional to the exponentially weighted distance of each agent from its local plane as shown in Equation (4). It is important to note that these planes or tiles are locally generated by each agent and are not subject to any global smoothness constraints. Vector quantities are expressed in bold.


With component interactions described by:

IjNet Resultant vector of interactions applied to particle i
IjIn-plane attr Resultant vector of in-plane attraction on particle i
IjIn-plane rep Resultant vector of in-plane repulsion on particle i
uij Unit vector along "line of sight" between particles i and j
dp ij In-plane distance of particle i to particle j
dp attr Crit Cut-off distance for in-plane attraction
dp rep Crit Cut-off distance for in-plane repulsion
dp ij In-plane distance of particle i to particle j
αrep Exponent for in-plane repulsion
ko Pre-exponent multiplier for in-plane repulsion
N Total number of neighbours for particle i
ωin-plane Weighting coefficient for in-plane attraction-repulsion
ωout-plane Weighting coefficient for out of plane attraction

IiOut of plane Resultant vector for interaction normal to plane
uni Unit vector normal to local plane at i
dni Distance of particle i normal to local plane
dn Crit Cut-off distance normal to plane
αn Exponent for attraction normal to plane

These navigational driving forces are vector quantities and each is computed as the vector sum of the contributions made by individual neighbour particles located within a given radius of interaction. Similarly, the overall resultant or net interaction is also a vector quantity computed from the weighted vector sum of the net attraction, repulsion and alignment components described in Equations (1) through (4). The desirability of individual candidate moves is then assigned in accordance with Equation 5. This value is based on the degree of alignment of each potential move with the unit directional vector corresponding to the overall resultant or net interaction being applied to a particle. Depending on the choice of component weights, this loose assembly of small local tiles can collectively generate very complex geometries.

Dm,iOut of plane Desirability of candidate move m for particle i
um Unit directional vector of candidate move m
uRes i Unit directional vector of resultant interaction IjNet
Denotes the vector dot product

Parallel computing

The development and initial testing of trial concepts was carried out by using MatLab1 as a rapid prototyping environment. Prototype code was then further developed to enable parallel computation by conducting a first set of staging trials in a parallel MatLab environment2 running under a master-slave hierarchical message passing scenario. This staging was carried out on an HP Netserver Lt6000 equipped with 4, 700 MHz processors and 2 GB of RAM, operating under Redhat Linux version 7.1. When the performance of a new component of the parallel code is verified, this code is then translated to C language using the MatLab Compiler toolkit. This translated C code is then linked to a main communication backbone, designed and written in native C, using IBM's mpcc MPI standard compiler. Contrary to the master-slave scenario used under MatLab, this native C backbone is based on a slave-slave communication scenario which effectively reduces communication requirements by a factor of two. Medium-scale experiments are conducted on an shared-memory IBM p690 server equipped with 16 Power 4 processors (1.1 GHz) and 32 GB of RAM, operating under AIX version 5.1. The status of the simulation space is reported in text file format. These data files describing the position and chemical state of each agent are then subjected to post-processing steps to extract relevant statistics and produce graphical representations. Graphical rendering was performed using scripts written in AVS /Express3, both running on a 2-processor IBM Intellistation. Animations of recent simulation results may be viewed at the Project CyberCell website4.


Proto-membrane stability

The first test of this formulation is the validation of its ability to maintain a stable concave geometry if initialized accordingly. Overall stability should not however come at the expense of agent mobility; for example a solution whereby all agents remain at their initial positions is a trivial and unacceptable solution since lateral diffusion of membrane particles must be maintained. To illustrate this property, different levels of internal pressure were created inside a proto-cell by filling the cell with small and highly energetic inert particles. The membrane was constructed assuming a radius of 0.4 nm for the particle agents, each a hard-sphere equivalent of a lipid bilayer component. This size of model particle is consistent with estimates based on values for cell surface area (Goodsell [14]) and lipid population size (Neidhardt and Umbarger [9]) reported for a typical E. coli cell. It should be noted that E. coli is used only as a convenient reference to establish scale and proportionality. The model membrane presented is in fact much simpler in structure than that of an E. coli cell and more closely resembles a mammalian cell membrane.

A proto-cell with a diameter of 60 nm, roughly the size of a typical transport vesicle, was constructed using 13,000 particle agents. This was found to be a convenient problem size for rapid testing. Chemically inert, rapidly diffusing, particles were used in the intra-cellular space as fill material. Three experiments were conducted with a first experiment using an empty intracellular space, a second using a fill level sufficient to balance the surface tension forces, and a third where the internal pressure exceeded the limit strength of the proto-membrane. In these experiments, groups of molecular agents were initially positioned to form the surface of a sphere and then released to pursue their individual trajectories. Each proto-cell is presented graphically at different time steps in Figure 1.

Figure 1: Stability of the overall geometry and surface dynamics of a 60 nm diameter proto-cell where (A) no internal pressure is exerted by fill particles (case I - 9,000 lipid particles and 2,000 diffusional markers), (B) where pressure is lower than elastic limit of the shell (case II - 9,000 lipid particles, 2,000 diffusional markers and 2,000 fill particles) and (C) where pressure exceeds the elastic limit (case III - 9,000 lipid particles, 2,000 diffusional markers and 21,000 fill particles), at time steps i = 1, 2501, and 5001 (1 time step = 0.2 μs). Lipid particles are colored grey, diffusional marker lipids are colored red and fill particles are colored yellow.

Results in Figure 1A show an empty cell initialized as a sphere. In the complete absence of internal pressure the proto-cell shrinks as membrane particles draw closer to one another. This increased level of spatial crowding results in a loss of mobility for individual particle agents. This is demonstrated by the gradual decrease in incremental diffusion path illustrated by curve I in Figure 2B. Nonetheless minimal lateral diffusion is maintained and the proto-cell settles into a new stable configuration. Shrinkage of the proto-cell is effectively arrested in the second experiment as shown in Figure 1B. After a short period of settling-in the membrane appears to reach a "relaxed" state where the agents form a relatively smooth monolayer defining a concave hollow sac. This layer is very fluid and local irregularities such as bulges appear and ripple across the membrane surface. Small transient pores form, alter their shape and disappear dynamically. Furthermore, a higher level of mobility is maintained by each agent with movement being directed primarily, but not exclusively, in the plane of the membrane. Allowing such excursions from the local plane enables the model membrane to evolve morphologically.

Similarly, results in Figure 1C show the expansion and eventual rupture of a model cell where the internal pressure generated by the fill particles exceeds the structural limitations of the membrane. The corresponding increase in surface area is documented by the results in curve III of Figure 2A. The initial rate of expansion is quite high. However, as the proto-cell expands beyond the membrane's limit of deformation, it ruptures and releases its contents into the surrounding media. As fill particles leave the cell and internal pressure is released, the rate of corresponding expansion slows and eventually halts. Results presented in curve III of Figure 2B, show the resulting incremental change in total diffusion path exhibited by the particle agents at the surface of this ruptured cell. Although high at the very outset as a result of rapid initial expansion of the cell, the rate of diffusion rapidly approaches values similar to those observed at the surface of the equilibrated cell. This diffusion rate remains slightly lower nonetheless as the discontinuities in the surface created by the tears offer fewer opportunities for in-membrane lateral diffusion to occur.

Figure 2: (A) Change in total external surface area (x10 nm2) for a simulated proto-cell, 60 nm in diameter, where no fill particles are used (case I - 9,000 lipid particles and 2,000 diffusional markers), where the internal pressure exerted by fill particles is contained by the membrane (case II - 9,000 lipid particles, 2,000 diffusional markers and 2,000 fill particles), and where the internal pressure exceeds the elastic limit of the shell (case III - 9,000 lipid particles, 2,000 diffusional markers and 21,000 fill particles), as it evolves over 5000 time steps (1 time step = 0.2 μs).
(B) Average Incremental change in diffusion path (1 grid unit = 0.1 nm) traveled by the molecular agents in the proto-membrane for cases where internal pressure exerted by fill particles is inexistent (I - 11,000 particles), is contained by the membrane (II - 13,000 particles), and exceeds the elastic limit of the shell (III - 32,000 particles) as it evolves over 5000 time steps (1 time step = 0.2 μs).

The diffusion patterns of membrane particles can be explored further by examining the frequency histograms of diffusion paths presented in Figure 3. In all three experiments, the histograms are Gaussian in shape. As particles follow a pseudo-random path or 'drunken walk' in their tangential diffusion within the membrane, the average length of the total diffusion path for a given experiment is significantly greater than the 'line of sight' net positional change for this same experiment. The effects of crowding on the membrane particles in the collapsing cell and the discontinuities in the surface of the ruptured cell can also be observed in these frequency histograms. Results presented in Figure 3A show that the shape of the distribution in net diffusion path is similar for all three cases leaving the mean value as the main distinguishing feature between experiments. The net path lengths for particles diffusing in a collapsing (curve I) or a ruptured membrane (curve III) are on average higher than those of an equilibrated cell (curve II). The increased bulk motion associated with both these cases will certainly contribute to producing a higher average net path length. This intuitive conclusion is further supported by results presented in Figure 3B where the average total diffusion path length for the collapsing cell as well as for the ruptured cell do not exceed the average total diffusion path traveled by membrane particles in the equilibrated cell. The fact that even with a higher average net path length, a lower or equivalent total displacement is obtained suggests that membrane particles in an equilibrated cell surface travel with a significantly greater degree of random mobility. In these experiments therefore, a higher average net displacement of the model lipid particles can be ascribed mainly to bulk motion of the entire membrane. This same bulk motion leads to more narrow distributions in displacement for both the collapsing and ruptured cell membrane particles since less random motion involved with this component.

Figure 3:: Distribution frequency (in number of particles) of net (A) and total (B) diffusion path length (1 grid unit = 0.1 nm) for autonomous agents forming a proto-membrane for cases where internal pressure exerted by fill particles is inexistent (I - 11,000 particles), is contained by the membrane (II - 13,000 particles), and exceeds the elastic limit of the shell (III - 32,000 particles) as it evolves over 5000 time steps (1 time step = 0.2 μs).

The lateral diffusion velocity of particles forming the membrane falls well within a range of biological relevance. Lipid particles traverse the surface of a bacterial cell, a distance of about 2 μm, in roughly 1 second at 37°C whereas the majority of trans-membrane proteins diffuse at a rate of between 100 to 500 times slower [15]. In general terms therefore, membrane components can be considered to diffuse at velocities in the nm/s to the μm/s order of magnitude range. In the following simulations, a grid increment of 0.1 nm and a simulation time step of 0.2 μs per iteration were used. Applying this scaling to the results obtained from the equilibrated cell simulation presented as curve II in Figure 3A, we obtain an average direct line-of-sight lateral diffusion velocity of approximately 3.27 μm/s. This is within two-fold of experimentally determined values. Furthermore, model parameters can be easily tuned within the framework presented here to more accurately reflect in vivo conditions.

With reasonable scaling adjustments, we have also found that the mechanical properties of the simulated membrane approximate experimentally determined values for pure lipid bilayers. The elastic properties of a lipid vesicle, as expressed by Young's modulus of elasticity, is typically determined by measuring vesicle swelling and the release of luminal solutes as a function of osmotic pressure. An example can be found in work by Hallet et al. [16]. This relative swelling or strain ε is linearly related to the applied stress s during elastic deformation by the modulus of elasticity or Young's modulus k as shown in Equation (6a). Values for Young's modulus k have been reported within the range of 3 to 9 x 107 N/m2. In the case of experiment III, lysis begins at roughly 400 iterations. This corresponds to a relative increase in surface area ΔA/A of about 20% at the point of membrane rupture as shown in Figure 2A. The solute concentration in the luminal volume of the proto-cell was roughly 0.3 mol/l. Combined with a system temperature of roughly 34°C, this solute concentration produces and an osmotic pressure gradient Δp of approximately 790 kPa across the model membrane. This system temperature was estimated based on Maxwell's relation linking particle velocity to temperature as applied to a typical 160 kDa protein solute, moving one grid point (0.1 nm) in 9 out of 10 iterations. Particle limit velocity was also scaled by a factor of 1 x 106, as estimated by Goodsell [14], to adjust between particles traveling in solution as opposed to particles traveling in a vacuum. The model cell in experiment III was designed with a radius r of 30 nm and a wall thickness t of 1.6 nm. Together with the applied osmotic pressure gradient, this results in a value of 7.2 x 107 N/m2 for the elastic modulus of the simulated membrane using Equation (6b), a value that corresponds closely to experimental measurements. Therefore, the settings used to produce a qualitatively life-like behaviour also result in mechanical properties that are similar to those obtained by experiment.

σ = k · ε(6a)
σ normal stress (N/m2)
k modulus of elasticity or Young's modulus (N/m2)
ε strain or relative deformation (m/m)

By rewriting, as in Hallett et al. [16], we obtain

Δp pressure gradient across the vesicle wall (N/m2)
r, t vesicle radius and wall thickness respectively (nm)
k modulus of elasticity or Young's modulus (N/m2)
ΔA change in vesicle surface area (m2)
A initial surface area of the vesicle (m2)

Growth and division

The membrane prototype presented here is a population of freely moving agents navigating under the influence of local incentives. As such, it is possible to have new particles co-opted into a spatially organized group. It should therefore be possible to grow an existing proto-membrane by adding membrane building blocks to the environment. The results of such an experiment are shown in Figures 4A and Figures 4B. Membrane precursor molecules were injected at regular intervals into the intra-cellular space according to a pre-programmed pulse train. Diffusing freely through the cavity or virtual cytoplasm, these precursor particles reacted upon collision with the existing membrane wall changing their chemical state from that of a precursor to that of a membrane particle as prescribed by Equation (7). Once integrated into the monolayer, these particles become subject to the same component interactions as other members of the membrane.

LipidMembrane + LipidPrecursor 2 LipidMembrane(7)
LipidMembrane Lipid bilayer component particle
LipidPrecursor Lipid bilayer precursor
Pf Forward reaction probability
Pr Reverse reaction probability

The rate of growth was controlled in two experiments by altering the delay between pulses and the amount of precursor material injected into the proto-cell at each pulse. Additional control can be exercised by further adjusting the threshold probability for the forward reaction occurring upon collision of precursor particles with the existing membrane. In the first experiment, the precursor injection rate was set proportional to the expected increase in membrane surface area. As demonstrated in Figures 4A, the resulting proto-cell will grow in a stable manner, gradually increasing in volume. The proto-cell retains its structural integrity and its basic hollow geometry. However, at high injection rates an unexpected result was obtained. In the experiment shown in Figures 4B, the proto-cell remains roughly the same size and the rapid injection of precursors leads to the formation of a network of walls in the luminal space. These walls span the diameter of the cavity and divide the proto-cell into separate chambers. It is emphasized that this behaviour was not explicitly programmed into the simulation but rather arose from the combination of local conditions and the interaction of membrane particles navigating under local behavioural rules. Because these simulations are based on stochastic methods, no two simulations evolve in exactly the same way. Notably, the size, shape and number of chambers will differ in experiments that have identical start points. The final structure is therefore ultimately determined by small fluctuations in local environment from experiment to experiment. This can be compared to bifurcation points occurring in high order systems that lead to multiple separate and stable system states.

Figure 4: Growth dynamics of a proto-cell, initially 60 nm in diameter, fuelled by the introduction of membrane building blocks in the intracellular space at a rate proportional to the increase in surface area (A) and at a significantly higher rate (B) at time steps i = 1, 14001 and 24001 respectively (1 time step = 0.2 μs). Lipid particles are coloured grey, new membrane particles are coloured red and fill particles are coloured yellow. At final state, the slow-growth proto-cell (A) was composed of 11,000 lipid particles, 3,500 new membrane particles and 4,000 inert fill particles, whereas final composition of the rapid-growth proto-cell (B) was 11,000 lipid particles, 9,000 new membrane particles and 3,000 inert fill particles.


One of the defining attributes of a living cell is its ability to carry out biochemical reactions in isolation from its surrounding environment, where conditions can be controlled locally to promote certain reactions and inhibit others. The membrane enables the cell to create and maintain a closed environment as well as providing a selective and directional barrier to molecules which can be handpicked for entry into the cell or transport out of the cell. Clearly the design and implementation of a biologically realistic model is a key element in any whole cell simulation. In current whole cell models the membrane is typically represented by a static surface. However membranes are dynamic structures which change in composition, shape and surface area as they adapt to changing conditions.

Much of the work directed at understanding membrane mechanics has involved classical as well as coarse grain molecular dynamics simulations conducted at the atom or atom cluster scale. For example, coarse grain MD and DPD have been used to study the self-assembly of bilayers [17] as well as the formation of transient pores in the membrane [18, 19]. Because of their high computational overhead these atomic-scale studies have been limited to the simulation of relatively small slices of lipid bilayer but have done so in great mechanistic detail. Recent work has also been directed at examining phenomena that occur over larger bilayer surfaces. Mesoscale undulations of cell-scale lipid bilayers have been studied by recruiting methods from the field of statistical mechanics and engineering. These methods however require mean-field continuum assumptions to provide the necessary computing speed. Recent extensions to molecular dynamics methods, such as those presented in work by Ayton et al. [20], have also been used to model membrane motion from basic mechanical properties. Although not a replacement for these rigorous mechanistic models, the work presented here offers an approximate solution which is strikingly similar to the microscopic and macroscopic properties of a biological membrane. A simple system composed of relatively unsophisticated particle agents therefore is capable of manifesting characteristic membrane properties such as lateral diffusion, response to osmotic pressure gradients, growth and division. Furthermore, simple qualitative tuning of the model can quantitatively align these properties surprisingly well with those of its biological counterpart.

Perhaps the main advantage of this framework is its flexibility and its predisposition to interface with other methods. Within a hybrid model scheme this discrete formalism can be paired with differential equation-based methods to add spatial resolution and numerical stability in situations where continuum assumptions do not hold true. For example metabolites and other small molecules which are abundant and adhere to a smooth concentration gradient could be managed by ordinary or partial differential equations. Whereas the larger, less abundant and more spatially disparate molecules which make up much of the proteome would be managed in parallel by the discrete particle formalism presented here. Both discrete particle and continuous field solutions could then be superimposed at each iteration to provide a complete picture of the system. Such a scheme would make maximum use of the numerical stability of a particle-based method and the speed of resolution associated with equation-based methods. Integration within a hybrid scheme is not however a requirement for the effective use of this formalism which scales well with problem size and available resources. The results of early benchmarking trials carried out with a prototype level code are presented in Figure 5. These results indicate that a close to ideal speed-up factor is maintained as the problem size is increased and additional computer resources are added. This suggests that with the communication and parallelization strategy currently implemented, additional processors are used quite effectively by the simulation code in enabling the rapid execution of larger and more complex systems.

Figure 5: Speedup factor achieved during benchmarking experiments show high usage rate of parallel resources and linear scale-up in performance with increases in the size of the processor pool.

As this method implicitly represents individual molecules and molecular assemblies, we can attempt to quantify the computational requirements associated with a complete cell. Even though this simple membrane model is more closely related to that of a mammalian cell in structure, we will use readily available measurements of the size and molecular composition of E. coli as a basis for estimating the complexity of a simple organism. In recent work, Goodsell [14] estimates that in a state of balanced growth, conducted in a glucose medium at 37°C, a typical E. coli membrane would encapsulate a volume of 0.88μm3. A spherical cell with an equivalent volume would have a radius of approximately 600nm and an external surface area of 4.5 μm2. Neidhardt and Umbarger [9] report that a typical E. coli cell will contain in the order of 22 million lipid molecules. Since most of these will make up the cell membrane, an idealized hard-sphere membrane particle would therefore have a radius in the order of 0.4 - 0.6 nm, which corresponds well with our estimates. Data from the same authors suggests that this typical cell will contain roughly 2.5 million protein molecules. Assuming that ions, cofactors and metabolites constitute 4.3% of dry cell weight [9], and that they are contained in a net cytoplasmic volume of approximately 0.6 μm3 [14], we obtain an additional 36 million biomolecules. This produces a model E. coli cell with a total population of approximately 50 million biomolecules at full scale. In comparison, Mycoplasma is roughly 1/50th the size of an E. coli cell in terms of molecular population. This spherical cell with a radius of 150 nm offers an attractive intermediate problem size for constructing a basic minimal cell. Using concentrations which are consistent with estimates for E. coli, a mycoplasma cell would include approximately 50,000 protein molecules, 500,000 metabolites, cofactors and ions, and be bounded by a membrane composed of roughly 600,000 lipid molecules, representing a total of approximately 1 million biomolecules.

The computational power required to model a whole cell can be estimated from scale-up experiments performed by our group on a variety of platforms. The simulation code is being designed to standards which make it fully compatible with IBM's BlueGene L supercomputing platform. This platform will incorporate roughly 65,000 processors. Our most current results indicate that 24 hours of computing time on BlueGene L would enable us to capture 60 seconds of mycoplasma life with a spatial resolution of 1 nm and time step of 1 μs. Accordingly, a complete cell cycle for mycoplasma could be simulated in roughly 3 weeks. Using Moore's law to loosely approximate the time horizon for achieving real-time performance, we expect it will be possible to execute real-time simulations representing every molecule of an E. coli. cell on a typical high-performance computing platforms in less than 10 years.


First and foremost the authors wish to thank the Alberta Science and Research Authority (ASRA) and the Canadian Foundation for Innovation (CFI) for their continued and generous support of Project CyberCell. The authors also extend special thanks to IBM for their significant contributions both through consultation with their technical experts and through the generous donation of high-performance computing hardware. Thank you to Dr. Peter Tieleman for his helpful comments regarding the performance of state of the art molecular dynamics software. Acknowledged also is the valued contribution of Brendan Cote to our code optimization efforts, Torin Huzel for his assistance with the graphical rendering, Einar Brandt for making available the parallel MatLab environment, as well as Ron Senda and his staff at the University of Alberta CNS Research Group.


  1. Tomita, M., Hashimoto, K., Takahashi, K., Shimuzu, T.S., Matsuzaki, Y., Miyoshi, F., Saito, K., Tanida, S., Yugi, K., Venter, J. C. and Hutchison, C. A. (1999). E-Cell: Software environment for whole cell simulation. Bioinformatics 15, 72-84.

  2. Reed, J. L. and Palsson, B. O. (2003). Minireview: Thirteen years of building constraint-based in silico models of Escherichia coli. J. Bacteriol. 185, 2692-2699.

  3. Schaff, J. C., Slepchenko, B. M., Choi, Y., Wagner, J., Resasco, D. and Loew, L. (2001). Analysis of nonlinear dynamics on arbitrary geometries with the Virtual Cell. Chaos 11, 115-131.

  4. Pastor, R. W. (1994). Molecular dynamics and Monte Carlo simulations of lipid bilayers. Curr. Opin. Struct. Biol. 4, 486-492.

  5. Groot, R. D. and Warren, P. B. (1997). Dissipative particle dynamics: bridging the gap between atomistic and mesoscopic simulation. J. Chem. Phys. 107, 4423-4435.

  6. Chopard, B. and Droz, M. (1988). Cellular automata approach to non-equilibrium phase transitions in a surface reaction model: static and dynamic properties. J. Phys. A: Math. Gen. 21, 205-211.

  7. Krafczyk, M., Schulz, M. and Rank, E. (1998). Lattice gas simulations of two-phase flow in porous media. Commun. Numer. Meth. Engng. 14, 709-717.

  8. Berry, H. (2002). Monte Carlo Simulations of Enzyme Reactions in Two Dimensions: Fractal Kinetics and Spatial Segregation. Biophys. J. 83, 1891-1901.

  9. Neidhart, F. C. and Umbarger, H. E. (1996). Chemical Composition of Escherichia coli. In: Escherichia coli and Salmonella. Cellular and Molecular Biology. Curtiss, R., Ingraham, J. L., Lin, E. C. C., Low, K. B., Magasanik, B., Rfznikopp, W. S., Riley, M., Schaechter, M. and Umbarger, H. E. (eds.), Vol. 1, 2nd edn., American Society for Microbiology Press, Washington, DC, pp. 13-16.

  10. Franks, K. M., Bartol, T. M. and Sejnowski, T. J. (2002). A Monte Carlo model reveals independent signal at central glutamatergic synapses. Biophys. J. 83, 2333-2348.

  11. Viscido, S. V., Miller, M. and Wethey, D. S. (2001). The response of a selfish herd to an attack from outside the group perimeter. J. theor. Biol. 208, 315-328.

  12. Inada, Y. and Kawachi, K. (2002). Order and flexibility in the motion of fish schools. J. theor. Biol. 214, 371-387.

  13. Gillespie, D. T. (1977). Exact stochastic simulation of coupled chemical reactions. J. Phys. Chem. 81, 2340-2361.

  14. Goodsell, D. S. (1991). Inside a living cell. Trends Biochem. Sci. 16, 203-206.

  15. Horton, H. R., Moran, L. A., Ochs, R. S., Rawn, J. D. and Scrimgeour, K. G. (2002). Lipid Bilayers and Membranes are Dynamic structures. In: Principles of Biochemistry, 3rd edn., Challice J. (ed.), Prentice Hall, Upper Saddle River, NJ, pp. 281-285.

  16. Hallett, F. R., Marsh, J., Nickel, B. G. and Wood, J. M. (1993). Mechanical properties of vesicles II: a model for osmotic swelling and lysis. Biophys. J. 64, 435-442.

  17. Goetz, R. and Lipowsky, R. (1998). Computer simulations of bilayer membranes: Self-assembly and interfacial tension. J. Chem. Phys. 108, 7397-7409.

  18. Müller, M. Kastov, K. and Schick, M. (2002). New mechanism of membrane fusion. J. Chem. Phys. 116, 2342-2345.

  19. Groot, R. D. and Rabone, K. L. (2001). Mesoscopic simulation of cell membrane damage, morphology change and rupture by nonionic surfactants. Biophys. J. 81, 725-736.

  20. Ayton, G., Smondyrev, A. M., Bardenhagen, S. G., McMurtry, P. and Voth, G. (2002). Interfacing molecular dynamics and macro-scale simulations for lipid bilayer vesicles. Biophys. J. 83, 1026-1038.


Footnote 1: MatLab is a registered product of The MathWorks, Inc. of Natick, MA (see

Footnote 2: MatLab Parallelisation Toolkit is made available by Einar Brandt Heiberg, Research Engineer at the Department of Medicine and Care, Clinical Physiology, Linkoping University, Linkoping, Sweden (

Footnote 3: AVS/Express is a registered product of Advanced Visual Systems Inc. of Waltham, MA. (see

Footnote 4: