Cellular automata simulation of topological effects on the dynamics of feed-forward motifs
© Apte et al; licensee BioMed Central Ltd. 2008
Received: 01 October 2007
Accepted: 27 February 2008
Published: 27 February 2008
Feed-forward motifs are important functional modules in biological and other complex networks. The functionality of feed-forward 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 feed-forward family of motifs, including double and triple feed-forward 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.
Feed-forward 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 feed-forward constructions according to the rate of dynamics they enable. At the same number of network nodes and constant other parameters, the bi-parallel and tri-parallel motifs provide higher network efficacy than single feed-forward motifs. Additionally, a topological property of isodynamicity was identified for feed-forward motifs where different network architectures resulted in the same overall rate of the target production.
It was shown for classes of structural motifs with feed-forward 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 whole-cell models.
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 . In particular, mathematical modeling has been applied to diverse areas of science , including chemistry  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 scale-free networks . Within the scale-free framework, genes, proteins, and metabolites are further organized into functional modules with specific structural motifs . 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 feed-forward (FF) loop. Feed-forward loops have been studied in biological systems and have been found to be involved in a number of different processes including regulatory mechanisms  and cell differentiation . Transcriptional regulation is found in most organisms and the dynamics of feed-forward motifs in gene regulatory networks has modeled in detail by Alon and co-workers [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 feed-forward control in the transcriptional network that promotes cell growth was recently elucidated by Palomero et al.  and inhibitory feed-forward effects in neural circuits have been studied by Klyachko and Stevens . Cordero and Hogeweg  have shown that gene evolution depends on the topology of gene regulatory networks. A recent review  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  and motif specific spatio-temporal 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 feed-forward motifs with different architecture, and the relationship between topology and overall process rate in feed-forward loops. A variety of stochastic simulation methods exist that can be used to model networks dynamics [29–32]. In our study, cellular automata  were selected as a promising method for dynamic modeling of all possible topologies of feed-forward 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 non-stochastic. Results of cellular automata topodynamic pattern simulations were verified with those obtained by parallel ODE and non-linear differential equations simulations.
Cellular automata modeling of biochemical networks
Space and time are discrete: There is a two-dimensional 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, PM, joining, PJ (XX) and PJ (XY), and breaking away, PB(XX) and PB(XY), are PM = PJ = PB = 1. This means that all cells may move, join, and break apart with equal probability.
The overall probability of a movement, PM = 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 substrate-enzyme complex (SE(i)) is formed. This complex has an assigned transitional probability, PT = 0.5, of changing to a new product-enzyme complex PE(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, PT 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 non-stochastic 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 four-cell 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 steady-state 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 constant-coefficient 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 feed-forward motifs, starting from the same initial conditions.
where k1 through k4 and through denote rate constants (which for extracting the topological effects on motif productivity are assumed equal) and [E1] through [E4] 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.
Selection of the major CA parameters
Cellular Automata simulation of the dynamics* of feed-forward series FFA shown in Fig. 1
Number of Nodes
Lattice Density %
Number of Runs
More detailed simulations were run to quantitatively capture the manner in which the feed-forward dynamics vary with different lattice density, D, and the number of nodes, V, in the FFA-series. 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 R2 = 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 R2 equal to 1.0000 and 0.9901, respectively.
Influence of the size of the feed-forward motif on the dynamics of 100% conversion of the source substrate into target product
469 – 1425
1603 – 4473
2966 – 8444
4435 – 13282
6397 – 20526
Modeling seven feed-forward series having different topology
Linear dependence of the overall rate of feed-forward motifs on the number of motive nodes
804 – 2160
219.00 N + 227.14
1735 – 3736
320.18 N + 483.46
2043 – 4530
402.36 N + 481.36
1375 – 3428
329.93 N + 148.50
1080 – 2247
202.93 N + 258.64
785 – 1574
133.36 N + 235.07
1955 – 4104
365.61 N + 395.89
2966 – 8444
920.75 N + 187.50
4187 – 9044
800.93 N + 998.07
4542 – 10565
993.71 N + 556.00
3291 – 8067
802.64 N + 89.357
3505 – 7599
683.57 N + 766.14
2408 – 4927
418.21 N + 705.07
2721 – 5086
397.61 N + 1121.8
The comparison of the efficacy of performance of the ten four-node 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 four-node 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 Feed-Forward Pattern 1 (DFFP1)
The shorter the graph distance d(S→T) between the source node and the target node in a feed-forward motif, the higher the overall conversion rate:
A(d = 3) < B, C (d = 2) < D, E, G, H, I (d = 1)
Note, that the bi-parallel motifs F and J do not obey this rate inequality.
Dynamic Feed-Forward Pattern 2 (DFFP2)
The shorter the average path length l(S→T) between the source node and the target node in a feed-forward 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 singles-out network I as the most efficient four-node structure, in agreement with the result obtained by the nonlinear ODE model. The bi-parallel motifs F and J do not obey this feed-forward topodynamic trend, which is more important in larger networks where the number of S → T paths increases rapidly.
Dynamic Feed-Forward Pattern 3 (Isodynamicity)
Some feed-forward 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 Feed-Forward Pattern 4 (DFFP4)
A < B < C < D
The ring-closures described by this pattern are shown in Fig. 7 with serial numbers 1, 2, and 3. The generality of this topology-dynamics 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 Feed-Forward Pattern 5 (DFFP5)
Adding a second feed-forward edge (double feed-forward 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 substrate-to-target conversion is higher when the second FF-link starts in a node located on the longer source-target 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 triple-feed-forward 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 feed-forward link does not necessarily result in acceleration and no stable trend exists.
Dynamic Feed-Forward Pattern 6 (DFFP6)
Feed-Forward < Bi-Parallel < Tri-Parallel
Three such conversions:
D < F, E < I, J < H
are shown in Fig. 6, where they are denoted by asterisks.
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 topological-dynamic (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 feed-forward 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 feed-forward loop. The topological hierarchy established in this study for four-node 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 input-to-output 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 topology-dynamics studies involving construction of networks from combinations of such structural blocks will aid in increasing our understanding of complex biological networks.
- CA – Cellular Automata:
ODE – Ordinary Differential Equations, FFP – Feed-Forward Pattern, FFA through FFG – Feed-Forward Series A through G
The authors are indebted to Dr. L.B. Kier, VCU, for introducing them to the field of CA modeling.
- Massoud TF, Hademenos GJ, Young WL, Gao E, Pile-Spellman J, Vinuela F: Principles and philosophy of modeling in biomedical research. FASEB Journal. 1998, 12: 275-285.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: 491-496.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: E-CELL: software environment for whole-cell simulation. Bioinformatics. 1999, 15: 72-84.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 high-throughput and computational data elucidates bacterial networks. Nature. 2004, 429: 92-96.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: 17480-17484.View ArticleGoogle Scholar
- Jeong H, Tombor B, Albert R, Oltvai ZN, Barabási AL: The large-scale organization of metabolic networks. Nature. 2000, 407: 651-654.View ArticleGoogle Scholar
- Milo R, Shen-Orr S, Itzkovitz S, Kashtan N, Chklovskii D, Alon U: Network motifs: Simple building blocks of complex networks. Science. 2002, 298: 824-827.View ArticleGoogle Scholar
- Harary F: Graph Theory. 1969, Reading, MA: Addison-WessleyGoogle 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: 13773-13778.View ArticleGoogle Scholar
- Mattick JS: A new paradigm for developmental biology. Journal of Experimental Biology. 2007, 210: 1526-1547.View ArticleGoogle Scholar
- Mangan S, Alon U: Structure and function of the feed-forward loop network motif. Proc Natl Acad Sci USA. 2003, 100: 11980-11985.View ArticleGoogle Scholar
- Milo R, Itzkovtiz S, Kashtan N, Levitt R, Shen-Orr S, Ayzenshtat I, Sheffer M, Alon U: Superfamilies of evolved and designed networks. Science. 2004, 303: 1538-1542.View ArticleGoogle Scholar
- Kashtan N, Itzkovitz S, Milo R, Alon U: Topological generalizations of network motifs. Phys Rev E. 2004, 70: 03-1909View 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: 450-461.View ArticleGoogle Scholar
- Dobrin R, Beg QK, Barabási A-L, Oltvai ZN: Aggregation of topological motifs in the E. coli transcriptional regulatory network. BMC Bioinformatics. 2004, 5: 10-View 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: 15682-15687.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 c-MYC and activates a feed-forward-loop transcriptional network promoting leukemic cell growth. PNAS. 2006, 103: 18261-18266. See also "Correction" in PNAS, 2007: 104, 4240View ArticleGoogle Scholar
- Klyachko VA, Stevens CF: Excitatory and feed-forward inhibitory hippocampal synapses work synergistically as an adaptive filter of natural spike trains. PLoS Biology. 2006, 4: e207-View ArticleGoogle Scholar
- Cordero OX, Hogeweg P: Feed-forward loop circuits as a side effect of genome evolution. Mol Biol Evol. 2006, 23: 1931-1936.View ArticleGoogle Scholar
- Yuan Qi, Hui Ge: Modularity and dynamics of cellular networks. PLoS Computational Biology. 2006, 2: e174-View ArticleGoogle Scholar
- Prill RJ, Iglesias PA, Levchenko A: Dynamic properties of network motifs contribute to biological network organization. PLoS Computational Biology. 2005, 3: e343-View ArticleGoogle Scholar
- Zhigulin VP: Dynamical motifs: building blocks of complex network dynamics. 13: arXiv:cond-mat/0311330
- Sardanyés J, Solé RV: Spatio-temporal dynamics in simple asymmetric hypercycles under weak parasitic coupling. Physica D. 2007, 231: 116-129.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: 403-434.View ArticleMathSciNetGoogle Scholar
- Rathinam M, Petzold LR, Cao Y, Gillespie DT: Stiffness in stochastic chemically reacting systems: The implicit tau-leaping method. J Chem Phys. 2003, 119 (24): 12784-12794.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 C-K, Testa B, Carrupt P-A: A cellular automata model of enzyme kinetics. J Molec Graphics. 1996, 14: 227-231.View ArticleGoogle Scholar
- Weimar JR: Cellular Automata. Edited by: Bandini S, Tomasini M, Chopard B. 2002, Berlin: Springer, 294-303.View ArticleGoogle Scholar
- Kier LB, Bonchev DG, Buck GA: Modeling biochemical networks: A cellular automata approach. Chem Biodiversity. 2005, 2: 233-243.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: 581-591.Google Scholar
- Kier LB, Cheng CK: A cellular automata model of an anticipatory system. J Molec Graphics. 2000, 18: 29-32.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: 529-534.View ArticleGoogle Scholar
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.