Cellular automata simulation of topological effects on the dynamics of feedforward motifs
 Advait A Apte^{1},
 John W Cain^{2, 3},
 Danail G Bonchev^{2, 3} and
 Stephen S Fong^{1, 3}Email author
DOI: 10.1186/1754161122
© Apte et al; licensee BioMed Central Ltd. 2008
Received: 01 October 2007
Accepted: 27 February 2008
Published: 27 February 2008
Abstract
Background
Feedforward motifs are important functional modules in biological and other complex networks. The functionality of feedforward motifs and other network motifs is largely dictated by the connectivity of the individual network components. While studies on the dynamics of motifs and networks are usually devoted to the temporal or spatial description of processes, this study focuses on the relationship between the specific architecture and the overall rate of the processes of the feedforward family of motifs, including double and triple feedforward loops. The search for the most efficient network architecture could be of particular interest for regulatory or signaling pathways in biology, as well as in computational and communication systems.
Results
Feedforward motif dynamics were studied using cellular automata and compared with differential equation modeling. The number of cellular automata iterations needed for a 100% conversion of a substrate into a target product was used as an inverse measure of the transformation rate. Several basic topological patterns were identified that order the specific feedforward constructions according to the rate of dynamics they enable. At the same number of network nodes and constant other parameters, the biparallel and triparallel motifs provide higher network efficacy than single feedforward motifs. Additionally, a topological property of isodynamicity was identified for feedforward motifs where different network architectures resulted in the same overall rate of the target production.
Conclusion
It was shown for classes of structural motifs with feedforward architecture that network topology affects the overall rate of a process in a quantitatively predictable manner. These fundamental results can be used as a basis for simulating larger networks as combinations of smaller network modules with implications on studying synthetic gene circuits, small regulatory systems, and eventually dynamic wholecell models.
Background
Modeling is a means of making predictions and testing our understanding. In some sense, our level of understanding of an entity can be measured by how well we can model that entity [1]. In particular, mathematical modeling has been applied to diverse areas of science [2], including chemistry [3] and biology [4–6]. The quantitative nature of mathematical modeling has the benefit of yielding detailed, objective descriptions and predictions of processes. An accurate mathematical model can help clarify the roles of individual components within a process and generate specific, testable hypotheses and predictions. The quantitative results of a mathematical model also provide an objective basis for evaluating the accuracy of a model when compared to experimental results and can enable iterative improvement of a model [7, 8].
One of the main benefits of quantitative modeling and analysis is the ability to identify emergent, general properties. An example of this is that complex networks from the internet to biological metabolism have been found to organizationally function and expand as scalefree networks [9]. Within the scalefree framework, genes, proteins, and metabolites are further organized into functional modules with specific structural motifs [10]. Motifs are defined in terms of graph theory [11, 12] as simple connected subgraphs, the abundance of which in a given network is very different from a random graph having the same number of vertices and edges.
One of the most prevalent topological motifs is the feedforward (FF) loop. Feedforward loops have been studied in biological systems and have been found to be involved in a number of different processes including regulatory mechanisms [13] and cell differentiation [14]. Transcriptional regulation is found in most organisms and the dynamics of feedforward motifs in gene regulatory networks has modeled in detail by Alon and coworkers [15–19] to successfully reproduce basic temporal dependencies. Detailed characterization of specific topological motifs should lead to new analyses such as predicting gene regulatory patterns resulting from the aggregation of different topological motifs [20, 21]. The mechanism of feedforward control in the transcriptional network that promotes cell growth was recently elucidated by Palomero et al. [22] and inhibitory feedforward effects in neural circuits have been studied by Klyachko and Stevens [23]. Cordero and Hogeweg [24] have shown that gene evolution depends on the topology of gene regulatory networks. A recent review [25] analyzed the relation between structural modules and dynamics of cellular networks, as a basis for engineering cells to produce desired properties. More abstractly, studies relating topology to function can lead to better understanding of biological network dynamics such as robust dynamical stability [26] and motif specific spatiotemporal dynamics [27, 28]. All of these examples illustrate the connection between biological function and the topology of biological networks.
The present study focuses on the qualitative and quantitative characterization of a variety of feedforward motifs with different architecture, and the relationship between topology and overall process rate in feedforward loops. A variety of stochastic simulation methods exist that can be used to model networks dynamics [29–32]. In our study, cellular automata [33] were selected as a promising method for dynamic modeling of all possible topologies of feedforward loops due to its flexibility, robustness and accuracy. Widely applied in a variety of areas of science and technology, the cellular automata method recently showed great promise for modeling dynamics of complex biological systems [34–38]. In modeling biochemical processes, we followed the general method used by Kier and Cheng [39, 40]. Specifically, this study is focused on studying the purely topological effects on motif productivity. Thus, the cellular automata simulations were conducted in a manner so as to keep all probabilistic parameters constant. This effectively "freezes" the process stochasticity, making our simulations de facto nonstochastic. Results of cellular automata topodynamic pattern simulations were verified with those obtained by parallel ODE and nonlinear differential equations simulations.
Methods
Cellular automata modeling of biochemical networks

Space and time are discrete: There is a twodimensional toroidal grid of cells that is viewed at subsequent equidistant time steps.

At each time step, each cell has a single state – empty or occupied. The cell may be occupied by an enzyme, a substrate, a product, a substrate/enzyme complex or a product/enzyme complex.

The state of a cell at a given time step depends only on its own state and the cell states in its neighborhood, all taken at the previous step.

Each cell has four sites (von Neumann neighborhood), on which interactions can be simulated.

The contents of a cell may break away from an occupied cell or move to join a cell that is occupied. The probabilities for moving, P_{M}, joining, P_{J} (XX) and P_{J} (XY), and breaking away, P_{B}(XX) and P_{B}(XY), are P_{M} = P_{J} = P_{B} = 1. This means that all cells may move, join, and break apart with equal probability.

The overall probability of a movement, P_{M} = 1, is divided into probabilities for movements onto k grid directions, where k = 1–4 is the number of unoccupied neighboring cells.

The only joining of two cells that has a consequence is that between the specific substrate S(i) and the specific enzyme E(i). When such an encounter occurs, a substrateenzyme complex (SE(i)) is formed. This complex has an assigned transitional probability, P_{T} = 0.5, of changing to a new productenzyme complex PE(i).$\text{S}(\text{i})+\text{E}(\text{i})\leftrightarrow \text{SE}(\text{i})\stackrel{{\text{P}}_{\text{T}}}{\to}\text{PE}(\text{i})\to \text{P}(\text{i})+\text{E}(\text{i})$(1)
In such reactions, the transitional probability is regarded as a measure of enzyme activity or propensity, and may be varied within the entire range of values between 0 and 1. In this study, P_{T} was uniformly assumed to be 0.5 to allow reversibility of the first reaction step and to isolate the purely topological effects on motif productivity. By assuming the movement, joining and breaking probabilities are all equal to 1, and by fixing the only essential probabilistic parameter, transitional probability, to a constant value (0.5), we "freeze" the model stochasticity, and have de facto a nonstochastic CA modeling.
The rules, applied at random to all cells, represent one iteration of the modeling procedure, which determines the cells new states and trajectories. A cell and its fourcell environment can acquire any of the five states defined above. The initial state of the system is random and, thus, does not determine subsequent configurations at any iteration. After many iterations, the system reaches a relatively constant configuration (analogous to a chemical steady state), characterized by counts of cells. The model is statistical; many runs are performed, and the number of iterations needed to attain a steadystate is averaged. The specific CA parameters: number of runs, number of cells for different species, grid density, and percent of conversion of the source substrate into the target product were optimized, as discussed in RESULTS.
Differential equations modeling of networks
Importantly, the resulting constantcoefficient linear systems of ODEs can be solved explicitly, yielding exact formulas for each state variable. In particular, the formula describing the evolution of the target state (T) can be used to determine how much time is needed to achieve a specified level of conversion. One may then study the effect of network topology on the dynamics by ranking the conversion times for various feedforward motifs, starting from the same initial conditions.
where k_{1} through k_{4} and ${\overline{k}}_{1}$ through ${\overline{k}}_{4}$ denote rate constants (which for extracting the topological effects on motif productivity are assumed equal) and [E_{1}] through [E_{4}] denote enzyme concentrations. Because nonlinearity prevents solving the ODEs explicitly, they were solved numerically by the forward Euler method with a time step of 0.001 units. There was no need to use a more sophisticated numerical method as we considered the special case in which all rate constants are equal, and the system is not stiff (the eigenvalues of the Jacobian matrices associated with these systems of ODEs are (i) all negative and (ii) are of the same order of magnitude).
where r_{ i }stands for the rate constant of the reversed reaction step i.
Results
Selection of the major CA parameters
Cellular Automata simulation of the dynamics* of feedforward series FFA shown in Fig. 1
Number of Nodes  Lattice Density %  Number of Runs  SD** average  SD %  

50  100  250  1000  
3  3.6  447  469  458  460  7  1.51 
10  143  145  144  144  8  5.56  
60  21  20  21  20  8  39.0  
6  3.6  939  895  893  906  8  0.89 
10  287  293  290  288  10  3.45  
60  41  42  41  41  9  21.6  
9  3.6  1430  1425  1426  1422  2  0.14 
10  458  460  456  458  2  0.44  
60  65  65  65  65  9  13.9 
More detailed simulations were run to quantitatively capture the manner in which the feedforward dynamics vary with different lattice density, D, and the number of nodes, V, in the FFAseries. Equation (12), relating these quantities, was derived from the data of Table 1:
I = 640.5VD^{1.0862}  8.281ln(D) + 36.08
The equation derived shows considerable acceleration of the FF processes with the increase in density in agreement with the law of mass action, since density is proportional to a substrate concentration modeled as number of cells per unit lattice area.
In deriving Eq (12) we used the linear dependence of the number of iterations on the number of FF loop nodes: I = aV + b, which for V = 3, 6, and 9 nodes were obtained with correlation coefficient R^{2} = 0.9961, 0.9988, and 0.9998, respectively. The regressions best expressing the dependence of a and b on the lattice density D were
a = 640.5D^{1.0862} and b = 8.281ln(D) + 36.08
with R^{2} equal to 1.0000 and 0.9901, respectively.
Influence of the size of the feedforward motif on the dynamics of 100% conversion of the source substrate into target product
N*  Iterations Range  R^{2}  SD % 

100  469 – 1425  0.9960  0.99 
300  1603 – 4473  0.9994  0.76 
500  2966 – 8444  0.9996  0.31 
700  4435 – 13282  0.9989  0.28 
900  6397 – 20526  0.9931  0.19 
Modeling seven feedforward series having different topology
Linear dependence of the overall rate of feedforward motifs on the number of motive nodes
Series  N  Iterations range  Regression  R^{2}  SD % 

FFA50  3–9  804 – 2160  219.00 N + 227.14  0.9796  4.14 
FFB50  4–10  1735 – 3736  320.18 N + 483.46  0.9966  1.43 
FFC50  4–10  2043 – 4530  402.36 N + 481.36  0.9980  1.71 
FFD50  4–10  1375 – 3428  329.93 N + 148.50  0.9939  2.38 
FFE50  4–10  1080 – 2247  202.93 N + 258.64  0.9874  2.62 
FFF50  4–10  785 – 1574  133.36 N + 235.07  0.9981  1.70 
FFG50  4–10  1955 – 4104  365.61 N + 395.89  0.9900  2.78 
FFA100  3–9  2966 – 8444  920.75 N + 187.50  0.9996  0.31 
FFB100  4–10  4187 – 9044  800.93 N + 998.07  0.9994  0.29 
FFC100  4–10  4542 – 10565  993.71 N + 556.00  0.9990  0.25 
FFD100  4–10  3291 – 8067  802.64 N + 89.357  0.9993  0.34 
FFE100  4–10  3505 – 7599  683.57 N + 766.14  0.9994  0.33 
FFF100  4–10  2408 – 4927  418.21 N + 705.07  0.9991  0.50 
FFG100  4–10  2721 – 5086  397.61 N + 1121.8  0.9991  0.42 
Discussion
The comparison of the efficacy of performance of the ten fournode networks (Fig. 6) shows that CA and linear ODE order the motifs in the same way with a minor exception. Namely, CA ranks structures H and I with the same highest conversion rate (2408 ± 13 and 2427 ± 15, respectively), since the iteration numbers overlap within the range of their standard deviations. The linear ODE also ranks H and I as the fastest fournode topologies, adding a third structure G, not only showing the same conversion time 2.17, but this time T(t) is given by the exact same formula in all three cases (eq. 17).
The nonlinear ODE models with reversible and irreversible first steps produce identical ordering of the ten structures. It coincides with the ordering of the first seven structures, described above by CA and linear ODEs, while suggesting that network I is the fastest, G and H having very close performance, H shown as slightly slower than G.
The topological analysis of the nine networks revealed some useful patterns of their dynamics. Although the networks analyzed here are relatively simple, they could be of use when analyzing local topology in large complex networks. Several of the observed topodynamic patterns are described below.
Dynamic FeedForward Pattern 1 (DFFP1)
The shorter the graph distance d(S→T) between the source node and the target node in a feedforward motif, the higher the overall conversion rate:
A(d = 3) < B, C (d = 2) < D, E, G, H, I (d = 1)
Note, that the biparallel motifs F and J do not obey this rate inequality.
Dynamic FeedForward Pattern 2 (DFFP2)
The shorter the average path length l(S→T) between the source node and the target node in a feedforward motif, the higher the overall conversion rate:
A ( l = 3) < B, C ( l = 2.5) < D, E, G, H ( l = 2) < I ( l = 5/3)
Accounting for all S → T paths is a slightly more sensitive pattern, which singlesout network I as the most efficient fournode structure, in agreement with the result obtained by the nonlinear ODE model. The biparallel motifs F and J do not obey this feedforward topodynamic trend, which is more important in larger networks where the number of S → T paths increases rapidly.
Dynamic FeedForward Pattern 3 (Isodynamicity)
Some feedforward motifs with different topology produce the same overall S → T conversion rate by the CA and linear ODE models:
CA: H (2408 ± 13) = I(2427 ± 15)
ODE: G = H = I = 2.169053700
The CA and the nonlinear ODE simulations showed these two motifs with different, although relatively close efficacy, the structure J being the slower one:
CA: J (3348 ± 18) < F (3291 ± 17)
Irreversible Nonlinear ODE: J (6.05) < F (5.54)
Reversible Nonlinear ODE: J (9.00) < F (8.51)
Linear ODE: J(2.679) = F(2.679)
The property of isodynamicity is a surprising novel network pattern, which could warrant further detailed studies.
Dynamic FeedForward Pattern 4 (DFFP4)
A < B < C < D
The ringclosures described by this pattern are shown in Fig. 7 with serial numbers 1, 2, and 3. The generality of this topologydynamics relationship was verified for the entire series FFA, FFB, and FFC (Fig. 5) having up to ten motif nodes. In all cases, the standard deviation the number of CA iterations was found to be more than two orders of magnitude smaller than that number. This pattern goes beyond the simple topological patterns 1 and 2 shown above, which cannot discriminate between FFB and FFC series.
Dynamic FeedForward Pattern 5 (DFFP5)
Adding a second feedforward edge (double feedforward motif), between any pair of nodes in the longer path of the FF loop, speeds up the dynamics of the source substrate conversion into the target product (Fig. 8).
D < E < H
These inequalities for the number of iterations, illustrated in Fig. 8 with four node motifs, were verified and found valid with no exceptions for all sizes of the three FF series examined (four to ten loop nodes). Comparing the FFF and FFE series, one may generalize that the acceleration of substratetotarget conversion is higher when the second FFlink starts in a node located on the longer sourcetarget path and ends into the target node, rather than to start in the source node and end in another node before the target one. Since the structures of the triplefeedforward motif FFG (see graph G in Fig. 8) combine the CA trends of the FFF and FFE series, the acceleration in this series is intermediate between the ones of FFE and FFF. However, the ODE models do not confirm this result with the linear model showing G and H to be isodynamic, whereas the two nonlinear models shows G as slightly more efficient than H. Therefore, adding a third feedforward link does not necessarily result in acceleration and no stable trend exists.
Dynamic FeedForward Pattern 6 (DFFP6)
FeedForward < BiParallel < TriParallel
Three such conversions:
D < F, E < I, J < H
are shown in Fig. 6, where they are denoted by asterisks.
Conclusion
The technique of dynamic modeling with cellular automata shows great promise in modeling complex biological systems. Such systems can be broken down to subsystems of smaller scale (to ease computational time) and simulated independently so as to shed light on the processes on a larger scale. The essential element in such applications is the extraction of useful topologicaldynamic (topodynamic) patterns, which identify specific effects of topological structure on the dynamics of network processes while keeping all kinetic parameters constant. The beauty of the topological approach in studies of dynamics is in the generality of the patterns found, which are independent on the nature of the processes, and may be applied to any process of chemical transformation, as well as to any process of mass, energy or information transfer down the forward direction of the motifs.
The dynamics of the feed forward motifs observed in this study revealed important aspects of networks with such components. Not only does any feedforward link added to a linear cascade of chemical/biochemical reactions accelerate the process, but the acceleration is further enhanced by adding a second forward link in the feedforward loop. The topological hierarchy established in this study for fournode motifs predicts that the acceleration of the overall process in such motifs continue increasing with the decrease in the distance (both along the shortest path and along all paths) between the input and output nodes, whereas at the same distance the cellular automata and differential equation simulations produce in a similar manner a further distinction between the motifs dynamic efficacy. The intriguing property of isodynamicity was identified showing motifs with the same number of nodes and different topology to have the same overall rate of inputtooutput transformation. If shown to be present in larger biological networks, the observed isodynamic property could indicate a level of biological robustness at a topological level. Further topologydynamics studies involving construction of networks from combinations of such structural blocks will aid in increasing our understanding of complex biological networks.
Abbreviations
 CA – Cellular Automata:

ODE – Ordinary Differential Equations, FFP – FeedForward Pattern, FFA through FFG – FeedForward Series A through G
Declarations
Acknowledgements
The authors are indebted to Dr. L.B. Kier, VCU, for introducing them to the field of CA modeling.
Authors’ Affiliations
References
 Massoud TF, Hademenos GJ, Young WL, Gao E, PileSpellman J, Vinuela F: Principles and philosophy of modeling in biomedical research. FASEB Journal. 1998, 12: 275285.Google Scholar
 Basmadjian D: The Art of Modeling in Science and Engineering. 1999, Boca Raton, FL: Chapman & Hall/CRCView ArticleMATHGoogle Scholar
 Dryer FL: The Phenomenology of Modeling Combustion Chemistry. 1991, New York: WileyGoogle Scholar
 Kauffman KJ, Prakash P, Edwards JS: Advances in flux balance analysis. Current Opinion in Biotechnology. 2003, 14: 491496.View ArticleGoogle Scholar
 Tomita M, Hashimoto K, Takahashi K, Shimizu TS, Matsuzaki Y, Miyoshi F, Saito K, Tanida S, Yugi K, Venter JC, Hutchison CA: ECELL: software environment for wholecell simulation. Bioinformatics. 1999, 15: 7284.View ArticleGoogle Scholar
 Bellouquid A, Delitala M: Mathematical Modeling of Complex Biological Systems. 2006, Berlin: SpringerMATHGoogle Scholar
 Covert MW, Knight EM, Reed JL, Herrgard MJ, Palsson BO: Integrating highthroughput and computational data elucidates bacterial networks. Nature. 2004, 429: 9296.View ArticleGoogle Scholar
 Reed JL, TR , Chen KH, Joyce AR, Applebee MK, Herring CD, Bui OT, Knight EM, Fong SS, Palsson BO: Systems approach to refining genome annotation. PNAS. 2006, 103: 1748017484.View ArticleGoogle Scholar
 Jeong H, Tombor B, Albert R, Oltvai ZN, Barabási AL: The largescale organization of metabolic networks. Nature. 2000, 407: 651654.View ArticleGoogle Scholar
 Milo R, ShenOrr S, Itzkovitz S, Kashtan N, Chklovskii D, Alon U: Network motifs: Simple building blocks of complex networks. Science. 2002, 298: 824827.View ArticleGoogle Scholar
 Harary F: Graph Theory. 1969, Reading, MA: AddisonWessleyGoogle Scholar
 JL , Yellen J: Handbook of Graph Theory. 2004, Boca Raton, FL: CRC PressGoogle Scholar
 Kashtan N, Alon U: Spontaneous evolution of modularity and network motifs. PNAS USA. 2005, 102: 1377313778.View ArticleGoogle Scholar
 Mattick JS: A new paradigm for developmental biology. Journal of Experimental Biology. 2007, 210: 15261547.View ArticleGoogle Scholar
 Mangan S, Alon U: Structure and function of the feedforward loop network motif. Proc Natl Acad Sci USA. 2003, 100: 1198011985.View ArticleGoogle Scholar
 Milo R, Itzkovtiz S, Kashtan N, Levitt R, ShenOrr S, Ayzenshtat I, Sheffer M, Alon U: Superfamilies of evolved and designed networks. Science. 2004, 303: 15381542.View ArticleGoogle Scholar
 Kashtan N, Itzkovitz S, Milo R, Alon U: Topological generalizations of network motifs. Phys Rev E. 2004, 70: 031909View ArticleGoogle Scholar
 Alon U: Introduction to Systems Biology – Design Principles of Biological Circuits. 2006, Boca Raton, FL: Chapman&Hall/CRCMATHGoogle Scholar
 Alon U: Network motifs: theory and experimental approaches. Nature Rev Genet. 2007, 8: 450461.View ArticleGoogle Scholar
 Dobrin R, Beg QK, Barabási AL, Oltvai ZN: Aggregation of topological motifs in the E. coli transcriptional regulatory network. BMC Bioinformatics. 2004, 5: 10View ArticleGoogle Scholar
 Wong SL, Zhang LV, Tong A, Li Z, Goldberg DS, King OD, Lesage G, Vidal M, Andrews B, Bussey H, Boone C, Roth FP: Combining biological networks to predict genetic interactions. PNAS. 2004, 101: 1568215687.View ArticleGoogle Scholar
 Palomero T, Lim WK, Odom DT, Sulis ML, Real PJ, Margolin A, Barnes KC, O'Neil J, Neuberg D, Weng AP, Aster JC, Sigaux F, Soulier J, Look AT, Young RA, Califano A, Ferrando AA: NOTCH1 directly regulates cMYC and activates a feedforwardloop transcriptional network promoting leukemic cell growth. PNAS. 2006, 103: 1826118266. See also "Correction" in PNAS, 2007: 104, 4240View ArticleGoogle Scholar
 Klyachko VA, Stevens CF: Excitatory and feedforward inhibitory hippocampal synapses work synergistically as an adaptive filter of natural spike trains. PLoS Biology. 2006, 4: e207View ArticleGoogle Scholar
 Cordero OX, Hogeweg P: Feedforward loop circuits as a side effect of genome evolution. Mol Biol Evol. 2006, 23: 19311936.View ArticleGoogle Scholar
 Yuan Qi, Hui Ge: Modularity and dynamics of cellular networks. PLoS Computational Biology. 2006, 2: e174View ArticleGoogle Scholar
 Prill RJ, Iglesias PA, Levchenko A: Dynamic properties of network motifs contribute to biological network organization. PLoS Computational Biology. 2005, 3: e343View ArticleGoogle Scholar
 Zhigulin VP: Dynamical motifs: building blocks of complex network dynamics. 13: arXiv:condmat/0311330
 Sardanyés J, Solé RV: Spatiotemporal dynamics in simple asymmetric hypercycles under weak parasitic coupling. Physica D. 2007, 231: 116129.View ArticleMathSciNetMATHGoogle Scholar
 Wilkinson D: Stochastic Modeling for Systems Biology. 2006, Boca Raton FL: Chapman & Hall/CRCGoogle Scholar
 Gillespie DT: A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J Comput Phys. 1976, 22: 403434.View ArticleMathSciNetGoogle Scholar
 Rathinam M, Petzold LR, Cao Y, Gillespie DT: Stiffness in stochastic chemically reacting systems: The implicit tauleaping method. J Chem Phys. 2003, 119 (24): 1278412794.View ArticleGoogle Scholar
 Bremaud P: Gibbs Fields, Monte Carlo Simulation, and Queues. 2001, Berlin, New York: Springer, 2nd printingGoogle Scholar
 Wolfram S: A New Kind of Science. 2002, Champaign, IL: Wolfram MediaMATHGoogle Scholar
 Deutsch A, Dormann S: Cellular Automaton Modeling of Biological Pattern Formation. 2005, Boston, MA: BirkhauserMATHGoogle Scholar
 Kier LB, Cheng CK, Testa B, Carrupt PA: A cellular automata model of enzyme kinetics. J Molec Graphics. 1996, 14: 227231.View ArticleGoogle Scholar
 Weimar JR: Cellular Automata. Edited by: Bandini S, Tomasini M, Chopard B. 2002, Berlin: Springer, 294303.View ArticleGoogle Scholar
 Kier LB, Bonchev DG, Buck GA: Modeling biochemical networks: A cellular automata approach. Chem Biodiversity. 2005, 2: 233243.View ArticleGoogle Scholar
 Bonchev DG, Kier LB, Cheng CK: Cellular automata (CA) as a basic method for studying network dynamics. Lecture Series on Computer and Computational Sciences. 2006, 6: 581591.Google Scholar
 Kier LB, Cheng CK: A cellular automata model of an anticipatory system. J Molec Graphics. 2000, 18: 2932.View ArticleGoogle Scholar
 Neuforth A, Seybold PG, Kier LB, Cheng CK: Cellular automata models of kinetically and thermodynamically controlled reactions. Int J Chem Kinet. 2000, 32: 529534.View ArticleGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.