How does pattern formation occur accurately when confronted with tissue growth and stochastic fluctuations (noise) in gene expression? Dorso-ventral (D-V) patterning of the mandibular arch specifies upper versus lower jaw skeletal elements through a combination of Bone morphogenetic protein (Bmp), Endothelin-1 (Edn1), and Notch signaling, and this system is highly robust. We combine NanoString experiments of early D-V gene expression with live imaging of arch development in zebrafish to construct a computational model of the D-V mandibular patterning network. The model recapitulates published genetic perturbations in arch development. Patterning is most sensitive to changes in Bmp signaling, and the temporal order of gene expression modulates the response of the patterning network to noise. Thus, our integrated systems biology approach reveals non-intuitive features of the complex signaling system crucial for craniofacial development, including novel insights into roles of gene expression timing and stochasticity in signaling and gene regulation.
A fundamental question in developmental biology is pattern formation, i. e. the acquisition of positional identity in cells resulting in spatially organized domains of gene expression. Computational analyses have long sought to address how patterning occurs in growing tissues that change their size and shape by modeling morphogen gradients, signaling between cells, geometric transformations and other mathematically-amenable aspects of development [1–4]. More recently, computational modeling has revealed how signaling networks integrate with one another and the importance of feedback loops in precise regulation of early developmental patterning systems [5–11]. Models for more complex developing structures such as vertebrate limb buds [12,13], hair follicles [14,15], pigment cells in the skin [16], the spinal cord [17,18] or the palate [19,20] require integrating multiple signals within rapidly expanding three-dimensional (3D) tissues. Pharyngeal arches are bilateral, segmentally-repeated structures that form in the ventral head of vertebrate embryos and give rise to skeletal, muscle and connective tissues of the face and neck, including the upper and lower jaws. Arches are complex both in their 3D morphologies and in their embryonic cellular origins. Streams of cranial neural crest (NC) cells migrate into each arch segment and surround cores of myogenic/vasculogenic mesoderm. The surrounding ectodermal and endodermal epithelia produce signals that subsequently pattern the arch along its dorso-ventral (D-V) axis, resulting in at least three early domains: ventral (V), intermediate (I), and dorsal (D) [21–23]. D-V arch patterning involves a highly-conserved signaling network consisting of the Bone morphogenetic protein 2/4/7 (Bmp) and Endothelin-1 (Edn1) signaling pathways, secreted by ventral arch epithelia [24–28], and dorsal Jagged1 (Jag1) /Notch signaling [29,30]. Errors in these signals can lead to craniofacial birth defects, such as auriculocondylar syndrome in humans, in which Edn1 signal transduction is disrupted leading to partial homeotic transformation of ventral skeletal elements to a dorsal fate [31,32]. Understanding how D-V domains arise in the midst of NC migration and arch growth is both an experimental and a computational challenge. Both Edn1 and Bmp are crucial for ventral and intermediate arch development, but with distinct effects on gene expression [24,27,30,33–35]. Bmp, which acts as a morphogen in many contexts [36–38], primarily induces and maintains genes expressed ventrally such as Hand2, a critical transcription factor for ventral mandibular identity. In contrast, while Edn1 also induces ventral genes initially in a concentration-dependent manner [26,39] it later becomes primarily required for expression of intermediate genes such as Dlx5/6, which are required for ventral/intermediate mandibular (lower jaw and jaw joint) development, and less dependent on Bmp [24,29,30]. Craniofacial patterning defects in edn1-/- mutants can be largely rescued by injection of Edn1 protein throughout the arch [27], and recent work suggests Edn1 plays a more permissive than instructive role [40]. Given their common targets, what are the advantages of having these two ventral morphogens acting in parallel during early D-V arch patterning? In addition, how does patterning occur robustly in the face of continuous cell divisions, rearrangements, and noise in both the signaling molecules and their downstream gene regulatory networks (GRNs)? To address these questions, we have developed the first computational model of arch D-V patterning that incorporates growth, migration, gene expression and different sources of noise. We represent the known components of the arch GRN with a system of ordinary differential equations (ODEs) that accurately reproduces published genetic perturbations of arch D-V patterning. We establish the model using measurements of spatiotemporal patterns of gene expression and 3D NC cell movements obtained from time-lapse movies of live zebrafish embryos. Quantitative temporal gene expression data reveal that intermediate domain genes are expressed before genes marking the ventral domain and dorsal genes are expressed last. The model confirms that this temporal order of intermediate-ventral-dorsal (IVD) patterning improves some aspects of robustness of D-V patterning (precision, referring to consistency across simulations/embryos), while making other aspects (accuracy, referring to closeness to the theoretical ideal pattern) more sensitive. The model further suggests that Bmp signaling primarily establishes the sizes and positions of patterning domains, while Edn1 plays a permissive role, and that noise in the GRN and each of the signaling pathways affects patterning differently. Our model reveals novel features of the early spatiotemporal dynamics of gene expression that are critical for patterning the complex 3D structure of the craniofacial skeleton during embryogenesis.
Institutional Animal Care and Use Committee protocol #2000–2149.
Neural crest (NC) -derived ectomesenchymal cells in pharyngeal arches 1 (mandibular) and 2 (hyoid) in zebrafish are patterned into three D-V domains between 14–36 hpf, which give rise to distinct skeletal elements in the adult (Fig 2A–2C). Arch D-V length roughly doubles over this period (from 30 to 60 μm). Previous in situ hybridization (ISH) studies have shown that dlx3/4/5/6 are expressed together in an early ventral-intermediate (V-I) domain that later separates into V and I [24,29,30,42]. hand2, the homolog of which is induced by Dlx5/6 in mice [34,60,61], marks the new V domain and represses Dlx genes ventrally, restricting their expression to the I domain [42,62]. However, the precise timing of gene expression between 14–20 hpf, when these domains first appear, remains unclear. To address this, we have measured transcript levels of seven D-V patterning genes in FAC-sorted arch cells using NanoString analysis (dorsal: jag1b, hey1; intermediate: dlx3b, dlx4b, dlx5a, dlx6a; ventral: hand2) (Fig 2D). Surprisingly, expression of dlx3b peaks early at 20 hpf, followed six hours later by three other intermediate genes (dlx4b, dlx5a, dlx6a), the ventral gene (hand2), and the dorsal gene jag1b at 26 hpf. Similarly with hybridization chain reaction (HCR) in situs, to facilitate co-localization and quantitation of expression (Fig 2E–2T) we find that within the domain of dlx2a expression, which marks NC cells in the entire arch, strong dlx3b expression is detected at ~17 hpf, at least an hour before dlx5a expression appears faintly at ~18 hpf. Meanwhile, hand2 expression is not detected until ~20 hpf. Interestingly expression of hand2 arises abruptly, while other genes such as dlx5a appear more slowly, yet both peak at a similar time point in the NanoString analysis. Thus, intermediate genes are the first to be expressed in the D-V sequence of arch patterning followed by ventral and finally dorsal genes. The 1D model (Fig 3A–3E) recapitulates the relative timing and sizes of D-V domains in the mandibular arch. Initially intermediate gene expression extends from the ventral end to approximately halfway up the D-V axis. Subsequently, dorsal genes are expressed at the dorsal end of the arch, leaving a section of “unpatterned” cells (white regions, in which gene expression is below the arbitrary 20% threshold) between I and D domains, which gradually diminishes. Initiation of ventral gene expression at 28 hpf creates a narrow V domain, which moves the I domain dorsally. The 2D model (Fig 3F–3J) captures these spatiotemporal dynamics of D-V domain formation. Here, individual cells are also defined as patterned if their gene expression exceeds 20%, and they are gradually colored correspondingly, such that grey indicates cells not yet expressing genes above the cut-off and more deeply colored cells indicate higher and higher levels of gene expression. Initially none of the cells express any of the genes above the 20% cut-off (grey cells). Between 22–35 hpf the arch elongates anteriorly and ventrally. During this tissue deformation the cells acquire D, I and V fates, (yellow, blue and pink, respectively) and form domains of the correct size and shape. The simulation results agree with live imaging of hand2: GFP: sox10: lyn-tdTomato double transgenics (Fig 3K–3O) or dlx5a: GFP; sox10: lyn-tdTomato double transgenics (Fig 3P–3T). We note that the boundaries of gene expression as shown by transgene reporter intensity are sharp (S2 Fig). To generate a 2D cell-center based model that reflects arch morphogenesis as accurately as possible we have measured mandibular arch deformation in images of 6 sox10: lyn-tdTomato transgenic zebrafish embryos, the average of which is used to generate an arch outline (Fig 3K–3T). By further analyzing time-lapsed images of sox10: nEOS transgenics we find that: 1) cell number roughly doubles between 22–36 hpf, 2) ~90% of this increase in cell number is due to cell division and 3) only ~10% is due to cells migrating into the arch dorsally. These parameters are incorporated into the model to compute the time intervals of cell division and cell migration. The 2D model also reproduces previously reported phenotypes of genetic or pharmacological perturbations that disrupt D-V patterning (Fig 4). V/I domains do not form and the D domain expands in embryos lacking Bmp or Edn1 signaling (reduced in the modeling simulations to 1% of wild-type expression [24,27] (Fig 4A and 4D). The I domain does not form and D expands ventrally in embryos overexpressing the dorsal factor Jag1 (500% of wild-type production) (Fig 4C). Conversely, the I domain expands to replace D in a jag1-/- mutant (1% of wild-type expression) or when Edn1 is overexpressed (500% of wild-type gradient maximum) (Fig 4B and 4F). Importantly, the model also recapitulates the normal D-V patterning observed experimentally with moderate, uniform Edn1 expression (50% of wild-type gradient maximum), achieved with Edn1 protein injections [24,27] (Fig 4E). Thus, despite the minimal GRN on which it was based, the model captures many aspects of patterning observed experimentally in vivo. For any three sets of D-V patterning factors, there are six possible different temporal orders of gene expression (VID, VDI, IVD, IDV, DIV, DVI). Our gene expression studies indicate that dlx3b and its associated I domain appear first, so we asked how this order fares in our model as compared with other possible orders. By varying the production and degradation parameters for each gene group we simulate all six temporal orders (S1 and S2 Tables) and examine if any one is more robust than another. Changing the temporal order does not affect the final pattern at 35 hpf (S3 Fig). To compare sensitivity between the three temporal orders of gene expression we compute si (t) (Eq 9), evaluated at 10 equidistant time points between 22–35 hpf. We summarize the 10 resulting data points in box plots (Fig 5, S4–S6 Figs), where the line in the middle denotes the median, the bottom and top edges of the box the 25th and 75th percentile, respectively, and the whiskers extend to the most extreme data points. We normalize the data with respect to the highest value of si (t) for both parameters p1 and p5 at all time points (for unprocessed, time-dependent data see S7 Fig). The distance between the median line and zero (dashed black lines) indicates how sensitive gene expression is to perturbations in the control parameters p1-p7, and the size of the box and length of the whiskers indicate how much the sensitivity changes with time. In general, all model simulations are relatively robust to parameter variations, indicating that the results are not due to the specific choice of parameters. When cells are exposed to morphogen concentrations typical for the V and I domains (high Bmp and Edn1 concentrations) the GRN is most sensitive to perturbations in the BMP signaling parameters p1 (activation of V genes) and p5 (activation of I genes) (Fig 5; S4 and S5 Figs), and less so to variations in Edn1 signaling parameters p2 and p7. This is particularly true in cases where I genes are expressed last (VDI and DVI), since median sensitivity levels of dorsal genes deviate further from zero (dashed black lines; Fig 5 and S4 and S5 Figs). However, when cells are exposed to low Bmp and Edn1 concentrations (S6 Fig), typical for the D domain, perturbations in the Bmp gradient only affect the ventral genes, while intermediate and dorsal genes are sensitive to perturbations in the Edn1 parameters p2 and p7. This agrees with published evidence that Edn1 plays a primarily permissive role in ventral and intermediate gene expression [27,40]. A permissive role for Edn1 is further evident when we compute the sensitivity to parameter variations in the GRN if only one of the two morphogen gradients is present (S8 Fig). The results are similar to those for two morphogens in all three V, I and D domains. In all cases, the most crucial parameter is the one controlling the effect of the morphogen on the I domain, p2. To simulate the stochasticity that may occur in vivo, we have added Gaussian noise to the morphogen gradients and to the GRN (Eqs 11–14, Figs 6 and 7, S9–S11 Figs and S6 and S7 Movies). To compute a large number of stochastic simulations for statistics, and for better visualization, we investigate the effects of noise first in the 1D model. We plot the mean (thick lines) and ±σ (shaded regions) over 100 simulations (Fig 6B–6E, 6A for comparison without noise). For the 2D model, we show end-states for single simulations (Fig 6B’–6E’, 6A’ for comparison without noise). The simulations reveal that noise in Bmp signaling is transmitted differently than noise in Edn1 into the GRN. While intermediate gene expression is affected by both signals, ventral gene expression is not affected by noise in the Edn1 gradient, since the only input into the V domain is the Bmp gradient. Dorsal genes are most robust to morphogen fluctuations (Fig 6C–6D’), since they are mostly controlled indirectly through the ventral genes, but they are the most sensitive to gene regulation noise (Fig 6B and 6B’). Since Edn1 is required for intermediate gene expression, but increasing Edn1 signaling does not expand the I domain (Fig 4E and 4F), fluctuations in Edn1 only act in one direction, meaning that lower Edn1 levels due to noise reduce the I domain but higher levels have no effect. As a result the I domain is reduced with fluctuations in Edn1 (Fig 6D and 6D’), when compared to the deterministic case (Fig 6A and 6A’, S9 Fig) in both 1D and 2D models. This means that Edn1 acts “unidirectionally” as a permissive factor. Gene expression noise appears to be the strongest driver of fluctuations in patterning, since even relatively small fluctuations (ν = 0. 05) in the GRN alter gene expression profiles substantially (Fig 6B and 6B’). The effects of individual fluctuations are additive when all sources of noise are combined (Fig 6E and 6E’, S11 Fig), i. e. both increased fluctuations in dorsal gene expression due to GRN noise, and expansion of the V domain due to noise in the Edn1 gradient. Precision refers to the level of variation in boundary positions in stochastic simulations. Similar to our analysis of sensitivity to different parameters in the GRN, we examine if the sequence of D-V domain formation influences precision, as a measure of robustness in response to noise. When fluctuations are limited to the Bmp gradient the genes expressed first absorb most of the noise (Fig 7). The exception is when dorsal genes are expressed earliest, which are in general the most robust to Bmp noise, such that the DIV and DVI sequences are the least susceptible to Bmp fluctuations (Fig 7E and 7F). This is in contrast to noise in the GRN (v = 0. 05), where genes expressed earliest are the most robust and genes expressed later are more susceptible to noise (S9 Fig). This is particularly true when genes expressed in the I domain are first in the D-V sequence. This indicates that patterning I first is beneficial since fluctuations in intermediate gene expression affect both the precision of the V-I and I-D boundaries. With fluctuations in the Edn1 gradient, the I domain is severely reduced when intermediate genes are expressed last, while there is still a distinct I domain in the case of either V or I being expressed first (S10 Fig), further indicating that an early expression of intermediate genes is beneficial for boundary accuracy. Noise in the Edn1 gradient has the most similar effect across the 6 possible temporal orders. Early expression of intermediate genes, however, leads to slightly stronger fluctuations in their expression as a result of Edn1 noise, while still preserving a distinct I domain. Due to the different effects of the three sources of noise in the context of different temporal orders, combining all sources of noise indicates that a later onset of intermediate gene expression leads to the strongest fluctuations and higher sensitivity (S11 Fig). When differences in the temporal order are enforced explicitly in the simulations (S12 Fig), the responses to the different sources of noise are similar to the case of intrinsic regulation, only with a more severe loss of the I domain due to noise in Edn1 (S13 Fig). However we do not see any distinction between the temporal orders in the extrinsic model with noise in the Bmp gradient (S14 Fig) in contrast to the intrinsic model with Bmp noise (Fig 7). Accuracy refers to boundary positions relative to the measured wild-type positions (Fig 8). We simulate 10 runs of the 2D model, compute statistics of the boundary error (E) in Eq 8, normalize according to the highest value and plot the mean (line) and ±σ (error bars). Accuracy depends strongly on which gene group is expressed last, especially for the I-D boundary. When either Bmp or Edn1 are noisy, the I-D boundary is positioned most accurately when I is last, slightly less accurate when D is last and inaccurate when V is last (Fig 8B and 8D). In contrast, when Bmp is noisy the V-I boundary shows no clustering of temporal orders of patterning (Fig 8A), but when Edn1 is noisy accuracy increases when I is last (Fig 8C). The boundary accuracy depends less distinctly on the temporal order of patterning when the gene regulation is noisy (Fig 8E and 8F). When we combine all sources of noise a late appearance of the I domain clearly improves positioning of the V-I boundary and slightly improves the I-D boundary (Fig 8G and 8H). In general the domain boundaries are more sensitive to fluctuations in Bmp than in Edn1, especially at the V-I boundary. Presumably the mutual inhibition between intermediate and dorsal genes makes the I-D boundary more robust to stochastic fluctuations. When the sequence of patterning is controlled by extrinsic factors not inherent in our minimal GRN (S12 Fig), the accuracy of boundary positioning is more sensitive to individual sources of noise in all cases except for positioning the V-I boundary with fluctuations in only Edn1 (S15 Fig). However, with the external control of gene expression timing the noise effects appear to be less additive and the extrinsic model positions the boundaries more accurately. From the above results, different temporal orders lead to different degrees of noise in gene expression profiles (Figs 6,7 and S9–S11), while there are observed differences in accuracy of domain boundary positioning (Fig 8). We now quantify these effects to investigate how accuracy and precision relate to each other at both domain boundaries, and with the different sources of noise (Fig 9). Precision is measured as the maximal standard deviation (Figs 7 and S9–S11), corresponding to the maximal width of the shaded regions. Accuracy is the mean value of the boundary error (Fig 8). Hence for both, a low value indicates higher accuracy or higher precision. When all sources of noise are included in the simulation, an anti-diagonal trend is observed, indicating a trade-off between precision and accuracy, where temporal patterns with more precise boundary positioning have less accuracy, and vice versa. The observed IVD pattern favors precision over accuracy at both boundaries, suggesting that the patterning system has evolved to maximize precision in gene expression domain boundaries.
We have analyzed spatiotemporal patterns of expression of D-V patterning genes during pharyngeal arch morphogenesis in zebrafish embryos and combined our experimental observations with published data to generate the first computational model of the developing mandibular arch. Previous efforts have compiled information from imaging and gene expression databases, such as the FaceBase consortium [63], and analyzed cellular behavior during craniofacial morphogenesis [20], but ours is the first to integrate spatiotemporal gene expression patterns with morphogen gradients, tissue measurements, and known mutant phenotypes. Our model captures many of the spatial and temporal features of arch development and recapitulates genetic perturbations. It also provides novel insights into developmental constraints on the system, including: 1) Bmp is responsible for providing positional information to the cells, while Edn1 is permissive, 2) the temporal order of patterning is important for the system’s capacity to account for noise, and 3) the temporal order favors precision over accuracy in boundary positioning. Many developmental processes that involve periodic patterning reflect underlying reaction-diffusion systems that deal efficiently with noise through their intrinsic feedback loops [12,16,19]. We find that the temporal order of gene expression provides a previously unappreciated factor in improving responses to noise. Surprisingly in our experiments it is the intermediate gene dlx3b that is the first detected at 16–17 hpf by HCR and slightly later (20 hpf) in our NanoString analyses, which may reflect the fact that sorted cells used for NanoString are derived from multiple arches (mandibular, hyoid, branchial) while with HCR we image the first arch directly. dlx3b might serve as an early response factor, possibly integrating input from multiple signals and priming other patterning genes in the GRN. The early appearance of the I domain means that the arch is not patterned consecutively from one end to the other (i. e. the VID or DIV orders), but rather in the peculiar sequence of establishing the middle first. Intuitively one can assume that patterning the central domain first already establishes both boundaries. However, that is not the case here either, since expression of the intermediate genes initially extends to the ventral end of the arch and is only later replaced by the ventral genes, such that both boundaries are established independently and non-simultaneously. We have investigated if this temporal order is beneficial for the response of the system to stochastic fluctuations. Our model simulations suggest that avoiding having intermediate genes expressed last improves the robustness to perturbations in Bmp signaling parameters and precision in the positioning of D-V domain boundaries. While dlx3b knockdown does not cause significant patterning defects, this could reflect compensation by dlx4b or dlx5a [42]. The unique function of dlx3b in the GRN is supported by its distinct spatial expression pattern from the dlx5/6 pair and also from its neighboring dlx4 cluster counterpart [42,64–66]. Our results suggest that altering the temporal order of D-V gene expression, (e. g. using optogenetic approaches to induce expression of hand2 in its normal spatial domain but prior to onset of intermediate gene expression and thus creating a VID sequence) will disrupt the accuracy and robustness of D-V domain boundaries. Future studies are also needed to determine if the temporal order of gene expression in this system is controlled by gene-intrinsic differences in sensitivity to signals (as is the case in our minimal model), or if trans-acting factors play a greater role. We also note that the precision and accuracy of D-V domain boundary formation are distinctly susceptible to noise, with the boundary between V and I domains especially sensitive to noise in Bmp signaling or in the downstream GRN. Precision and accuracy in the boundaries of gene expression domains are both goals for patterning systems, but our model reveals divergent paths to reach each of those goals in the mandibular arch. We find that the most precise gene expression profiles in response to Bmp signaling are achieved when dorsal genes are expressed first (Fig 7), while the most accurate boundaries were obtained when I genes were expressed last (Fig 8). Instead, our data in zebrafish indicate that I genes are expressed first, reflecting a trade-off that seems to favor precision over accuracy (Fig 9). It is important to note here that the model does not achieve high boundary positioning accuracy for the I-D boundary with any of the temporal patterns even in the absence of noise (S1A Fig). This might be a result of measurement inaccuracy and due to the projection of a 3D geometry into 2D, or that other factors not included in our minimal GRN are responsible for pushing the I-D boundary more towards the ventral end. Such trade-offs between precision and accuracy in other systems typically involve negative feedback [12,18,67], while it is temporal gene expression control that navigates this tradeoff in our model. Future analysis will reveal if feedback loops also operate at arch domain boundaries. Our modeling results for precision and accuracy derive from combinations of 10–100 simulations, somewhat analogous to variability in patterning that can occur between individual embryos. We speculate that developmentally this suggests that patterning will favor precision (i. e. show a narrow range of variability between individuals), even if that range differs from the species ideal. Evolutionarily, this might permit rapid variation in craniofacial morphology during radiation events while still maintaining intra-population characteristics. Changes in craniofacial structures are one of the most striking adaptations in rapidly evolving vertebrate populations such as African rift lake cichlids, and it would be interesting to examine the parameters of patterning in these species [68,69]. Bmp and Edn1 signaling have overlapping but distinct functions in D-V patterning of the mandibular arch, with Edn1 primarily maintaining intermediate gene expression and playing a permissive rather than instructive role [24,27,30,40]. Our modeling results are consistent with Bmp signaling acting as the principal instructive ventralizing signal, at least for inducing V and I domains. They also suggest that moderate perturbations or noise in Bmp signaling, at levels that do not completely eliminate expression of various factors or formation of later cartilage structures, will still have effects on D-V arch patterning. In several reports, changes in levels of patterning signals result in a moderate-to-severe spectrum of phenotypes. Reduction in Bmp signaling by heat-shock induction of a dominant-negative Bmp receptor can lead to varying degrees of loss of ventral structures, depending on the timing of heat-shock and the dose of the dominant-negative transgene [24,30,70]. Increasing Bmp levels moderately by overexpression of Bmp ligand or injection/bead insertion of Bmp proteins causes homeotic transformations, such as the more dorsal palatoquadrate cartilage acquiring characteristics of the ventral Meckel’s cartilage, while more severe changes in signaling levels leads to significant losses throughout the jaw cartilage [24,30,71]. Modulation of Edn1 signaling in either direction similarly results in varying degrees of patterning changes, although consistent with its permissive role even very high levels of Edn1 do not expand the I domain, while reducing the D domain [24,26,27,30,39]. These divergent phenotypes, along with our analyses of parameter sensitivity and the effects of noise, suggest a region of stability where the patterning network can compensate for changes in signals, with moderate phenotypes appearing at the edges of this region, and severe phenotypes occurring when signaling falls outside the stable region altogether. Our computational model suggests that further experiments up- and down-regulating signaling will reveal transition points where stability breaks down, and those data will in turn improve the quality of our model. We can also extend our framework in the future to incorporate other signals that have been implicated in D-V patterning and growth, including Wnts [72] and Fgfs [73] as well as other GRN components [40,42]. Boundaries between different D-V domains in the mandibular arch are sharp, which is especially prominent in our live imaging (S2 Fig). How is this sharpness achieved? Our model simulates sharp domain boundaries in the absence of noise. However, when stochastic fluctuations are present the domain boundaries are no longer sharp since all cells are responding independently to the noisy signals. Computationally, we can extend our model to investigate if cell-cell communication and cohesion between cells of the same identity will improve boundary sharpening. Experimentally, we can examine if changes in signaling too small to disrupt overall patterning can nevertheless degrade boundary sharpness, which would suggest that this sharpness comes from the patterning pathways themselves. Although our current analysis of live cell dynamics does not suggest a significant degree of NC cell rearrangements once they have reached the pharyngeal arches, cell sorting contributes to domain boundary sharpening in some contexts, such as the neural tube [74], and automated tracking of large numbers of arch NC cells in the future might reveal subtle but important cell movements in arches as they are patterned. D-V arch patterning of the mandibular arch and its associated GRN are largely conserved among vertebrates [21,23,75]. However, differences in the resulting anatomy of the jaw and skull in the adult, as well as the size, shape, and growth of the mandibular arch primordium make it difficult in some cases to draw clear homologies across species e. g. zebrafish and mouse or human [21,23]. In addition, differences in the timing of gene expression and our ability to identify clear functional homologues across species have made comparisons of the genetic pathways involved challenging. Having made a computational model for the zebrafish jaw that incorporates these spatiotemporal features we can now extrapolate to other species to determine if similar constraints apply. In addition, we can begin to ask questions at a more integrated, systems biology level, about how such changes in size, shape and timing may have arisen and how the GRN has adjusted to compensate for changes in such features as signal propagation and noise.
Proper development of the body requires boundaries to form between regions in which cells will form different structures, and these boundaries need to be properly organized in space. This must occur accurately even in moving, dividing cells and in the presence of the noise that is inherent in all biochemical processes. We use development of the upper and lower jaw as a model to study boundary formation. In this work, we combine detailed experimental measurements with computational modeling to investigate the role the timing of gene expression plays in organizing spatial boundaries, and find that the different orders of gene expression navigate a tradeoff between precision and accuracy in boundary positioning.
fish
medicine and health sciences
engineering and technology
signal processing
vertebrates
animals
noise reduction
simulation and modeling
animal models
organisms
osteichthyes
developmental biology
model organisms
experimental organism systems
molecular development
embryos
research and analysis methods
embryology
animal studies
gene expression
head
morphogens
signal transduction
zebrafish
eukaryota
cell biology
anatomy
jaw
genetics
biology and life sciences
cell signaling
bmp signaling
Modeling craniofacial development reveals spatiotemporal constraints on robust patterning of the mandibular arch