Abstract
Background
The behaviors of cells in metazoans are context dependent, thus largescale multicellular modeling is often necessary, for which cellular automata are natural candidates. Two related issues are involved in cellular automata based multicellular modeling: how to introduce differential equation based quantitative computing to precisely describe cellular activity, and upon it, how to solve the heavy time consumption issue in simulation.
Results
Based on a modified, language based cellular automata system we extended that allows ordinary differential equations in models, we introduce a method implementing asynchronous adaptive time step in simulation that can considerably improve efficiency yet without a significant sacrifice of accuracy. An average speedup rate of 4–5 is achieved in the given example.
Conclusions
Strategies for reducing time consumption in simulation are indispensable for largescale, quantitative multicellular models, because even a small 100 × 100 × 100 tissue slab contains one million cells. Distributed and adaptive time step is a practical solution in cellular automata environment.
Background
Quantitative multicellular modeling, used, for example, in the electrophysiological modeling of the heart, in neural networks, and in the metabolic modeling of cells in tissues/organs [1,2], is an important area of computational biology [3]. Computational efficiency is an important issue for such in silico research, because, first, many ordinary differential equations (ODE) employed to describe intracellular activity are nonlinear; and second, such models often consist of a large number of homo or heterogeneous cells (Fig 1). Electrophysiological modeling of the heart, in which the description of electrical activity in each of a large number of heterogeneous cells is based on HodgkinHuxley (HH) type membrane equations, is a typical case. Effective reduction of time required is essential to make simulation affordable on available computer resources.
Figure 1. Quantitative cellular automata with different neighborhoods (radius = 1). (A) Moor neighborhood. (B) A userdefined neighborhood.
Figure 2. Simulation results of the heterogeneous 2D model locating in a 256 × 256 cell space. Different colors stand for different values of the dynamically computed membrane potential of cells. 36675 and 36606 are running steps; each step corresponds to 0.01 ms (mini second). The color strips between the inner surface and the outer surface show the normal endocardiumtoepicardium repolarization. At 366 ms most cells have undergone depolarization and are at the end of repolarization, the high accordance between picture A and picture B therefore indicates the validity and safety of the adaptive time step method. A. The distribution of membrane potential of cells in epicardiumtoendocardium repolarization without the adaptive time step. B. The distribution of membrane potential of cells with the adaptive time step. Pictures are captured with GIMP, a screen capture tool under Linux.
Several methods, including the very popular adaptive time step method, have been proposed to enhance the computational efficiency of largescale systems. Obviously, the effectiveness of a method is chiefly dependent on the details of implementation as well as the scale and nature of models. For example, the classic Rush & Larson method for runtime adjustment of time step, though very straightforward for singlecell modeling, faces such challenges as parallelism, asynchronism and heterogeneity when used in multicellular biological models [4]. Very often, a method coded and workable in one programming environment cannot be directly adopted in another. Thus, new methods, accompanied by new modeling strategies and tools, are reported from time to time.
While partial differential equations (PDE) are often used to model propagation phenomena in multicellular situations [5], there are cases where they are not applicable, and propagation must be simulated at a celltocell level. This occurs, for example, when a system is heterogeneous consisting of different kinds of components. When a system evolves dynamically and shows emergent behaviors, centralized, PDE based descriptions also become unfitting. Modeling arrhythmias with a heart model is such a case in which the electrical properties of assorted cardiac cells and gap junctions may change substantially in simulation. For such systems, decentralized modeling is needed, which calls for relevant optimal numerical methods to conquer the time consumption problem.
Modeling with cellular automata has gained much prevalence recently [69], especially for tissue/organ modeling, not only for its ability to simulate discrete intercellular communication, but also for its implicit largescale parallelism. In such modeling, a natural mapping between each biological cell and each automaton cell is often assumed. Although most reported applications are qualitative in nature [912], cellular automata are not necessarily restricted to discrete modeling. Various extensions can be made to existing systems to facilitate various types of cellular automata computation [13]. Effective methods to improve the efficiency of quantitative cellular automata computation are therefore of much interest. The purpose of this short paper is to introduce a novel and widely applicable method to implement distributed and asynchronous adaptive time step in quantitative cellular automata modeling. A heterogeneous electrophysiological model of the heart built with five groups of membrane equations of cardiac cells is described as the test model, and performance evaluation based on it is made. As in Rush & Larson method, cyclic activity of cells, either synchronous or asynchronous, is the basis of the distributed asynchronous adaptive time step method.
Results
The test model is a two dimensional one consisting of six kinds, 4300 cardiac cells whose activity is described by relevant HH type membrane equations (see Fig 2). Electrical activity is cyclically created in the sinoatrial node cells by the membrane equations in these cells. Then, the potential difference between repolarized cells and resting cells triggers the membrane equations in resting cells, driving the propagation of cardioelectrical activity in the heart. The long repolarization and resting stages in each cell in each cycle makes adaptive time step method useful. Simulations were carried out on a PC with a 1.7 GHz Pentium4 CPU. First, we used MatLab to perform simulation with two ventricular cells connected by a gap junction (the first cell is stimulated by an external current; the second one is triggered by the transjunctional current) to find out appropriate thresholds for Ina and factor. A workable (but possibly not optimal) set of Ina value for ventricular cells is 0.01, 0.001, 0.002 (mA), and the corresponding values of factor are 1, 2, 5, and 7. The second cell executes 10071 integrations in a 600 ms simulation. A speedup rate of 5.96 ( = 60000/10071) is acquired.
To confirm the validity and stability of this method in a realistic heterogeneous model, we carried out simulation with the sparse 2D test model. It shows that the method worked reliably and produced an action potential distribution at repolarization stage that agrees quite well with that produced by simulation without adaptive time step (Fig 2). At the time of 100,000 steps (1000 ms), we checked five ventricular cells to see how many skips occurred in each. A group of extra fields, num1 to num7, were set in the Cellang program to record the times of integration under different skips (see Fig 3). For these five cells, an average speedup rate of 4.456 was reached. As in the twocell simulation, small fluctuations of transmembrane current and gating variable m occur in the phase 4 of action potential caused by larger time steps, but they did not affect action potential and transjunctional propagation.
Figure 3. The number of integrations in a 1000 ms simulation. Five cells locating at <93,75>, <18,38>, <96,42>, <54,37>, <54,72> are chosen to display the number of integration happened in each cell in a simulation with the adaptive time step. numX gives the number of integrations; the real time step is Δt*X.
Figure 4. The twocell simulation results (the second cell) using MatLab. A. The action potential. B. The state of gating variable m that controls the sodium channel. In phase 4, m fluctuates within a range, but the action potential keeps correct.
Since the heterogeneous model does not occupy the whole cell space, the noncardiac automaton cells should also spend a certain amount of time on housekeeping operations. To further check the efficiency of the method, we made 400 ms simulations with a 128 × 128 homogeneous 2D model fully occupied by ventricular cells except 1 or 4 sinoatrial node cells locating at 1 or 4 corners to initiate excitation. With the adaptive time step method, the onesinoatrialnodecell version consumed 56 m 20 s to run 400 ms, the foursinoatrialnodecell version spent 56 m 58 s. However, without the adaptive time step method, the onesinoatrialnodecell version needs 142 m 56 s to make a 400 ms simulation. We remark that the fact that the agreement of the speedup rate for both propagation stimulated by one corner and propagation stimulated by four corners is significant for the simulation of abnormal (arrhythmic) propagation, as it indicates that highly asynchronous cellular activities, such as ectopic beats and reentrant movement, do not degrade the performance as they do in some methods [20].
Discussion
We described a novel method to implement distributed and asynchronous adaptive time step in a quantitative cellular automaton environment and showed that this method is numerically efficient. The variation of cellular activity that admits different scale of time step of integration is the precondition to employ the method. Simulation results acquired from two electrophysiological models prove the effectiveness of this method, in which a speedup rate of about 5 can be reached. Since each cell determines its time step independently based on the value of Ina, the dynamic adjustment of time step of all cells runs in a highly asynchronous way. Simulation results of the heterogeneous and sparse 2D model show that the method can work effectively in realistic models. Also, as we have indicated above, the method is effective in the presence of highly asynchronous cellular activities. We believe this is because the method itself works in a fully asynchronous way and a finegrained and implicit domain decomposition exists naturally (each cell is a subdomain). When the global Δt is small enough and the threshold values of factor and Ina, the indicator of cellular activity, are set within suitable bounds determined by singlecell simulations, the method is safe and stable. Simulation shows that the computational gain always surpasses the overhead of using two extra fields to check the value of Ina and to keep tract of skip.
In comparison with Quan et al's method [20], the advantage of this method is that there is neither an explicit domain decomposition nor a priority queue. Actually, a finegrained domain decomposition is implicit in cellular automata computation as each automaton cell is both a computing unit and a subdomain. Thus, there is no extra cost of domain decomposition. Besides, in a complex and highly heterogeneous model, an algorithm for explicit domain decomposition may be very complicated and dependent on many details of the model such as dimension and geometrical properties. In contrast, the model based on quantitative cellular automata described here is completely general.
We remark that Quan et al obtained a reduction of computation as high as 17 times (from 3 to 17), with the highest speedup rate attained in the resting phase (phase 4) of the action potential [20]. This is, however, not supported by our results. Both the twocell straightforward simulation and the sparse heterogeneous simulation show that the dynamics of the gating variable m is unstable in phase 4, especially when time step is large (Fig 4). By contrast, we obtained the highest stable speedup rate in the plateau phase. We found the overall 5times speedup rate we obtained is the same as that reported by Cherry et al [21], who used an explicit mesh refinement method. How large a speedup rate can be obtained and when the largest speedup rate is attained depend on the dynamic properties of the ODEs involved. Our data show that the activity of ventricular cells described with the LuoRudy phase I model is more stable in phase 2 than in phase 4, in which a 7times speedup rate cannot be sustained.
In a specific model built with quantitative cellular automata, how much time can be saved depends on several factors. The first is the structure of a model. The sparser a model is, like the heart with four chambers and an irregular surface, the less we benefit from the method, because of the housekeeping operations spent on nonbiological cells. The second factor is the time spent on cell communication when istim is also used to jointly control the skip. In a 3D model with the Moore neighborhood, each cell has 26 neighboring cells. Thus, the time spent on computing istim will be tripled compared to a 2D model in which a cell has only 8 neighbors. In electrophysiological modeling, when a complex algorithm is used to compute transjunctional current [22], the time spent on cell communication can be considerable. Thirdly, the complexity of differential equations has a significant impact on the performance. If newly published action potential models, which use many more ODEs to describe more channels and the concentration fluctuation of ions, are used, the suspension of integration can significantly improve the performance of the simulation. If, on the other hand, the action potential model is very simple, for example, of the FitzhughNagumo type, the extra cost arising from the use of two additional fields might neutralize the benefit of adaptive time step. Finally, the proper definition of thresholds of Ina and factor is crucial.
Conclusion
Modeling with cellular automaton has become very popular recently, especially in the modeling of tissues/organs in computational biology. Our work demonstrates that quantitative modeling, which requires an extension of classical cellular automata, is feasible. Furthermore, the simple yet effective method described in this paper shows that distributed and asynchronous adaptive time step can be readily implemented and the performance of simulation can be significantly improved without a significant loss of accuracy of simulation. Besides electrophysiological modeling, this method can be applied to other biological modeling. For example, in cell cycle control and signal transduction models of tissue growth and morphogenesis [23,24], or in multicellular model of hostpathogen interaction [12], even a smallscale 100 × 100 × 100 3D model contains 1 million cells, whose complex, local activities described with ODEs will lead to much more complex, global behavior through cell to cell communication. The adaptive time step method in quantitative cellular automata environment described here may render the simulation practicable.
Methods
Quantitative Cellular Automata
Celluar is a cellular automata system with the programming language Cellang [14]. A Cellang program is shared by all automaton cells and describes the computation within and between cells. A data file, which can be edited manually or created by a program, specifies the location of cells in an ndimensional cell array and the initial value of cell field(s) in each cell. The data file functions as input to the Cellang program. A predefined variable time, which increases by 1 after each step, provides synchronization for all cells. Though common arithmetical and logical operators are provided, Cellang lacks the necessary facilities for numerical computation.
Discrete in value should not be compulsory for modeling with cellular automata in many cases. Exploiting the fact that a Cellang program uses C files as intermediate codes, we built numerical computation facilities into it by adding floatingpoint data type, function calls (including the mathematical functions in the C library), and other quantitative facilities. With these extensions, numerical solutions of ODE can be coded and largescale parallel solution of ODEs (over a large number of cells) can be realized in a simple and straightforward way. We have built several models with the extended system, including an electrophysiological model of the whole heart for which the proposed adaptive time step method is very important.
Models of the Heart
To validate the efficiency of the asynchronous adaptive time step method, two prototypical models were designed and tested. The first one is a 128 × 128 2D homogeneous ventricular sheet, with cells of sinoatrial node located at one or four corners to initiate excitation. The second one, also with a 128 × 128 resolution, is a 2D heterogeneous cardiac sheet consisting of ventricular and atrial cells, cells of atrioventricular and sinoatrial nodes, and those of conduction fibers. The electrical activities in these cells are described by the corresponding HH type action potential models [1519]. We note that a 128 × 128 2D model involves 16384 action potential models, i.e., 16384 systems of ODEs. Numerical solution of ODEs is carried out using the explicit Euler method. In both models, a standard Moore neighborhood is adopted, i.e., each cell has 8 neighbors. A simple and static gap junction model, which uses the potential difference between two cells and the resistance of the gap junction to determine the transjunctional current, is adopted. The resistance of gap junctions does not change with membrane potentials.
Among the various quantitative cellular automata models in which intracellular activity is described by ODE, the heterogeneous electrophysiological model may be the best one to illustrate the parallel solution of nonlinear (HH type) equations and the implementation of asynchronous adaptive time steps. In this model, the activity of different kinds of cardiac cells, which we assume to be in onetoone correspondence with automaton cells, is rendered as follows: A field type is defined in the Cellang program, with initial value set in the data file. According to the different value of type, the Cellang program is divided into several sections with the ifthen statement. Different cells execute codes in different sections. A set of cell fields are declared to store the current value of the membrane potential, transjunctional currents and gating variables. Windows are created to display a selected field, such as the transmembrane potential, of all cells (Fig 2), as well as various cellular and channel electrical activities of any selected cells (Fig 3).
Asynchronous Adaptive Time Step
Although there are different methods to implement adaptive time step in biological models built with ODEs, they actually share the same biological basis, i.e., the fluctuation of cell activity such as cell cycle activity and electrical activity in which the change rate of key variables is different in different times. What we did is to implement the method in a quantitative and heterogeneous cellular automata environment. To use the quantitative cellular automata to simulate a system with ODEs in a large number of cells, besides the predefined time that is used to iterate integration in cells, a userdefined floating type time step Δt is defined, which is automatically shared by all cells, for the numerical solution of ODEs. The value of Δt depends on the nature of ODEs. In our two illustrative models it is selected to be 0.01 ms. Δt is changeable in runtime. However, to describe asynchronous cellular activities caused and connected by normal and abnormal cell communication (such as excitation propagation in this case), it is not feasible to dynamically change Δt. For a fixed Δt, to simulate 1000 ms of electrical activity, every cell in the model should run 100,000 times integration. Yet, there is another way to make usage of the cyclic property of cellular activity to reduce computation in simulation. In the case of the electrophysiological model described here, the integration of membrane equation can be skipped when the cellular activity is not very active.
As in other methods [20], for the modeling of ventricular cells we can use Ina, the ionic current of sodium, to reflect the activity of each ventricular cell. We use an extra cell field to store Ina. At any step of the computation, if a ventricular cell's Ina value is smaller than a threshold, the integration of all membrane equations at this step can be skipped. Based on results acquired from singlecell simulation, we can specify a set of Ina thresholds to determine how many steps of integration can be safely skipped in different situations. In our examples, four Ina thresholds are set to determine whether integration of the HH equations should be executed immediately, or skipped 1, 4, or 6 times respectively (specified by a temporary variable factor whose value is 1, 2, 5, 7). At the same time, we also record how many steps have been skipped since the last integration with a cell field skip. Thus, the real time step used in each integration for each cell is Δt*skip.
We note, however, that a static skip mechanism (a specific value of Ina determines a specific value of factor) cannot accurately deal with emergent depolarization that can be triggered by either normal or abnormal transjunctional currents. The latter situation frequently occurs in various arrhythmias. The Ina determined, unmodifiable factor makes the final Δt*skip often too large, leading to errors or even an overflow. Hence, we make the skip dynamically adjustable as follows:
1. In each step, initially set the default value of factor to 1, meaning to do integration immediately;
2. Read the field Ina, according to its value and the preset thresholds of Ina, determine the value of factor;
3. If factor <= skip, do integration and reset skip to 1;
4. Otherwise, skip integration and increase the skip number by skip = skip + 1.
In this way, for every cell, the real time step used in each integration is dynamically adjusted at runtime according to the strength of its activity.
For cells where Ina is not the major inward depolarizing current, dvdt, the change rate of membrane potential, can be the indicator of cellular activity. Although the algorithm can be put before computing the transjunctional currents to skip even transjunctional currents computation, a safer choice is to put it after the computation and use istim, the sum of transjunctional currents that functions as the stimulating current to a cell, to jointly control the suspension of integration.
Authors' contributions
HZ invented the method and drafted the paper. PP provided critical comments on mathematical issues and polished the paper. YS was involved in implementing the method. PD was participated in and in charge of the project. All authors prepared, read and approved the final manuscript.
Acknowledgements
This work was supported by ASTAR, the Agency of Science, Technology And Research, of Singapore.
References

Thakor NV Jr, Ferrero JM, Saiz J, Gramatikov BI, Ferrero JM Sr: Electrophysiologic models of heart cells and cell networks.
IEEE Eng Med Biol Mag 1998, 17(5):7383. PubMed Abstract  Publisher Full Text

Segre D, Vitkup D, Church GM: Analysis of optimality in natural and perturbed metabolic networks.
Proc Natl Acad Sci 2002, 99:1511215117. PubMed Abstract  Publisher Full Text

Rush S, Larson H: A practical algorithm for solving dynamic membrane equations.
IEEE Trans Biomed Eng 1978, 25:389392. PubMed Abstract

Qu Z, Garfinkel A: An advanced algorithm for solving partial differential equation in cardiac conduction.
IEEE Trans Biomed Eng 1999, 46:11661168. PubMed Abstract  Publisher Full Text

Ermentrout GB, EdelsteinKeshet L: Cellular automata approach to biological modeling.
J Theor Biol 1993, 160:97133. PubMed Abstract  Publisher Full Text

Siregar P, Sinteff JP, Julen N, Le Beux P: An interactive 3D anisotropic cellular automata model of the heart.
Comput Biomed Res 1998, 31:323347. PubMed Abstract  Publisher Full Text

Patel AA, Gawlinski ET, Lemieux KL, Gatenby RA: A cellular automata model of early tumor growth and invasion: the effects of native tissue vascularity and increased anaerobic tumor metabolism.
J Theor Biol 2001, 213:315331. PubMed Abstract  Publisher Full Text

Benyoussef A, El HafidAllah N, El Kenz A, EzZahraouy H, Loulidi M: Dynamics of HIV infection on 2D cellular automata.
Physical A 2003, 322:506520. Publisher Full Text

Griffeath D, Moore C: New Constructions in Cellular Automata. Oxford University Press; 2003.

Eckart JD: A cellular automata simulation system: version 2.0.

Yanagihara K, Noma A, Irisawa H: Reconstruction of sinoatrial node pacemaker potential based on the voltage clamp experiments.

Liu Y, Zeng W, Delmar M, Jalife J: Ionic mechanisms of electronic inhibition and concealed conduction in rabbit atrioventricular nodal myocytes.
Circulation 1993, 88:16341646. PubMed Abstract

Nygren A, Fiset C, Firek L, Clark JW, Lindblad DS, Clark RB, Giles WR: Mathematical model of an adult human atrial Cell: the role of K+ currents in repolarization.
Circ Res 1998, 82:6381. PubMed Abstract  Publisher Full Text

Luo CH, Rudy Y: A model of the ventricular cardiac action potential, depolarization, repolarization, and their interaction.
Circ Res 1991, 68:15011526. PubMed Abstract

McAllister RE, Noble D, Tsien RW: Reconstruction of the electrical activity of cardiac purkinje fibres.
J Physiol 1975, 251:159. PubMed Abstract

Quan W, Evans SJ, Hastings HM: Efficient integration of a realistic twodimensional cardiac tissue model by domain decomposition.
IEEE Trans Biomed Eng 1998, 45:372384. PubMed Abstract  Publisher Full Text

Cherry EM, Greenside HS, Henriquez CS: A spacetime adaptive method for simulating complex cardiac dynamics.
Phys Rev Lett 2000, 84:13431346. PubMed Abstract  Publisher Full Text

Vogel R, Weigart R: Mathematical model of vertebrate gap junctions derived from electrical measurements on homotypic and heterotypic channels.
J Physiol 1998, 510:177189. PubMed Abstract

Chen KC, CsikaszNagy A, Gyorffy B, Val J, Novak B, Tyson JJ: Kinetic analysis of a molecular model of the budding yeast cell cycle.
Mol Biol Cell 2000, 11:369391. PubMed Abstract  Publisher Full Text

Schoeberl B, EichlerJonsson C, Gilles ED, Muller G: Computational modeling of the dynamics of the MAP kinase cascade activated by surface and internalized EGF receptors.
Nat Biotechnol 2002, 20:370375. PubMed Abstract  Publisher Full Text