Computational analysis of BMP gradients in dorsal-ventral patterning of the zebrafish embryo.

The genetic network controlling early dorsal-ventral (DV) patterning has been extensively studied and modeled in the fruit ﬂy Drosophila . This patterning is driven by signals coming from bone morphogenetic proteins (BMPs), and regulated by interactions of BMPs with secreted factors such as the antagonist short gastrulation (Sog). Experimental studies suggest that the DV patterning of vertebrates is controlled by a similar network of BMPs and antagonists (such as Chordin, a homologue of Sog), but differences exist in how the two systems are organized, and a quantitative comparison of pattern formation in them has not been made. Here, we develop a computational model in three dimensions of the zebraﬁsh embryo and use it to study molecular interactions in the formation of BMP morphogen gradients in early DV patterning. Simulation results are presented on the dynamics BMP gradient formation, the cooperative action of two feedback loops from BMP signaling to BMP and Chordin synthesis, and pattern sensitivity with respect to BMP and Chordin dosage. Computational analysis shows that, unlike the case in Drosophila , synergy of the two feedback loops in the zygotic control of BMP and Chordin expression, along with early initiation of localized Chordin expression, is critical for establishment and maintenance of a stable and appropriate BMP gradient in the zebraﬁsh embryo.


Introduction
The proper functioning of tissues and organs requires that each cell differentiates appropriately for its position. In many cases, the positional information that instructs cells about their prospective fates is conveyed by concentration gradients of morphogens bound to cellular receptors. Morphogens are signaling molecules that, when bound to receptors, assign different cell fates at different concentrations (Teleman et al., 2001;Wolpert et al., 2002).
Morphogen action is of special importance in understanding development, as it is a highly efficient way for a population of uncommitted cells in an embryo to create complex patterns of gene expression in space (Gurdon and Bourillot, 2001;Lander et al., 2002;Lander, 2007).
One of the well known systems that provides positional information is the genetic network controlling early dorsalventral (DV) patterning. This patterning, studied in organisms as diverse as fruit flies (Drosophila) and amphibians (Xenopus), is orchestrated by bone morphogenetic proteins (BMP). The binding of BMPs to receptors leads to the phosphorylation and activation of Smad1/5/8 proteins, which translocate to the nucleus to regulate gene expression (Massague and Chen, 2000). BMP activity is regulated by several secreted factors including the antagonist chordin or its invertebrate homologue, short gastrulation (Sog). Recently, mathematical models have provided useful insights into the establishment and function of the BMP activity gradient in Drosophila embryo and its imaginal discs (Eldar et al., 2002;Lander et al., 2002;Eldar et al., 2003;Lander et al., 2006;Lou et al., 2004;Lander et al., 2005Lander et al., , 2007Lou et al., 2005;Mizutani et al., 2005;Shimmi et al., 2005). In the embryo, the functions of BMPs (Dpp, Scw and their heterodimers) are inhibited by the binding of Sog (together with the secreted protein twisted gastrulation (Tsg)) and this inhibition is relieved by the cleavage of BMP-bound Sog by the Tolloid (Tld) protease, which releases free Dpp. Localized expression of Dpp and Sog, together with the ability of Sog binding and cleavage to drive facilitated transport of BMPs up their own concentration gradients, allows for the rapid and reliable formation of a steep BMP activity gradient. In this system, mathematical modeling has been especially valuable in explaining a variety of experimental results.
In vertebrate embryos, experiments suggest that BMPs and Chordin also work together in patterning the DV axis, and probably act in similar ways (Holley et al., 1995;Piccolo et al., 1996;Blader et al., 1997;Holley and Ferguson, 1997;Piccolo et al., 1997;Connors et al., 1999), albeit with an inversion of organism orientation, i.e. the invertebrate ventral-to-dorsal direction needs to be understood as homologous to the vertebrate dorsal-toventral one (Kishimoto et al., 1997;Neave et al., 1997;Nguyen et al., 1998;Nikaido et al., 1999). The experimental study of vertebrate DV patterning has utilized a variety of systems, including amphibians, mammals, and fish.
Of late, the zebrafish, Danio rerio, has seen increasing use as a model for early vertebrate development. Its advantages include a short life cycle of approximately 12 weeks, which facilitates genetic analysis; and transparency of the embryo, which enables the fate of individual cells during development to be observed easily (Kimmel et al., 1995). As in other vertebrates, homologues of Drosophila Dpp and Scw (BMPs 2b, 7 and others), Sog (chordin/dino) Tsg, and Tolloid (minifin/Xolloid/BMP1) are present in the zebrafish (Mullins, 1998;Connors et al., 1999), and experiments show that BMPs are required to impose ventral fates; chordin and other inhibitors are required to antagonize BMPs and establish dorsal fates; and Tolloidlike proteases mediate the proteolytic cleavage of chordin and relief of chordin-mediated inhibition (Hammerschmidt and Mullins, 2002).
Despite fundamental similarities between DV patterning in invertebrates and vertebrates, there are substantial differences in the pace of development (Lander, 2007), the geometry over which patterning occurs, and the presence of other factors that interact with the BMP/ chordin system (e.g. Reversade and De Robertis, 2005;Rentzsch et al., 2006). In addition, there are important differences in when and where BMPs and BMP inhibitors are expressed, and how their expression is controlled.
For example, in the Drosophila embryo, the shape of the BMP gradient is determined by fixed locations of zygotic production of Dpp (only in the dorsal region) and Sog (only in the ventral region). In vertebrate embryos (e.g. Xenopus and zebrafish), expression of BMP and chordin is more dynamic and flexible: An initial domain of chordin expression is highly localized (Hammerschmidt et al., 1996a;Hibi et al., 2002), but BMP expression is relatively uniform throughout the embryo (Hemmati-Brivanlou and Thomsen, 1995;Schmidt et al., 1995;Hammerschmidt et al., 1996b;Thomsen, 1997;Nishimatsu and Thomsen, 1998;Hammerschmidt and Mullins, 2002;Wolpert et al., 2002). As the maternal cues responsible for setting initial expression domains decay away, patterns of BMP and chordin expression come under the influence of zygoticallyacting transcriptional positive feedback loops (e.g., BMP signaling upregulates BMP expression and downregulates chordin expression) (Schulte-Merker et al., 1997). Such processes ultimately operate in an embryo in which all cells seem to have the potential to express either chordin or BMP, and sharp expression domains therefore need to be maintained actively.
Is the transient, localized expression of a BMP inhibitor such as chordin essential for producing a stable BMP gradient in vertebrate embryos? Are synergistic feedback loops required for maintaining the gradient? What are the specific roles of those feedbacks? How do the geometry, size, and developmental pace of typical vertebrate embryos interact with forming BMP gradients? In this paper, we study these questions by computational analysis of a mathematical model for zebrafish embryos between the end of blastula stage and the beginning of gastrulation, when a DV gradient of bmp transcripts is most prominent (Hammerschmidt and Mullins, 2002). This model is based on known biochemical interactions among extracellular diffusing ligands BMP and Chordin; a non-diffusing cell surface receptor; and an enzyme, Tolloid, which can cleave and destroy chordin. The feedback of BMP signaling on BMP and Chordin expression is incorporated in the system, as well as an initial unstable pattern of transient BMP and Chordin expression. Additional regulatory components and interactions identified through recent experimental studies (e.g. Reversade and De Robertis, 2005;Rentzsch et al., 2006) are omitted from the analysis, but the potential effects of some of them are discussed.
Through computational analysis, we investigate quantitatively the role of Chordin and Tolloid in the BMP activity and, importantly, observe that synergistic feedback loops in zygotic gene expression cooperate with initial asymmetries (dependent on maternal factors) to regulate the gene network and create a stable BMP morphogen gradient pattern. The model is studied on a three-dimensional shell that mimics the geometry of the early embryo of zebrafish during 30%-epiboly $shield stage (Wolpert et al., 2002). The inclusion of three-dimensional geometry is important for realistic modeling of a complex embryonic shape. Extensive three-dimensional numerical simulations are conducted, and they agree well with experiments, including those in which various components are removed by genetic mutation. We also examine the robustness of BMP gradients with respect to changes in bmp and chordin dosage. Finally, we discuss the relevance of these findings to other well-studied vertebrate systems, such as Xenopus.

Models, equations and numerical methods
Although the zebrafish embryo undergoes a series of cell movement and shape change during its life cycle (Wolpert et al., 2002), from the end of blastula period to the middle of gastrula period (approximately 4.5-7.5 h after fertilization) it roughly has a shape of a three-dimensional open spherical ring, perched upon a large mass of yolk ( Fig. 1(a)). In spherical polar coordinates, this ring, O, can be defined as where R 1 and R 2 are the radii for inner sphere and outer sphere, respectively, and y O is the angel from the vertical plane. Zebrafish DV pattern formation has been divided into three different phases: establishment of the Spemann-Mangold organizer, establishment of the morphogenetic BMP gradient, and morphogenetic interpretation of the gradient by target cells (Hammerschmidt and Mullins, 2002). During Spemann-Mangold organizer formation, around the start of zygotic transcription by the embryo, a coarse DV pattern is set up, dictated by factors generated by the mother and deposited into the egg (maternal control) (Hibi et al., 2002). The zebrafish Spemann-Mangold dorsal-organizer region corresponds to the location of the dorsal embryonic shield and adjacent paraxial mesodermal and neuroectodermal tissues (Ho, 1992;Shih and Fraser, 1996;Saude et al., 2000;Hibi et al., 2002). In the second phase, the initial DV pattern is refined through the actions of zygotically expressed gene products. Ventralizing signals promote the maintenance of initial ventral specification and factors generated by the organizer, for example, Chordin, antagonize ventral specification and promote dorsal fate (Hammerschmidt et al., 1996a, b). An interaction between ventralizing BMP signals and dorsalizing BMP antagonistic factors generates a BMP activity gradient that leads to the specification of a diverse set of cell types dependent on the interpretation of the graded BMP signal (the third phase) (Gurdon and Bourillot, 2001;Green, 2002;Hammerschmidt and Mullins, 2002).
The Spemann-Mangold organizer region, where Chordin is produced substantially in the first (''maternal control'') period, is at the corner of the dorsal region, and it is roughly one-quarter the size of an embryo in both y and f directions (Shih and Fraser 1996;Hibi et al., 2002). We model this region as To represent the maternally controlled expression of chordin in the region O o , we use the expression V C org e Àat , which implies production proportional to an exponentially declining maternal signal. During the same period, BMP is also expressed under the control of maternal factors, throughout almost the entire embryo (Kishimoto et al., 1997;Nguyen et al., 1998;Koos and Ho, 1999;Nikaido et al., 1999;Hammerschmidt and Mullins, 2002), which we model as V L mat e Àbt defined everywhere in the embryo O. The production of BMP and Chordin under the influence of zygotic factors we represent as V L and V C , respectively. V L is taken to be a positive feedback function of the BMP signal, and V C to be a negative feedback function of the same signal (Schulte-Merker et al., 1997).
In Fig. 1(b), a schematic diagram for the biochemical interactions among BMP, BMP receptors, Chordin, and Tolloid is presented. With the concentration of BMP denoted by [L], the concentration of BMP-receptor complexes (to which BMP signaling is assumed to be proportional) by [LR], total receptor concentration by R 0 , the concentration of free Chordin by [C], and the concentration of BMP-Chordin complex by [LC]. D L , D C , and D LC represent for the three diffusion coefficients for BMP, Chordin, BMP-Chordin complexes. k on , k off , k deg , j on , j off , and t are the binding and degradation rates for BMP, Chordin, and their complexes. With this notation, the reaction-diffusion model is formulated as following:

ARTICLE IN PRESS
where The outer sphere surface represents the boundary between the embryo and the outside environment (corresponding in biological terms either to the enveloping layer, or to the egg shell); hence the boundary condition on r ¼ R 2 is modeled as a reflective boundary condition. The inner sphere surface r ¼ R 1 and lateral surface y ¼ y O perch upon the yolk, and we assume that molecules may diffuse into yolk. Hence leaky boundary conditions are applied on r ¼ R 1 as well as at y ¼ y O . In summary, the boundary conditions are: (1) On r ¼ R 2 (reflective boundary conditions): (2) On r ¼ R 1 (leaky boundary condition): (3) On y ¼ y O (leaky boundary condition): Actually, the simulations in the next section indicate concentrations of all molecules have no significant changes if these constants are smaller than 0.01 (not shown). Larger values for those constants imply a strong absorptive boundary effect on r ¼ R 1 and y ¼ y O . The morphogen gradients from numerical simulations using such boundary conditions are inconsistent with the experimental observations (Hammerschmidt and Mullins 2002). It suggests that the boundary effect should be more reflective than absorptive.
In this paper, the parameters related to the geometry of embryo are chosen to be R 1 ¼ 280 mm, R 2 ¼ 350 mm, and y O ¼ 801 (Wolpert et al., 2002). The diffusion coefficients are D L ¼ D C ¼ D LC ¼ 85 mm 2 s À1 , the reaction constants take on the values: and all other type of molecules are assumed with a zero initial concentration (Hammerschmidt and Mullins, 2002;Mizutani et al., 2005). We solve the system (3) numerically by a finite difference scheme (Gustafsson et al., 1995). The diffusion terms are approximated by the second order central difference. The adaptive Runge-Kutta-Fehlberg-2-3 method (Stoer and Bulirsch 1993) is used for the temporal discretization. For the results presented in this paper, 2601 (17 Â 9 Â 17) spatial grid points are used. Convergence of the calculations and better resolution are observed when the spatial meshes are refined. The overall accuracy of numerical simulations is second order in space and third order in time.

Results
Experiments show that in wild-type embryos, BMP expression is progressively lost in dorsal regions, leading to a DV gradient of BMP transcripts, which comprises both the presumptive mesoderm in the marginal region and the presumptive ectoderm in the animal region of the embryo (Hammerschmidt et al., 1996a, b;Schulte-Merker et al., 1997;Hammerschmidt and Mullins, 2002).
In Fig. 2, numerical simulations are shown for the BMP signals and other molecule concentrations using contour flood plots of the three-dimensional embryo. Additionally, the temporal evolution of the pattern of BMP expression is plotted with respect to y, i.e. along the DV axis, with fixed It is noted that the initial local production of Chordin around the dorsal organizer region leads to a coarse gradient of occupied BMP receptors. This gradually develops into a steep gradient along the DV direction and evolves to a reasonable approximation of the experimentally observed BMP signal distribution within a time scale consistent with the establishment of BMP activity gradient from 30%-epiboly to the shield stage (Hammerschmidt and Mullins, 2002). As shown in Fig. 3, the simulation results of Chordin concentration and BMP signal distribution (Figs. 3(d) and (e)) agree well with the experimental observations (Figs. 3(a) and (b)) for the wild type embryo.

Chordin sharpens the gradient
Absence of chordin is known to result in a loss of dorsal fates (Hammerschmidt et al., 1996a;Miller-Bertoglio, 1997). Fig. 2. A wild type case using the parameters:  In the ventralized Chordin mutant, the retraction of the BMP expression domain is also affected, resulting in the maintenance of high levels of BMP transcripts in dorsal regions (Hammerschmidt et al., 1996a, b;Miller-Bertoglio et al., 1997;Hammerschmidt and Mullins, 2002;Wagner and Mullins, 2002). This behavior (Fig. 3(c)) is observed in numerical simulations of the model with Chordin removed (Fig. 3(f)). High levels of BMP-receptor complexes are now distributed over most of the embryo including the dorsal organizer.

ARTICLE IN PRESS
In the absence of Tolloid metalloprotease, Chordin should be able to act only as a BMP inhibitor. Experimentally, BMP signaling is found to be reduced, and animals are dorsalized (Blader et al., 1997;Mullins, 1998;Connors et al., 1999). This behavior is observed in the model as well (not shown). When Tolloid is removed, the BMP signal is much lower and broader than wild type (Fig. 2), and free Chordin and BMP-Chordin complex dominate the whole embryo. As in Drosophila (Shimell et al., 1991;Ferguson and Anderson, 1992;Francois et al., 1994;Holley et al., 1996;Yu et al., 1996;Marques et al., 1997;Ashe and Levine, 1999;Yu et al., 2000;Decotto and Ferguson, 2001;Ross et al., 2001;Scott et al., 2001;Eldar et al., 2002;Eldar et al., 2003;Kao et al., 2003;Shimmi and O'Connor, 2003;Lou et al., 2005;Mizutani et al., 2005;Shimmi et al., 2005), a balance between Chordin and Tolloid is required to create a sharp BMP signaling gradient along the DV axis of the embryo.

Two synergistic feedbacks drive and regulate the BMP gradient
Unlike the situation in the Drosophila embryo, in which fixed expression domains of BMP (i.e. Dpp) and chordin (Sog) generate a gradient of BMP protein and activity, in vertebrate embryos gradient formation reflects the remodeling of these expression domains over time. Because the capacity to produce BMP or Sog appears to reside everywhere, long-term patterns of expression should be a function of both the initial conditions and factors that determine how the system evolves over time.
In simulations of the model, we observe that one of these factors is the balance between the influence of BMP signaling on BMP vs. Chordin expression. When these two feedback loops are not well tuned, it is found that the BMP-receptor concentration approaches the following constant (in space) steady state away from the boundaries of the embryo: This expression happens to be the steady state solution for BMP-receptor concentration when diffusion is set to zero and the boundary conditions are no-flux in the system (3). The choice of the two feedback parameters g L and g C is thus critical for creation of non-uniform morphogen gradients. We define the quanitity G max y ð½LRÞ À min y ð½LRÞ (4) on the one-dimensional curve as an approximate measure of the highest value of the slope of the gradient of BMP signaling in the embryo. The presence of a sharp BMP gradient within the embryo requires G to be of a magnitude comparable to the initial total receptor concentration R 0 . Fig. 4(a) illustrates the sensitivity of G to the two feedback parameters. Note that G reaches a maximum around g L ¼ 10 and g C ¼ 10. When the feedback parameter g C is fixed, halving g L (corresponding to a stronger positive feedback from the signal to the BMP expressions) leads to a uniform BMP signals because increased BMP leads to receptor saturation ( Fig. 4(b)); on the other hand, doubling g L results in a relatively flat BMP signaling gradient, since the weaker BMP expression does not provide enough BMP to overcome inhibition by Chordin. As illustrated in Fig. 4(c), a larger g C corresponding to stronger negative feedback on Chordin production, leads the gradient of BMP signaling to extend further into the organizer region, and also to become shallower. The bi-stable nature of this two-feedback system, in which high BMP signals induce more BMP and low BMP signals induce more Chordin, which inhibits BMP, drives Fig. 4. (a) G (defined in (4)) as a function of the feedback parameters g L and g C at time ¼ 3 h. (b) The BMP activity gradient at time ¼ 3 h along y axis for four different g L with g C ¼ 10. (c) The BMP activity gradient at time ¼ 3 h along y axis for four different g C with g L ¼ 10.

ARTICLE IN PRESS
the separation of two morphogen states and creates a gradient. This also suggests that an initial morphogen gradient is needed. In fact, the interaction between the maternally-driven phase of initial Chordin expression, which is localized in the dorsal organizer, and the initial phase of expression of BMP, which is taken to be everywhere in the whole embryo initially, may result in only a shallow BMP-receptor gradient. As the initial effects quickly disappear, the two zygotically-driven feedbacks take over and dominate the dynamics.
However, it is found that the strength and the decay rate of the initial expression of BMP must be chosen in an appropriate manner to generate a stable desired gradient. As illustrated in Fig. 5(a), initial BMP production, V L mat , should be large when the decay rate, b, is fast and small when the decay is slow. The observed relationship has an approximate form: V L mat =b % 3 ¼ R 0 . Interestingly, solving the equation for BMP assuming zero diffusion and V L mat bV L max , we obtain ½L ¼ ðV L mat =bÞð1 À e Àbt Þ. Therefore, the relationship in Fig. 5(a) indicates that initial BMP expression and its decay rate must be balanced such that free ½L is comparable to the total receptor concentration. Only when the free BMP concentration is large enough to generate a substantial initial BMP gradient during the maternal control period, the feedback loops can be effective, and drive the initial pattern to a stable steady BMP gradient. In simulating the model, a decay time is generally chosen that leads such expression to finish in approximately 10 min (600 s); consequently, we choose V L mat ¼ R 0 b ¼ 0:0501 mM=s as given in Fig. 2. Interestingly, the morphogen pattern generally is much less sensitive to the change of decay rate and strength for the initial expression of Chordin (Fig. 5(b)).

Robustness of the BMP activity gradient with respect to changes in chordin or BMP gene dosage
One natural test of the sensitivity of the BMP activity gradient with respect to chordin or bmp gene dosage is to model zebrafish embryos heterozygous for mutations in chordin and bmp. It is not surprising that the BMP activity gradient is very sensitive to bmp gene dosage ( Fig. 6(a)). In this case all BMP synthesis parameters were taken to be half of their wild-type values (Fig. 2), and the BMP activity gradient was observed to disappear quickly to a uniform state along the whole embryo. The BMP pattern is less sensitive to chordin dosage, as illustrated in Fig. 6(b). Here, halving chordin synthesis parameters not only increased BMP activity near the organizer region, but also decreased it at the ventral end of the embryo. Fig. 6(c) shows how the maximum steepness of the BMP activity gradient changes as a function of simultaneous changes of both chordin and bmp dosages. One can see how the gradient becomes shallower with decreasing chordin dosage, and how it collapses altogether for BMP dosages about two fold higher or lower than the optimum.
These results share many similarities with the situation in Drosophila embryos. There, a heterozygous mutation in Dpp is lethal, whereas heterozygous mutation of Sog is viable, albeit with an abnormal BMP activity gradient. Specifically, reduction in Sog dosage decreases BMP signaling at the dorsal midline (where it is normally highest), and increases it away from the midline (Mizutani et al., 2005).

Effect of organizer size
Another interesting question is how the size of the organizer region affects the development of BMP gradients. To study this, we vary the organizer region in both DV and lateral directions, and compare the new dynamics with the wild type case shown in Fig. 2. During early development of the BMP gradient for such a case (Fig. 7), the transient pattern generated by initial, maternallydetermined expression of BMP and Chordin, is quite different from the wild type in Fig. 2. The embryo with the larger dorsal organizer (compared to the wild type) initially produces a much larger region of low BMP signaling ( Fig. 7(a)), and a smaller dorsal organizer generates a much smaller initial region of low BMP signaling (Fig. 7(c)). Howerver, when both systems are close to their steady states, their BMP patterns (Figs. 7(b) and (d)) then become very similar to the wild type ( Fig. 2(b)).

ARTICLE IN PRESS
Next, we mimic the shield removal experiments (Saude et al., 2000) by excluding a shield from the organizer region in the model such that the organizer contains two isolated regions in the lateral direction. As seen in Fig. 7(e) in the animal view, the shield region in the dorsal corner still has high BMP signaling during the transient dynamics. This is because there is no maternally-determined expression of Chordin in that region. As the maternal loading fades away, the system recovers and a normal pattern (Fig. 7(f)) emerges when the regulation of the two zygotic feedbacks, which are independent of the size and shape of the dorsal organizer, start to dominate the dynamics. The recovery of the system after shield removal is consistent with the experimental observation (Saude et al., 2000).
All the simulations in Fig. 7 were carried up to 20 h, and the solutions using different organizers all converge to a pattern similar to the wild-type. Mathematically, this is expected because as time becomes large enough, the Chordin synthesis rate in (3) is independent of the organizer region. This suggests that the patterning can proceed normally when organizer region is changed.

Conclusion and discussion
In zebrafish DV patterning, establishment of a morphogenetic BMP activity gradient is driven by a multi-protein network that depends on spatially localized initial conditions of synthesis, a potential for facilitated BMP transport, and two synergistic transcriptional feedback loops. The goal of the present study was to use mathematical modeling and computation to elucidate how these processes work together to create and maintain the BMP activity gradient.  We began by establishing a model for DV patterning of the early zebrafish embryo during the end of blastula to the middle gastrula period. The shape of the embryo is modeled as a three-dimensional open spherical ring. In the model, a ligand (BMP) diffuses freely, binds with its receptors, and degrades through the ligand-receptor complex. BMP also binds to another diffusing molecule, Chordin, to form a complex. This complex then can be cleaved by Tolloid to free BMP and destroy Chordin. The BMP-receptor complex, regarded as the signal, in turn regulates the synthesis of BMP and Chordin. The resulting reaction-diffusion system in three spatial dimensions was solved using a finite-difference method.

ARTICLE IN PRESS
We found that key experimental (e.g. genetic) observations are well captured by the proposed model. The simulations and analysis reveal that Chordin, along with Tolloid, are required to form a sharp gradient, and the two synergistic feedbacks on BMP and Chordin synthesis are critical for maintaining and refining the coarse pattern initiated by the maternally-dependent control of early zygotic BMP and Chordin. The BMP gradient is sensitive to the strength of the positive feedback on BMP expression itself, and the negative feedback on Chordin expression limits the ability of Chordin to inhibit BMP. Thus, it confines the scope of dorsalization by Chordin.
The results bear both similarities with and differences from models of DV patterning in Drosophila. In both organisms, the action of Tolloid on Chordin/Sog drives transport of BMP up its own concentration gradient, leading to gradient sharpening. In both organisms, the pattern is predicted to be highly sensitive to BMP dosage, and less so to Chordin/Sog dosage. A major difference is that Drosophila DV patterning is traditionally conceived of as an open-loop process, wherein domains of BMP and Sog expression are assumed to be fixed during patterning. More recently, some feedback has been discovered during late stages of Drosophila DV patterning (Wang and Ferguson, 2005), although it appears to act at the level of cellular responses to BMPs, not BMP expression.
As a result of the presence of transcriptional feedback loops in zebrafish DV patterning, steady state gradient steepness is relatively insensitive to initial chordin expression level or size of the domain of chordin expression. This result may help explain experimental manipulations in which the embryonic shield (a major portion of the dorsal organizer) was extirpated from zebrafish embryos, yet normal DV patterning ensued (Saude et al., 2000).

Another significant difference between zebrafish and
Drosophila is in the three dimensional geometry of the embryo. In Drosophila, BMP diffusion occurs within a thin, bounded compartment surrounding a roughly ellipsoidal embryo; DV gradient formation is readily simulated as a reaction-diffusion problem in one dimension, with little loss of accuracy (Eldar et al., 2002;Lander et al., 2002;Eldar et al., 2003;Lander et al., 2006;Lou et al., 2004;Lander et al., 2005Lander et al., , 2007Lou et al., 2005;Mizutani et al., 2005;Shimmi et al., 2005). In view of the more complex arrangement of the zebrafish embryo, we chose here to carry out fully three-dimensional simulations, which are computationally much more costly. One may also attempt a one-dimensional approximation (in the y direction) of the zebrafish embryo (e.g. at fixed r ¼ 0.5(R 1 +R 2 )), as shown in Fig. 8. Here Eq. (3) were solved in one-dimension using the same parameters as those in Fig. 2. The results are very similar to the corresponding one-dimensional cut of the three-dimensional simulation in Fig. 2(c). However, when the onedimensional cut is made closer to the lateral margins of the embryo, significant differences are observed, due to boundary effects (not shown). The fact that one-dimensional simulations capture many (albeit not all) of the salient features of DV patterning in the zebrafish, supports the view that the basics of patterning are relatively insensitive to changes in embryonic shape. This, in turn, suggests how essentially the same molecular system that works in the fish embryo (a spherically curved sheet of cells atop a ball of yolk) can also control DV patterning in the embryos of the frog (an asymmetrically hollowed out ball of cells), the chick (a flat sheet of cells atop yolk), and the mouse (an irregular, hollow cylinder of cells). More direct confirmation of this conjecture is worthwhile, especially in the amphibian, Xenopus, in which DV patterning has been well studied, and many experimental manipulations have been done. Xenopus differs from zebrafish not just in geometry, but also in some of the detailed mechanisms of how and where the organizer is formed (Wolpert et al., 2002).
It is generally believed that embryonic patterning systems employ complex networks of gene regulation and molecular interaction in order to generate patterns that are highly robust to environmental and genetic perturbations (Lander, 2007). The mechanisms explored in the present study form BMP gradients with shapes that are reasonably robust to some perturbations (e.g. changing chordin dosage, or organizer size; Figs. 5-7), however it is ARTICLE IN PRESS noteworthy that the precise location along the DV axis at which the BMP gradient undergoes its major decline, is much less robust than the local shape of the gradient at that position. This suggests that additional regulatory mechanisms may be needed to help control the position of the BMP gradient. Recent studies in both zebrafish and Xenopus have identified new feedback processes that may play such roles.
In zebrafish, for example, it was recently shown that crossveinless-2 (cv2), a diffusible molecule structurally related to Chordin, is expressed in response to BMP signaling, and competitively inhibits the binding of chordin to BMP, thus creating an additional positive feedback loop in BMP signaling (Rentzsch et al., 2006). Making matters more complicated, cv2 apparently undergoes proteolytic cleavage that affects its association with the extracellular matrix, altering its ability to promote BMP function.
In both zebrafish and Xenopus, it has been found that a single BMP family member, known as ADMP, is expressed in a pattern opposite to that of classical BMPs, i.e. it is produced in the organizer region and not elsewhere, where it seems to play a negative feedback role in controlling organizer size (Lele et al., 2001). Recent experiments in Xenopus indicate that ADMP is required for embryos to exhibit a remarkable feat of self-regulation: when early Xenopus embryos are bisected, dorsal halves often develop normally, initiating BMP signaling on their ventral sides (Reversade and De Robertis, 2005). Other negative feedback factors in the embryo include molecules such Sizzled and BAMBI, which are expressed on the ventral side (Reversade and De Robertis, 2005). Incorporation of any of these additional molecular interactions into reactiondiffusion schemes outlined in Fig. 1(b) is straightforward. An important goal for future work will be to determine whether such interactions can confer and increase the robustness of both the location and shape of the BMP gradients that pattern the DV axis of vertebrate embryos.
One should notice that in our model the two diffusion coefficients of the two morphogens are the same, and there is no long-ranging inhibition associated with diffusion, which is a critical component required for pattern-forming in Meinhardt-Turing system (Meinhardt, 1982). In our model, it is the two zygotic feedbacks that refine and amplify initially shallow pattern by maternally-driven to a stable morphogen gradient that to some degree is independent of the initial morphogen concentrations. This limited gradient forming capability unlike the Meinhardt-Turing patterning, is consistent with the limited regeneration capabilities in the zebrafish and in Xenopus.
In summary, the results presented here of modeling and simulating the establishment of BMP gradient in DV patterning of the early zebrafish embryo indicate the indispensability of the interactions of BMP with nonsignaling secreted proteins Chordin and Tolloid, the regulation of two synergistic feedback loops, and an initial localization of Chordin in the dorsal organizer during the earliest stages of patterning. In addition, other feedback processes are likely to be required to ensure robust and regulative pattern formation.