## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# Morphogenesis at criticality

Contributed by William Bialek, January 17, 2014 (sent for review September 28, 2013)

## Significance

Biological networks are described by many parameters, and the behavior of a network is qualitatively different (monostable, bistable, oscillating, etc.) in different parts of parameter space. Critical points and surfaces are the borders between such qualitatively different regimes, as with phase transitions in equilibrium thermodynamics. We argue that, as expected from the thermodynamic case, genetic regulatory networks should exhibit behaviors near criticality that are independent of most molecular details. Analyzing recent experiments on the gap gene network in the early *Drosophila* embryo, we find that these signatures of criticality can be seen, quantitatively. This raises the question of why evolution has tuned this network to such a special point in its parameter space.

## Abstract

Spatial patterns in the early fruit fly embryo emerge from a network of interactions among transcription factors, the gap genes, driven by maternal inputs. Such networks can exhibit many qualitatively different behaviors, separated by critical surfaces. At criticality, we should observe strong correlations in the fluctuations of different genes around their mean expression levels, a slowing of the dynamics along some but not all directions in the space of possible expression levels, correlations of expression fluctuations over long distances in the embryo, and departures from a Gaussian distribution of these fluctuations. Analysis of recent experiments on the gap gene network shows that all these signatures are observed, and that the different signatures are related in ways predicted by theory. Although there might be other explanations for these individual phenomena, the confluence of evidence suggests that this genetic network is tuned to criticality.

Genetic regulatory networks are described by many parameters: the rate constants for binding and unbinding of transcription factors to their target sites along the genome, the interactions between these binding events and the rate of transcription, the lifetimes of mRNA and protein molecules, and more. Even with just two genes, each encoding a transcription factor that represses the other, changing parameters allows for several qualitatively different behaviors (1). Strong mutual repression between two genes can lead to bistability, or switch-like behavior, such that the system can be stable in both “on–off” and “off–on” states (2). On the other hand, if interactions are weak, the two interacting genes have just one stable state, and the expression levels in this state are controlled primarily by the inputs. Importantly, if we imagine smooth changes in the strength of the repressive interactions, the transition from graded response to switch-like behavior is not smooth: the behavior is qualitatively different depending on whether the relevant interactions are stronger or weaker than a critical value. Here we explore the possibility that the gap gene network in the *Drosophila* embryo might be tuned to such a critical point.

Early events in the fruit fly embryo provide an experimentally accessible example of many questions about genetic networks (3⇓–5). Along the anterior–posterior axis, for example, information about the position of nuclei flows from primary maternal morphogens to the gap genes, shown in Fig. 1 (6⇓–8), to the pair rule and segment polarity genes. Although the structure of the gap gene network is not completely known, there is considerable evidence that the transcription factors encoded by these genes are mutually repressive (6, 9⇓–11). If we focus on a small region near the midpoint of the embryo (near ), then just two gap genes, hunchback (Hb) and Krüppel (Kr), are expressed at significant levels, and this is repeated at a succession of crossing points or expression boundaries: Hb–Kr, Kr–Kni (knirps; ), Kni–Gt (giant; ), and Gt–Hb , as we move from anterior to posterior. In each crossing region, it is plausible that the dynamics of the network is dominated by the interactions among just the pair of genes whose expression levels are crossing.

We argue that criticality in a system of two mutually repressive genes generates several clear, experimentally observable signatures. First, there should be nearly perfect anticorrelations between the fluctuations in the two expression levels. As a result, there are two linear combinations of the expression levels, or “modes,” that have very different variances. Second, fluctuations in the large variance mode should have a significantly non-Gaussian distribution, whereas the small variance mode is nearly Gaussian. Third, there should be a dramatic slowing down of the dynamics along one direction in the space of possible expression levels. Finally, there should be correlations among fluctuations at distant points in the embryo. These signatures are related: The small variance mode will be the direction of fast dynamics, and under some conditions the large variance mode will be the direction of slow dynamics; the fast fluctuations should be nearly Gaussian, whereas the slow modes are non-Gaussian; and only the slow mode should exhibit long-ranged spatial correlations. We will see that all of these effects are found in the gap gene network.

We begin by analyzing a two-gene network to make precise what we mean by criticality. Crucially, this analysis shows that there are quantitative, experimentally observable signatures of criticality that do not depend on the molecular details of the network. We then analyze experiments on the gap gene network, searching for these experimental signatures.

## Criticality in a Two-Gene Network

One of the most important results in statistical mechanics and dynamical systems theory is that the behavior of systems near a critical point is classifiable, and independent of many details (12, 13). In this spirit, we want to identify behaviors of a genetic network that would emerge near criticality. We focus on a network of two genes, because this is the simplest case and because, as emphasized above, the gap gene network is dominated by two genes at each of the crossing regions identified in Fig. 1. The analysis that follows is rather standard, but we make it explicit to clarify how our predictions follow from the idea of criticality itself, rather than depending on the many (largely unknown) molecular parameters.

We consider a broad class of models for a genetic regulatory circuit in which the rate at which gene products are synthesized depends on the concentration of all of the relevant transcription factors, and the gene products are degraded. To simplify, we ignore delays, so that the rate at which the protein encoded by a gene is synthesized depends instantaneously on the other protein (transcription factor) concentrations, and we assume that degradation obeys first-order kinetics. We also focus on a single cell, leaving aside (for the moment) the role of diffusion. Then, by choosing our units correctly we can write the dynamics for the expression levels of two interacting genes as

where and are the normalized expression levels of the two genes, and are the lifetimes of the proteins, and *c* represents the external (maternal) inputs. The functions and are the “regulation functions” that express how the transcriptional activity of each gene depends on the expression level of all of the other genes; with our choice of units, the regulation function runs between zero (gene off) and 1 (full induction). All of the molecular details of transcriptional regulation are hidden in the form of these regulation functions (14), which we will not need to specify. Finally, the random functions and model the effects of noise in the system.

If the interactions are weak, then for any value of the external inputs *c* there is a single steady-state response, defined by expression levels and . We can check whether this hypothesis is consistent by asking what happens to small changes in the expression levels around this steady state. We write , and similarly for , and then expand Eqs. **1** and **2** assuming that and are small. The result isHere and are effective decay rates for the two proteins, which must be positive if the steady state we have identified is stable. The parameter reflects the incremental effect of gene 2 on gene 1: means that the protein encoded by gene 2 is a repressor of gene 1, whereas means that the protein encoded by gene 2 is an activator of gene 1, and similarly for . The noise terms and play the same role as and , but have different normalization.

If the steady state that we have identified is stable, then the matrixmust have two eigenvalues with negative real parts. This is guaranteed if the interactions are weak , but as the interactions become stronger it is possible for one of the eigenvalues to vanish. This is the critical point. Notice that we can define the critical point without giving a microscopic description of all of the interactions that determine the form of the regulation functions.

The linearized Eq. **3** predicts that the relaxation of average expression levels to their steady states can be written as combinations of two exponential decays,where and are the slow and fast eigenvalues of . Thus, although we measure the two expression levels, there are combinations of these expression levels—directions (, ) in the plane—that provide more natural coordinates for the dynamics, such that motion along each direction is a single exponential function of time. As we approach criticality, the dynamics along the slow direction becomes very slow, so that .

The linearized Eq. **3** also predicts the fluctuations around the steady state. As we approach criticality, things simplify, and we find the covariance matrixwhere is the variance in the expression level of the first gene. As with the dynamics, there are two “natural” directions in the plane corresponding to eigenvectors of this covariance matrix (principal components). In this linear approximation, the critical point is the point where we “lose” one of the dimensions, and the fluctuations in the two expression levels become perfectly correlated or anticorrelated. In addition, the direction with small fluctuations is the direction of fast relaxation.

Eqs. **1** and **2** refer to a single cell, or, in the case of the syncitial *Drosophila* embryo, to a single nucleus. If we allow that the proteins encoded by genes 1 and 2 can move among neighboring nuclei, we should write coupled equations for the dynamics in these neighbors. A simple approximation is to imagine that nuclei are sufficiently dense that a plot of protein concentration vs. nuclear position becomes continuous, as with the profiles in Fig. 1. If the dynamics of molecular motion along the length of the embryo is diffusive, then the natural generalization of Eqs. **1** and **2** becomes

where for simplicity we assume that the diffusion constants of the two species are the same. These equations have the same qualitative form as those considered by Turing, but we consider a system that does not cross the transition to spontaneous pattern formation considered in Turing’s work (15).

The linearized version of Eqs. **7** and **8**, analogous to Eq. **3**, iswhere is the same matrix as in Eq. **4**. In general this matrix depends on position *x*, and it is not easy to exhibit the solution of these equations in closed form. However, if we focus on small regions of the embryo, such as the crossings where we expect our two-gene model is most useful, we can make the approximation that is constant. At the end of our discussion, we will try to stitch together these local approximations into a more global view.

If the noise sources driving the randomness in gene expression arise from local processes, as seems plausible, then the “forces” should be uncorrelated at different positions, although the output expression levels can become correlated via diffusion. These correlations again are expressed most naturally in terms of the slow and fast directions in the plane. Concretely, we can solve Eq. **9** in the approximation of constant to givewhere the “correlation length” ; similarly, for the fast mode we have . Thus, connected to the slowing of the dynamics at criticality , we should see spatial correlations extending over long distances .

To summarize, the approach to criticality generates slow dynamics ( in Eq. **5**), strong local (anti-)correlations (Eq. **6**), and long-ranged correlations in space ( in Eq. **10**). These predictions are connected, so that the anticorrelated combination of expression levels is the mode that exhibits long-ranged correlations, and this mode is plausibly the same mode that has slow dynamics. As we will see below, there is one more somewhat subtle prediction, namely that the distribution of fluctuations along the slow mode will be non-Gaussian.

## Local Correlations

Testing the predictions of criticality requires measuring the time dependence of gap gene expression levels, with an accuracy better than the intrinsic noise levels of the system. Absent live movies of the expression levels, the progress of cellularization provides a clock that can be used to mark the time during nuclear cycle 14 at which an embryo was fixed (16), accurate to within 1 min (7). Fixed embryos, with immunofluorescent staining of the relevant proteins, thus provide a sequence of snapshots that can be placed accurately along the time axis of development. Immunofluorescent staining itself provides a measurement of relative protein concentrations that is accurate to within of the maximum expression levels in the embryo (7).

In Fig. 2*A* we show the correlations between fluctuations in pairs of gap genes at each position. Gap gene expression levels plateau at into nuclear cycle 14 (7), and the mean expression levels are shown as a function of anterior–posterior position in Fig. 1. At each position we can look across the many embryos in our sample, and analyze the fluctuations around the mean, as in ref. 8. We see that, precisely in the “crossing region” where Hb and Kr are the only genes with significant expression (marked “A” in Fig. 2*A*), the correlation coefficient approaches , perfect (anti-)correlation, as expected at criticality. This pattern repeats at the crossing between Kr and Kni (B), at the crossing between Kni and Gt (C), and, perhaps less perfectly, at the crossing between Gt and Hb (D). [Because we observe substantial anticorrelations both in the Gt–Hb pair and in the Hb–Kni pair, it is likely that the system is more nearly 3D in the neighborhood of crossing D, so that no single pair can achieve perfect correlation.]

These strong anticorrelations are shown explicitly in Fig. 2*B*, where we plot the two relevant gene expression levels against one another at each crossing point. In all cases, the direction of small fluctuations is along the positive diagonal, whereas the large fluctuations are along the negative diagonal.

It is important that the strong anticorrelations tell us something about the underlying network, rather than being a necessary (perhaps even artifactual) corollary of the mean expression profiles. A notable feature of Fig. 2*A* thus is what happens away from the major crossing points. There is a Hb–Kni crossing at (E), but this does not have any signature in the correlations, perhaps because spatial variations in expression levels at this point are dominated by maternal inputs rather than being intrinsic to the gap gene network (17, 18). This is evidence that we can have crossings without correlations, and we can also have correlations without crossings, as with Hb and Kr at point H; interestingly, H marks the point where an additional posterior Kr stripe appears during gastrulation (19, 20). We also note that strong correlations can appear when expression levels are very small, as with Hb and Kni at points F and G; there also are extended regions of positive Kr–Kni and Kni–Gt correlations in parts of the embryo where the expression levels of Kr and Kni both are very low. Taken together, these data indicate that the pattern of correlations is not simply a reflection of the mean spatial profiles, but an independent measure of network behavior.

It should be admitted that strong anticorrelations can be explained without criticality. As an example, if *A* represses *B*, and the dominant source of noise in the expression of *B* is the fluctuating concentration of *A*, then *A* and *B* will be perfectly anticorrelated. Alternatively, if both *A* and *B* respond to a common input (e.g., as the gap genes respond to maternal morphogens), then if the dominant source of noise in the output is the response to variations in this input, again one can have strong (anti-)correlations. Although we could try to eliminate particular versions of these models as explanations for what we see in the gap genes, it seems more useful to note that none of these models has a natural way of generating the other signatures of criticality, to which we now turn.

## Non-Gaussianity

If we transform Eq. **3** to a description in terms of the fast and slow modes and , then precisely at criticality there is no “restoring force” for fluctuations in and formally the variance in Eq. **6** should diverge. This divergence is stopped by higher-order terms in the expansion of the regulation functions around the steady state, and this leads to a non-Gaussian distribution of fluctuations in . Although the data are limited, we do find, as shown in Fig. 2*C*, that fluctuations in the small variance (fast) direction are almost perfectly Gaussian, whereas the large variance (slow) direction shows significant departures from Gaussianity, in the expected direction.

## Dynamics

The time dependence of Hb and Kr expression levels during nuclear cycle 14 is shown, at the crossing point , in Fig. 3. Criticality predicts that if we take a linear combination of these expression levels corresponding to the direction of small fluctuations in Fig. 2*B* (cyan), then we will see relatively fast dynamics, and this is what we observe. In contrast, if we project onto the direction of large fluctuations (magenta), we see only very slow variations over nearly 1 h. Indeed, the expression level along this slow direction seems almost to diffuse freely, with growing variance rather than systematic evolution. Thus, strong (anti-)correlations are accompanied by a dramatic slowing of the dynamics along one direction in the space of possible expression levels, and a similar pattern is found at each of the crossings, Kr–Kni, Kni–Gt, and Gt–Hb. Again, this is consistent with what we expect for two-gene systems at criticality.

If we move along the anterior–posterior axis in the vicinity of the crossing point, the sum of expression levels of two genes, which is proportional to the fast mode, remains approximately constant, whereas the difference, which is proportional to the slow mode, changes. Therefore, the dynamics of the slow mode, shown in Fig. 3, will generate motion of the pattern along the anterior–posterior axis. This slow shift is well known (7, 21).

## Spatial Correlations

The slow dynamics associated with criticality implies that correlations should extend over long distance in space. In particular, as we approach criticality, the relaxation rate vanishes and the associated correlation length in Eq. **10** becomes infinitely long, in practice being limited by the size *L* of the embryo itself. Searching for these long-ranged correlations is complicated by the fact that the system is inhomogeneous, but we have a built-in control, because we should see the long-ranged correlations only in the slow, large variance mode , and not in . This also helps us discriminate against systematic errors that might have generated spurious correlations.

In Fig. 4*A* we show the normalized correlation functionwith *x* held fixed at the Hb–Kr crossing and *y* allowed to vary. We see that this correlation function is essentially constant throughout the crossing region. In contrast, the same correlation computed for the fast mode decays rapidly, with a length constant , just a few nuclear spacings along the anterior–posterior axis.

The dominant slow mode corresponds to different combinations of expression levels in different regions of the embryo. Generally, we can write the slow mode as a weighted sum of the different expression levels,where we label the gap genes for Hb, for Kr, for Kni, and for Gt. Near the Hb–Kr crossing, labeled “A” in Fig. 2*A*, we have , where are the weights that give us the antisymmetric combination of Hb and Kr, as drawn in Fig. 2*B*: , , and . Similarly, near the Kr–Kni crossing, labeled “B” in Fig. 2*A*, we have , with , , and , and this generalizes to crossing regions C and D. Using these weights, we obtain approximations to the slow mode,and we expect that these approximations are accurate in their respective crossing regions. Now we can test for correlations over longer distances by computing, for example,holding in the crossing region *A* while letting *y* vary through the crossing region *B*, and similarly for and . The results of this analysis are shown in Fig. 4*B*; note that is the correlation we have plotted in Fig. 4*A*.

Fig. 4*B* shows that the slow mode is correlated over very long distances. We can see, for example, in , correlations between fluctuations in expression level at the Hb–Kr crossing region and at the Kni–Gt crossing region, despite the fact that these regions are separated by of the length of the embryo and have no significantly expressed genes in common. These peaks in the correlation functions appear also at points anterior to the crossing regions, presumably at places where our approximations in Eq. **13** come close to some underlying slow mode in the network. The pattern of correlations has an envelope corresponding to an exponential decay with correlation length (dashed line in Fig. 4*B*), and similar results are obtained for the correlation functions , , , etc.. This means that fluctuations in expression level are correlated along essentially the entire length of the embryo, as expected at criticality.

## Discussion

To summarize, the patterns of gap gene expression in the early *Drosophila* embryo exhibit several signatures of criticality: near-perfect anticorrelations of fluctuations in the expression levels of different genes at the same point, non-Gaussian distributions of the fluctuations in the large variance modes, slowing down of the dynamics of these modes, and spatial correlations of the slow modes that extend over a large fraction of the embryo. Although each of these observations could have other explanations, the confluence of results strikes us as highly suggestive. Note that we have focused on aspects of the data that are connected to the hypothesis of criticality in a very general way, independent of other assumptions, rather than trying to build a model for the entire network or fit parameters.

The possibility that biological systems might be poised near critical points, often discussed in the past (22), has been reinvigorated by new data and analyses on systems ranging from ensembles of amino acid sequences to networks of neurons to flocks of birds (23). Although we would like to avoid too much speculation, it does seem natural to ask how being tuned to a critical point could be functionally advantageous for the embryo.

Criticality may have a special meaning in the context of pattern formation. Usually, if we have a network of biochemical reactions and diffusion, there are length scales set by combinations of rate constants and diffusion constants, and these length scales set the size of features in the pattern. In contrast, the patterns formed in embryonic development are thought to scale in proportion to the size of the embryo (24, 25). One way of expressing the problem of scaling is that moving the boundary set by the posterior end of the egg must have an effect on the patterns of gene expression even near the anterior end, whereas conventional models predict that perturbations will penetrate only a short distance ξ from the point where they are applied, comparable to the distance over which fluctuations in expression are correlated. At criticality, however, this correlation length becomes the size of the embryo itself, as in Fig. 4. Whereas this does not provide an explicit mechanism for scaling, it may provide a different point of view on the problem.

At a critical point, fluctuations are large, and it might seem that this is a disadvantage. However, criticality also implies a very large gain in response to small changes in the input signals (e.g., the maternal morphogens), and the long time scales can serve to average out components of the noise in the system. Analysis of simple models shows that there is a regime in which transcriptional regulatory systems can maximize the flow of information from inputs to outputs by operating near a critical point (26). There is evidence that noise levels in the gap gene system are nearly as low as possible given the finite concentrations of the relevant molecules and the time available for decisions to be made (27), and that the representation of positional information by the gap genes is optimized by matching the distribution of input signals to the noise characteristics of the network (8, 28). Tuning to a critical point may provide another layer of optimization along the same functional axis. If correct, this would place criticality in the gap gene network alongside other examples of optimizing signal-to-noise or information flow in biological systems, ranging from molecule counting in bacterial chemotaxis (29) to photon counting in human vision (30), among many other examples (25).

Even leaving aside the possibility of criticality, the aspects of the data that we have described here are not at all what we would see if the gap gene network is described by generic parameter values. There must be something about the system that is finely tuned to generate such large differences in the time scales for variation along different dimensions in the space of expression levels, or to ensure that correlations are so nearly perfect and extend over such long distances.

## Acknowledgments

We thank S. Blythe, K. Doubrovinski, O. Grimm, J. J. Hopfield, S. Little, M. Osterfield, A. M. Polyakov, M. Tikhonov, G. Tkačik, and A. Walczak for helpful discussions. We are especially grateful to E. F. Wieschaus for discussions both of the ideas and their presentation. This work was supported in part by National Science Foundation Grants PHY-0957573 and CCF-0939370, by National Institutes of Health Grants P50 GM071508 and R01 GM097275, by the Howard Hughes Medical Institute, by the W. M. Keck Foundation, and by Searle Scholar Award 10-SSP-274 (to T.G.).

## Footnotes

- ↵
^{1}To whom correspondence may be addressed. E-mail: wbialek{at}princeton.edu or dkrotov{at}princeton.edu.

Author contributions: This is a collaboration between theorists (D.K. and W.B.) and experimentalists (J.O.D. and T.G.), involving the analysis of previously published data. All authors contributed to all aspects of the paper, to the best of their abilities.

The authors declare no conflict of interest.

Freely available online through the PNAS open access option.

## References

- ↵
- ↵
- ↵
- ↵
- Lawrence PA

- ↵
- Gerhardt J,
- Kirschner M

- ↵
- ↵
- ↵
- Dubuis JO,
- Tkačik G,
- Wieschaus EF,
- Gregor T,
- Bialek W

- ↵
- ↵
- ↵
- ↵
- Guckenheimer J,
- Holmes P

- ↵
- Parisi G

- ↵
- ↵
- Turing AM

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Bak P

- ↵
- Mora T,
- Bialek W

- ↵
- ↵
- Bialek W

- ↵
- ↵
- ↵
- Tkačik G,
- Callan CG Jr.,
- Bialek W

- ↵
- ↵

## Citation Manager Formats

## Sign up for Article Alerts

## Article Classifications

- Physical Sciences
- Physics

- Biological Sciences
- Biophysics and Computational Biology