Critical Behavior and Axis Defining Symmetry Breaking in Hydra Embryonic Development

The formation of a hollow cellular sphere is often one of the first steps of multicellular embryonic development. In the case of Hydra, the sphere breaks its initial symmetry to form a foot-head axis. During this process a gene, ks1, is increasingly expressed in localized cell domains whose size distribution becomes scale-free at the axis-locking moment. We show that a physical model based solely on the production and exchange of ks1-promoting factors among neighboring cells robustly reproduces the scaling behavior as well as the experimentally observed spontaneous and temperature-directed symmetry breaking. Multicellular self-organization is a most intriguing phenomenon in the quest for understanding biological complexity and development [1]. One of the model organisms for the study of multicellular development is the freshwater polyp Hydra [2]. This $1 cm long animal has a mainly cylindrical shape bounded by a head (including tentacles and a mouth opening) at one end and a foot on the other. Hydra can regenerate lost body parts. The entire animal can even reform from a cell aggregate as small as 10 4 cells, made from dissociated and randomly condensed Hydra cells [3,4]. The reformation starts with the cells forming an inflating, hollow, spherical cell bilayer [5,6]. The cells then define a new foot-head axis, and further development takes place according to its direction. Axis locking occurs when a few cells start the genetic Hydra Wnt cascade, which is responsible for head formation [7]. In the adult, the Hydra specific ks1 gene is strongly expressed immediately before the cells terminally differentiate to become part of the head. Accordingly, ks1 has been described as a marker of ''cell head forming potential'' [8,9]. We have previously studied the expression of ks1 on the reforming cell sphere by the technique of in situ hybrid-ization [10]. We have observed that ks1 is expressed in fractal domains on the surface of the cell sphere. The domain size distribution exhibits a power-law behavior at the axis-locking moment. Our observations highlight a critical state of the Hydra molecular network consisting of collective fluctuations of gene expression patterns in coincidence with the axis-definition process [10,11]. The size of the fluctuations increases with time until a scale-free distribution of ks1 domain sizes is reached and axis-locking occurs. Critical dynamics have been observed in biological systems as diverse as ant colonies [12], gene expression [13], brain states [14], and neuronal networks [15], where it is seen as a …

Multicellular self-organization is a most intriguing phenomenon in the quest for understanding biological complexity and development [1].One of the model organisms for the study of multicellular development is the freshwater polyp Hydra [2].This $1 cm long animal has a mainly cylindrical shape bounded by a head (including tentacles and a mouth opening) at one end and a foot on the other.Hydra can regenerate lost body parts.The entire animal can even reform from a cell aggregate as small as 10 4 cells, made from dissociated and randomly condensed Hydra cells [3,4].The reformation starts with the cells forming an inflating, hollow, spherical cell bilayer [5,6].The cells then define a new foot-head axis, and further development takes place according to its direction.
Axis locking occurs when a few cells start the genetic Hydra Wnt cascade, which is responsible for head formation [7].In the adult, the Hydra specific ks1 gene is strongly expressed immediately before the cells terminally differentiate to become part of the head.Accordingly, ks1 has been described as a marker of ''cell head forming potential'' [8,9].
We have previously studied the expression of ks1 on the reforming cell sphere by the technique of in situ hybridization [10].We have observed that ks1 is expressed in fractal domains on the surface of the cell sphere.The domain size distribution exhibits a power-law behavior at the axis-locking moment.Our observations highlight a critical state of the Hydra molecular network consisting of collective fluctuations of gene expression patterns in coincidence with the axis-definition process [10,11].The size of the fluctuations increases with time until a scalefree distribution of ks1 domain sizes is reached and axislocking occurs.
Critical dynamics have been observed in biological systems as diverse as ant colonies [12], gene expression [13], brain states [14], and neuronal networks [15], where it is seen as a strategy to build up a coordinated response to slight stimuli or perturbations.
In this Letter we provide new data on Hydra head formation and we propose a simple model of its selforganization.While biochemical knowledge of the determinants of Hydra development is still largely incomplete [16] and many more details can be involved, here we consider a simplification which focuses on ks1.We take into account that intercellular communication is a ubiquitous phenomenon during embryogenesis [17] and that exchange of signaling molecules between neighboring cells, through e.g., gap junctions, is a prerequisite for Hydra patterning [18].A threshold mechanism in the production of ks1 can be indirectly assumed since in the experiments isolated, intermittent ks1 clusters are observed rather than a uniform, constant activation background.Thus, we assume that one or several ks1 (head) promoting agents are produced by the cells and spread from one cell to its neighbors when a threshold is reached and ks1 is expressed.In such a scenario, nearest-neighbor cell-cell coupling leads to a slow build up of correlations and, eventually, a large spatial fluctuation defines the head position and the symmetry breaking axis.
Model.-A minimal model incorporating such observations is (a) a ks1-promoting factor X is produced in each cell and (b) when a threshold X c is reached, ks1 is expressed and X is released to the surrounding cells, where further production or accumulation of X is induced.The value of X describes in abstract terms a molecular state of the cell that can be transferred to neighboring cells when a threshold is reached.
We model the Hydra ball as a spherical lattice of N ¼ 10 4 hexagonal sites representing Hydra cells [Fig.1(a)].Each cell is endowed with a quantity X of the ks1-promoting factor, which is initially randomly assigned at below-threshold dosage.The X factor is then produced at a rate .When the X level in the ith cell X i reaches the threshold value X c ¼ 1, a fraction C i of X i is released in equal parts to its nearest-neighbor cells, while X i in the discharging cell is reset to X i ¼ 0. We assume that the relaxation process is much faster than the time scale À1 ; i.e., we consider the !0 limit [19].The coupling between cells induces chain effects (avalanches) when stimulated neighboring cells also exceed the threshold value and propagate the discharge process.Interestingly, such a model can be mapped to a variant of a well-established model of self-organized criticality [19], originally proposed by Olami, Feder, and Christensen (OFC) for earthquake occurrence [20].
An important parameter in the model is the average conservation level C. Real cells are not identical, and differences in gene expression and metabolic activity exist: to account for these and similar stochastic effects, in the model we add a small, random cell-dependent term C i to the conservation level C, i.e., C i ¼ C þ C i .We find that for C i & 5% the results are insensitive to the noise level.Similarly, we have verified that our results are unchanged if a <5% noise component is added to the activation threshold or, equivalently, to the reset value (see Supplemental Material [21]).
Real Hydra is formed by a double layer of ectodermal and endodermal cells, but ks1 is expressed only in the ectodermal layer [8]; therefore, in our simulations we model the Hydra ball as a single-layered structure.We analyzed, however, the influence of a second layer and we observed that the overall model dynamics was essentially unchanged, with identical critical exponents.
Analysis of the model.-Wecarried out numerical simulations of the dynamics of the X factor in the model, to monitor avalanche spatial patterns corresponding to ks1 expression, and the conditions for the emergence of criticality.We explored the entire range C 2 ½0; 1 and found the best agreement with the experimental data for a conservation level C ' 0:95, treated as a reference below.The dependence of the critical behavior on C and noise level C is discussed in the Supplemental Material [21].Our numerical results are averaged over 100 random realizations.
At each time point we compute the probability distribution Pðs !Þ to observe a domain, i.e., an avalanche of X-discharging cells of size equal or larger than s (see [19]).As shown in Fig. 2, Pðs !Þ at short times falls off exponentially, as expected from the initial random distribution of X i across cells.However, as the dynamics evolves, a new regime emerges characterized by a powerlaw scaling behavior Pðs !Þ $ s À ( is the critical exponent of the cumulative distribution).The simulation time T, associated with the emergence of scale invariance in Pðs !Þ, corresponds in our case to N $ 2000 cascades.The general properties of the distribution Pðs !Þ are robust when C is varied, but the exponent , as much as T, depends on the conservation level C: consistently with the original OFC model, decreases with C (see Fig. 1 of the Supplemental Material [21]).We also recorded the fractal structure, described by the cluster area D cl and perimeter dimension D per of large avalanches (see [21]).In particular, for C ' 0:95 we found that ¼ 0:8 AE 0:05, and D cl ¼ 1:75 AE 0:05, D per ¼ 1:28 AE 0:05.Examples of avalanches in the model are shown in Fig. 1(b).At short times only small cascades of X-releasing cells are observed, which appear randomly located on the ball.As such a transient evolves in time, the typical spatial scale of domains increases, and after approximately a time T, a roughly steady regime is reached where domains of all sizes are found.
Comparison to experimental data.-Figure1(b) shows examples of experimental ks1 expression patterns (see [21] for details), corresponding to the early stages of regeneration (t ' 16 h), axis formation (t ' 25 h), and early head formation (t ' 40 h) in Hydra.The same figure also shows clusters of cells from our model simulations, in the case C ¼ 0:95, corresponding to three time points, T=2, T, 2T, where T is the time of appearance of criticality.Since in the model the duration of an avalanche is taken to be vanishingly small, only one avalanche per time frame is present by definition.By extending the persistence time of the ks1 activation cluster induced by each X avalanche, many a ks1 activation cluster could be seen simultaneously.It is worth noting here that in the experiments just a few small ks1 clusters could be seen at short times, while at longer times, i.e., when large clusters are present, most Hydra spheres contain only one ks1 domain.This observation supports the assumption that ks1 expression patterns have a comparatively short lifetime.
The shape and relative sizes of experimental and model domains look qualitatively similar.The model avalanche area and perimeter dimensions D cl ¼ 1:75 AE 0:05, D per ¼ 1:28 AE 0:05, agree within numerical errors with the experimental values D cl ' 1:7, D per ' 1:3 [10].The good agreement suggests to associate T with the time required for Hydra axis formation, i.e., about 25 h.In turn, that can be used to estimate the average production rate as ¼ X T =T, where X T ' 0:2X c is the average amount of X produced in each cell before time T, as measured in the simulations.That gives $ 10 À6 s À1 .With the plausible value X c $ 10 3 molecules [22], this is a rate X c $ 10 À3 molecule=s, compatible with average mRNA and protein production rates [22] and with our assumption of using !0.
In Fig. 2 we show the model distribution Pðs !Þ at the three time points, T=2, T, 2T, for the case C ¼ 0:95, along with the corresponding experimental distributions.Considering the level of simplifications of the model, the simulations reproduce well the experimental behavior: at criticality, a power-law decay is observed over two decades, with an identical value (within experimental errors) of the critical exponent ' 0:8 and similar cutoffs [23].The evolution in time of the experimental distribution is also reproduced, suggesting that the model is well capturing the general mechanism underlying the emergence of longrange correlations.The differences between the model at T=2 and the experiment at T ¼ 16 h may be due to time scale separation not holding in the experimental system at early times.
Temperature gradient.-Ourexperiments of Hydra regeneration in a temperature gradient [10,11,21] show that the metabolic rate difference between the cold and hot regions triggers a developmental asymmetry where heads form preferably at 0 or 180 with respect to the gradient direction [Fig.3(a)].In the case of a difference of ÁT ¼ 0:6 C, the head forms at the hot ( ' 0 ) or cold ( ' 180 ) sides with similar probability, while a stronger gradient ÁT ¼ 0:9 C results in a slight tendency for the head to form at the cold side.In the absence of the gradient or if the gradient is applied after the axis-locking moment, the axis direction is randomly oriented.The probability distribution functions for the angle formed between the direction of the gradient and the orientation of the center of mass of clusters of size larger than 20% of the system (i.e., s !20%N), from computer simulations.Four gradient amplitude Á are considered.For nonzero gradients, the angle distribution is bimodal with a peak at 0 and at 180 .The peak at 180 , i.e., in the direction opposite to the gradient, grows with Á, indicating that the largest clusters tend to be in antiphase with the gradient.To test the model against such a nontrivial observation, we introduced a spatial gradient Á in the X-factor production rate, (that is expected to have an Arrhenius-type dependence on temperature [24]).We set the angular coordinate with respect to the gradient direction to be ¼ 0 for the region of highest production rate, and ¼ 180 for the one with the lowest rate [Fig.3(b)].Head formation is associated with the appearance of ks1-expressing clusters of size comparable to the entire Hydra ball.Thus, we record the angular distribution of simulated clusters with a size s larger than 20% of the entire system (s > s 0 ¼ 20%N).We observe that the presence of a gradient substantially modifies the spatial distribution of the avalanches [Fig.3(b)]: avalanches larger than s 0 form preferably at 0 and at 180 .This bimodal distribution is in qualitative agreement with the experimental data [Fig.3(a)].More pronounced peaks correspond to larger gradients Á , with a higher peak at the cold side [Fig.3(b)], in ''antiphase'' with the gradient.A possible mechanism for such an asymmetry is the phase locking [25] induced by regions with slower firing dynamics (in the present case, regions with lower X production rates).This simple mechanism could explain the experimentally observed tendency of the head to emerge in antiphase with the gradient as the temperature difference increases [Fig.3(a)].
Considering the simplicity of the model the overall correspondence with experiments is surprising.Based on the insight provided by the model, we propose the following scenario for head formation in Hydra: single ks1-expressing cells have a small probability to trigger Wnt expression; as development advances, larger ks1 domains of cells will appear and, accordingly, the probability to observe a Wnt head spot will also increase.This way, a direct correlation emerges between the location of head formation and large ks1 domains.In such a scenario, a Poissonian probability exists to observe more than a single head, in agreement with observations performed on large Hydra aggregates [7,10], where several heads can appear.
Discussion.-Our simple model illustrates, from minimal assumptions, how a local coupling between production and exchange of ks1-promoting factors can naturally induce the scaling distribution of ks1-expressing domains observed during Hydra development, and the peculiar response of the Hydra ball to temperature gradients.The rationale behind the model is that the ks1 phenomenology in Hydra embryogenesis could derive from purely local exchange of intercellular signals in a tissue of otherwise identical cells, triggering large activation avalanches, similar to those observed in OFC systems.This self-organized scenario is particularly appealing since, despite a long search, long-range signaling mechanisms have not been identified in Hydra development, and the precise nature of the ks1-promoting factor X is still unclear.
Our results are unchanged in the presence of & 5% noise in the X factor conservation level, activation threshold, or reset value.Higher levels of noise increasingly destroy spatial correlations (cf. the Supplemental Material [21] and Refs.[25][26][27][28]); therefore, noise in Hydra must be controlled accordingly.Indeed, in biological development several mechanisms of noise control exist, which could make the system adequately resilient to fluctuations [29].
Self-organized critical behaviors have been observed in a variety of physical systems [19].They entail the development of large scale correlations through nearest-neighbor interactions, allowing coordinated macroscopic functions, efficient transmission of information, and enhanced sensitivity to external stimuli [30,31] without the need for any external tuning.This might have been a crucial mechanism in the transition from a disorganized set of unicellular organisms to an organized multicellular community.

FIG. 1 .
FIG.1.Avalanches of ks1 activation in Hydra.(a) In the model, all cells are initially assigned a random value X i < X c of the ks1-promoting factor (the darker the cell, the higher the concentration of X).At each time step, the concentration of X in each cell increases as X þ until the threshold level X c (black) is reached.The ith cell then discharges a fraction c i ¼ C i =6 of X i to the neighboring cells and is subsequently reset to zero (white).Avalanche-like dynamics occurs when neighboring cells also reach the threshold level and discharge.(b) Comparison of ks1 expression patterns (black spots) from simulations and experiments, at three different developmental times.For simulations, T is the time of appearance of criticality.

FIG. 3 (
FIG. 3 (color online).Symmetry breaking in a temperature gradient.(a) Experimental results for the distribution of head orientations .Three conditions are shown: ÁT ¼ 0 (dotted line), ÁT ¼ 0:6 C (solid green line), and ÁT 0:9 C (dashed red line).(b) The probability distribution functions for the angle formed between the direction of the gradient and the orientation of the center of mass of clusters of size larger than 20% of the system (i.e., s !20%N), from computer simulations.Four gradient amplitude Á are considered.For nonzero gradients, the angle distribution is bimodal with a peak at 0 and at 180 .The peak at 180 , i.e., in the direction opposite to the gradient, grows with Á, indicating that the largest clusters tend to be in antiphase with the gradient.