## Abstract

Morphogen-mediated patterning is the predominant mechanism by which positional information is established during animal development. In the classical view, the interpretation of positional signals depends on the equilibrium distribution of a morphogen, regardless of the dynamics of gradient formation. The problem of whether or not morphogen dynamics contribute to developmental patterning has not been explored in detail, partly because genetic experiments, which selectively affect signalling dynamics while maintaining unchanged the steady-state morphogen profile, are difficult to design and interpret. Here, I present a modelling-based approach to identify genetic mutations in developmental patterning that may affect the transient, but leave invariant the steady-state signalling gradient. As a case study, this approach is used to explore the dynamic properties of Hedgehog (Hh) signalling in the developing wing of the fruitfly, *Drosophila melanogaster*. This analysis provides insights into how different properties of the Hh gradient dynamics, such as the duration of exposure to the signal or the maximum width of the transient gradient, can be genetically perturbed without affecting the steady-state distribution of the Hh concentration profile. I propose that this method can be used as an experimental design tool to investigate the role of transient morphogen gradients in developmental patterning and discuss the generality of these ideas in other problems.

## 1. Introduction

A classical paradigm in developmental biology is that cells in a developing embryo or tissue acquire information about their relative spatial location as a function of the local concentration of chemical signals in their environment called morphogens. The key idea is that the interpretation of positional information results in the establishment of gene expression patterns whose boundaries correspond to concentration thresholds of the morphogen gradient [1–3]. A common assumption of the classical morphogen model is that concentration gradients are interpreted approximately at the steady state in order to produce stable patterns of gene expression [4,5]. However, recent studies in some systems have proposed that morphogens are dynamic during much of the patterning process and have raised criticisms to the classical morphogen concept [6–9]. For instance, if gradient formation evolves in time, how and when are concentration thresholds interpreted? Or, more generally, how do morphogen dynamics contribute to positional information? One possibility is that patterns evolve as the gradient develops, giving rise to transient patterns that converge to stable gene expression domains as the gradient approaches the equilibrium [10]. This scenario supports a ‘no-role’ model for transient gradients (figure 1*a*), but this seems not to be always the case. For example, the identity of the most posterior digits in the vertebrate limb is not solely determined by the steady-state concentration of a morphogen, but it also depends on the duration of the signal [11–13]. Furthermore, a similar interplay between signal strength and duration of exposure has been recently shown to operate in the specification of distinct neural fates in the vertebrate spinal cord [14]. These studies propose a model of patterning in which cell fate decisions are determined by integrating the strength of the morphogen signal over time (figure 1*b*). Alternatively, different patterns of gene expression may result from the dynamic properties of a morphogen gradient. For example, we recently proposed that the specification of different spatial domains of expression in response to Hedgehog (Hh) signalling in the *Drosophila* wing disc depend on a spatial shift of the Hh gradient prior to reaching the equilibrium [15]. In particular, we showed that cells exposed only transiently to Hh express a different combination of genes compared with cells constantly receiving the signal or cells that never receive it at all (figure 1*c*). Thus, in contrast to the classical morphogen model (figure 1*a*), patterning is often not exclusively dependent on the steady-state readout of the morphogen concentration, but also on the full history of exposure to the signal (figure 1*b*,*c*). Nonetheless, the problem of how the dynamics of gradient formation contribute to morphogen interpretation has been difficult to demonstrate by classical genetic experiments because mutants that affect gradient dynamics are also likely to perturb the steady-state morphogen distribution. Despite the recent introduction of techniques to quantify and measure the dynamic properties of morphogen gradients in living embryos or tissues [16–18], the identification of mutants that selectively affect gradient dynamics with little or no effect on the equilibrium profile remains challenging.

For over 50 years, mathematical modelling in developmental biology has provided powerful tools to identify potential mechanisms and properties of morphogen gradient formation and their interpretation in developmental patterning [19–24]. However, much of this literature has been built on the assumption that morphogens are interpreted at the steady state and transient gradients are often ignored. Reducing a dynamical system to its steady state is usually mathematically convenient as the equations of the model often simplify considerably and theoretical tools are readily available to study the stability of solutions. In contrast, less analytical tools exist to investigate the temporal evolution of a gradient and studies that have taken morphogen dynamics into consideration usually rely on numerical simulations [25,26]. In this paper, I present a general theoretical approach to study the role of the dynamics of morphogen gradients in developmental patterning. Given a mathematical model of a morphogen-based pattern formation system, we consider parameter perturbations, referred to here as *steady-state invariant perturbations*, that may change the transient shape of a morphogen without effectively affecting its steady-state distribution. We introduce a theoretical framework to study these perturbations and discuss how they could be used as a design tool for genetic experiments that permit us to assay the role of gradient dynamics in the specification of patterns of gene expression. As a case study, the dynamic properties of the Hh gradient in the *Drosophila* wing disc are systematically investigated to illustrate how this approach can provide insights into the roles of signalling dynamics during development.

## 2. Results

### 2.1. Theoretical framework: definitions and examples

In order to introduce the concept of steady-state invariant perturbations, consider a simple example, namely a naive model of a one-dimensional morphogen produced at a point source (*x* = 0) that establishes a concentration gradient by diffusion and linear degradation,
2.1where [A] denotes the concentration of morphogen A, *δ*(*x*) is the Dirac delta distribution, and *α*, *β* and *D* are the production, degradation and diffusion rates of A, respectively. At the steady state, the shape of the [A] gradient is exponential^{1} and given by
2.2

In equation (2.2), *A*_{0} and *λ* represent the amplitude and characteristic length of the steady-state gradient, respectively, and their values determine uniquely this solution. As *A*_{0} and *λ* are defined in terms of the parameters of the system, a perturbation on the wild-type parameter values will maintain the steady-state solution invariant if and only if the values of *A*_{0} and *λ* remain unchanged. To formalize, let be the wild-type parameter values and consider a mutant that perturbs the system such that the effective parameter values of the system change to (*α*′, *β*′, *D*′). We say that the mutant preserves the steady-state gradient or is *steady-state invariant* if the following equations hold:
2.3

The set of all steady-state invariant perturbations can be represented geometrically in parameter space and will be referred to as the *steady-state invariant set* of the system. In general, a steady-state invariant perturbation Δ can be denoted as a vector in parameter space that is based on the point of wild-type parameter values and ends on another point of the steady-state invariant set (figure 2*a*). In this example, it follows that every perturbed parameter vector (*α*′, *β*′, *D*′) satisfying equation (2.3) must be of the form
2.4

Equation (2.4) represents a straight line in parameter space that crosses through the origin and contains the wild-type parameter vector (figure 2*b*). Thus, the line defined by equation (2.4) is the steady-state invariant set of this system. Importantly, a steady-state invariant perturbation may affect the transient evolution of the gradient but always maintains the steady-state morphogen gradient unchanged. In this simple example, the effects on morphogen dynamics along the steady-state invariant set are simple to interpret; for *δ* < 1, gradient formation is slower than in the wild-type case, while for *δ* > 1 the steady-state gradient is approached faster than in the wild type (figure 2*c*–*e*). In addition, the rates in which the steady state is approached are space dependent; cells adjacent to the source approach equilibrium faster than cells away from it, and this property also holds after steady-state invariant perturbations (see the electronic supplementary material, figure S1). Despite the simplicity of this example, the definition of steady-state invariant sets and the method to compute them can be generalized (box 1).

Consider the following general model of developmental patterning. Let **G** = ([g_{1}], [g_{2}], … , [g_{k}]) a vector denoting concentrations of gene products g_{1}, … , g_{k} and assume that the dynamics of gene concentrations are described by a reaction–diffusion equation of the form
B 1

In this equation, ∇^{2} denotes the Laplacian operator (), **D** is a vector of diffusion coefficients, ** f** is the reaction function that describes the interactions of gene products, and

*μ*is a vector of kinetic rates or parameters of the system. Note that the systems considered in the text (equations (2.1) and (2.5)) are particular cases of equation (B 1). At the steady state, we set the time derivatives to zero so that the steady-state concentrations,

**G**

^{ss}, obey the following equation: B 2

Assume that the solution of (B 2) exists and let us denote it as . A parameter perturbation (e.g. a genetic perturbation) of the system is a function Δ of the form Δ(*μ*) = (*δ*_{1}*μ*_{1}, *δ*_{1}*μ*_{2}, … , *δ*_{r}*μ*_{r}) ≡ (*μ*′_{1}, *μ*′_{2}, … , *μ*′_{r}) for some positive constants *δ*_{1}, … , *δ*_{r}. We are interested in the study of perturbations that leave the steady-state solution unchanged at least in a region of space **S**. A *steady-state invariant perturbation* is a parameter perturbation Δ that satisfies the following property:
B 3

Our goal is to find the set of steady-state invariant perturbations satisfying equation (B 3). For this purpose, it is useful to consider the following geometric representation. Note that there is a one-to-one correspondence between parameter perturbations and points in parameter space (i.e. points of the form (*μ*′_{1}, *μ*′_{2}, … , *μ*′_{r}) with *μ*′_{i} ≥ 0, i = 1, … , r). For each **x** (fixed), and, therefore, we can consider the points (*μ*′_{1}, *μ*′_{2}, … , *μ*′_{r}) in parameter space that satisfy
B 4

Denote by *Ω*_{x} the set of points in parameter space satisfying (B 4) for a given **x**, and define the *steady-state invariant set*, , of the system (in **S**) as the set that results from the intersection of *Ω*_{x} for all , i.e.
B 5(see the electronic supplementary material, figure S2). Hence, the steady-state invariant set is a geometric representation of the set of mutants that leave invariant the equilibrium gene concentrations in cells located within the region **S**.

### 2.2 Hh signalling in the *Drosophila* wing disc: a case study

With the aim of motivating the applications of the concept of steady-state invariant sets to cases of biological interest, I considered a model of Hh signalling in the *Drosophila* wing imaginal disc (figure 3*a*). In the developing wing disc, Hh is secreted by cells in the posterior compartment and forms a concentration gradient into the anterior compartment. In the absence of Hh, the Hh receptor, Patched (Ptc), represses signalling activity in anterior cells. Upon Hh reception, inhibition of Ptc on the signal transduction component Smoothened (Smo) is alleviated. Although the molecular details by which Ptc inhibits Smo and the mechanisms downstream of Smo that result in the interpretation of Hh signalling activity remain largely unclear, experimental data support a phenomenological model which suggests that Hh signal activity is a function of the ratio of bound to unbound Ptc [27].

A general property of the Hh signalling pathway is the transcriptional upregulation of *ptc* in response to Hh activity that results in increased Ptc expression in cells abutting the anterior–posterior (AP) boundary. This feedback property of the Hh signalling pathway is necessary to limit the range of Hh signalling activity and has been implicated in the robustness of the gradient to fluctuations in Hh dosage [23,28]. Furthermore, in a recent study, we provided experimental evidence for a transient expansion (or ‘overshoot’) of the Hh gradient that results from the Hh-dependent Ptc upregulation feedback in the system [15]. We proposed that the Hh overshoot is required to establish three distinct domains of gene expression (figure 3*b*), supporting the importance of morphogen gradient dynamics in this system. However, the problem of how specific properties of this transient overshoot (e.g. its duration) contribute to patterning has not been studied in detail. Here, I use the concept of steady-state invariant perturbations to make predictions about which genetic perturbations of the system affect the dynamic properties of this overshoot without modifying the equilibrium distribution of the Hh gradient. We consider the mathematical model of the Hh signalling pathway (figure 3*a*) as originally presented in our previous work [15]:
2.5where [Hh], [*ptc*], [Ptc] and [Hh_Ptc] are the concentrations of Hh, *ptc* (mRNA), Ptc (protein) and the Hh–Ptc complex, respectively. The coefficients *α*, *β*, *γ* and *T* represent the rates of synthesis, degradation, complex formation and translation, respectively. *S*^{+}(*x*) (or *S*^{−} (*x*)) is a step function of the form *S*^{+}(*x*) = 1 if *x* > 0 (or *S*^{−}(*x*) = 1 if *x* < 0), and zero otherwise, that is used to represent compartment-specific reactions. For example, Hh is produced exclusively in the posterior compartment (*x* < 0), but the signal can only be activated in the anterior compartment (*x* > 0). The variable [Signal] denotes the concentration of a factor that models the activity of the pathway at the intracellular level.

At the steady state, algebraic substitutions reduce the set of equations (2.5) to the following second-order equation, which is valid exclusively in the anterior compartment (see the electronic supplementary material, appendix):
2.6with , and . This equation is nonlinear and an analytical solution may be difficult (or impossible) to obtain. However, equation (2.6) can be linearized to obtain approximate solutions in regions where Hh is present at sufficiently high () or sufficiently low () levels, i.e. near and far from the AP boundary, respectively (see the electronic supplementary material, appendix, for details). In each of these regions, equation (2.6) reduces to the form,
2.7with *B* and *C* coefficients that result from the linearization. Using appropriate boundary conditions in the regions where the approximations are valid (see §5), it is possible to obtain analytical solutions of equation (2.7). A summary of the explicit solutions of equation (2.7) near or far from the AP boundary is included in table 1 (see the electronic supplementary material, appendix, for details).

In a similar manner that *A*_{0} and *λ* in equation (2.2) represent the amplitude and characteristic length of the gradient, respectively, it is useful to provide a physical interpretation of the parameters in the solutions of equation (2.7). We note that, assuming [Ptc]_{SS} is approximately constant (), the steady-state version of the first equation in the set of equations (2.5) can be written in the form of equation (2.7) with and *C* = 0, where represents the effective degradation rate of Hh. Therefore, our approximate solutions near and far from the AP boundary are equivalent to assuming that [Ptc]_{SS} is constant in each of these regions and provide an estimate of these constants as a function of the kinetic parameters of the model (table 1). Furthermore, these approximations permit analytical expressions to be obtained that define the steady-state invariant set of the system. In particular, parameter perturbations that maintain each of the solutions of equation (2.7) unchanged will maintain the global steady-state distribution of the Hh morphogen approximately invariant in much of the anterior compartment. Thus, the steady-state invariant set of the linearized model is given explicitly by
2.8where the primed symbols represent perturbed parameters and the WT subscripts denote constants obtained by substitution of the wild-type parameter values into the coefficient expressions in table 1. Therefore, any set of perturbed parameter values that satisfy the set of equations (2.8) will maintain the approximated steady-state gradient unchanged. However, the steady-state invariant set as described by equations (2.8) is not very useful in practice because it involves a very large number of parameters. While steady-state invariant sets in high-dimensional parameter spaces are difficult to analyse, it is possible to study particular parameter subsets that satisfy equations (2.8) in lower dimensions. As a first case, we consider the following perturbation:
2.9with denoting the vector of wild-type parameter values and all other parameters kept unperturbed (i.e. at their wild-type values). By substitution into equation (2.8), it is clear that equation (2.9) represents a steady-state invariant subset of the linearized system. In fact, equation (2.9) also leaves invariant the original steady-state equation of the system (equation (2.6)) and, therefore, it is an exact steady-state invariant subset of the system. Geometrically, equation (2.9) represents a line in the *α*_{ptc}_{0} − *α*_{ptc} − *β*_{ptc} parameter subspace (figure 3*c*). Numerical simulations suggest that a few-fold changes in *δ* affect the duration but have little effect on the amplitude of the overshoot (figure 3*d,e*). In particular, we noted that the duration of the overshoot increases rapidly as *δ* decreases (*δ* < 1), but decreases slowly as *δ* increases (*δ* >1). Therefore, perturbations along this steady-state invariant subset provide the opportunity to study the role of signal duration in further detail.

An even simpler case of a steady-state invariant subset of the linearized system can be obtained by noticing that the system of equations (2.8) does not involve the degradation rate of the Hh–Ptc complex, *β*_{Hh_Ptc}. Therefore, under these approximations, the steady-state solution does not depend on the value of *β*_{Hh_Ptc}. For simplicity, consider the approximate steady-state invariant subset given by
2.10with all other parameters held constant at their wild-type values (figure 3*f*). The constraint 0 < *δ* < 1 was enforced so that the assumption remains valid near the AP boundary. In contrast to the subset defined by equation (2.9), perturbations along the subset described by equation (2.10) will affect the amplitude of the overshoot, but its duration will be maintained relatively constant (figure 3*g*,*h*). The simplicity of this subset permits the design of genetic steady-state invariant perturbations in the system.^{2} The *ptc*^{14} allele has been characterized as a mutant that is defective in endocytosis-dependent internalization and degradation of the Hh–Ptc complex, suggesting that the stability of the Hh–Ptc complex is increased in Ptc^{14} mutants [29]. Importantly, other properties of the Ptc protein appear unaffected in *ptc*^{14} mutant clones. For example, Ptc^{14} proteins bind Hh and repress Hh signal transduction in a similar manner to wild-type Ptc [29]. Furthermore, Ptc^{14} is not overexpressed in anterior clones away from the AP boundary, suggesting that unbound Ptc^{14} degradation is normal. This genetic evidence suggests that *β*_{Hh_Ptc} decreases in *ptc*^{14} mutant mosaics, while other parameters approximately maintain their wild-type values. Therefore, we predict that *ptc*^{14} mutants satisfy equation (2.10) and can be considered as a steady-state invariant perturbation of the system. Based on our modelling results, perturbations along this steady-state invariant subset affect the amplitude of the overshoot, but not the overall shape and duration of the transient response, suggesting that more anterior cells would be transiently exposed to the signal (figure 3*g*,*h*). Experiments using *ptc*^{14} mutant clones abutting the AP show that the expression domains of Hh target genes expand in the region of these clones [29], consistent with the predictions of this work (figure 3*g*).

Other subsets can be similarly obtained from equations (2.8) and their effects on transient dynamics can be systematically analysed. However, sets involving many parameters are difficult to visualize and perturbations along these subsets may not be easily achieved by experimental manipulations (see below).

## 3. Discussion

### 3.1. Analysis and computation of steady-state invariant sets

The problem of computing the steady-state invariant set from a given model is in general non-trivial and, with the exception of very simple models (e.g. equation (2.1)), cannot be obtained analytically. Nonetheless, there are different manners to overcome this problem and get a partial or approximate set of steady-state invariant perturbations. The simplest way to obtain at least a steady-state invariant subset of a system is to inspect the steady-state equation and ask which parameter perturbations leave the equation invariant. For example, perturbations that satisfy equation (2.9) leave the steady-state problem (equation (2.6)) unchanged. In these cases, there is no need to obtain the solution of the steady-state equation, but this approach will only provide a subset of steady-state invariant perturbations. We will say that a steady-state invariant set is *direct* if it can be deduced from the steady-state equation. Note that equation (2.4) is a direct subset of equation (2.1), and, from the steady-state solution (2.2), we can show that it contains all the steady-state invariant perturbations of the system. However, this is not necessarily true in general. Therefore, an interesting theoretical challenge for the future will be to find conditions to determine in which cases the full steady-state invariant set is direct, i.e. when we can obtain all steady-state invariant perturbations of the system without needing to solve the equations. Alternatively, steady-state invariant sets can be obtained approximately by linearization (as we did for equation (2.6)), or by spatially discretizing the system, i.e. turning the problem into a system of ordinary differential equations (as described in box 1). Although these approaches do not always result in exact steady-state invariant perturbations, they are generally applicable and often provide perturbations that are apt for experimental design (e.g. equation (2.10)). In addition, in these approaches, the property of steady-state invariance needs not to be global, i.e. it can be considered in a desired subset of the whole developmental field, and this might often be convenient in practice (e.g. when studying local patterning effects using genetic mosaics).

Another technical challenge that deserves further discussion is how to analyse and visualize steady-state invariant sets. This problem is common in practice because steady-state invariant sets are usually contained in high-dimensional spaces. For example, the dimension of the parameter space of our highly simplified model of Hh signalling (equation (2.5)) is 16. More realistic models may involve an even larger number of parameters, making the resulting steady-state invariant set very difficult to analyse and visualize. Our analysis in this paper reduces to the study of subsets (equations (2.9) and (2.10)) that are contained in particular parameter subspaces and result in tractable manipulations, but an interesting future direction is to use theoretical tools (e.g. nonlinear extensions of principal component analysis; see [30] for a review) to reduce the dimensions of steady-state invariant sets in a way that simplifies their use for experimental design.

### 3.2. Steady-state invariant perturbations as a tool to investigate the roles of morphogen gradient dynamics

The main goal of this work is to use it as a tool to design genetic experiments that separate transient from equilibrium effects in developmental patterning. One problem that arises in practice is that steady-state invariant perturbations often impose constraints on many independent kinetic parameters of the system. For instance, in our example of a free morphogen established by diffusion and linear degradation (equation (2.1)), only those perturbations in which all three parameters are altered in the same proportion remain on the steady-state invariant set (equation (2.4)). Despite the simplicity of this example, a genetic perturbation subject to these constraints (i.e. such that equation (2.4) approx. holds) is difficult to achieve experimentally because kinetic rates, in general, cannot be precisely tuned using standard genetic techniques. This example illustrates a fundamental property of this approach as a tool for modelling-based experimental design: simple models that result in relatively simple steady-state invariant sets do not necessarily predict simple genetic experiments.

In practice, it is convenient to know, based on the steady-state invariant set, how difficult it is to design a genetic experiment that causes a steady-state invariant perturbation in the system. In other words, how likely is it for a random parameter perturbation to be steady-state invariant? A geometric notion that provides insights into this question is the concept of co-dimension.^{3} For the practical purposes of this study, we define the co-dimension of a set as follows. Let *M* be a steady-state invariant set that affects *s* parameters and can be parametrized by *r* variables; then the co-dimension of *M* is defined as codim(*M*) = *s* − *r*. For instance, the steady-state invariant subsets defined by equations (2.9) and (2.10) depend on three (*α*_{ptc}_{0}, *α*_{ptc} and *β*_{ptc}) and one (*β*_{Hh_Ptc}) parameters, respectively, and both are parametrized by the single variable *δ*. Therefore, their co-dimensions are two and zero, respectively. Intuitively, the co-dimension of a set is inversely related to the concept of *degrees of freedom*, i.e. the number of parameters that can be arbitrarily varied without abandoning the set. In general, steady-state invariant subsets of small co-dimensions are expected to be more useful in practice because, in these cases, perturbations have the property that independent kinetic parameters are loosely coupled (codim(*M*) > 0) or not coupled at all (codim(*M*) = 0; e.g. equation (2.10)).

As a modelling-based experimental design tool, the approach introduced in this paper also exhibits some technical limitations. One consideration is that steady-state invariant sets depend on the details of the mathematical models. For example, the model of Hh signalling described by equations (2.5) assumes that the signal is activated as a function of the ratio of bound to unbound Ptc [27]. However, this model has been controversial in the literature as it fails to explain experimental data from a study that strongly support that Ptc inhibits Smo catalytically [31]. In order to test whether the prediction that Ptc^{14} is approximately a steady-state invariant perturbation depends on the choice of the model; an alternative model that assumes the inhibition of the signal by the catalytic action of Ptc was considered (see the electronic supplementary material, appendix equation (A 8)). Interestingly, using similar assumptions near and far from the AP boundary, it is possible to show that Ptc^{14} is also a steady-state invariant perturbation in this alternative model (see the electronic supplementary material, appendix, for details). However, this is not necessarily true in general, and these situations need to be considered on a case-by-case basis.

Another technical limitation of these tools is that mathematical models are always based on highly simplified representations of the biological system, so that steady-state invariant perturbations predicted by the theory may not be so experimentally. Therefore, these tools should be used as predictors of steady-state invariant genetic perturbations but appropriate control experiments should be performed to show that the steady-state distribution of the signal is in fact unchanged by the proposed perturbation.

## 4. Concluding remarks

There is much recent interest in the problem of how the dynamics of morphogen gradients contribute to developmental patterning (reviewed in [32]). This work presented a theoretical approach that may help to identify mutations that affect the transient establishment of morphogen gradients without changing their equilibrium distribution. Although the experimental implementation of these tools to study the role of morphogen gradient dynamics in developmental pattern formation is left to future studies, the details of how they can be used in practical cases were provided. As an example, our analysis of steady-state invariant subsets in a model of Hh signalling in the *Drosophila* wing provided the opportunity to decouple two properties of the Hh gradient overshoot, namely duration and amplitude. The duration of the overshoot can be modulated in a steady-state invariant manner by varying the rates of *ptc* production and degradation (figure 3*c*–*e*). Furthermore, the amplitude of the overshoot can be affected by reducing the rate of Ptc-dependent Hh internalization and degradation. In fact, previous genetic experiments in which a mutant form of Ptc (Ptc^{14}) that is defective in Ptc-mediated Hh endocytosis, but otherwise appears to function normally, provide an example of an experimental steady-state invariant perturbation in this system.

The theory of steady-state invariant sets introduced here was motivated by the problem of whether or not morphogen dynamics contribute to developmental patterning. Nonetheless, the applications of the tools presented here are not limited to the design of experimental perturbations and may be broadly applicable to other problems. For example, in engineering, it is often desirable to design systems in which the steady state is robust to a certain class of perturbations. Therefore, the converse problem, namely to construct a dynamical system that will maintain a desired set of perturbations invariant, poses an interesting theoretical challenge for the future. The context in which this converse problem has been studied in some detail is robust control theory, and there are several examples in which imposing invariance into the steady-state values of some variables constrains the dynamic properties of the system. For instance, a well-known result in control theory is that the steady-state error in some desired signal remains negligible to external perturbations, or variations in internal components, if and only if the error dynamics are governed by integral feedback, and this has been shown to apply to perfect adaptation in chemotaxis [33]. Finally, it is likely that evolution has faced and solved this converse problem by selecting for optimal mechanisms of developmental patterning; in some systems, natural selection may have minimized the effects of certain classes of genetic perturbations to maintain invariant some essential gene expression patterns. Alternatively, animal evolution could have employed the diversity of morphogen dynamics to enrich the mechanisms underlying developmental patterning.

## 5. Methods

### 5.1. Wing disc immunostaining

In figure 3*b*, a third instar wing disc of genotype *dpp10 638/CyO* was immunostained with mouse anti-Ptc (Hybridoma Bank Developmental Studies at the University of Iowa) and rabbit anti-β-galactosidase (Invitrogen) antisera following standard techniques. *dpp10 638* is a lacZ enhancer trap reporter that expresses nuclear β-galactosidase under the control of the *dpp* enhancer.

### 5.2. Analytical solutions

The analytical solution of equation (2.1) was derived in Bergmann *et al.* [6] and is given by
with erfc, the complementary error function, . This solution was used to generate the gradient profiles in figure 2*c*–*e*. Approximate analytical solutions of equation (2.6) near and far from the AP boundary (equation (2.7)) were obtained by linearization, and the full details of these approximations are provided in the electronic supplementary material, appendix.

### 5.3. Numerical simulations

Our model of Hh-dependent patterning of the wing disc is assumed one-dimensional along the AP axis. We used a system of coordinates centred at the AP boundary (*x* = 0), where the posterior and anterior ends of the disc were assumed at and , respectively (with *L*_{P} = *L*_{A} = 100 µm, the length of each compartment). Equations (2.5) were numerically solved with the parameter values as shown in the electronic supplementary material, table S1, and using a Forward-in-Time-Centred-in-Space algorithm implemented in Matlab. As in previous work [15], we use zero initial conditions, except for [*ptc*] and [Ptc], which were taken as

Furthermore, we imposed zero-flux boundary conditions at the disc extremes (; ), and assumed continuity of the [Hh]_{SS} profile and its derivative at the AP boundary (*x* = 0; see the electronic supplementary material, appendix, for further details).

## Acknowledgements

I thank John Doyle for his advice and insightful discussions about this work. I also thank Angelike Stathopoulos and Arthur Lander for comments on the manuscript. This work was supported by a grant from the National Institutes of Health (R01-GM078992) to John Doyle.

## Footnotes

↵1 In this example, we ignore boundary conditions by assuming that the system is infinitely long and diffusion occurs in both directions from the source. These assumptions allow us to solve the full time-dependent problem (2.1) exactly using the method of Fourier transforms (see §5) and solution (2.2) arises as a limit of the time-dependent solution when t→∞.

↵2 Strictly, our assumptions hold in each region where the steady-state levels of Ptc are nearly uniform (i.e. near the AP boundary, where Ptc levels are upregulated; and away from the Hh signalling domain, where Ptc is present at basal levels). Thus, steady-state invariant perturbations satisfying equation (2.11) should be valid in much of the anterior compartment; however, there may be a small domain separating these regions where Ptc levels cannot be considered constant and where steady-state invariance cannot be guaranteed.

↵3 Steady-state invariant sets are defined by algebraic equations and therefore the concepts of dimension and co-dimension can be rigorously defined (see the electronic supplementary material, appendix).

- Received December 25, 2010.
- Accepted February 18, 2011.

- This Journal is © 2011 The Royal Society