BMC Bioinformatics BioMed Central Methodology article Asynchronous adaptive time step in quantitative cellular automata

Background The behaviors of cells in metazoans are context dependent, thus large-scale multi-cellular modeling is often necessary, for which cellular automata are natural candidates. Two related issues are involved in cellular automata based multi-cellular 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 large-scale, quantitative multi-cellular 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 multi-cellular 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 Hodgkin-Huxley (HH) type membrane equations, is a typical case. Effective reduction of time required is essential to make simulation affordable on available computer resources.
Several methods, including the very popular adaptive time step method, have been proposed to enhance the computational efficiency of large-scale 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 single-cell modeling, faces such challenges as parallelism, asynchronism and heterogeneity when used in multi-cellular 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 multi-cellular situations [5], there are cases where they are not applicable, and propagation must be simulated at a cell-to-cell 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 [6][7][8][9], especially for tissue/organ modeling, not only for its ability to simulate discrete intercellular communication, but also for its implicit large-scale parallelism. In such modeling, a natural mapping between each biological cell and each automaton cell is often assumed. Although most reported applications are quali-tative in nature [9][10][11][12], 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 Quantitative cellular automata with different neighborhoods (radius = 1) 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 two-cell 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.
Since the heterogeneous model does not occupy the whole cell space, the non-cardiac 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 homoge-Simulation results of the heterogeneous 2D model locating in a 256 × 256 cell space neous 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 one-sinoatrial-node-cell version consumed 56 m 20 s to run 400 ms, the four-sinoatrial-node-cell version spent 56 m 58 s. However, without the adaptive time step method, the one-sinoatrial-node-cell 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 asyn- The two-cell simulation results (the second cell) using MatLab Figure 4 The two-cell 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.

The number of integrations in a 1000 ms simulation
chronous 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 fine-grained 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 single-cell 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 fine-grained 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 two-cell 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 5-times 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 Luo-Rudy 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 non-biological 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 Fitzhugh-Nagumo 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 multi-cellular model of host-pathogen interaction [12], even a small-scale 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.

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 n-dimensional 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 floating-point 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 large-scale 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 [15][16][17][18][19]. 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 transjunc-tional 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 one-to-one 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 if-then 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 user-defined 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 ven-tricular 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 single-cell 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.