Signal dynamics in Sonic hedgehog tissue patterning

For decades the morphogen hypothesis has helped to explain a wide range of tissue patterning processes (Crick,1970; Turing,1952). The hypothesis states that chemical signals, termed`morphogens', are secreted from signaling centers, and that the resulting static extracellular morphogen concentration gradient emanating from the center spatially organizes and patterns tissue architecture. That is, rapid morphogen transport creates a concentration gradient invariant over the time window of tissue patterning (i.e. at steady state). Sonic hedgehog (Shh),which forms a concentration gradient to pattern the limb bud(Riddle et al., 1993),midbrain (Britto et al., 2002),forebrain (Ericson et al.,1995) and spinal cord (Roelink et al., 1995) during vertebrate development, is ostensibly a canonical example of the morphogen hypothesis. However, very recent studies suggest that Shh concentration gradient dynamics play a crucial role in tissue patterning. Both the time of exposure of a cell to a given Shh concentration(Ahn and Joyner, 2004; Harfe et al., 2004; Kohtz et al., 1998; Park et al., 2004; Wolff et al., 2003; Yang et al., 1997) and the timing of Shh source secretion (Ericson et al., 1996) are crucial determinants of Shh tissue patterning. The classical morphogen hypothesis does not account for such dynamics in gradient formation and cellular response.

Beyond passive diffusion, morphogen systems can have a number of`accessory' mechanisms that modulate ligand transport. Specifically, each transport mechanism potentially modifies not only the steady state concentration gradient, but also the rate of morphogen transport at various times in the patterning process. For example, studies in Drosophilahave identified numerous genes essential for actively transporting morphogens upon their release from secreting cells, and mutating such transporters disrupts tissue patterning (Chen et al.,2004; Han et al.,2004; Takei et al.,2004). Also, high affinity interactions of morphogens with cell surface (Chuang and McMahon,1999) and extracellular matrix (ECM) components(The et al., 1999) serve the putative roles of depleting or immobilizing extracellular-diffusing morphogens to limit long-range signaling. Any of these mechanisms can affect the temporal evolution of concentration gradients in developing tissue, which is overlooked by the focus of the morphogen hypothesis on the steady state concentration gradient.

To investigate temporal effects of transport and signaling, we model Shh regulation of dorsoventral spinal cord patterning in chick embryonic development stages 10-26 [∼33-116 hours after egg laying(Ricklefs and Starck, 1998)]. Shh, secreted from the floorplate, diffuses into the neural tube(Roelink et al., 1995), and,as its concentration decreases within the tissue from approximately 15 nM at the floorplate to 0.5 nM at dorsal edge of the neural tube, target cells switch at threshold values from mature ventral to dorsal phenotypes [e.g. from V3→MN→V2→V1 in Fig. 1A, as has been previously reviewed(Persson et al., 2002)]. Our model begins after the neural fold appears and essentially as the neural tube closes at stage 10, when Shh is first secreted from the floorplate. We subsequently track the cell fate switch between V3 interneurons (V3) and motoneurons (MN) that occurs through stage 26.

The structure of Shh, as well as its various interacting proteins, has complicated a simple understanding of how its transport establishes a gradient during this process. Shh is covalently modified by hydrophobic moieties,including a C-terminal cholesterol and a N-terminal palmitic acid(Pepinsky et al., 1998; Porter et al., 1996), which may anchor the ligand to cell membranes and thereby significantly reduce its diffusivity. However, Shh is still capable of signaling at a large distance,up to 20 cell diameters, away from its source. In addition, the Shh receptor Patched (Ptc) is upregulated by Shh signaling, and its subsequent binding and receptor-mediated internalization of Shh depletes the ligand from the extracellular space (Chen and Struhl,1996; Marigo and Tabin,1996). Furthermore, Shh and its Drosophila homolog Hedgehog can form multimers, and the transmembrane protein Dispatched (Dis) is likely to be involved in regulating their assembly and intercellular transport(Kawakami et al., 2002). Moreover, a membrane glycoprotein, Hedgehog-interacting protein (Hip), binds Shh with high affinity to modulate its signaling activity(Chuang and McMahon, 1999). ECM proteins also regulate Shh transport, as high-affinity binding of Shh to vitronectin in the neural tube has been suggested to aid in the proper presentation of Shh to differentiating motoneurons(Pons and Marti, 2000). Finally, the effective transport of Drosophila Hedgehog depends upon the activity of heparan sulfate proteoglycans (HSPG)(Bornemann et al., 2004; Takei et al., 2004; The et al., 1999), and Shh has also been shown to bind HSPG (Rubin et al., 2002). The individual and synergistic contributions of each of these highly complex elements to the ability of Shh to pattern tissue are unclear. Shh transport via diffusion was previously modeled in the vertebrate limb bud, using a simple signal transduction mechanism without consideration of these accessory transport mechanisms(Dillon et al., 2003). Therefore, to complement, synthesize, and guide experimental work, we have applied a systems biology analysis to explore the effects of diffusion,receptor-ligand dynamics and gene regulation dynamics on Shh gradient formation and tissue patterning.

We build upon a previous single cell model(Lai et al., 2004) to analyze Shh patterning of the developing spinal cord and find that the concentration gradient initially established by diffusion can be modified by several long timescale mechanisms, including morphogen binding to ECM and gene expression. Modeling results are able to reproduce the experimental profiles and clarify the potential roles of six modular mechanisms involved in Shh gradient formation and cellular signaling (Fig. 1B). This investigation suggests that different components can be assembled in a modular fashion to dynamically pattern a morphogen gradient according to the needs of specific tissues.

We consider dorsoventral patterning of the chick neural tube arising from Shh transport in one dimension. The spatial axis of our model tracks the diffusion of Shh away from its floorplate source and through the neural tube,and partitions the neural tube into a mesh of discrete 10 μm cubic elements, each containing one cell (Fig. 1A). Inside each 10 μm cube, the extracellular space consists of interconnected channels of unspecified geometry, but that in sum occupy 20%of the volume based on empirically measured void fractions of neural tissue(Incardona et al., 2002). The diffusion coefficient was modified to account for tortuosity(Lander et al., 2002). A Shh molecule moving through the continuous extracellular space of the mesh(Fig. 1A) can smoothly diffuse through the extracellular regions in cubes A, B and then C. By contrast, all cell surface and intracellular species (e.g. receptors and receptor-ligand complexes) are completely restricted within their cellular volume compartment,which is surrounded by a small 10 nm plasma membrane element/barrier.

At the ventral-most region of the chick neural tube, high level Shh expression is initiated exclusively in the floorplate(Fig. 1A) during stages 10-12(∼34 hours after egg laying). At this time (t=0), all cells in the tube have the same initial gene expression profile. As time progresses,Shh diffuses dorsally from its floorplate source through the mesh and binds to receptors or other components, and high Shh signal levels induce a cell phenotype switch (Lai et al.,2004). The position of the mature phenotypes seen in wild-type embryos after stages 26 (>80 hours after laying) is shown(Fig. 1A). Recent work in mouse embryos indicates that early Shh secretion from the notochord may diffuse far into the neural tube to affect MN commitment(Jeong and McMahon, 2005), and such scenarios in chick embryos could readily be incorporated by adding extra elements for the notochord cells (at x<0), with appropriate Shh secretion dynamics.

The Shh signaling network is represented as a set of differential equations that track the rates of change in the concentrations of network constituents,and whose individual terms represent rates of diffusion, protein synthesis and degradation (Fig. 1B, Fig. 2A). At the single cell level, we build upon the Shh signaling network derived by previously(Lai et al., 2004) to include cellular internalization effects (see Figs S3, S4 in the supplementary material). As Shh is increased above a threshold concentration, it stimulates Gli production to the point where Gli positively feeds back upon its own expression and rapidly switches the state of the network to `on'(Fig. 2B). The activities of Gli2, which overlap with those of both Gli1 and Gli3, are highly context dependent, and its molecular interactions in the neural tube progenitor cells require further characterization (Bai et al., 2004; Ruiz i Altaba,1999). As a result, we have effectively parsed the effects of Gli2 into two types: either a pure transcriptional activator, the `Gli1' type, or a transcription factor of both repressor and activator functions, the `Gli3'type. Thus, in the model, the effects of Gli2 are effectively lumped into a Gli1 term and a Gli3 term. We report the Gli1 concentration as the important system output, as the on/off gli1 expression interface demarcates the V3/MN boundary.

Kinetic, diffusive and binding parameter values were either directly taken from the literature or estimated based upon analogous biological systems(Fig. 2 legend, see also Table S1 in the supplementary material). Parameter estimates were chosen to meet the three following experimental observations: switching `on' of homeodomain nkx2.2 (which serves as the gli1 domain marker in the case of the model) expression at a ∼3 nM Shh threshold at steady state(Ericson et al., 1995; Ericson et al., 1997b);∼50-hour kinetic timescale of Shh secretion from the floorplate [based upon timescales for MN specification from figure 3D in Ericson et al.(Ericson et al., 1997b)]; and a Nkx2.2 protein fluorescence intensity spatial profile in a wild-type chick embryo [see figure 3B in Ericson et al.(Ericson et al., 1997b)]. To satisfy the last criterion, because the Shh secretion rate from the floorplate has not been quantitatively determined, we chose it such that our pattern matched the Nkx2.2 switching interface seen at 70 μm from the floorplate(Ericson et al., 1997b). For the ventral-most cells (close to x=0), floorplate induction, marked by an increase in hnf3β, occurs above a 10 nM Shh concentration(Briscoe et al., 2000). Such floorplate induction, probably due to additional downstream targets of gli1 or other signals not included in this model, accounts for the decrease in nkx2.2 expression seen experimentally in and near the floorplate (Fig. 3B). For each parameter in the core-signaling pathway, we conducted sensitivity analysis for parameter values over four orders of magnitude to observe whether the single cell response to Shh varied (Fig. 2C). Model behavior was investigated for parameter values that changed the 3 nM Shh switching threshold at steady state from 1 to 10 nM and pattern evolution time from 30 to 150 hours Shh secretion, and all conclusions and trends discussed below remained qualitatively the same.

The set of equations in Fig. 2A are presented in dimensional form. However, relationships between groups of variables can be intuitively easier to interpret than individual parameters. In addition, grouping variables reduces the number of independent parameters that are necessary to describe the model. As a result,the following types of variables were non-dimensionalized by the corresponding parameters: concentrations by K_Gli3, space by 1 cm, and time by 1/k_deg. The corresponding non-dimensional equations are shown in Figs S3, S4 (in the supplementary material) and were coded into the FEMLAB software.

Dynamic expression of transcription factors (e.g. Gli family) governs the V3/MN/V2/V1 pattern of the neural tube(Fig. 1A), a process experimentally analyzed by staining for numerous homeodomain proteins(Briscoe et al., 2000; Jessell, 2000), among the earliest of which are the markers Nkx2.2 in the V3 region and Pax6 in the remainder of the tube. In this work, we focus on the ventral-most binary cell fate switch delineating the V3/MN boundary, which results in the Nkx2.2/Pax6 histological demarcation. Homeodomain protein expression patterns are also regulated by BMP, FGF and retinoid signals (reviewed by Jessell, 2000), events that can be incorporated into the model in future work to evolve the model behavior from a binary switch into a ladder of cell fates. Given that the V3 domain and Nkx2.2 expression are entirely lost in mutants with compromised Shh signaling[e.g. Shh^-/- (Litingtung and Chiang, 2000), Smo^-/-(Wijgerde et al., 2002),Shh^-/-/Gli3^-/- (Bai et al., 2004; Litingtung and Chiang, 2000), Gli3^-/-/Smo^-/-(Wijgerde et al., 2002), and Gli1^-/-/Gli2^-/-(Park et al., 2000) mice], we use Nkx2.2 purely as a robust marker of cells where Shh signaling is active,and where Gli is presumably expressed. Also, because Nkx2.2 is rapidly upregulated following Shh signaling, we assume that it can be used as a marker of Shh signaling activity, without requiring knowledge of the mechanistic and molecular interactions that connect Gli to Nkx2.2.

To analyze the dynamics of the V3/MN binary cell fate switch, we first built upon a single cell Shh signaling network model to include cellular internalization of the Shh receptor Ptc and Ptc-Shh complexes(Fig. 2A, see also Figs S3, S4 in the supplementary material). As a result of a positive-feedback loop,arising from Gli1-binding sites within its own promoter, and a negative-feedback loop via Ptc upregulation, this signaling network exhibits a robust on/off switch in Gli1 expression at steady state(Lai et al., 2004). However,we found that the Shh signaling network also undergoes robust switching behavior before reaching steady state, particularly within the 3-day time window for V3/MN pattern completion in chick embryos(Fig. 2D). In contrast to the canonical morphogen hypothesis, the tissue may not need to wait for the morphogen concentration to reach a steady state in order to undergo patterning. That is, as soon as morphogen signaling increases the expression of a target transcription factor above a key threshold value, it can induce a cell phenotype switch. Quantitative results from other systems indicate that as little as a twofold increase in the concentration of a transcription factor can switch cell fate (Niwa et al.,2000; Shimizu and Gurdon,1999), indicating that it is possible or even likely that cells begin to interpret the Shh gradient information before the gradient reaches steady state. In support of this hypothesis, numerous recent studies have shown that not just the level but also the duration of Shh exposure is an important determinant of cell response (Ahn and Joyner, 2004; Harfe et al., 2004; Kohtz et al.,1998; Park et al.,2004; Wolff et al.,2003; Yang et al.,1997). The model supports these studies. If we conservatively assume that a cell can switch phenotype once it achieves a sevenfold increase in Gli1 concentration from its initial basal concentration, the model predicts that a cell exposed to different Shh concentrations will switch fate after different durations (Fig. 2D,results were qualitatively similar for other threshold increases such as threefold to tenfold, data not shown). In reality, as the gradient evolves,cells in the neural tube will be exposed to a complex sequence of Shh concentrations that will drive cell fate changes (i.e. the V3/MN boundary) at different times after exposure, and probably before Shh reaches a steady state gradient.

For all subsequent analysis, we analyze the Shh gradient at 83 hours, the time at which the neural tube has been experimentally shown to have a mature V3/MN demarcation, but before the morphogen gradient has achieved steady state. Note that a higher Shh concentration threshold must be present to switch cell fate within early time windows than if the tissue were allowed to proceed to steady state (5.7 nM for a switch at 43 hours, 4.5 nM for 83 hours,and 3.5 nM for steady state, respectively, in Fig. 2D). Accessory transport mechanisms will modulate cell phenotype patterning by dynamically varying the Shh concentration a cell is exposed to during this developmental time window.

To gain comprehensive spatial as well as temporal insights into how the Shh signal propagates in the neural tube, we dynamically tracked the concentrations of Shh network constituents in a multicellular, spatial model. For the core Shh signaling pathway, which includes intracellular signal transduction and ligand internalization (within the dashed line in Fig. 1B), the simulation correctly reproduced a wild-type tissue pattern with distinct regions of gli1 `on' and `off' cells (Fig. 2; see also Movies 1-8 in the supplementary material).

During this patterning process, extracellular Shh rapidly built to a high concentration within 5 hours, but then fell to approach steady state levels at 60-80 hours (Fig. 2E). Note that Shh rapidly (t<5 hours) rose well above the static levels found to induce a MN to V3 phenotype switch [∼3 nM(Ericson et al., 1997b)]. Because protein expression from Shh target genes does not build appreciably within this relatively short timescale, receptor-ligand internalization and passive diffusion alone governed the early evolution of the Shh concentration profile. Shh levels then more rapidly declined as Ptc, a direct Shh transcriptional target, increased in the ventral-most portions of the embryo and began to mediate Shh degradation via receptor-mediated endocytosis(Fig. 2F). In contrast to the more continuous Shh profile that is smoothened by transport, Ptc profiles exhibited discrete on/off regions due to induced ptc expression above the Shh switching threshold (Fig. 2D,E). Near the floorplate, free Ptc initially decreased as extracellular Shh levels rapidly elevated and bound free Ptc. The resulting Ptc-Shh complexes were rapidly internalized but degraded more slowly, so that the predominant form of Ptc was an internal, complexed form(Fig. 2G), consistent with previous observations (Incardona et al.,2002; Incardona et al.,2000). After 20 hours, ptc was highly upregulated due to high Shh signaling near the floorplate.

Shh also drove a dynamic Gli expression pattern. Cells at the interface exhibited a transient increase followed by a decrease in Gli1 concentration,and after 83 hours of Shh secretion, an 8- to 10-fold higher Gli1 protein level was seen for the ventral `on' cells versus the dorsal `off' cells(Fig. 2H). All concentrations in the neural tube approach but do not reach steady-state levels in the developmental time window of 3 days for wild-type chick V3/MN patterning(Fig. 2E-H). The resulting intricate concentration profiles of all the core pathway components at various times during patterning would be difficult to predict intuitively due to the multiple mechanisms occurring at many different timescales: ligand diffusion over a cell diameter over tenths of seconds [(10×10^-4cm)²/1×10^-7 cm²/s], ligand binding and internalization on the order of minutes (1/K_shhk_on and 1/k_Cin; see Table S1 in the supplementary material), and gene expression and protein degradation over a number of hours (1/k_deg,1/k_Pout and 1/k_Cdeg).

As described above (Materials and methods/Results), sensitivity analysis of parameters was conducted, and values were chosen to satisfy Shh threshold switching levels and gene expression spatial profiles seen experimentally. We next tested the ability of the model to predict experimental results by simulating a previous study in which the neural tube was transfected to express a dominant-negative Ptc mutant that lacks the capacity to bind Shh yet still inhibits Smo signaling (Briscoe et al., 2001). In this study, embryos at stages 10-12 were transfected on one side of the neural tube and then imaged at stages 20-24(t=39-63 hours). To simulate the transfection in Fig. 3A, cells between 40-60μm and 80-110 μm were supplied with a signaling-defective Ptc by making Smo activity insensitive to Ptc in the differential equations governing intracellular signaling [i.e. by setting the ratio(K_Gli3/K_ptc)→∞ in Fig. S2A in the supplementary material]. Briscoe et al. found that transfected cells (green in Fig. 3A) lacked nkx2.2expression (Briscoe et al.,2001), whereas the contralateral, untransfected side exhibited the previously described wild-type pattern(Ericson et al., 1997b). The model captured the effects of the mutant Ptc on pattern formation, as the Shh signaling range increased in both the experimental and modeling results for the transfected side (transfected range ∼90 μm>∼70 μm for wild type in Fig. 3B). The transfected cells cannot sense Shh signal, and therefore do not upregulate ptc, which on the untransfected side serves as a barrier to Shh transport.

After the model had successfully reproduced patterning perturbed by dominant-negative Ptc, it was used to analyze the effects of several accessory mechanisms likely to have a strong effect on Shh transport and tissue patterning (I-VII in Fig. 1A). These mechanisms operate over a wide range of time scales, ranging from transport-hindering mechanisms that exert rapid effects on the developing Shh gradient to negative-feedback loops operating at long timescales. First, the Hip negative-feedback loop is shown to restrict the range of Shh patterning by acting on both short and long timescales. Then, we analyze various mechanisms that retard or promote Shh diffusion, which leads us to classify patterning into three distinct regimes.

Hip is a transmembrane glycoprotein that functions as an inducible antagonist of Shh signaling, because it is a non-signaling transcriptional target of Shh signaling that binds and sequesters Shh(Chuang and McMahon, 1999; Jeong and McMahon, 2005). We simulated Gli upregulation of Hip and allowed it to bind Shh reversibly, as well as to undergo internalization to accelerate Shh degradation(Fig. 1B, II). Ventral shifts in the wild-type pattern occur when Hip is added to the model(Fig. 4). By sequestering extracellular Shh, Hip acts as a `shunt' to remove free Shh from the tissue. Therefore, as the initial (t=0) Hip concentration(Fig. 4A) or the maximal Hip synthesis rate (Fig. 4B) was increased, the extracellular Shh concentration progressively decreased and shifted the interface ventrally, consistent with recent experimental results(Stamataki et al., 2005).

The most unanticipated result of our simulations was that mechanisms that slow ligand transport could actually induce phenotype switching deeper into the tissue. For example, at very low Shh diffusion coefficients(Fig. 5), the signal did not diffuse far from the floorplate, yielding only a small region in which the cell fate switched. However, as the diffusivity was gradually increased, the morphogen concentration rapidly increased within a large region near its source, leading to a deeper patterning(Fig. 5A,B). However, once the Shh diffusion coefficient reached a high value [>10^-6cm²/s], the morphogen rapidly diffused away from its source to establish a very shallow morphogen gradient where the concentration at most points was too low to change cell phenotype(Fig. 5A). The result is biphasic behavior, where the pattern interface position (i.e. Gli1 level)undergoes a maximum value as a function of a parameter such as Shh diffusivity(Fig. 5B,C). Patterning at low diffusion constants constitutes a `signal accumulation regime', whereas high diffusion constants yield a `signal dispersal regime'. Note that small changes in diffusivity do not influence the interface position near 10^-7cm²/s, an estimate for the diffusion constant of Shh in solution. The Shh signaling system therefore may have evolved to be robust in this regime and to consistently pattern an interface at 70 μm. This biphasic behavior is observed at t=83 hours as well as at steady state.

The transmembrane protein Dis interacts with cholesterol-modified Hh(Kawakami et al., 2002) and has been hypothesized to be involved in packaging lipid-modified Shh into clusters (Zeng et al., 2001). Dis is incorporated in our model as a catalyst for complexing membrane-bound Shh into freely diffusing aggregates [six Shh units large(Zeng et al., 2001)], which can then more readily diffuse dorsally and induce signaling(Fig. 1B, III). To mimic membrane-bound Shh, we reduced the diffusivity of monomeric Shh to 10^-10 cm²/s and observed the effects of Dis concentration and aggregate diffusivity on patterning. Similar to the Shh diffusivity results, varying both of these parameters yielded a biphasic response (Fig. 6A,B). At low Dis concentration (10^-4-10^-1 nM), a sufficient aggregate of high diffusivity was formed to signal deeper into the tube. However, at high Dis concentration (>10¹ nM), most Shh was absorbed into the hexameric form. The resulting net decrease in the number of signaling molecules therefore counteracted the benefits of the enhanced diffusion of these aggregates and shifted the interface ventrally. In addition, the biphasic response to Dis concentration(Fig. 6A) and aggregate diffusivity (Fig. 6B) occurs for the same reasons as for monomeric Shh diffusivity(Fig. 5).

As Shh diffuses, it can encounter various components of the ECM. Because such interactions have been proposed to modify the neural tube pattern, we analyzed the effects of reversible Shh binding to one constitutively expressed ECM component, HSPG [mechanism IV in Fig. 1B (Gould et al.,1995)]. A biphasic response was observed as a function of HSPG extracellular concentration (Fig. 7A,B). High HSPG concentrations act as a high-capacity morphogen`sponge' that prevents free Shh from building up to levels sufficiently high to effectively signal. By contrast, at very low HSPG concentrations, Shh diffusion is unhindered, and the morphogen rapidly spreads to large distances,leading to a dilution of the factor to levels too low to effectively signal(as in Fig. 5). However, at intermediate ranges of HSPG, this ECM component can concentrate Shh within a relatively broad region near its source, leading to a dorsal shift in the interface. Adding immobilized HSPG to the system [5-30 μM] therefore functions analogously to reduce Shh diffusivity and, paradoxically, increases the range of effective Shh signaling by hindering its transport. However, for an estimated Shh affinity for HSPG of 350 nM(Loo et al., 2001), relatively high levels of HSPG are required to modify the pattern.

We have simulated Shh transport and signaling in the neural tube to investigate the relative impacts of free diffusion, binding to cell surface and ECM components, intracellular trafficking, and aggregate formation on developmental pattern formation. The influence of a number of parameters and mechanisms associated with the Shh signaling system on patterning dynamics were analyzed during a crucial ∼3-day developmental window (stages 10-26)(Ericson et al., 1996; Ericson et al., 1997b; Roelink et al., 1995).

With the core-signaling pathway (dashed circle in Fig. 1B), the simulation predicts a highly dynamic profile for both intracellular and extracellular Shh network constituents. Full knowledge of such profiles can help to both interpret and guide experimentation. Initially, Shh levels rapidly escalate and extend deep into the tissue for t<20 hours, followed by substantial reductions at longer times(Fig. 2E). Therefore,experimental analysis of the Shh concentration profile only at later times could miss the early Shh build-up and thereby misinterpret the ligand signaling range. The steady state Shh profile (approximately shown in Fig. 2E for t>60 hrs) has a rapid decay close to the morphogen source and a slow decay further away from the source, consistent with previous theoretical studies of steady state morphogen profiles (Eldar et al.,2003). Note that early in our patterning results (0-24 hours),there is a smooth Gli1, Ptc, Gli3R and Gli3A gradient(Fig. 2H,F, and Movies 1-8 in the supplementary material), but this profile begins to sharpen into a discrete interface at 72 hours (stage 16-26). Quantitative assays of the expression of Gli1, Gli3 or Ptc at a broad range of times from stage 12 through 26 can further test directly whether our framework accurately predicts their expression patterns, and many experimental stains for Gli or Ptc expression at times before steady-state are consistent with smooth concentration gradients in the MN domain(Lei et al., 2004; Stamataki et al., 2005). In particular, snapshots of the Ptc profile in the neural tube indicate a highly dynamic, initially graded profile consistent with our results (see Fig. S5 in the supplementary material), and at later times processes not included in the model (e.g. Shh-independent signals influencing the relatively uncharacterized Ptc promoter) are likely to modulate ptc expression, especially in the more dorsal sections. Recent work supports a role for a graded Gli3 profile in early patterning (Stamataki et al., 2005), although several experimental details preclude a direct, quantitative comparison between these results and our simulations. In particular, endogenous Gli1, Gli2 or Gli3 expression is not directly measured,and exogenously introduced Gli3 is expressed at levels that vary over time. However, this important work provides strong evidence for the role of gene expression dynamics in tissue patterning.

The mechanism by which a single cell interprets a morphogen gradient can occur at the transcription factor level(Niwa et al., 2000; Shimizu and Gurdon, 1999). Quantitative differences in Gli have not been experimentally tested for V3/MN patterning in the vertebrate neural tube; however, given the Shh threshold and kinetic data from chick V3/MN patterning used for our parameter estimates (see Materials and methods), the model predicts that V3 cell fate specification could be established at a sevenfold increase in Gli1 from its initial basal concentration [i.e. a Gli1 concentration at the time of V3 specification of>11 nM (or 7×Gli1 concentration at t=0)]. We tested model behavior for V3 cell fate specification occurring at a range of thresholds from three- to tenfold Gli (4.9-16.3 nM) increases, and all of our conclusions remain qualitatively the same (data not shown). The sevenfold increase above basal/initial Gli expression levels corresponds to a two- to threefold Gli1 difference across the V3/MN interface position at approximately t=50 hours (Fig. 2H). Therefore, the model behavior is consistent with the twofold increase in Oct3/Oct4 expression in embryonic stem cells (Niwa et al.,2000) and the threefold increase in SMAD complexes in a Xenopus blastula cell (Shimizu and Gurdon, 1999) that have been found to trigger cell fate switching. It is notable that certain borderline cells(Fig. 2H) experience transient two- or threefold increases in Gli1, and a deeper investigation of the induction kinetics for the next generations of transcription factors downstream of Gli may reveal whether these cells transiently express MN markers (Ericson et al., 1996)or permanently commit to an MN fate. As we have previously discussed in a single cell model, stochastic effects, which can in the future be incorporated into this spatial model, may account for the transient co-expression of markers of different cell fates (Lai et al., 2004). Finally, future incorporation of more detailed mechanisms of interaction between the transcription factors Gli1-Gli3, Nkx2.2,Pax6 and others, as these interactions are further elucidated, would help update the model to match or predict future patterning results.

Vertebrate Shh patterning can be further complicated by the fact that Ptc is not the only receptor that mediates ligand degradation. Like Ptc, Hip also binds Shh with high affinity and is a Shh transcriptional target, yet Hip mediates Shh endocytosis and degradation without transducing a signal. A shunt in an electrical circuit is an alternate pathway that diverts current away from the remainder of the circuit, analogous to the receptor-mediated endocytosis and ensuing degradation that divert a morphogen from signaling. Intracellular shunting by Ptc and Hip, as seen in the modeling results(Fig. 4A,B), attenuates Shh signaling, consistent with several studies in various organisms(Chuang and McMahon, 1999; Jeong and McMahon, 2005). Negative-feedback loops, which establish shunts via molecules like Hip, limit Shh penetration and can `stabilize' patterning, a mechanism previously proposed for morphogen gradients (Eldar et al., 2003).

Although Hip expression has been detected near all Shh signaling centers,its basal concentration and extent of upregulation upon Shh signaling are both parameters that vertebrates may use to regulate Hedgehog signaling with great spatial precision (Chiang et al.,1999; Chuang et al.,2003; Tojo et al.,2002). For example, Hip prevents the spread of excess Shh ligand beyond odontogenic mesenchyme in tooth development, thus restricting the Shh signaling to specific regions of the oral axis(Coulombe et al., 2004). Other than basal concentration and extent of Hip upregulation, other rates in the Hip pathway may be modulated in different organisms, as in the mouse neural tube where Hip-Shh complex internalization appears to be slow(Chuang and McMahon, 1999; Jeong and McMahon, 2005). The ventralization of the tube observed when Hip is overexpressed(Fig. 4), and the non-cell-autonomous nature of this expansion, is very consistent with recent experimental work (Stamataki et al.,2005). Interestingly, soluble, diffusible forms of Hip have recently been found in the mature brain(Coulombe et al., 2004). Our modeling results suggest that this new mechanism may potentially extend the Shh signaling range by protecting Shh from binding to the cell surface Hip variant, Ptc, or even HSPG (data not shown).

Although receptor binding can restrict the morphogen signaling range, other mechanisms may unexpectedly enhance it. Two ostensibly opposing experimental observations have been difficult to reconcile: the long-range signaling ability of Shh and membrane anchorage of the ligand by hydrophobic modification. Shh associates with the membrane through the addition of two lipid tethers during its synthesis, a N-terminal palmitic acid(Pepinsky et al., 1998) and a C-terminal cholesterol (Porter et al.,1996). Despite the reduced diffusivity accompanying membrane association, many studies have demonstrated the long range signaling ability of Shh in the neural tube over 20 cell diameters (>200 μm)(Briscoe et al., 2001; Ericson et al., 1997a; Gritli-Linde et al.,2001).

The simulation counterintuitively predicts that lower diffusion constants can actually concentrate the signal in the ventral neural tube and thereby extend its signaling range (Fig. 5A). Specifically, a typical 20 kDa protein in solution has a diffusivity of order ∼1×10^-7 cm²/s, and hydrophobic modification is likely to decrease the diffusivity to a range between 10^-8 and 10^-10 cm²/s(Creighton, 1992). This reduction may actually help Shh to extend its signaling range two additional cell layers further from the floorplate(Fig. 5C). In one study in the vertebrate limb, Shh was detected ∼200 μm from the source, whereas knock-in of a non-lipid modified Shh led to ligand detection only at lower levels and much closer to the source (see Lewis et al., 2001). The interpretation was that lipid modification was necessary for long-range transport. By contrast, our results indicate that the rapidly diffusing,non-lipid-modified form may be rapidly diluted within the tissue to fall below the experimental threshold of detection, whereas the lipid modification concentrates the ligand.

Via independent mechanisms, both Dis and ECM can also extend the signaling range of Shh. First, although experiments have yet to quantify a spatial profile for Shh aggregates in the neural tube, our modeling results show that Dis-catalyzed aggregation of highly diffusible Shh aggregates can propagate a Shh signal (Fig. 6A,B). Analogous to the monomeric Shh results(Fig. 5), there is a biphasic response in the V3/MN pattern to both the diffusivity of the aggregate and the Dis-catalyzed rate of aggregate generation. In the neural tube of dismutant mice, ventral fates are not properly specified, and Shh immunoreactivity is detected only in Shh-producing cells in the ventral neural tube, somite and limb (Kawakami et al.,2002), indicating that this mechanism exerts significant control over the range of signaling. Although Dis may have additional functions in the Hedgehog pathway, Dis-catalyzed Shh aggregation may be a general transport modulating mechanism operating in Drosophila(Burke et al., 1999) and zebrafish (Nakano et al.,2004).

Distinguishing among the potential roles of ECM in ligand presentation,stabilization and accumulation has been difficult(Bornemann et al., 2004; Giráldez et al., 2002; Pons and Marti, 2000; Rubin et al., 2002; The et al., 1999). Our results indicate that a very simple mechanism, the reversible binding of Shh to HSPG,can either lengthen or restrict the signaling range depending on the HSPG concentration (Fig. 7A). Both HSPGs and the EXT genes involved in their synthesis are abundantly expressed in a developmentally regulated manner in the mammalian central nervous system,suggesting their functional roles in neural tube patterning(Gould et al., 1995; Inatani et al., 2003; Yamaguchi, 2001). As HSPG is actively remodeled by proteases in many tissues, and the ligand affinity for HSPG can be correspondingly modulated, HSPG concentration and affinity can serve as robust, tunable parameters to modulate ligand signaling in particular tissues. Our results also indicate that vitronectin, a direct Shh transcriptional target, can also modulate ligand transport (see Fig. S2 in the supplementary material).

Simulations of both the Dis-catalyzed ligand aggregation mechanism(Fig. 6) and ECM function(Fig. 7) indicate a similar`biphasic' behavior (as seen in Fig. 5 with Shh diffusivity), where the effective signaling range undergoes a maximum as a function of a key parameter. This behavior occurs because aggregating Shh into freely diffusible compounds allows it to overcome diffusion barriers, whereas ECM components function as additional diffusion barriers to concentrate the signal.

This analysis has led to the identification of three major regimes in which accessory Shh transport mechanisms modulate patterning(Fig. 8). Signal accumulation mechanisms limit Shh diffusion by restricting morphogen accumulation to near the source. By contrast, signal dispersal mechanisms promote high Shh diffusion to provide a shallow signal gradient over the developing field, such that the concentration is sufficiently high to change cell fate only near the source. Signaling range is maximized between these two extremes. Furthermore,shunting mechanisms promote intracellular degradation of the Shh signal across the entire profile to decrease the spatial Shh concentration profile. Finally,caution should be taken in interpreting experimental results, as the intersecting profiles indicate that experimental interpretations are highly dependent on assay sensitivity. For example, immunostaining for Shh with a 2.5 nM sensitivity would indicate that signal accumulation mechanisms result in the deepest Shh transport, whereas an assay with a 0.5 nM sensitivity would indicate that signal dispersal mechanisms give rise to the farthest transport.

Correlating patterning with morphogen gradients only at steady state can overlook many of the processes that contribute to tissue patterning, a hypothesis that emerging evidence in a variety of morphogen systems supports. First, in the vertebrate limb, sensitivity to the time a cell is exposed to a given Shh concentration (Ahn and Joyner,2004; Harfe et al.,2004; Kohtz et al.,1998; Park et al.,2004; Wolff et al.,2003; Yang et al.,1997) indicates that specification clearly does not only occur at or after steady state but throughout the entire morphogen transport process. Second, waves of Shh source secretion in the neural tube(Ericson et al., 1996) were crucial determinants of MN specification, suggesting that the sequence of Shh concentration levels may be important for specification. Third, morphogen degradation during transport can determine its signaling range, as the stability and therefore the signaling range of Nodal are regulated by its hydrophobic modifications, which are intriguingly similar to those of Shh(Le Good et al., 2005). Lastly, because cellular competence to a morphogen, a phenomenon not yet considered in our model, can be transient and regulated by several signals,steady state morphogen gradients established after a cell loses competence are irrelevant to tissue patterning. For all of these reasons, accurate monitoring of cell fate as a function of time is important in studying tissue pattering by morphogens.

Mechanistic models serve several roles in developmental biology: they test which experimental results can be accounted for with all that is currently known about a system, they can guide/suggest future experimentation, and they can test whether several hypothesized mechanisms can account for a novel phenotype. For example, high throughput siRNA approaches have recently identified novel components of the Hh pathway(Lum et al., 2003), and modeling can be used to test their potential mechanisms of action. Other molecules can also be included, such as Ptc2, Megalin, and the You class of proteins, whose precise mechanism in the Hedgehog pathway is not yet known(Carpenter et al., 1998; Ding et al., 1998). For example, our preliminary modeling work does not support a function of Megalin as an active Shh transporter through the cytoplasm(McCarthy et al., 2002), as long-range signaling is not observed, even with extremely fast intracellular transport rates of Shh-Megalin complexes (data not shown). Furthermore, a unique advantage of this finite element numerical approach is that it can incorporate the effects of cell division and death during patterning. Future work can expand the current static model geometry to a `living mesh', so that elements are added (or subtracted) as they arise (or die) in irregular geometries and at various times. Finally, the simulation indicates that the rate of Shh secretion from the floorplate, which feeds the patterning process through the activity of a highly complex promoter not yet modeled(Epstein et al., 1999), is a highly important parameter that should be experimentally measured.

The model guides and predicts numerous additional experiments to analyze Shh patterning in the neural tube, and potentially other tissues where Shh transport is involved. First, knock-in of various forms of Shh with varying extents of lipid modification (Feng et al., 2004) is predicted to pattern to differing depths, as a function of Shh diffusivity (Fig. 5). Likewise, Shh diffusivity can be tuned by injecting soluble HSPG. Similar to experiments in the brain where the spread of molecules with an HSPG interaction domain have been increased(Nguyen et al., 2001), high levels of soluble HSPG would induce Shh signal dispersal(Fig. 8). Second, the secretion rate of Shh can be varied in the floorplate through the use of a regulable promoter, and the model can be used to predict the precise position of the V3/MN interface as more MN are produced at the expense of V3 neurons with increasing secretion. Finally, neural tube explants can be incubated in blocking antibodies against Hip, Dis, HSPG or vitronectin. The concentration of such antibodies would tune the levels of active Hip, Dis, HSPG and vitronectin to the levels shown in Figs 4, 6 and 7, and Fig. S1 (in the supplementary material), respectively, to test the predicted biphasic responses.

In conclusion, we have analyzed how complex mechanisms acting at various times regulate morphogen transport and modulate tissue patterning. The results guide and suggest further experiments on how these mechanisms work in concert to provide robust Shh neural tube patterning. Future modeling work can explore the effects of numerous modular mechanisms in the Hedgehog pathway, and may intriguingly suggest further experiments on how these modules have been evolutionary `plugged' into particular tissues to restrict or propagate the Hedgehog signal when required.

This work was supported by a National Science Foundation Graduate Fellowship to K.S. and a University of California Cancer Research Coordinating Committee Award and NIH NS048248 to D.V.S.