- Email: [email protected]
Cell divisions as a mechanism for selection in stable steady states of multi-stationary gene circuits
Physica D 218 (2006) 70–76 www.elsevier.com/locate/physd
Cell divisions as a mechanism for selection in stable steady states of multi-stationary gene circuits Vitaly V. Gursky a,∗ , Konstantin N. Kozlov b , Alexander M. Samsonov a , John Reinitz c a Theoretical Department, The Ioffe Physico-Technical Institute of the Russian Academy of Sciences, 26 Politekhnicheskaya st., St. Petersburg, 194021, Russia b Department of Computational Biology, Center for Advanced Studies, St. Petersburg State Polytechnic University,
29 Politekhnicheskaya st., St. Petersburg, 195259, Russia c Department of Applied Mathematics and Statistics, and Center for Developmental Genetics, Stony Brook University, Stony Brook, NY 11794-3600, USA
Received 9 December 2005; received in revised form 23 March 2006; accepted 12 April 2006 Available online 15 May 2006 Communicated by Y. Nishiura
Abstract We study the properties of attractors and their basins of attraction in gene circuits whose dynamics is accompanied by cell divisions (mitosis). The dynamical equations for such circuits are formulated as hybrid systems with a phase space of variable dimension. As an example, we consider a two-gene circuit which mimics cross-repressing gap genes in the segmentation system of the early Drosophila embryo. We show that the circuit exhibits multi-stationarity and that cell divisions in the circuit substantially reduce the number of steady states reachable by the system, playing, in this way, the specific dynamical role. We clarify two factors underlying the mechanism of selection in steady states and quantify their effects. Because multi-stationarity is a general property underlying such biological phenomena as differentiation and memory, we believe the dynamical role of cell divisions which we report on may be an important element in the correct description of these processes. c 2006 Elsevier B.V. All rights reserved.
Keywords: Multi-stationarity; Regulatory gene networks; Attractors; Cell division; Drosophila segmentation
1. Introduction The mathematical concepts of attractors and their basins of attraction are of fundamental importance for understanding biology. The phenomenon of cell determination demonstrates that morphogenetic fields in embryos have basins of attraction. These exhibit themselves through the phenomenon of error correction (called regulation in classical developmental biology), in which different initial conditions give rise to the same cellular state. In general many cellular states must be produced in a precise spatial pattern, so the overall dynamics involves many basins in a state space of high dimension. A problem arises in explicitly representing such basins mathematically, because even in a simple dynamical model these basins may have very complicated structures. One ∗ Corresponding author. Tel.: +7 812 2929352; fax: +7 812 2971017.
E-mail addresses: [email protected] (V.V. Gursky), [email protected] (K.N. Kozlov), [email protected] (A.M. Samsonov), [email protected] (J. Reinitz). c 2006 Elsevier B.V. All rights reserved. 0167-2789/$ - see front matter doi:10.1016/j.physd.2006.04.009
way to overcome this difficulty is to directly look for related attractors [1–3]. It is frequently the case in developmental biology that these attractors are stationary, and thus it is important to characterize multi-stationarity. In order to do so, it is necessary to understand the effects of certain general biological processes on the dynamics, one of which is cell division, or mitosis. In this paper, we characterize the effects of mitosis on the dynamics of a simple system chosen to capture key elements of realistic models of embryogenesis. The maximum number of possible stationary attractors is an important quantity which characterizes a multi-stationary system. If the biological system is multi-cellular and its state is treated as a direct product of cellular states, it is evident that the true state space of the system will be astronomically huge [4,5]. However, some mechanisms reduce this complexity, which is the topic of this paper. These mechanisms control the variety of reachable steady states. We show that cell divisions, together with passive cell–cell communication, provide such a mechanism in multi-cellular organisms. As an example, we
V.V. Gursky et al. / Physica D 218 (2006) 70–76
consider a model of segmentation gene expression in the early development of Drosophila [6–9]. This process is accompanied by a series of essentially synchronous divisions of nuclei, which populate the surface of the embryo during the first hours of its life and get surrounded by the cell membranes at the onset of gastrulation [10]. In subsequent sections, we show how nuclear cleavages incorporated in the dynamical equations lead to an effective reduction of the set of reachable steady states. The features of the model related to mitosis and responsible for this reduction are clarified, and their relative roles are quantified in a two-gene example. 2. Model Segmentation genes in the fruit fly Drosophila control the development of segments, which are repeating units forming the body of the fly [10]. Immediately following the deposition of a Drosophila egg, a rapid series of 13 almost synchronous nuclear divisions takes place, without the formation of cells. The period between two subsequent nuclear divisions is called the cleavage cycle. In this work it is convenient to denote each cleavage cycle by an integer n and associate cycle n = 0 with the actual cycle 14A, at the end of which segment determination is complete and gastrulation begins. In this terminology, the previous cycles have negative indices n = −1, n = −2, etc. The expression of segmentation genes is to a very good level of approximation a function only of distance along the anterior–posterior (A–P) axis (the long axis of the embryo ellipsoid). This allows us to use models containing only a onedimensional array of nuclei along the A–P axis. Let us denote by M(n) the number of nuclei under consideration in cleavage cycle n. M(n) has the simple recurrence relation M(n) = 2M(n − 1). Denoting the concentration of the ath protein gene product in a nucleus i at time t by via (t), we write a set of ordinary differential equations for N zygotic genes as: ! N X dvia (t) ab b a a = R g T vi (t) + m i − λa via (t) + dt b=1 a + D a (n) (1 − δi1 )(vi−1 (t) − via (t)) a + (1 − δi M(n) )(vi+1 (t) − via (t)) , (1) where a = 1, . . . , N ; i = 1, . . . , M(n); δi1 and δi M(n) are the Kronecker symbols, used to treat the boundary nuclei. The first term on the right-hand side of (1) describes the regulated rate of synthesis of protein from the ath gene. The function g(·) is a monotonicsigmoid rangingfrom zero to one, p and we use g(y) = (1/2) 1 + y/ y 2 + 1 . The regulation of gene a by gene b is characterized by the regulatory matrix element T ab . The term m ia describes the aggregate regulatory effect of various maternal transcription factors on gene a in nucleus i, which is constant in time in the Drosophila embryo. For the small system considered here, we also take m ia to be independent of i. The maximum rate of synthesis for protein a is given by the constant R a . The second term on the right-hand
71
side of Eq. (1) describes the degradation of the ath protein, which is modeled as first order decay with constant rate λa . Diffusion of protein from nucleus i to two adjacent nuclei is described in the third term. The Kronecker symbols in the diffusion term are needed to take care of the fact that there is only one adjacent nucleus for each of the two boundary nuclei. The diffusion constants are assumed to be inversely proportional to the squared distance between adjacent nuclei, a quantity which does not appear explicitly. We instead suppose that after each division the distance is halved, and hence the diffusion constants vary with the cleavage cycle number according to the rule D(n) = 4D(n − 1). The Eq. (1) are augmented with initial conditions, whose values depend on the precise biological situation being modeled. The number of equations in model (1) is not constant but changes with time due to mitosis. Namely, the model includes 1 , . . . , t s at which the maternally defined schedule of times tdiv div nuclear divisions take place. These points in time delimit the k , the number of nuclei is doubled, cleavage cycles. At each tdiv and hence the number of state variables via and the number of equations for them are doubled as well. We apply a rule k for symmetric cell division to the new state variables at tdiv as follows. The concentrations of proteins in the two daughter nuclei are equal to the concentration of the same protein in the mother nucleus [6]: k k a k v aj (tdiv ) = via (tdiv ) = vi+1 (tdiv ),
(2)
where j indexes nuclei before the division and i indexes the nuclei after the cleavage (nucleus j has divided into nuclei i and i + 1), and k denotes the times at which divisions take place. To finish the statement of the model, we introduce the set of time periods τ 1 , . . . , τ s which precede the division times and represent the finite duration of chromosome condensation, separation, and decondensation which take place during k − mitosis. The model equations during the time interval (tdiv k ) has the same form as in (1), but without the synthesis τ k , tdiv term (R a should be set to zero). This is to incorporate the experimental fact that during mitosis protein synthesis shuts down. The model was successfully used in [7] to describe the formation of stripes by the pair-rule gene even-skipped as the result of regulation from gap and maternal genes. The gap gene system was analyzed in depth in a further study using this model [8,9]. Here we address the general question of the collapse of multi-stationarity using a simple two-gene system (N = 2) with easily understood behavior. To elucidate the influence of nuclear divisions on variability of attractors, we take the system starting at various numbers of nuclei and undergoing various numbers of nuclear divisions. The specific parameter values in (1), shown in Table 1, were chosen in order to give the system the key properties of a real gap gene circuit. In particular, the elements of the regulation matrix were taken as to mimic the gap network topology, which exhibits auto-activation and cross-repression [8,9]; thus, we have T aa > 0 and T ab < 0 for a 6= b. It can be shown that such a mutually repressing circuit in one nucleus
72
V.V. Gursky et al. / Physica D 218 (2006) 70–76
Table 1 Parameter values in Eq. (1) with N = 2 T ab Ra ma λa D a (0) 0.12 −0.17 10.37 −3.82 0.04 0.04 −0.16 0.25 11.20 5.50 0.05 0.02 The first line of numbers corresponds to a = 1, and the second to a = 2. As mentioned in the text, we assume m ia = m a . The column with diffusion constants contains their values in cycle n = 0; for earlier cycles the constants should be calculated by using the rule D(n) = 4D(n − 1).
can produce only two stationary attractors, characterized by a maximal expression level of one of the genes and the minimal expression level for the other. The orders of magnitude for all parameter values in Table 1, k and τ k ), were taken as well as the mitosis schedule (times tdiv to be the same as in the gap gene circuit reconstructed from data in [8,9]. A program analyzing attractors and their basins in system (1) was implemented in C and Perl. For each of many initial conditions, the program solves the equations until a preset large time. It then verifies that the state at large times is constant up to a given accuracy and hence that the established state can be regarded as an attractor. (Limit cycles can also be detected but they did not arise in the present application.) The identified attractor is then compared to others already found. If the attractor is new, the attractor and its initial conditions are written to a file; otherwise the initial conditions are added to a list of initial conditions belonging to the basin of a previous identified attractor. The equations were solved using a Bulirsch–Stoer adaptive step size method [11] with the solver numerical accuracy set to 10−8 . The stable steady states of the equations were found with the help of a script written in Mathematica package. 3. Number of attractors in a two-gene system It can be seen from the model equations that, because of nuclear divisions, the system propagates from a phase space of smaller dimension to one of larger dimension. The dimension of the current phase space is N × M(n), which is doubled in each subsequent cleavage cycle. In this section, we show how this fact influences the number of attractors in the system. The final phase space is related to cleavage cycle n = 0. The stationary attractors in this phase space should be linearly stable fixed points, which are solutions to the functional equations obtained by setting the right-hand side of (1) to zero with n = 0. However, due to mitosis, few such stable steady states are reachable by the system whatever initial conditions we take in the previous cleavage cycles. To show this, we have calculated the stable steady states and the actual attractors that occur in the two-gene system with various numbers of nuclear cleavages and initial nuclei. The initial conditions were varied in a rescaled phase box which was constrained by zero concentration from below and by 1 from above. The results of these calculations are shown in Table 2. It can be seen from the table that the presence of cell divisions in the system leads to an essential decrease in the number of reachable attractors: the more nuclear cleavages the system undergoes on its way to cycle 0, the fewer attractors are accessible to it.
Table 2 Numbers of attractors in the two-gene system starting at various sets of initial nuclei and undergoing various amounts of nuclear divisions, in comparison with the corresponding numbers of stable steady states in the system M(0)
S
A1
A2
A3
1 2 4 8 16
2 4 12 116 6076∗
– 2 4 16 256
– – 2 4 16
– – – 2 4
In the table, M(0) is the number of nuclei at cycle n = 0, S is the number of stable steady states, A1 , A2 , and A3 are the numbers of attractors when the start is at cycle n = −1 (1 division), cycle n = −2 (2 divisions), and cycle n = −3 (3 divisions), respectively. By “−”, we mark situations when the calculations are impossible because we cannot have fewer than 2d nuclei after d consecutive divisions. Each number of attractors in the table was calculated by using 105 different initial conditions. The number in the table with mark “∗” is only a current lower bound of the actual number of steady states. It was obtained by the Mathematica script starting from 4 377 057 random initial approximation points in the Newton method and took several weeks of Intel P4 3.7 GHz computer work. However, at the point of paper submission the script was still adding new steady states in a file quite regularly; therefore, it is most probable that the system has more steady states.
In Fig. 1, we show examples of attractors that can be reached by the system and stable steady states forbidden for the dynamics due to mitosis. The graph represents calculations in the two-gene system which undergoes one division and finishes at 8 nuclei. It can be seen from Table 2 that there are in total 116 stable steady states in such system and only 16 of them occur as attractors when initial conditions are varied in the phase box. 4. Unfolding the mechanism of attractor selection There are two factors responsible for the mechanism of attractor selection due to cell divisions. Both factors represent different constraints on the system, which prevent it from moving to certain steady states. Here, we describe and quantify these factors in the particular example of the two-gene system starting in cycle n = −1 at 4 nuclei and undergoing 1 nuclear division (116 stable steady states and 16 attractors; see Table 2). 4.1. The constraint by the division rule The first constraint is clearly set by the division rule (2). This rule forbids certain states in the phase space for the dynamical system in cycle 0 immediately following the cleavage, whatever initial conditions it takes in cycle −1. Specifically, it follows from (2) that, at this time, the system may occupy only those states via which obey the constraint: a via (tdiv ) = vi+1 (tdiv ),
i = 1, 3, 5, . . . , M(0) − 1,
where tdiv is the cleavage time delimiting cycle −1 and cycle 0, considered as belonging to the latter cycle. For example, there exist no initial conditions in cycle −1 that permit the concentrations v1a and v2a to be different at the moment of entrance into the final phase space. It is evident that such forbidden states are very numerous and are apt to lie in the basins of attraction for certain stable steady states, which may become unreachable for the system dynamics.
V.V. Gursky et al. / Physica D 218 (2006) 70–76
73
Fig. 1. Examples of attractors (left column of panels) reached by the system from different initial conditions, and some of the stable steady states (right column of panels) forbidden for the dynamics due to mitosis. Protein concentrations for gene one (solid curves) and gene two (dashed curves) are shown in conventional units as functions of nucleus number. The model solved is the two-gene system starting at 4 nuclei at cycle n = −1 and finishing at 8 nuclei at cycle n = 0, undergoing 1 nuclear cleavage. The attractors in the figure are solutions at t = 10 000. The steady states are linearly stable roots of the right-hand side of model equations in the final phase space. The concentration values calculated in each nucleus are connected on the picture by the straight lines in between the nuclei to ease observation.
To quantify how strong this factor is, we have calculated the number of attractors in the case where protein concentrations in the daughter nuclei immediately after the cleavage are not exactly the same as in the mother one but are perturbed with some noise. This assumption is more realistic than the rule (2), because these values for the concentrations may differ significantly in data [12]. Therefore, we consider a modified division rule of the form: via (tdiv ) = v aj (tdiv ) + ξia ,
a a vi+1 (tdiv ) = v aj (tdiv ) + ξi+1 , (3)
where, as above, nucleus j ( j = 1, . . . , 4) from cycle −1 divides into nuclei i and i + 1 in cycle 0; ξia is an addition due to the noise effects. We have calculated the number of attractors in the model with the modified division rule (3), taking ξia as random quantity uniformly distributed in the interval [−ξmax , ξmax ], where the noise amplitude ξmax > 0 can be varied at will.
During the calculation, the value of ξia (for each i and a) was different for each different initial condition. Therefore, as we have used 105 different initial conditions to calculate attractors, we encountered the same number of random values of ξia in the calculation. We performed several such runs with various values ξmax , obtaining various numbers of reachable attractors. The idea of the numerical experiment was to try to cover with the noise a larger and larger part of the final phase space at the division time, in order to let the system fall into more and more basins of attraction for attractors which are forbidden in the unperturbed case. The result of the calculations is shown in Fig. 2. The figure shows that, at a noise value up to 0.08, the number of attractors stays close to the original value of 16, given in Table 2. As the amplitude of noise rises, the number of attractors grows, since the system has a larger range of starting points in cycle 0. Finally, with ξmax close to 0.35, the number of attractors
74
V.V. Gursky et al. / Physica D 218 (2006) 70–76
Fig. 2. Numbers of attractors calculated in the same system as for Fig. 1, except that the modified division rule (3) is used with various noise amplitudes ξmax . The values of ξmax are put on the horizontal axis, and the corresponding numbers of attractors are on the vertical one. For a given amplitude ξmax , the attractors are calculated by using 105 different initial conditions.
approaches 116, the number of all possible stable steady states. 4.2. Constraint by history The second constraint on the system’s ability to choose between stable steady states is provided by the history of the system, i.e. by its dynamics during the time preceding the entrance into the final phase space. It is evident that the solution trajectories come into cycle 0 already arranged into domains corresponding to the basins of attraction of attracting states in the last-but-final phase space, which might be actual attractors if the division were not to occur. This effect is shown in Fig. 3. It is clear from Fig. 3A that in the particular nucleus the solution trajectories gather in two
domains. Consequently, the whole blank area of the picture will not be represented in cycle 0, and attracting states from that area will not be reached. The effect of clustering of different attracting states in Fig. 3B is related to the weakness of coupling between adjacent nuclei, as values of the diffusion constants are small [7,9]. This effect means that, despite the large number of potential attractors, the related final protein concentrations in a given nucleus would accumulate in isolated islands inside the two-dimensional projection of the phase space. To quantify the constraint by history, we calculated attractors, starting from many uniformly distributed random initial conditions, in a series of models which was obtained from the basic model under consideration by gradual reduction of the duration of cycle −1. It is clear that as cycle −1 becomes shorter, there is less time to canalize trajectories prior to cycle 0, and hence the system will be free to occupy more states during this cycle. The results are shown in Fig. 4, where the numbers of attractors are plotted for 10 durations of cycle −1, ranging from 0.1 to 21.1 min, which is the actual duration of the cycle prototype in the real embryo. It is seen from the figure that the number of attractors rises only when cycle −1 is shorter than about 9 min, i.e. about half its natural length. Another interesting feature is that 57 attractors occur with cycle −1 set to the shortest duration of 0.1 min, compared to 116 stable steady states total in the model. About half of the steady states still lie in forbidden basins because of the constraint (2), since the division rule is used in all calculations. 5. Conclusions During the early development of multicellular organisms, many molecular processes, such as gene expression, occur in cells which undergo numerous divisions. Segmentation gene
Fig. 3. Two snapshots of a two-dimensional projection of the phase space in the same system as for the previous figures. The axes stand for the concentrations, in conventional units, of proteins from the two genes in the first nucleus (variables v11 and v12 ). In panel A, the snapshot is made at the end of cycle n = −1, right before the division. The dots in the panel represent the two-dimensional parts of solutions of the model at the end of cycle −1, with 104 different initial conditions taken at the beginning of cycle −1. Panel B corresponds to the snapshot at infinite time in the first nucleus, which is one of the two progenies of nucleus from panel A (the first nucleus in cycle −1 divides into nuclei 1 and 2 in cycle 0). Panel B shows 16 actual attractors (triangles) which the system goes to from the same initial conditions as for panel A. The panel also contains those of 116 stable steady states (daggers) which can never be reached by the system from any initial condition. Of the 16 triangles in B, 8 are hardly seen because they stay very close to the other 8. However, this does not mean that these two sets of attractors are identical: the triangles represent only the two-dimensional part of the attractors, and the two attractor sets are not similar in other dimensions. The same remark is valid for closely standing daggers of non-reachable stable steady states.
V.V. Gursky et al. / Physica D 218 (2006) 70–76
Fig. 4. Numbers of attractors in the same system as for previous figures, but with 10 different durations (in minutes) of cycle n = −1. The smallest duration is 0.1 min, and the largest is 21.1 min. Each number of attractors was calculated by using 105 different initial conditions.
expression in the early Drosophila embryo is an example of such a process. We have shown that nuclear divisions, which accompany gene expression in the system, play the dynamical role in this process. As each nucleus divides, adding new state variables to the system, the dimension of the state space increases from cycle to cycle. The complexity of the final phase space allows the model to have many stable steady states, which can be identified as the potential final expression patterns. However, because the system starts in a space of smaller dimension, only a limited number of these steady states are actually accessible to the system. Therefore, cell divisions provide a mechanism for the selection of possible attractors. We have identified two factors underlying the mechanism of selection and numerically illustrated their effect in a small two-gene system. The first factor that renders some basins of attraction inaccessible is provided by symmetric nuclear divisions. The second constraining factor is the dynamics itself at early times, which arranges the solution trajectories in narrow domains in the phase space, such that the attracting sets of some steady states are not captured by these domains. Inspection of Table 2 provides a rough estimate of the magnitude of the reduction in the number of attractors because of cell divisions. It is seen that the number of stable steady states grows nonlinearly and very rapidly with the final dimension of the phase space. At the same time, we find that the number Ad of attractors in the two-gene system with the chosen parameter values obeys the binary rule Ad = 2 Min , where d is the total number of divisions that occur in the model, and Min is the number of initial nuclei; this rule is valid only for d > 0. (The binary nature of this formula stems from the fact that the two genes have mutually exclusive states of expression because of mutual repression.) Thus, whatever the final dimension of the phase space, the number Ad is determined as soon as the numbers Min and d are chosen. This means that the relative fraction of reachable to total attractors in the model decreases rapidly with increasing d and/or Min . In particular, it follows from Table 2 that, for Min rising from 1 to 8, the fraction of attractors in the one-cleavage system, measured as a percentage of the number of all stable steady states, decreases from 50%
75
to, at maximum, 4%. Similarly, in the two-cleavage system with Min varying from 1 to 4 this fraction ranges from 17% to, at maximum, 0.3%. If we consider the two-gene systems finishing in the same phase space of dimension 32 (16 nuclei in the last cleavage cycle) but having different number d of nuclear divisions, then, for d varying from 1 to 3, the fraction decreases from, at maximum, 4% to, at maximum, 0.07%. It should be noted though that we do not have evidence that the described binary rule for the number of attractors, together with the precise percent calculated, can be quantitatively extended to more general systems. However, it is likely that the observed tendency of decrease in the portion of accessible attractors is a generic mechanism in the presence of approximately symmetric mitoses. Molecular transport between cells provides another mechanism for reducing the number of attractors in multi-cellular systems. In model (1), the transport of material takes place by the exchange of proteins between adjacent nuclei by diffusion. The effect of diffusion can be seen in the S-column of Table 2. Each entry in this column contains the number of stable steady states which are attractors in a system that has no divisions. As described above, there can be only 2 attracting states in one nucleus and, therefore, 2 Min of them if the multi-nuclear system is treated as the direct product of nuclei. The S-column of the table shows that this rule breaks down if the system has no divisions and more than 2 nuclei. The numbers of attractors are less than 2 Min in such cases, and this fact is a consequence of the diffusive coupling between nuclei. Acknowledgements The support of this work by NIH Grants RUB1-1578 and CRDF GAP award RBO1286, and by RFBR-NWO grant 047.011.2004.013 is gratefully acknowledged. References [1] R. Thomas, M. Kaufman, Multistationarity, the basis of cell differentiation and memory. II. Logical analysis of regulatory networks in terms of feedback circuits, Chaos 11 (1) (2001) 180–195. [2] R. Thomas, R. D’Ari, Biological Feedback, CRC, Boca Raton, FL, 1990. [3] S.A. Kauffman, Origins of Order: Self-Organization and Selection in Evolution, Oxford University Press, Oxford, 1993. [4] L. S´anchez, D. Thieffry, A logical analysis of the Drosophila gap-gene system, J. Theoret. Biol. 211 (2001) 115–141. [5] L. S´anchez, D. Thieffry, Segmenting the fly embryo: A logical analysis of the pair-rule cross-regulatory module, J. Theoret. Biol. 224 (2003) 517–537. [6] E. Mjolsness, D.H. Sharp, J. Reinitz, A connectionist model of development, J. Theoret. Biol. 152 (1991) 429–453. [7] J. Reinitz, D.H. Sharp, Mechanism of eve stripe formation, Mech. Dev. 49 (1995) 133–158. [8] J. Jaeger, S. Surkova, M. Blagov, H. Janssens, D. Kosman, K.N. Kozlov, Manu, E. Myasnikova, C.E. Vanario-Alonso, M. Samsonova, D.H. Sharp, J. Reinitz, Dynamic control of positional information in the early Drosophila embryo, Nature 430 (2004) 368–371. [9] J. Jaeger, M. Blagov, D. Kosman, K.N. Kozlov, Manu, E. Myasnikova, S. Surkova, C.E. Vanario-Alonso, M. Samsonova, O.H. Sharp, J. Reinitz, Dynamical analysis of regulatory interactions in the gap gene system of Drosophila melanogaster, Genetics 167 (2004) 1721–1737.
76
V.V. Gursky et al. / Physica D 218 (2006) 70–76
[10] P.A. Lawrence, The Making of a Fly, Blackwell Sci. Publ., Oxford, 1992. [11] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in C, Cambridge University Press, Cambridge, UK, 1992.
[12] E. Poustelnikova, A. Pisarev, M. Blagov, M. Samsonova, J. Reinitz, A database for management of gene expression data /in situ/, Bioinformatics 20 (2004) 2212–2221.
Cell divisions as a mechanism for selection in stable steady states of multi-stationary gene circuits Vitaly V. Gursky a,∗ , Konstantin N. Kozlov b , Alexander M. Samsonov a , John Reinitz c a Theoretical Department, The Ioffe Physico-Technical Institute of the Russian Academy of Sciences, 26 Politekhnicheskaya st., St. Petersburg, 194021, Russia b Department of Computational Biology, Center for Advanced Studies, St. Petersburg State Polytechnic University,
29 Politekhnicheskaya st., St. Petersburg, 195259, Russia c Department of Applied Mathematics and Statistics, and Center for Developmental Genetics, Stony Brook University, Stony Brook, NY 11794-3600, USA
Received 9 December 2005; received in revised form 23 March 2006; accepted 12 April 2006 Available online 15 May 2006 Communicated by Y. Nishiura
Abstract We study the properties of attractors and their basins of attraction in gene circuits whose dynamics is accompanied by cell divisions (mitosis). The dynamical equations for such circuits are formulated as hybrid systems with a phase space of variable dimension. As an example, we consider a two-gene circuit which mimics cross-repressing gap genes in the segmentation system of the early Drosophila embryo. We show that the circuit exhibits multi-stationarity and that cell divisions in the circuit substantially reduce the number of steady states reachable by the system, playing, in this way, the specific dynamical role. We clarify two factors underlying the mechanism of selection in steady states and quantify their effects. Because multi-stationarity is a general property underlying such biological phenomena as differentiation and memory, we believe the dynamical role of cell divisions which we report on may be an important element in the correct description of these processes. c 2006 Elsevier B.V. All rights reserved.
Keywords: Multi-stationarity; Regulatory gene networks; Attractors; Cell division; Drosophila segmentation
1. Introduction The mathematical concepts of attractors and their basins of attraction are of fundamental importance for understanding biology. The phenomenon of cell determination demonstrates that morphogenetic fields in embryos have basins of attraction. These exhibit themselves through the phenomenon of error correction (called regulation in classical developmental biology), in which different initial conditions give rise to the same cellular state. In general many cellular states must be produced in a precise spatial pattern, so the overall dynamics involves many basins in a state space of high dimension. A problem arises in explicitly representing such basins mathematically, because even in a simple dynamical model these basins may have very complicated structures. One ∗ Corresponding author. Tel.: +7 812 2929352; fax: +7 812 2971017.
E-mail addresses: [email protected] (V.V. Gursky), [email protected] (K.N. Kozlov), [email protected] (A.M. Samsonov), [email protected] (J. Reinitz). c 2006 Elsevier B.V. All rights reserved. 0167-2789/$ - see front matter doi:10.1016/j.physd.2006.04.009
way to overcome this difficulty is to directly look for related attractors [1–3]. It is frequently the case in developmental biology that these attractors are stationary, and thus it is important to characterize multi-stationarity. In order to do so, it is necessary to understand the effects of certain general biological processes on the dynamics, one of which is cell division, or mitosis. In this paper, we characterize the effects of mitosis on the dynamics of a simple system chosen to capture key elements of realistic models of embryogenesis. The maximum number of possible stationary attractors is an important quantity which characterizes a multi-stationary system. If the biological system is multi-cellular and its state is treated as a direct product of cellular states, it is evident that the true state space of the system will be astronomically huge [4,5]. However, some mechanisms reduce this complexity, which is the topic of this paper. These mechanisms control the variety of reachable steady states. We show that cell divisions, together with passive cell–cell communication, provide such a mechanism in multi-cellular organisms. As an example, we
V.V. Gursky et al. / Physica D 218 (2006) 70–76
consider a model of segmentation gene expression in the early development of Drosophila [6–9]. This process is accompanied by a series of essentially synchronous divisions of nuclei, which populate the surface of the embryo during the first hours of its life and get surrounded by the cell membranes at the onset of gastrulation [10]. In subsequent sections, we show how nuclear cleavages incorporated in the dynamical equations lead to an effective reduction of the set of reachable steady states. The features of the model related to mitosis and responsible for this reduction are clarified, and their relative roles are quantified in a two-gene example. 2. Model Segmentation genes in the fruit fly Drosophila control the development of segments, which are repeating units forming the body of the fly [10]. Immediately following the deposition of a Drosophila egg, a rapid series of 13 almost synchronous nuclear divisions takes place, without the formation of cells. The period between two subsequent nuclear divisions is called the cleavage cycle. In this work it is convenient to denote each cleavage cycle by an integer n and associate cycle n = 0 with the actual cycle 14A, at the end of which segment determination is complete and gastrulation begins. In this terminology, the previous cycles have negative indices n = −1, n = −2, etc. The expression of segmentation genes is to a very good level of approximation a function only of distance along the anterior–posterior (A–P) axis (the long axis of the embryo ellipsoid). This allows us to use models containing only a onedimensional array of nuclei along the A–P axis. Let us denote by M(n) the number of nuclei under consideration in cleavage cycle n. M(n) has the simple recurrence relation M(n) = 2M(n − 1). Denoting the concentration of the ath protein gene product in a nucleus i at time t by via (t), we write a set of ordinary differential equations for N zygotic genes as: ! N X dvia (t) ab b a a = R g T vi (t) + m i − λa via (t) + dt b=1 a + D a (n) (1 − δi1 )(vi−1 (t) − via (t)) a + (1 − δi M(n) )(vi+1 (t) − via (t)) , (1) where a = 1, . . . , N ; i = 1, . . . , M(n); δi1 and δi M(n) are the Kronecker symbols, used to treat the boundary nuclei. The first term on the right-hand side of (1) describes the regulated rate of synthesis of protein from the ath gene. The function g(·) is a monotonicsigmoid rangingfrom zero to one, p and we use g(y) = (1/2) 1 + y/ y 2 + 1 . The regulation of gene a by gene b is characterized by the regulatory matrix element T ab . The term m ia describes the aggregate regulatory effect of various maternal transcription factors on gene a in nucleus i, which is constant in time in the Drosophila embryo. For the small system considered here, we also take m ia to be independent of i. The maximum rate of synthesis for protein a is given by the constant R a . The second term on the right-hand
71
side of Eq. (1) describes the degradation of the ath protein, which is modeled as first order decay with constant rate λa . Diffusion of protein from nucleus i to two adjacent nuclei is described in the third term. The Kronecker symbols in the diffusion term are needed to take care of the fact that there is only one adjacent nucleus for each of the two boundary nuclei. The diffusion constants are assumed to be inversely proportional to the squared distance between adjacent nuclei, a quantity which does not appear explicitly. We instead suppose that after each division the distance is halved, and hence the diffusion constants vary with the cleavage cycle number according to the rule D(n) = 4D(n − 1). The Eq. (1) are augmented with initial conditions, whose values depend on the precise biological situation being modeled. The number of equations in model (1) is not constant but changes with time due to mitosis. Namely, the model includes 1 , . . . , t s at which the maternally defined schedule of times tdiv div nuclear divisions take place. These points in time delimit the k , the number of nuclei is doubled, cleavage cycles. At each tdiv and hence the number of state variables via and the number of equations for them are doubled as well. We apply a rule k for symmetric cell division to the new state variables at tdiv as follows. The concentrations of proteins in the two daughter nuclei are equal to the concentration of the same protein in the mother nucleus [6]: k k a k v aj (tdiv ) = via (tdiv ) = vi+1 (tdiv ),
(2)
where j indexes nuclei before the division and i indexes the nuclei after the cleavage (nucleus j has divided into nuclei i and i + 1), and k denotes the times at which divisions take place. To finish the statement of the model, we introduce the set of time periods τ 1 , . . . , τ s which precede the division times and represent the finite duration of chromosome condensation, separation, and decondensation which take place during k − mitosis. The model equations during the time interval (tdiv k ) has the same form as in (1), but without the synthesis τ k , tdiv term (R a should be set to zero). This is to incorporate the experimental fact that during mitosis protein synthesis shuts down. The model was successfully used in [7] to describe the formation of stripes by the pair-rule gene even-skipped as the result of regulation from gap and maternal genes. The gap gene system was analyzed in depth in a further study using this model [8,9]. Here we address the general question of the collapse of multi-stationarity using a simple two-gene system (N = 2) with easily understood behavior. To elucidate the influence of nuclear divisions on variability of attractors, we take the system starting at various numbers of nuclei and undergoing various numbers of nuclear divisions. The specific parameter values in (1), shown in Table 1, were chosen in order to give the system the key properties of a real gap gene circuit. In particular, the elements of the regulation matrix were taken as to mimic the gap network topology, which exhibits auto-activation and cross-repression [8,9]; thus, we have T aa > 0 and T ab < 0 for a 6= b. It can be shown that such a mutually repressing circuit in one nucleus
72
V.V. Gursky et al. / Physica D 218 (2006) 70–76
Table 1 Parameter values in Eq. (1) with N = 2 T ab Ra ma λa D a (0) 0.12 −0.17 10.37 −3.82 0.04 0.04 −0.16 0.25 11.20 5.50 0.05 0.02 The first line of numbers corresponds to a = 1, and the second to a = 2. As mentioned in the text, we assume m ia = m a . The column with diffusion constants contains their values in cycle n = 0; for earlier cycles the constants should be calculated by using the rule D(n) = 4D(n − 1).
can produce only two stationary attractors, characterized by a maximal expression level of one of the genes and the minimal expression level for the other. The orders of magnitude for all parameter values in Table 1, k and τ k ), were taken as well as the mitosis schedule (times tdiv to be the same as in the gap gene circuit reconstructed from data in [8,9]. A program analyzing attractors and their basins in system (1) was implemented in C and Perl. For each of many initial conditions, the program solves the equations until a preset large time. It then verifies that the state at large times is constant up to a given accuracy and hence that the established state can be regarded as an attractor. (Limit cycles can also be detected but they did not arise in the present application.) The identified attractor is then compared to others already found. If the attractor is new, the attractor and its initial conditions are written to a file; otherwise the initial conditions are added to a list of initial conditions belonging to the basin of a previous identified attractor. The equations were solved using a Bulirsch–Stoer adaptive step size method [11] with the solver numerical accuracy set to 10−8 . The stable steady states of the equations were found with the help of a script written in Mathematica package. 3. Number of attractors in a two-gene system It can be seen from the model equations that, because of nuclear divisions, the system propagates from a phase space of smaller dimension to one of larger dimension. The dimension of the current phase space is N × M(n), which is doubled in each subsequent cleavage cycle. In this section, we show how this fact influences the number of attractors in the system. The final phase space is related to cleavage cycle n = 0. The stationary attractors in this phase space should be linearly stable fixed points, which are solutions to the functional equations obtained by setting the right-hand side of (1) to zero with n = 0. However, due to mitosis, few such stable steady states are reachable by the system whatever initial conditions we take in the previous cleavage cycles. To show this, we have calculated the stable steady states and the actual attractors that occur in the two-gene system with various numbers of nuclear cleavages and initial nuclei. The initial conditions were varied in a rescaled phase box which was constrained by zero concentration from below and by 1 from above. The results of these calculations are shown in Table 2. It can be seen from the table that the presence of cell divisions in the system leads to an essential decrease in the number of reachable attractors: the more nuclear cleavages the system undergoes on its way to cycle 0, the fewer attractors are accessible to it.
Table 2 Numbers of attractors in the two-gene system starting at various sets of initial nuclei and undergoing various amounts of nuclear divisions, in comparison with the corresponding numbers of stable steady states in the system M(0)
S
A1
A2
A3
1 2 4 8 16
2 4 12 116 6076∗
– 2 4 16 256
– – 2 4 16
– – – 2 4
In the table, M(0) is the number of nuclei at cycle n = 0, S is the number of stable steady states, A1 , A2 , and A3 are the numbers of attractors when the start is at cycle n = −1 (1 division), cycle n = −2 (2 divisions), and cycle n = −3 (3 divisions), respectively. By “−”, we mark situations when the calculations are impossible because we cannot have fewer than 2d nuclei after d consecutive divisions. Each number of attractors in the table was calculated by using 105 different initial conditions. The number in the table with mark “∗” is only a current lower bound of the actual number of steady states. It was obtained by the Mathematica script starting from 4 377 057 random initial approximation points in the Newton method and took several weeks of Intel P4 3.7 GHz computer work. However, at the point of paper submission the script was still adding new steady states in a file quite regularly; therefore, it is most probable that the system has more steady states.
In Fig. 1, we show examples of attractors that can be reached by the system and stable steady states forbidden for the dynamics due to mitosis. The graph represents calculations in the two-gene system which undergoes one division and finishes at 8 nuclei. It can be seen from Table 2 that there are in total 116 stable steady states in such system and only 16 of them occur as attractors when initial conditions are varied in the phase box. 4. Unfolding the mechanism of attractor selection There are two factors responsible for the mechanism of attractor selection due to cell divisions. Both factors represent different constraints on the system, which prevent it from moving to certain steady states. Here, we describe and quantify these factors in the particular example of the two-gene system starting in cycle n = −1 at 4 nuclei and undergoing 1 nuclear division (116 stable steady states and 16 attractors; see Table 2). 4.1. The constraint by the division rule The first constraint is clearly set by the division rule (2). This rule forbids certain states in the phase space for the dynamical system in cycle 0 immediately following the cleavage, whatever initial conditions it takes in cycle −1. Specifically, it follows from (2) that, at this time, the system may occupy only those states via which obey the constraint: a via (tdiv ) = vi+1 (tdiv ),
i = 1, 3, 5, . . . , M(0) − 1,
where tdiv is the cleavage time delimiting cycle −1 and cycle 0, considered as belonging to the latter cycle. For example, there exist no initial conditions in cycle −1 that permit the concentrations v1a and v2a to be different at the moment of entrance into the final phase space. It is evident that such forbidden states are very numerous and are apt to lie in the basins of attraction for certain stable steady states, which may become unreachable for the system dynamics.
V.V. Gursky et al. / Physica D 218 (2006) 70–76
73
Fig. 1. Examples of attractors (left column of panels) reached by the system from different initial conditions, and some of the stable steady states (right column of panels) forbidden for the dynamics due to mitosis. Protein concentrations for gene one (solid curves) and gene two (dashed curves) are shown in conventional units as functions of nucleus number. The model solved is the two-gene system starting at 4 nuclei at cycle n = −1 and finishing at 8 nuclei at cycle n = 0, undergoing 1 nuclear cleavage. The attractors in the figure are solutions at t = 10 000. The steady states are linearly stable roots of the right-hand side of model equations in the final phase space. The concentration values calculated in each nucleus are connected on the picture by the straight lines in between the nuclei to ease observation.
To quantify how strong this factor is, we have calculated the number of attractors in the case where protein concentrations in the daughter nuclei immediately after the cleavage are not exactly the same as in the mother one but are perturbed with some noise. This assumption is more realistic than the rule (2), because these values for the concentrations may differ significantly in data [12]. Therefore, we consider a modified division rule of the form: via (tdiv ) = v aj (tdiv ) + ξia ,
a a vi+1 (tdiv ) = v aj (tdiv ) + ξi+1 , (3)
where, as above, nucleus j ( j = 1, . . . , 4) from cycle −1 divides into nuclei i and i + 1 in cycle 0; ξia is an addition due to the noise effects. We have calculated the number of attractors in the model with the modified division rule (3), taking ξia as random quantity uniformly distributed in the interval [−ξmax , ξmax ], where the noise amplitude ξmax > 0 can be varied at will.
During the calculation, the value of ξia (for each i and a) was different for each different initial condition. Therefore, as we have used 105 different initial conditions to calculate attractors, we encountered the same number of random values of ξia in the calculation. We performed several such runs with various values ξmax , obtaining various numbers of reachable attractors. The idea of the numerical experiment was to try to cover with the noise a larger and larger part of the final phase space at the division time, in order to let the system fall into more and more basins of attraction for attractors which are forbidden in the unperturbed case. The result of the calculations is shown in Fig. 2. The figure shows that, at a noise value up to 0.08, the number of attractors stays close to the original value of 16, given in Table 2. As the amplitude of noise rises, the number of attractors grows, since the system has a larger range of starting points in cycle 0. Finally, with ξmax close to 0.35, the number of attractors
74
V.V. Gursky et al. / Physica D 218 (2006) 70–76
Fig. 2. Numbers of attractors calculated in the same system as for Fig. 1, except that the modified division rule (3) is used with various noise amplitudes ξmax . The values of ξmax are put on the horizontal axis, and the corresponding numbers of attractors are on the vertical one. For a given amplitude ξmax , the attractors are calculated by using 105 different initial conditions.
approaches 116, the number of all possible stable steady states. 4.2. Constraint by history The second constraint on the system’s ability to choose between stable steady states is provided by the history of the system, i.e. by its dynamics during the time preceding the entrance into the final phase space. It is evident that the solution trajectories come into cycle 0 already arranged into domains corresponding to the basins of attraction of attracting states in the last-but-final phase space, which might be actual attractors if the division were not to occur. This effect is shown in Fig. 3. It is clear from Fig. 3A that in the particular nucleus the solution trajectories gather in two
domains. Consequently, the whole blank area of the picture will not be represented in cycle 0, and attracting states from that area will not be reached. The effect of clustering of different attracting states in Fig. 3B is related to the weakness of coupling between adjacent nuclei, as values of the diffusion constants are small [7,9]. This effect means that, despite the large number of potential attractors, the related final protein concentrations in a given nucleus would accumulate in isolated islands inside the two-dimensional projection of the phase space. To quantify the constraint by history, we calculated attractors, starting from many uniformly distributed random initial conditions, in a series of models which was obtained from the basic model under consideration by gradual reduction of the duration of cycle −1. It is clear that as cycle −1 becomes shorter, there is less time to canalize trajectories prior to cycle 0, and hence the system will be free to occupy more states during this cycle. The results are shown in Fig. 4, where the numbers of attractors are plotted for 10 durations of cycle −1, ranging from 0.1 to 21.1 min, which is the actual duration of the cycle prototype in the real embryo. It is seen from the figure that the number of attractors rises only when cycle −1 is shorter than about 9 min, i.e. about half its natural length. Another interesting feature is that 57 attractors occur with cycle −1 set to the shortest duration of 0.1 min, compared to 116 stable steady states total in the model. About half of the steady states still lie in forbidden basins because of the constraint (2), since the division rule is used in all calculations. 5. Conclusions During the early development of multicellular organisms, many molecular processes, such as gene expression, occur in cells which undergo numerous divisions. Segmentation gene
Fig. 3. Two snapshots of a two-dimensional projection of the phase space in the same system as for the previous figures. The axes stand for the concentrations, in conventional units, of proteins from the two genes in the first nucleus (variables v11 and v12 ). In panel A, the snapshot is made at the end of cycle n = −1, right before the division. The dots in the panel represent the two-dimensional parts of solutions of the model at the end of cycle −1, with 104 different initial conditions taken at the beginning of cycle −1. Panel B corresponds to the snapshot at infinite time in the first nucleus, which is one of the two progenies of nucleus from panel A (the first nucleus in cycle −1 divides into nuclei 1 and 2 in cycle 0). Panel B shows 16 actual attractors (triangles) which the system goes to from the same initial conditions as for panel A. The panel also contains those of 116 stable steady states (daggers) which can never be reached by the system from any initial condition. Of the 16 triangles in B, 8 are hardly seen because they stay very close to the other 8. However, this does not mean that these two sets of attractors are identical: the triangles represent only the two-dimensional part of the attractors, and the two attractor sets are not similar in other dimensions. The same remark is valid for closely standing daggers of non-reachable stable steady states.
V.V. Gursky et al. / Physica D 218 (2006) 70–76
Fig. 4. Numbers of attractors in the same system as for previous figures, but with 10 different durations (in minutes) of cycle n = −1. The smallest duration is 0.1 min, and the largest is 21.1 min. Each number of attractors was calculated by using 105 different initial conditions.
expression in the early Drosophila embryo is an example of such a process. We have shown that nuclear divisions, which accompany gene expression in the system, play the dynamical role in this process. As each nucleus divides, adding new state variables to the system, the dimension of the state space increases from cycle to cycle. The complexity of the final phase space allows the model to have many stable steady states, which can be identified as the potential final expression patterns. However, because the system starts in a space of smaller dimension, only a limited number of these steady states are actually accessible to the system. Therefore, cell divisions provide a mechanism for the selection of possible attractors. We have identified two factors underlying the mechanism of selection and numerically illustrated their effect in a small two-gene system. The first factor that renders some basins of attraction inaccessible is provided by symmetric nuclear divisions. The second constraining factor is the dynamics itself at early times, which arranges the solution trajectories in narrow domains in the phase space, such that the attracting sets of some steady states are not captured by these domains. Inspection of Table 2 provides a rough estimate of the magnitude of the reduction in the number of attractors because of cell divisions. It is seen that the number of stable steady states grows nonlinearly and very rapidly with the final dimension of the phase space. At the same time, we find that the number Ad of attractors in the two-gene system with the chosen parameter values obeys the binary rule Ad = 2 Min , where d is the total number of divisions that occur in the model, and Min is the number of initial nuclei; this rule is valid only for d > 0. (The binary nature of this formula stems from the fact that the two genes have mutually exclusive states of expression because of mutual repression.) Thus, whatever the final dimension of the phase space, the number Ad is determined as soon as the numbers Min and d are chosen. This means that the relative fraction of reachable to total attractors in the model decreases rapidly with increasing d and/or Min . In particular, it follows from Table 2 that, for Min rising from 1 to 8, the fraction of attractors in the one-cleavage system, measured as a percentage of the number of all stable steady states, decreases from 50%
75
to, at maximum, 4%. Similarly, in the two-cleavage system with Min varying from 1 to 4 this fraction ranges from 17% to, at maximum, 0.3%. If we consider the two-gene systems finishing in the same phase space of dimension 32 (16 nuclei in the last cleavage cycle) but having different number d of nuclear divisions, then, for d varying from 1 to 3, the fraction decreases from, at maximum, 4% to, at maximum, 0.07%. It should be noted though that we do not have evidence that the described binary rule for the number of attractors, together with the precise percent calculated, can be quantitatively extended to more general systems. However, it is likely that the observed tendency of decrease in the portion of accessible attractors is a generic mechanism in the presence of approximately symmetric mitoses. Molecular transport between cells provides another mechanism for reducing the number of attractors in multi-cellular systems. In model (1), the transport of material takes place by the exchange of proteins between adjacent nuclei by diffusion. The effect of diffusion can be seen in the S-column of Table 2. Each entry in this column contains the number of stable steady states which are attractors in a system that has no divisions. As described above, there can be only 2 attracting states in one nucleus and, therefore, 2 Min of them if the multi-nuclear system is treated as the direct product of nuclei. The S-column of the table shows that this rule breaks down if the system has no divisions and more than 2 nuclei. The numbers of attractors are less than 2 Min in such cases, and this fact is a consequence of the diffusive coupling between nuclei. Acknowledgements The support of this work by NIH Grants RUB1-1578 and CRDF GAP award RBO1286, and by RFBR-NWO grant 047.011.2004.013 is gratefully acknowledged. References [1] R. Thomas, M. Kaufman, Multistationarity, the basis of cell differentiation and memory. II. Logical analysis of regulatory networks in terms of feedback circuits, Chaos 11 (1) (2001) 180–195. [2] R. Thomas, R. D’Ari, Biological Feedback, CRC, Boca Raton, FL, 1990. [3] S.A. Kauffman, Origins of Order: Self-Organization and Selection in Evolution, Oxford University Press, Oxford, 1993. [4] L. S´anchez, D. Thieffry, A logical analysis of the Drosophila gap-gene system, J. Theoret. Biol. 211 (2001) 115–141. [5] L. S´anchez, D. Thieffry, Segmenting the fly embryo: A logical analysis of the pair-rule cross-regulatory module, J. Theoret. Biol. 224 (2003) 517–537. [6] E. Mjolsness, D.H. Sharp, J. Reinitz, A connectionist model of development, J. Theoret. Biol. 152 (1991) 429–453. [7] J. Reinitz, D.H. Sharp, Mechanism of eve stripe formation, Mech. Dev. 49 (1995) 133–158. [8] J. Jaeger, S. Surkova, M. Blagov, H. Janssens, D. Kosman, K.N. Kozlov, Manu, E. Myasnikova, C.E. Vanario-Alonso, M. Samsonova, D.H. Sharp, J. Reinitz, Dynamic control of positional information in the early Drosophila embryo, Nature 430 (2004) 368–371. [9] J. Jaeger, M. Blagov, D. Kosman, K.N. Kozlov, Manu, E. Myasnikova, S. Surkova, C.E. Vanario-Alonso, M. Samsonova, O.H. Sharp, J. Reinitz, Dynamical analysis of regulatory interactions in the gap gene system of Drosophila melanogaster, Genetics 167 (2004) 1721–1737.
76
V.V. Gursky et al. / Physica D 218 (2006) 70–76
[10] P.A. Lawrence, The Making of a Fly, Blackwell Sci. Publ., Oxford, 1992. [11] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in C, Cambridge University Press, Cambridge, UK, 1992.
[12] E. Poustelnikova, A. Pisarev, M. Blagov, M. Samsonova, J. Reinitz, A database for management of gene expression data /in situ/, Bioinformatics 20 (2004) 2212–2221.