Canalization of Gene Expression and Domain Shifts in the Drosophila Blastoderm by Dynamical Attractors
Canalization refers to the constancy of the wild type phenotype under varying developmental conditions 1–4 ., In order to explain canalization , C . H . Waddington hypothesized that there must only be a finite number of distinct developmental trajectories possible , since cells make discrete fate decisions , and that each such trajectory , called a chreod , must be stable against small perturbations 5 ., One aspect of canalization , the buffering of phenotypic variability against genotypic variability in wild type , has received considerable experimental 2 , 6–10 and theoretical 11–13 attention ., The phenomenon of canalization of genotypic and environmental variation was seen by Waddington as a consequence of the underlying stability of developmental trajectories , an idea supported by theoretical analysis 13 ., But this central idea of Waddingtons has heretofore received little attention in real developmental systems because of a lack of relevant quantitative molecular data ., The further investigation of Waddingtons hypothesis is of great importance because it provides a scientific connection between the reliability and invariance of the formation of cell types and tissues in the face of underlying molecular variability , as we now explain ., Quantitative molecular data permitting the study of developmental canalization are now available for the segment determination process in Drosophila 14 ., The segmented body plan of the fruit fly Drosophila melanogaster is determined when the embryo is a blastoderm 15 by the segmentation genes 16 ., Quantitative spatiotemporal gene expression data show that the maternal protein gradients and the early expression patterns of the zygotic gap and pair-rule genes vary a great deal from embryo to embryo 14 , 17 ., The variation of the expression patterns of the gap and pair-rule genes decreases over time so that it is significantly lowered by the onset of gastrulation at the end of cellularization ( 14 , Fig . 1 ) ., The observed reduction of variability over time in the segmentation gene system suggests that the developmental trajectory of the Drosophila embryo is stable against perturbation ., The characterization of the stability properties of the developmental trajectory is central to our understanding of the mechanisms that underlie canalization 3 ., In the case of the gap genes , we have shown elsewhere 18 that variation reduction relative to the maternal gradient Bicoid ( Bcd ) occurs because of gap gene cross regulation ., Using a gene circuit model of the gap gene network 18–22 we identified specific regulatory interactions responsible for variation reduction in silico and verified their role in canalization experimentally ., Importantly , the model reproduces the observed low variation of gap gene expression patterns 18 , which provides an opportunity to analyze the properties of the system that give rise to developmental stability ., These results raise two generic problems that occur in the analysis of complex numerical models ., First , even if the model describes a natural phenomenon faithfully , understanding the natural phenomenon is only achieved when the models behavior can be understood as well ., The complexity of the model , unsurprising in terms of the underlying complexity of the biological system itself , poses a significant challenge to understanding model function ., Second , any model is an approximation to the actual mechanisms operating in an organism ., The models behavior must be robust to perturbation , since organisms develop and function reliably even though the underlying mechanisms are subject to a wide variety of perturbations and stresses ., There is extensive molecular variability among cells and embryos ( 14 , 17 , 23–26; reviewed in 27 ) and yet there is functional identity between equivalent cell types or conspecific individuals ., René Thom tried to resolve this apparent contradiction between the constancy of biological function and the variability in biological substructure by proposing a qualitative topological view of the trajectories of dynamical models 28 ., The term “topology” is used here to refer to properties of developmental trajectories that are invariant under continuous deformation ., The preservation of these properties ensures the robustness of model behavior , while a qualitative view often leads to an intuitive understanding of complex mechanisms ., One such robust property is an attractor state , or a stable steady state of a dynamical system , that attracts all trajectories in some neighborhood of itself ., Attractor states are locally stable under small perturbations of the dynamical model 29 , and for this reason it has been proposed that cell fates are attractors 11 , 13 , 30–33 ., The presence of an attractor state in the phase space of a system implies that there exists a region of phase space , called the basin of attraction , in which all trajectories approach the attractor asymptotically 34 , 35 ., This suggests that an attractor is the kind of qualitative robust property that could explain the stability of trajectories , and hence canalization ., There are , however , three important considerations to keep in mind when using attractors to describe the Drosophila blastoderm ., First , the reduction of variation due to attractors is only guaranteed at late times , but the reduction in the variation of the gap gene expression patterns takes place over about 100 minutes prior to gastrulation ., The reduction of variation before gastrulation is biologically essential as the expression patterns of engrailed and wingless , which form the segmentation prepattern , have a resolution of one nucleus and are created by the precise overlap of pair-rule and gap domains 14 , 36 ., Furthermore , at about the time of gastrulation the embryo undergoes the midblastula transition 37 , 38 at which time a qualitative change occurs in the genetic control of the embryo ., Second , in general there can be more than one attractor in the phase space 39–43 ., Thus , the basins of attraction need to correspond to biological initial conditions and be large enough to ensure robustness ., Finally , the set of attractors found must succeed not only in explaining canalization but also the morphogenetic properties of the system ., One such property is the anterior shift of gap gene domains located in the posterior region 14 , 21 , 44 ., These shifts are biologically significant and are difficult to reconcile with stable point attractors ., In this paper we show that the variation reduction of gap gene expression patterns is a consequence of the action of robust attracting states ., We further show that the complex patterning of the gap gene system reduces to the three qualitative dynamical mechanisms of ( 1 ) movement of attractors , ( 2 ) selection of attractors , and ( 3 ) selection of states on a one dimensional manifold ., The last of the three mechanisms also causes the domain shifts of the gap genes , providing a simple geometric explanation of a transient phenomenon ., In the Gap Gene Circuits section we briefly describe the gene circuit model; see 18 for a full description ., For each nucleus in the modeled anteroposterior ( A–P ) region , we identified the attractors in the gap gene phase space , calculated the trajectories , the basins of attraction and other invariant sets such as one dimensional attracting manifolds ( Stability Analysis of the Trajectories of the Gap Gene System section ) ., The stability of the trajectories was tested by varying the initial conditions within a biological range , based on gene expression data , that represents the variability of early gap gene expression ., We plotted the attractors and several trajectories corresponding to different initial conditions to make phase portraits that show the global qualitative behavior of the system ., Finally , we studied how the phase portraits changed as A–P position was varied to infer qualitative pattern formation mechanisms ., The biological conclusions about canalization and pattern formation arising from the dynamical characterization are presented in the Mechanisms of Canalization and Pattern Formation section ., The gene circuit used in this study models the spatiotemporal dynamics of the protein expression of the gap genes hunchback ( hb ) , Krüppel ( Kr ) , giant ( gt ) , and knirps ( kni ) during the last two cleavage cycles ( 13 and 14A ) before gastrulation 37 in the Drosophila blastoderm ., The protein products of these genes localize to nuclei 45–48 so that the state variables are the concentrations of the proteins in a one dimensional row of nuclei along the A–P axis of the blastoderm ., The concentration of the protein in the nucleus at time is denoted by ., In the model we considered a region , from 35% to 92% egg lenth ( EL ) along the A–P axis , which corresponds approximately to the region of the blastoderm fated to form the segmented part of the adult body 49 , 50 ., The gap genes are expressed in broad domains ( Fig . 2A , B; 14 ) under the control of maternal cues ., The anterior maternal system acts primarily through the protein gradient Bcd 51–53 which is essentially stationary and has an exponential profile ( Fig . 2C; 14 , 51 , 54 ) during the modeled time period ., The posterior maternal system is represented by the maternal Hb gradient ( Fig . 2C; 55–57 ) ., The terminal system regulates gap gene expression by activating tailless ( tll ) and huckebein ( hkb ) 58–61 ., The terminal system is represented in the model by the Tll gradient , which is expressed posterior to 80% EL in the modeled region during cycles 13 and 14 ( 14 and Fig . S1B ) ., tll is considered upstream of the gap genes since its expression pattern is unchanged in gap gene mutants 62 ., The concentration of Bcd in nucleus is denoted by and was determined using Bcd data from a representative cycle 13 embryo by an exponential fit , so that ( see 18 for details ) ., The concentrations of Tll and another upstream regulator , Caudal ( Cad ) 63 , 64 , were determined by interpolating average data in time 18 ., The concentrations of Tll and Cad are denoted by respectively , with an explicit dependence on time , since these gradients are not stationary ( Fig . S1 ) ., The dynamical equations governing are given by ( 1 ) where in a gene circuit with genes and nuclei ., The first term on the right hand side of Eq ., ( 1 ) represents protein synthesis , the second one represents protein transport through Fickian diffusion and the last term represents first-order protein degradation ., The diffusion constant , varies inversely with squared internuclear distance , and is the degradation rate ., The synthesis term is set to zero during the mitosis preceding the thirteenth nuclear division as synthesis shuts down 65 ., Following this mitosis , the nuclei are divided and daughter nuclei are given the same state as the mother nucleus ., is the maximum synthesis rate , and is a sigmoidal regulation-expression function ., The first term in the argument of represents the transcriptional cross regulation between the gap genes and the genetic interconnectivity is specified by the matrix ., Positive elements of imply activation while negative ones imply repression ., The regulation of the gap genes by Bcd is represented in the second term and is the regulatory strength ., The regulation of the gap genes by upstream time-varying inputs is represented in the third term and is the number of such inputs ., There are two such inputs in this model , Cad and Tll , and the elements of the matrix have the same meaning as those of ., The last term , , represents the effect of ubiquitous transcription factors and sets the threshold of activation ., The initial conditions for Hb are specified using cleavage cycle 12 data ., Cycle 12 data are a good approximation to the maternal Hb gradient since the zygotic expression of hb appears to begin in cleavage cycle 13 17 ., The initial conditions for Kr , Gt , and Kni are taken to be zero , since their protein expression is first detected in cycle 13 14 , 61 , 66–68 ., The gene circuits parameters were determined by performing a least-squares fit to a time series of averaged gap gene data 14 using the Parallel Lam Simulated Annealing algorithm ( see Methods ) ., This time series has nine points ( time classes; see Table S1 ) , one in cycle 13 and the rest in cycle 14A ., The output of the gene circuit ( Fig . 2E , F ) fits the data ( Fig . 2D ) well and its network topology ( Fig . 2K ) is consistent with previous results ( see 18 for discussion and parameters ) ., In order to characterize the stability of the trajectories of the gap gene system in terms of qualitatively robust features like attractors , we apply the tools of dynamical systems theory 34 , 69 ., Since the gene circuit has variables ( Gap Gene Circuits section ) its state is represented as a point in an -dimensional concentration space , or phase space ., In general the concentrations of gap proteins change with time , and hence , a solution of the gene circuit is a curve in this phase space ., The gene circuit can also have solutions which do not change with time ., Such a solution , called an equilibrium or steady state solution , is represented as a single point in phase space ., The positions of the equilibrium solutions in phase space and their stability properties determine the stability of a general time varying solution of the gene circuit ., The reader not familiar with linear analysis near an equilibrium point should see Protocol S2 for a pedagogical description of equilibria and their stability in two dimensions ., Based on the analysis in the previous section the region of interest , from 35% EL to 71% EL , can be divided into an anterior and a posterior region ( Fig . 2J ) having distinct modes of canalization and pattern formation ., The two regions are separated by a saddle-node bifurcation that occurs at 53% EL ( Fig . 3A ) , that is , at the peak of the central Kr domain ., We next demonstrate that in the anterior region ( Fig . 2J ) , which extends from the peak of the third anterior gt domain to the peak of the central Kr domain , the state of a nucleus at gastrulation is close to a point attractor ., The trajectories are stable by virtue of being in the basin of attraction of the nucleuss attractor state and hence canalize ., Pattern formation occurs by the selection of one state from many in a multistable phase space ., The concentrations of the Bcd and Cad gradients control pattern formation in the anterior by determining the sizes of the basins and the positions of the attractors , while maternal Hb concentration selects a particular attractor by setting the starting point in its basin ., Previous experimental 73 and theoretical 19 work suggested that Bcd and maternal Hb patterned the anterior of the embryo synergistically; our results identify specific roles for Bcd and Hb in anterior patterning ., The posterior region extends from the peak of the central Kr domain to the peak of the posterior gt domain ( Fig . 2J ) and its nuclei have phase spaces with very different properties ., In this region , the state of the nucleus is far from any attractor state at gastrulation ., Instead the state of a nucleus is close to a one-dimensional manifold and canalization is achieved due to attraction by this manifold ., Even though the phase space is multistable , the biological range of maternal Hb concentrations in the posterior region place all nuclear trajectories in one basin of attraction ., As a consequence , the modes of pattern formation operative in the anterior cannot function in the posterior ., Maternal Hb patterns the posterior by determining the position on the attracting manifold which a particular trajectory reaches by the time of gastrulation ., These results reveal the mechanism by which maternal Hb acts as a morphogen in the posterior 73–75 and also explain the dynamical shifts of gap gene domains 14 , 21 , a significant biological property of the posterior region ., We begin the presentation of detailed results by describing the phase spaces of typical nuclei in the two regions , highlighting mechanisms for canalization and pattern formation ., An equilibrium is labeled by either ( point attractor ) or ( saddle equilibrium ) , denoted by a superscript , with subscripts denoting the number of eigenvalues having positive or negative real parts ., For example denotes the second saddle equilibrium in the modeled region which has one eigenvalue with positive real part and three with negative real parts ., Equilibria are also given descriptive names based on which proteins are at high levels ( on ) ignoring the proteins that are at low levels ., For example , if a point attractor is at hb-on , Kr-off , gt-on , and kni-off , it is referred to as the “hb , gt-on” attractor ., A discrete 6 , 7 and buffered response to perturbations is the hallmark of a canalized developmental system ., Without recourse to molecular data , Waddington sought to explain these two properties of the response by postulating certain favored stable developmental trajectories which he called chreods ., Our results ( see Fig . 7 for summary ) show that dynamical systems with multiple attracting states possess both of these properties ., Small perturbations are damped because of phase space volume contraction driven by attractors ., A discrete response to larger perturbations is a consequence of the discontinuous boundaries between the basins of attraction of a multistable system or of bifurcations ., Using a model based on gene expression data , we can conclude that the trajectory of the gap gene system is a chreod ., The initial high variation of gap gene expression may arise from early events governed by stochastic laws ., Previous observations indicate that the first nuclei in which gap gene transcription is activated are selected probabilistically 68 , 76 ., Moreover , the gradients of Hb and Cad proteins are formed by translational repression from the Nanos and Bcd gradients respectively 57 , 77 , 78 under conditions of relatively low molecular number 23 , which is likely to lead to intrinsic fluctuations 79 ., Our results show that a deterministic description of gap gene dynamics is sufficient to account for the reduction of initial variation regardless of its source ., It is evident however that there are at least two other types of variation that the system might be subject to ., First , a natural population will have genotypic variation which , in the framework of the model , would be reflected in the variation of its parameters ., Second , gap gene expression itself is likely to be a stochastic process rather than a deterministic one ., Notwithstanding this fact , there is no evidence in Drosophila for the coupling of molecular fluctuations to phenotypic fluctuations as seen in prokaryotes 80 , suggesting that molecular fluctuations are buffered in some sense ., We emphasize that an attractor is stable against small perturbations of the model itself 29 , and hence is a model property that is preserved to an extent if there is genotypic variation in a population or if errors are introduced by stochastic gene expression ., However , further study of both of these aspects of canalization is required in order to more fully understand their role ., With regard to pattern formation in the blastoderm , the prevailing theory is that the border positions of downstream genes are determined at fixed values or thresholds of the Bcd gradient 23 , 52 ., This idea cannot , however , account for either the low variability of downstream gene border positions 14 , 17 , 18 , or the dynamical shifts of domains in the posterior 14 , 21 ., Fixed threshold specification also cannot explain precise placement of the borders in the posterior since the low molecular number of Bcd in the nuclei implies a high level of molecular noise 23 , 81 ., In the dynamical picture ( Fig . 7 ) , contrary to the threshold view , Bcd ceases to have a role in positional specification posterior to the peak of the Kr domain since , posterior to this position , the geometry of the phase space does not change qualitatively with A–P position ., Instead , maternal Hb acts as a morphogen , obviating the problems arising from a low molecular number of Bcd ., Maternal Hb has long been recognized as a morphogen 74 , 75 for the posterior region but the mechanism with which it specifies the posterior region pattern was not clear ., As is the case with Bcd , a threshold-based theory for positional specification by Hb 82 is incomplete and requires the postulation of thresholds that can be modified by their targets ., The qualitative dynamics provides a viable mechanism for posterior patterning ., The attracting manifold is the geometric manifestation of asymmetric repression between the gap genes in reverse order of gap gene domains , ., The initial Hb concentration determines which neighborhood of the manifold the trajectory traverses as it is reaches the manifold: Kr-on , kni-on , or gt-on ., In other words , posterior patterning works by triggering particular feedback loops in the gap gene network based on maternal Hb concentration ., This mechanism also accounts for domain shifts , a property particular to the posterior region , since the trajectories mimic the geometry of the manifold as they approach it ., The dynamical analysis of the gap gene system provides a simple and integrative view of pattern formation in the blastoderm ( Fig . 7 ) ., The existence of distinct anterior and posterior patterning systems was inferred from the effect of maternal mutations on larval cuticle phenotype and was subsequently characterized in terms of the effects of the Bcd 52 , 73 , 83 , 84 and maternal Hb gradients 56 , 57 , 77 ., But where and how is the control of patterning transferred from Bcd to maternal Hb ?, Our analysis shows that the hand-off occurs at the A–P position where the Kr-on attractor is annihilated through a saddle-node bifurcation , implying a sharp rather than gradual transfer ., With knowledge of the two dynamical regimes , the complex spatiotemporal dynamics of the gap gene system can be understood in the simple terms of three mechanisms: movement of attractors through phase space , selection of attractors by initial conditions , and the selection of states on an attracting manifold ( Fig . 7 ) ., Finally , we mention the advantage of having the unexpected mechanism of a one dimensional manifold for canalization and patterning ., The Bcd concentration is a bifurcation parameter of the dynamical equations ., If there were specific attractors corresponding to each gap gene state , with bifurcations creating and annihilating them successively as the Bcd concentration is varied , the molecular noise in Bcd 23 would give rise to “jitter” or rapid switching between attractors ., The manifold with its smooth dependence on maternal Hb is qualitatively robust to such fluctuations ., In a connectionist model of cognition 85 , one dimensional unstable manifolds connecting a sequence of saddle points have been proposed as a means of representing transient brain dynamics ., The gap gene phase space is a low dimensional projection of the high dimensional phase space of all the molecular determinants in the blastoderm ., It may well be that the attractors found in our analysis are actually saddle points in the high dimensional phase space and are way points , with manifolds connecting them , rather than final end points ., The methods used to obtain and characterize the quantitative data are as described in earlier work 14 ., All gene expression levels are on a scale of 0–255 chosen to maximize dynamic range without saturation ., The numerical implementation of the gene circuit equations is as described 18 , 21 ., The gap gene circuit was fit to integrated gap gene data 14 using Parallel Lam Simulated Annealing ( PLSA ) 86 , 87 ., PLSA minimizes the root mean squared ( RMS ) difference between model output and data ., For each nucleus , data were available at nine time points ( Table S1 ) ., Search spaces , penalty function , and other annealing parameters were as described 22 , 88 ., The circuit analyzed in detail had an RMS score of 10 . 76 , corresponding to a proportional error in expression residuals of about 4–5% ., Equilibria were determined by the Newton-Raphson method as described in Protocol S3 ., One-dimensional unstable manifolds of hyperbolic equilibria were calculated by solving the ODEs using the Bulirsch-Stoer 71 method with starting points in the unstable eigenspace of the equilibria 72 ., The basin boundaries on the Hb axis were calculated by finding starting points for trajectories that reach saddle points with one positive eigenvalue ( Protocol S3 ) ., The time evolution of volume phase space was calculated as described ( Protocol S8 ) ., The methods used to calculate the equilibria branches and to determine the type of bifurcations are described in Protocol S4 .
The variation in the expression patterns of the gap genes in the blastoderm of the fruit fly Drosophila melanogaster reduces over time as a result of cross regulation between these genes , a fact that we have demonstrated in an accompanying article in PLoS Biology ( see Manu et al . , doi:10 . 1371/journal . pbio . 1000049 ) ., This biologically essential process is an example of the phenomenon known as canalization ., It has been suggested that the developmental trajectory of a wild-type organism is inherently stable , and that canalization is a manifestation of this property ., Although the role of gap genes in the canalization process was established by correctly predicting the response of the system to particular perturbations , the stability of the developmental trajectory remains to be investigated ., For many years , it has been speculated that stability against perturbations during development can be described by dynamical systems having attracting sets that drive reductions of volume in phase space ., In this paper , we show that both the reduction in variability of gap gene expression as well as shifts in the position of posterior gap gene domains are the result of the actions of attractors in the gap gene dynamical system ., Two biologically distinct dynamical regions exist in the early embryo , separated by a bifurcation at 53% egg length ., In the anterior region , reduction in variation occurs because of stability induced by point attractors , while in the posterior , the stability of the developmental trajectory arises from a one-dimensional attracting manifold ., This manifold also controls a previously characterized anterior shift of posterior region gap domains ., Our analysis shows that the complex phenomena of canalization and pattern formation in the Drosophila blastoderm can be understood in terms of the qualitative features of the dynamical system ., The result confirms the idea that attractors are important for developmental stability and shows a richer variety of dynamical attractors in developmental systems than has been previously recognized .
C . H . Waddington predicted in 1942 that networks of chemical reactions in embryos can counteract the effects of variable developmental conditions to produce reliable outcomes ., The experimental signature of this process , called “canalization , ” is the reduction of the variation of the concentrations of molecular determinants between individuals over time ., Recently , Waddingtons prediction was confirmed in embryos of the fruit fly Drosophila by observing the expression of a network of genes involved in generating the basic segmented body plan of this animal ., Nevertheless , the details of how interactions within this genetic network reduced variation were still not understood ., We use an accurate mathematical model of a part of this genetic network to demonstrate how canalization comes about ., Our results show that coupled chemical reactions having multiple steady states , or attractors , can account for the reduction of variation in development ., The variation reduction process can be driven not only by chemical steady states , but also by special pathways of motion through chemical concentration space to which neighboring pathways converge ., These results constitute a precise mathematical characterization of a healing process in the fruit fly embryo .
developmental biology, computational biology/transcriptional regulation, mathematics, developmental biology/pattern formation, computational biology, evolutionary biology, computational biology/systems biology