Morphogenetic action through flux-limited spreading

Available online at www.sciencedirect.com

Physics of Life Reviews 10 (2013) 457–475 www.elsevier.com/locate/plrev

Review

Morphogenetic action through flux-limited spreading M. Verbeni a , O. Sánchez a , E. Mollica b , I. Siegl-Cachedenier c , A. Carleton d , I. Guerrero b , A. Ruiz i Altaba c , J. Soler a,∗ a Departamento de Matemática Aplicada, Universidad de Granada, 18071-Granada, Spain b Centro de Biología Molecular (CSIC-UAM), Cantoblanco, E-28049 Madrid, Spain c Department of Genetic Medicine and Development, University of Geneva Medical School, 1 rue Miquel Servet, CH-1211 Geneva, Switzerland d Department of Neuroscience, University of Geneva Medical School, 1 rue Miquel Servet, CH-1211 Geneva, Switzerland

Received 13 June 2013; accepted 17 June 2013 Available online 25 June 2013 Communicated by J. Fontanari

Abstract A central question in biology is how secreted morphogens act to induce different cellular responses within a group of cells in a concentration-dependent manner. Modeling morphogenetic output in multicellular systems has so far employed linear diffusion, which is the normal type of diffusion associated with Brownian processes. However, there is evidence that at least some morphogens, such as Hedgehog (Hh) molecules, may not freely diffuse. Moreover, the mathematical analysis of such models necessarily implies unrealistic instantaneous spreading of morphogen molecules, which are derived from the assumptions of Brownian motion in its continuous formulation. A strict mathematical model considering Fick’s diffusion law predicts morphogen exposure of the whole tissue at the same time. Such a strict model thus does not describe true biological patterns, even if similar and attractive patterns appear as results of applying such simple model. To eliminate non-biological behaviors from diffusion models we introduce flux-limited spreading (FLS), which implies a restricted velocity for morphogen propagation and a nonlinear mechanism of transport. Using FLS and focusing on intercellular Hh-Gli signaling, we model a morphogen gradient and highlight the propagation velocity of morphogen particles as a new key biological parameter. This model is then applied to the formation and action of the Sonic Hh (Shh) gradient in the vertebrate embryonic neural tube using our experimental data on Hh spreading in heterologous systems together with published data. Unlike linear diffusion models, FLS modeling predicts concentration fronts and the evolution of gradient dynamics and responses over time. In addition to spreading restrictions by extracellular binding partners, we suggest that the constraints imposed by direct bridges of information transfer such as nanotubes or cytonemes underlie FLS. Indeed, we detect and measure morphogen particle velocity in such cell extensions in different systems. © 2013 Elsevier B.V. All rights reserved. Keywords: Morphogenetic gradients; Hedgehog; Flux-limited; Filopodia-like structures

* Corresponding author.

E-mail addresses: [email protected] (A. Carleton), [email protected] (I. Guerrero), [email protected] (A. Ruiz i Altaba), [email protected] (J. Soler). 1571-0645/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.plrev.2013.06.004

458

M. Verbeni et al. / Physics of Life Reviews 10 (2013) 457–475

1. Introduction Patterning of animal tissues by signals exerting concentration-dependent effects is a central mechanism for the generation of cellular diversity. Such morphogen gradients have been thought to be formed and act through passive linear or normal diffusion, oozing from the source. All previous morphogen models assume this behavior, which, from a microscopic point of view, necessarily implies Brownian motion (as the weak limit of random walks), and are based on the pioneering work of Turing, Crick and Meinhardt [1–3]. In these models, randomly diffusible ligands create a spatial instructive gradient, as in the French Flag model of Wolpert [4], and have been applied to Hedgehog (Hh) [5–7], FGFs, TGFβ, Nodal an other morphogens [8–13]. In general, the attraction of such models relies on their simple reaction–diffusion equations and on the fact that they can strikingly recapitulate natural patterns, such as spots on snail shells and stripes on zebra skins [14] (see also the historical review [15]). Depending on the research area, the equation stemming from a continuous limit of molecular Brownian motion has different names: Fourier law of heat equation – in thermodynamics and in mathematical physics, Fick’s law – in biological processes, normal diffusion – in physical chemistry, and mostly frictionless Fokker–Planck or linear diffusion – in physics and mathematics. It is a mathematically provable fact that for initially confined sources of diffusing substances, the corresponding solution to linear diffusion is strictly positive for any point and for any time independently on the value of the diffusion coefficient (see for instance the classical reference [16]). This property is known as infinite speed of propagation or instantaneous spreading and was already pointed out in the pioneer work of Einstein over a century ago [17]. He discussed the possibility of considering the variance linear in time, which implies linear diffusion and infinite speed of propagation, which Einstein noted as impossible from the physical point of view. It was also explicitly stated that depending on the context of applicability this property has its drawbacks. Modeling morphogen gradients with linear diffusion necessarily implies that every cell in the responding field receives a low level of morphogen instantaneously and, hence, predicts the same time of exposure for every cell in the whole tissue [18]. This is at odds with the observation that the development of morphogenetic responses requires time. Indeed, low level Hh signaling has cumulative effects in the neural tube and limb buds [18,19]. Linear diffusion suggests a response in cells at the opposite end to the source within the morphogenetic field, as these cells receive low signal levels from the onset due to instantaneous spreading. However, this is not commonly apparent. For example, in the chick limb bud the Shh gradient appears to be shorter than expected, with the proliferation of affected cell groups that leave the zone of influence giving the late appearance of a large gradient [20]. The cells that escape retain their early determination driven by a given dose of Shh signaling but cease to receive/respond to the signal when they leave the zone of influence. Under linear (normal) diffusion, they would be predicted to continue to receive Shh due to its instantaneous spreading. Thus, although linear diffusion can suggest a plausible quantitative approach, modeling Hh morphogenetic output by linear diffusion [7] poses so far unsolvable problems from a qualitative point of view. In addition, it has even been suggested that linear diffusion cannot correctly model morphogen dispersion [21]. The vertebrate neural tube and the insect wing imaginal disk are prime examples of morphogenetic patterning. In the former, embryonic pseudostratified ventral epithelial cells are instructed to acquire specific fates in response to Sonic Hedgehog (Shh) morphogenetic signals derived from the ventral floor plate, and earlier the notochord [22]. In the latter, Hh signal from the posterior compartment establishes pattern in the anterior compartment of the developing wing epithelium. Cells closest to the source receive higher doses of a morphogen and for longer times than those further apart from it, but, importantly, distal-most cells do not respond. Shh/Hh binds the membrane receptor Patched 1 (Ptc1) and inhibits its repression of the transmembrane protein Smoothened (Smo). In the presence of Shh/Hh signaling, Ptc1 is inhibited, Smo is active and this modifies the Gli code, the combinatorial function of activators and repressors of the Gli/Ci zinc finger transcription factors that lead to the regulation of target genes [23]. Shh/Hh morphogenetic action and interpretation thus require time and have been shown to involve feedforward and feedback responses, time-dependent gradient build-up, including a temporal accommodation to the morphogen and differential sensitivity to varying concentrations [18,24–28]. Experimentally, a number of observations appear incompatible with the use of linear diffusion to model Hh morphogen action: First, Hh morphogens do not behave as very small particles in large spaces, thus invalidating one of the assumptions of Brownian motion: Hh has been detected in Drosophila as visible aggregates of ∼20–300 nm in diameter [29,30]. Such particles have been shown to be multimers [31], membrane vesicles [32], oligomers [30], and/or lipoprotein particles [33,34]. Second, the dynamics of morphogen-induced pattern formation are inconsistent with the action of linear diffusion. For example, Islet1, a marker of motor neurons that mature at a distance from the ventral

M. Verbeni et al. / Physics of Life Reviews 10 (2013) 457–475

459

midline of the neural tube, is observed at early stages in ventral midline cells intermixed and even co-expressed with FoxA1 (HNF3beta), a ventral-most floor plate marker [35]. Indeed, the time of exposure, and not only the steady state morphogen concentration, is critical to specify cell fate [18]. Similarly, the induction of mesodermal markers in frog animal caps by a source of Activin occurs in a dynamic fashion: markers of distal positions appear only proximally at early stages, which are then replaced by proximal markers [36]. Thus, while current linear diffusion models can strikingly reproduce natural final patterns [37] they cannot account for how these are formed. An exception appears to be the Nodal–Lefty system, where differential dispersivity of morphogen and inhibitor can support modeling by linear diffusion [11]. This shortcoming of linear diffusion models cannot be simply resolved by imposing an arbitrary minimum threshold concentration of non-instructive signaling or noise [7]. Similarly, resolving the inherent problems of linear diffusion by appealing to the smallness of the diffusion tails in graphs plotting morphogen concentration versus distance from the source is incorrect since the accuracy of the computers used to develop the numerical simulations imposes its own threshold, which is independent of any biological reality. A further shortcoming of linear diffusion is that it cannot yield moving concentrations fronts, which are essential for biologically-relevant morphogenetic action over time, precisely because low Shh levels are read cumulatively over time [18,19]. In this context it is important to note that the existence of a front and the velocity of morphogen transport are critical elements. However, the sharpness of the front is not, as it depends on the quantity of morphogen being sent and the kinetics of morphogen take up by the responding cell. For example, the front of a second morphogen wave over a field of cells that has already responded to a previous wave may be less sharp than the front of the first wave produced at the onset of morphogen secretion, but both are governed by the same velocity of morphogen transport. If simple linear diffusion models cannot describe Hh morphogenetic signaling, how can it be modeled? Two recent findings suggest a radically different perspective for modeling Hh-Gli signaling. The first finding is related to advances in flux-limited nonlinear diffusion (in connection with optimal mass transportation theory) [42–51], which allows the formulation of a model with non-artificial time of exposure through the incorporation of a finite speed of signal propagation that depends on the properties of each system. In the case of Hh, this velocity can be measured experimentally but is context-dependent. These considerations contrast sharply with the use of Fick’s law since the propagation speed inherently included in linear diffusion models is always infinite, independently of the biological setting. Flux limited diffusion phenomena appear not only in the class of biological systems, but also in the case of interacting living entities when these have the ability of perceiving gradients. For instance this is the case of the dynamics of traffic, swarms and crowds (see for example [52,53]). Second, Hh/Shh particles do not move freely or randomly in the extracellular matrix but rather they are found associated with and move on cell extensions [21,38–41]. In this context, direct delivery modality of morphogens to target cells involves long filopodia or cytonemes (also known as cytoplasmic threads [38]) that extend between morphogen producing and receiving cells, thus functioning as conduits for morphogen traffic: Cytonemes can act as tracks to convey vesicles containing signaling molecules. This modality involves cell-to-cell contact and it might be the most precise spatial and temporal mechanism to transfer information between cells. Importantly, the dynamic instability of cytonemes and the existence of cytonemes of different lengths can generate an overall graded morphogen distribution. Here we sought to modify linear diffusion models of Hh signaling by incorporating FLS in order to realistically develop the correct Gli-expression patterns in responding cells at different distances from the source of Hh morphogen. In addition, we propose that cytonemes and associated Hh macroparticles are at least one of the biological bases of FLS. Our proposed model is an extension to Fick’s law and of the pioneering ideas of Turing, but overcomes the qualitative weaknesses associated with linear diffusion. Our present analysis to develop a general model for Shh-Gli signaling in the vertebrate neural tube is based on a continuous feedback between mathematical modeling, numerical analysis and data collection from experiments on Hh particle movement in the insect wing imaginal disk epithelium and in mammalian cells. We further use our model with published results on Shh gradient formation and action in a generic vertebrate neural tube in order to directly compare and contrast our work with previous results using linear diffusion. 2. Results The problems highlighted above with linear diffusion models raised the possibility that restricted spreading could better describe and underlie morphogen action. We thus opted to analyze the limitation of particle spreading by

460

M. Verbeni et al. / Physics of Life Reviews 10 (2013) 457–475

Fig. 1. Fluorescence recovery after photobleaching (FRAP) experiments in the 3rd instar Drosophila wing imaginal disk using a Hh-GFP transgenic line, at time t = 0 (A) and t = 275 s (B). (C) Plot of the values of the average pixel intensity at various times after photobleaching as a function of the distance from the A/P border in the imaginal wing disk. The dotted line and the arrow in (B) indicate the apparent position of the signaling front. Normalized fluorescence intensities are shown in (C), determined on very narrow strips with widths ranging between 0.68 and 1.68 µm perpendicular to the front propagation, while in (D) the bins to calculate mean velocities where of 3 µm. The gray arrows in (C) point to the recovery fronts. (D) Data derives from 4 FRAPs experiments and is shown as a function of the distance from the floor plate. The finite velocity of propagation of the fluorescent Hh front, ranges from 0.07 to 0.02 µm/s.

regulated flux (FLS) as this could provide a conceptually novel physical and mathematical approach to morphogen modeling. FLS suggests constraints to limit morphogen particle spreading, perhaps by the retention of secreted morphogens [34,54,55], and, therefore, that the formation of a gradient requires time. Moreover, it also suggests that at the front of the morphogenetic traffic wave, there will be a naturally occurring discontinuity that could perfectly fit the observed cellular behavior [56,57], without the necessity to impose arbitrary noise thresholds (see Fig. 8 in [7]). 2.1. Measuring Hedgehog gradient dynamics in the Drosophila imaginal wing disk Our approach may apply, in principle, to various Hh morphogenetic gradients, including that in the fly wing imaginal disk. This idea is supported by the general conservation of the Hh signaling pathway and the finding that mammalian Shh can functionally replace insect Hh [58]. A new key parameter for FLS modeling (see below) is the average value of the velocity of morphogen macroparticles, that we hereafter denote by c, which determines the concentration front. To measure the speed of morphogen propagation we used fluorescence recovery after photobleaching (FRAP) in the Drosophila wing imaginal disk, where functional Hh-GFP is produced in the posterior (P) compartment and spreads anteriorly forming a concentration gradient (see [59]; (Fig. 1A–B)). The receiving cells located in a stripe of anterior tissue abutting the antero-posterior

M. Verbeni et al. / Physics of Life Reviews 10 (2013) 457–475

461

(A/P) compartment border interpret the morphogen signal so that the expression of different Hh target genes occurs according to the concentration of Hh experimented by the responding cells. Given the difficulty to image the first wave of Hh particles in the earliest imaginal disk primordium, we have used explanted 3rd instar imaginal disks at times to photobleach an experimental square area (eliminating GFP fluorescence but not the presence of Hh particles) in the anterior compartment adjacent to the antero-posterior boundary (Fig. 1A). Hh-GFP is recovered over a very short time (2 minutes) after photobleaching, suggesting that for our model we can ignore tissue growth. We also note that the time to set up the Hh gradient, estimated by the photobleaching experiments, is only of a few minutes. During recovery of fluorescence after photobleaching (Fig. 1B), we measured the values of the average GFP pixel intensities in four independent experiments. Multiple measurements per experiments were obtained at increasing distances from the antero-posterior border, and the results plotted as a function of the distance from this border (Fig. 1C). The data indicate that the concentration of Hh, as measured by fluorescent intensity, increases over time with Hh fronts propagating with finite speed: Cells located at increasing distances from the antero-posterior border gain florescence at increasingly longer times, as compared with those closer to the border. A mean value of 0.07 µm/s was found for c for cells closest to the A/P border in the first few seconds of measured fluorescent recovery, followed by an average value of 0.02 µm/s throughout the photobleached area (Fig. 1D). Our analyses of limited data previously reported for fluorescent intensities of Shh/Hh-GFP in other systems [56,57] yield approximate average values of 0.0013 µm/s, suggesting the general consistency but variable nature of c. 2.2. Flux-limited modeling of the Shh morphogen gradient Although the value of c has been determined in the fly disk (see above), we opted to build our model using the vertebrate neural tube paradigm since Hh behavior appears to be conserved, and there is a previous linear diffusion model for Shh signaling using this system [7]. Given the short time over which we measure Hh propagation in the imaginal disk, or model Shh propagation within the neural tube, we consider the effect of tissue growth negligible for our study, much as in a previous model [7]. After the pioneering work of Rosenau [42,43], recent advances in deducing flux-limited processes from a microscopic point of view derive from optimal mass transportation theory [45], see also Appendix C, and from the hyperbolic limit of kinetic multiscale biological tissue models for multicellular growing systems [47] (see also [49, 48] and references therein), where, in contrast to linear diffusion, secreted molecules are not considered to move randomly. In order to provide insights on the genesis of flux-limiters, we viewed linear diffusion from the perspective of kinetic transport equations, which leads us to consider the interactions among particles from a microscopic point of view using the Fick equation (linear or normal diffusion) in one spatial dimension,   ∂u ∂ 2 u ∂ ∂ ∂ = 2= u ln u = [u v], ∂t ∂x ∂x ∂x ∂x where v is the microscopic field velocity associated to the molecule expressed as ∂(S(u)/u)/∂x, and S(u) = u ln u is a function related to the relative entropy of the system. We considereda new microscopic velocity as a consequence of optimizing the motion of the particle with respect to the line element 1 + |v|2 dx leading to the flux-limited equation   ∂u |u| ∂x ∂  ∂t u = ν . 2 ∂u ∂x |u|2 + νc2 | ∂x |2 c represents the maximum allowed macroscopic speed of propagation of the front and ν the viscosity. It has been proved [49] that in the context of this system the signal propagation takes place precisely at speed c. When the velocity c goes to infinity we recover the classical Fick’s law   ∂u |u| ∂x ∂u ∂ 2u ∂  =ν −→ ν 2 . 2 ∂u c→∞ ∂x ∂t ∂x |u|2 + νc2 | ∂x |2 The study of the existence and qualitative properties of traveling waves associated with a nonlinear flux limited partial differential equation coupled to classical reaction terms has allowed the deduction of the existence and uniqueness of finite speed moving fronts of classical regularity, but also the existence of discontinuous entropy traveling wave

462

M. Verbeni et al. / Physics of Life Reviews 10 (2013) 457–475

Fig. 2. Schematic diagram of the Shh-Gli pathway with key interactions as detailed in the text. The temporal time line follows the flow of positive arrows and negative T bars from the action of the secreted Shh morphogen on the left. Note that the Gli code includes the function of Gli2 and Gli3 activators as well as Gli2 repressors. For simplicity we incorporate all repressor function in the term Gli3Rep . For activators, we include Gli3Act but consider Gli1Act as the major Gli activator.

solutions [51]. The biological implications of such mathematical structures allow each cell to receive the signal according to the position of the propagation front in the tissue and to give different responses in the absence or presence of morphogens. Importantly, here the movement of morphogens is directed rather than random as in linear diffusion. We modeled the propagation of Shh in a generic neural tube by introducing a flux-limited reaction–diffusion equation. This embryonic structure is polarized along the dorsal-ventral axis [0, L], representing a natural privileged direction for the description of morphogen transport, thus justifying the use of a one-dimensional model:  ∂[Shh] = ∂x a [Shh], ∂x [Shh] + koff [Ptc1Shhmem ] ∂t − kon [Shh][Ptc1mem ], where

 a [Shh], ∂x [Shh] = ν 

[Shh]∂x [Shh] [Shh]2 +

ν2 (∂ [Shh])2 c2 x

is the flux-limiter, the ligand–receptor complex Shh·Ptc1 concentration is denoted by [Ptc1Shhmem ] and [Ptc1mem ] indicate the concentration of Ptc1 at the cellular membrane. The second and third term in the above equation derive from applying the law of mass action to the reversible reaction in which Shh binds to its receptor Ptc1. Boundary conditions prescribe the value of the flux at the floor plate (x = 0) of the neural tube a([Shh], ∂x [Shh])(0) = β, and at the roof plate x = L, where we impose a weak Dirichlet boundary condition a([Shh], ∂x [Shh])(L) = −c[Shh](L) [48], being L  300 µm for a generic vertebrate neural tube. The use of a flux-limited reaction–diffusion equation involves the appearance of a morphogen front at the distal end of the gradient in the neural tube. This situation contrasts with the predictions of linear diffusion models [7], in which a([Shh], ∂x [Shh]) = ∂x [Shh] and where the signal results to reach all neural tube cells immediately after secretion. 2.3. Modeling morphogenetic responses Throughout this section the equations governing the signal transduction inside a responding cell make sense only after Shh particles have reached the cell in question. Given that the support of the function representing the Shh concentration grows like ct, since it is bounded by the solutions of the flux-limited equation [49], the equations for the proteins inside the cells start to become operative only in the space interval [0, ct]. Then, in the absence of Hh/Shh ligands, the expression of positive Gli targets is repressed: there is no production of Gli1Act or enhancement of Ptc1cyt . Moreover, the transcription of gli3 is not inhibited. Hence, we consider the concentrations of Gli3Act , Gli3Rep and Ptc1cyt as constants. To begin to model reactions triggered inside the cell by the activation of a simplified pathway (Fig. 2) with the law of mass action, we describe the concentration of the ligand/receptor as ∂[Ptc1Shhmem ] = kon [Shh][Ptc1mem ] − koff [Ptc1Shhmem ] + kCout [Ptc1Shhcyt ] − kCin [Ptc1Shhmem ], ∂t where [Ptc1Shhcyt ] denotes the concentration of the complex Shh·Ptc1 internalized via receptor mediated endocytosis [59], which we describe as ∂[Ptc1Shhcyt ] = kCin [Ptc1Shhmem ] − kCout [Ptc1Shhcyt ] − kCdeg [Ptc1Shhcyt ]. ∂t

M. Verbeni et al. / Physics of Life Reviews 10 (2013) 457–475

463

The evolution equation for the concentration of Ptc1, [Ptc1mem ], is ∂[Ptc1mem ] = koff [Ptc1Shhmem ] − kon [Shh][Ptc1mem ] + kcyt [Ptc1cyt ], ∂t where [Ptc1cyt ], represents the concentration of Ptc1 produced by the target gene ptc1, once the Shh signaling pathway becomes activated [60]. Receptors are degraded in the internalized ligand–receptor complex Ptc1Shhcyt [59,61]. After synthesis, [Ptc1cyt ] molecules are exported to the surface with a rate kcyt . To distinguish numerically between the two regions defined by the front (before and after arrival), we have considered the fractional occupancy: Y=

[Ptc1Shhmem ] . [Ptc1mem ] + [Ptc1Shhmem ]

In the region where Y > 0, Shh molecules have bound transmembrane receptor, thus triggering the activation of the signaling pathway inside responding cells. In the absence of Shh signaling Y = 0. The synthesis of the receptor Ptc1 upon activation of Shh signaling is described by 





 ∂[Ptc1cyt ] = −kcyt [Ptc1cyt ] + kP Ptr Gli1Act (t − τ ), Gli3Act (t), Gli3Rep (t) ΦPtc , ∂t where Ptr is a weight term evaluated with the statistical technique BEWARE (see SM1, [7] and references therein), used to give quantitative treatment of transcriptional events occurring through the interactions of regulatory proteins (in this case Gli1Act , Gli3Act and Gli3Rep ) with specific DNA sequences. The function Ptr contains the information of the fact that gli1 positively regulates ptc1 (see SM1) and gli1 itself. kP is the rate of Ptc1 synthesis. τ is the delay time elapsed between reception of the Shh morphogen and the synthesis of Gli1Act /Ptc1 inside the cell, which can amount to hours. Finally, the factor ΦPtc is defined by ΦPtc =

[Ptc10 ] , [Ptc10 ] + [Ptc1mem ]

where [Ptc10 ] is the initial value of [Ptc1mem ]. Multiplying by this factor, we take into account the effect of the negative feedback by the increase of Ptc1 expression due to pathway activity. In addition to c, we introduce a second important parameter in FLS modeling to account for the time delay between reception of the morphogen and the genomic response. For Shh/Hh signaling we equate it to the delay in the synthesis of Gli1Act and Ptc1. Cellular internalization is hereby limited to the participation of Ptc1, which is also a negative regulator of Smo and thus of Gli activators acting downstream [23]. Since Gli proteins have different activator levels (generally Gli1Act being the strongest) and repressor activities (generally Gli3Rep being the strongest) in different species and contexts, we consider pure activator Gli1Act and repressor Gli3Rep functions, as in Drosophila CiAct and CiRep [62], as the opposite poles in a graded series of states, the combination of which comprises all other intermediate values, including Gli2Act and Gli2Rep . To further simplify the model, we also assume that the concentration of Gli3Act remains low and constant (for instance driven by Wnt activity [63]) without the influence of Shh signaling. In terms of Gli regulation, Shh activity induces transcription of gli1, the only reliable biomarker of response to Shh signaling [64], and the levels of Gli1Act (see Fig. 3E) can be considered as the final Shh pathway output whose evolution is

 





 ∂[Gli1Act ] = −kdeg Gli1Act + kG Ptr Gli1Act (t − τ ), Gli3Act (t), Gli3Rep (t) ΦPtc , ∂t where kG is the rate of Gli1Act synthesis, and the first right-hand side term gives Gli1Act degradation. Evolution equations for the concentrations of the activator form of Gli3, Gli3Act , and its repressor form, Gli3Rep (see Fig. 3G), are given by  kg3r



γg3 ∂[Gli3Act ] − Gli3Act − kdeg Gli3Act = ∂t 1 + RPtc 1 + RPtc 

 k ∂[Gli3Rep ]

g3r − kdeg Gli3Rep , = Gli3Act ∂t 1 + RPtc

464

M. Verbeni et al. / Physics of Life Reviews 10 (2013) 457–475

Fig. 3. Plots of Shh, PtcShhmem , Gli1Act , and Gli3Rep concentrations versus distance from the floor plate at various times. The plots have been obtained numerically solving our FLS (A, C, E, G) and the Saha and Schaffer [7] linear diffusion model (B, D, F, G). Note that in linear diffusion modeling (B) there are no natural fronts and that an artificial threshold had to be imposed at 2.5 nM in order to achieve them (see Fig. 8 in [7]). However, note that even with such arbitrary noise threshold the fronts retrocede in a non-biologically-relevant manner. In contrast, correctly advancing or retreating fronts appear naturally in FLS (A, C, E, G).

M. Verbeni et al. / Physics of Life Reviews 10 (2013) 457–475

465

Fig. 4. (A, B) Evolution of Gli1Act over time at a distance of 1 µm from the floor plate resulting from signal accumulation and subsequent desensitization in our FLS model (A). This behavior is not observed in the linear diffusion model of Saha and Schaffer [7] (B). The FLS curve reproduces a temporal adaptation mechanism: after 12 hours cells become desensitized to Shh signal and the response decreases. (C) Numerical simulation of FLS with the value of c obtained in our Drosophila FRAP experiments, where c = 0.03 µm/s. Note the sharpness of the gradients and the advancing concentration fronts.

where RPtc =

[Ptc1Shhmem ] [Ptc1mem ]

represents the ratio of liganded to unliganded Ptc1 at the cellular membrane. Gli3 and Shh/Gli1 function repress each other [65]. Here γg3 represents a source term for Gli3 transcription factor, while kg3r is the rate of Gli3 proteolytic cleavage into the smaller repressor form Gli3Rep [38,66,67]. To specify repression by Shh, we divided these two terms by the factor 1 + RPtc . The use of this quantity is suggested by the finding that a cell’s response to Hh signaling is determined by the ratio of liganded to unliganded Ptc1mem [68]. Finally, the last terms in the above equations describe degradation. 2.4. Numerical experiments The values of the constants in the FLS coupled equations model (koff , kon , . . . , kg3r , kdeg ) are taken from previous studies [5,7,69], with emphasis on the values of the linear diffusion model of Saha and Schaffer (6), or derived from estimations (see SM2). To start the analysis of our FLS neural tube model we have first used the average value of c = 0.0013 µm/s. Consistent with the existence of morphogen fronts, the qualitative behavior of our FLS involves a singularity front from the origin that is propagated away from it over time (Figs. 3, 6). Given that the aim of the numerical experiments is to test the effects of such fronts on morphogenetic responses, our FLS system must therefore be treated by specific numerical methods such as the WENO interpolation scheme [70], which takes into account the fact that the characteristics of the proposed model are more hyperbolic than parabolic (see Material and Methods). Numerical simulation of FLS yielded curves for morphogen ligand, receptor, and transcriptional responses that, unlike linear diffusion, fit the biological data (Figs. 3, 6, 7). Critically, numerical simulation yielded dynamic morphogen distribution curves over distance from the source that have moving fronts, with Gli1Act being enhanced shortly after Shh morphogen reception and Gli3Rep being repressed as the Shh wave moves forward. One key consequence of these two outcomes is the tight regulation of the Gli code: the increase in the level of Gli1Act is progressive (Fig. 3E, 4A) and starts straight after the Shh/Hh signal front reaches responding cells. Moreover, the opposite evolution of Gli3Rep expression (Fig. 3G) inversely fits with the evolution of Shh/Hh signaling (Fig. 3A). These behaviors are absent in linear diffusion models, see (Figs. 3, 6) and [7]. Numerical simulations of FLS show that cells loose sensitivity to ongoing Shh signaling (Figs. 4A and 3E), an unpredicted result that is in agreement with experimental data (see [18,71]). Such desensitization process may be related to the negative feedback by Ptc1 – signaling output decreases in cells exposed to concentrations of Shh not sufficient to bind all Ptc1 receptors – and is considered by the negative action of Ptc1mem through ΦPtc in the equations for the concentrations of Ptc1cyt and Gli1Act . FLS, unlike linear diffusion, thus accounts for the ability of the morphogen to

466

M. Verbeni et al. / Physics of Life Reviews 10 (2013) 457–475

upregulate factors (e.g. Ptc1) that feedback to control the interpretation of the gradient (Figs. 3E, 4A, and 7). Such temporal adaptation mechanism predicts that morphogenetic action takes place through a dynamic interplay between the graded signal and the target cells [18]. To test if the value of c could alter gradient dynamics, further simulations with the value of 0.03 µm/s obtained in our Drosophila FRAP experiments (Fig. 1) resulted in gradient fronts with sharp slopes over short distances and times (Fig. 4C). Interestingly, this result nicely fits the biological reality of a short-distance morphogenetic field directly under the control of Hh secretion in the insect wing. These findings suggest that the value of c is a biologically-relevant variable that can yield gradients of different shapes and sharpness. 2.5. Are cell extensions or cytonemes flux limiters? Given the complexity of the system, our current FLS model is one-dimensional and does not take into account specialized cell structures. However, FLS predicts the existence of restricted modes of information transfer between cells, which sharply contrasts with the Brownian propagation which is the basis of linear diffusion. Extracellular binding partners have been suggested to alter linear diffusion of morphogens [11,55] and could underlie FLS. In addition, however, cell extensions or filopodia-like structures called cytonemes [38,39,21] or nanotubes [72,74] could also underlie spatial restriction of morphogenetic action. Cytonemes have been proposed to behave as intercellular signaling routes for multiple morphogens in Drosophila [39,21], including Hh signaling (see [73,75]), and nanotubes can transport cargo from one cell to another [72,74]. In the same context, very recently it has been shown that specialized filopodia direct long-range transport of Shh during vertebrate tissue patterning [41]. To test if we could detect such cell extensions and measure the velocities of putative macroparticle cargo, we have used fluorescent GFP and tdTomato (TOM) molecules tagged to the cell membrane in a non-biased manner, as well as Hh-GFP, Ihog-YFP and CDO-GFP protein fusions. Ihog in Drosophila and CDO in mammals are related Shh/Hh binding Ig-like and FNNIII domain-containing membrane proteins (reviewed in [76]) that form part of Shh/Hh macroparticles. Expression of Ihog-YFP, and Hh-GFP separately or together driven by a HhGal4 driver in the posterior compartment of the Drosophila imaginal wing disk resulted in the visualization of yellow (for Ihog-YFP) or double yellow/green YFP+ /GFP+ (for their co-expression) filopodia or cytonemes growing into the anterior compartment (Fig. 5A and not shown). In addition, spectral microscope imaging revealed their decoration by multiple YFP+ /GFP+ Ihog/Hh-containing macroparticles (Fig. 5B). Sequential imaging of the Hh-containing macroparticles on filopodia or cytonemes showed an average velocity of movement towards the A compartment of 0.03 µm/s with instantaneous velocities ranging from 0.025 to 0.07 µm/s (Fig. 1D). To extend these results to a divergent system, we tested for the presence of cell extensions and associated macroparticles in human glioblastoma multiforme (GBM) cells (Fig. 5C–K), which are known to express Shh and require Shh-GLI signaling for tumorigenicity and growth [77,78]. Culture of GFPmem , Ihog-YFP- or CDO-GFP-lentivector transduced human U87 GBM cells revealed the presence of long filopodia, likely nanotubes or cytonemes, often between cells (Fig. 5C–G). These human cell extensions were highly dynamic, appeared to show anchorage points at the tips and bifurcations, remodeled quickly and displayed associated particles. Similar structures were also seen in primary GBM cells and in U87 cells grown on a monolayer of growth-arrested mouse embryonic fibroblasts (data not shown), ruling out cell line or growth on plastic artefacts. Live recording of GBM Ihog-YFP and CDP-GFP macroparticle movement using 2-photon videomicroscopy revealed that while there was no difference between CDO-GFP+ or Ihog-YFP+ markings, and that individual particles, even on the same cell extension (Fig. 5G), displayed different velocities. The average velocity was 0.33 µm/s, with instantaneous speeds ranging from 0.005 to 0.9 µm/s (Fig. 5H, I). These velocities are similar to those for rapid axonal transport [79], and are compatible with those of the transport of Hh or Ihog macroparticles in Drosophila (Fig. 5B), which range between 0.036 µm/s and 1.333 µm/s. To ascertain that cell extensions could sustain information transfer, we co-cultured U87 GBM cells transduced independently with membrane-bound RFP mem- or GFP mem-expressing replication-incompetent lentivectors. In this assay, red cells mix with and touch green cells (Fig. 5J), thus allowing us to test if particles of one color can move on cell extensions of the other color. Indeed, we could detect green particles moving along red extensions (Fig. 5K), with varying speeds and an apparent external association to the substrate as the particles were seen bobbing to either sides of the cell extension while still attached and advancing. This behavior strongly suggests that human cells can exchange generic membrane cargo via cell extensions.

M. Verbeni et al. / Physics of Life Reviews 10 (2013) 457–475

467

Fig. 5. Detection of cell extensions and trafficking dynamics of membrane macroparticles in the Drosophila imaginal wing (A, B) and in human glioblastoma (GBM) cells (C–K). (A–B) Images of the Drosophila wing imaginal disk showing cell extensions, possibly cytonemes, labeled with Ihog-YFP-GFP, extending from the posterior compartment into the anterior compartment. Puncta (macroparticles) can be observed on the extensions (arrows in B). (C) Pseudo-colored fluorescence image of cell extensions observed between U87 cells labeled with membrane targeted GFPmem. The enlarged box (magenta stripped line) shows a high magnification of the corresponding detail, showing cell extensions between two cells and associated macroparticles (arrows). (D–E) Pseudo-colored images of thin filopodia – possibly nanotubes or cytonemes – observed between U87 cells labeled with CDO-YFP. (D, E) show that same image except that the look up table in (E) has been saturated as compared with that in (D) in order to see more clearly the cell extensions. (F) Pseudo-colored image of cell extensions observed between U87 cells labeled with Ihog-YFP. The magenta box shows that detail zone analyzed in (G). (G) Analysis of Ihog-YFP + particle trafficking dynamics on cell extensions. Two particles (colored arrows) move independently and at fluctuating velocities on the same extension. Time in seconds (s) between two frames is indicated on each image. (H) Distribution of the frame-to-frame velocity computed for multiple particles under all cultures and staining methods as in (C–G). (I) Average and individual particle velocity in cell cultures labeled with different constructs as indicated. (J) Low magnification image of co-cultured U87 cells labeled either with membrane targeted GFPmem or membrane targeted RFPmem. (K) Red (RFP+ ) vesicle (arrows) moving onto a green cell extensions over time. Time in seconds (s) between frames is indicated on each image. The dotted white box shows a high magnification image of the red particle seen on the cell extension. In this case it ‘hangs’ over the right, while in the frames at 0, 132 and 264 s it hangs over the left side.

Fig. 6. Evolution of Ptcmem concentration versus distance from the floor plate at various times. The plots have been obtained numerically solving the Saha and Schaffer [7] system of coupled equations with linear diffusion (B) and our FLS model (A).

468

M. Verbeni et al. / Physics of Life Reviews 10 (2013) 457–475

Fig. 7. Evolution of Ptccyt concentration versus distance from the floor plate at various times in our model.

3. Discussion Here we use novel flux-limiters in mathematical modeling and determine the average velocities of Hh macroparticles to formulate their morphogenetic action without invoking linear diffusion. FLS resolves so far unaddressed fundamental problems with linear diffusion and, importantly, it predicts restricted movement of morphogen particles. The most important outcomes of our FLS model, which are not possible with linear diffusion, are: (i) The elimination of the unphysical infinite propagation speed of morphogen particles leading to the presence of biologically-critical concentration fronts (Fig. 4C and 3A). (ii) The assurance that different cell groups will be exposed to different morphogen concentrations at different times (Fig. 4C and 3A, E). (iii) The desensitization of the morphogen response in relation to time, as seen with the level of Gli1Act (Fig. 4), which was unpredicted in previous linear diffusion models of Shh action [7,18], but suggested – for different reasons – in a previous model by Meinhardt [80]. (iv) The correct regulation of the GLI code [81], with inverse regulation of Gli activators and Gli repressors as shown for Gli1Act and Gli3Rep , over time and space. (v) The notion of spatially restricted information transfer and the possibility that filopodia-like structures such as cytonemes or nanotubes could represent the flux-limiters predicted by FLS. We suggest that FLS modeling provides a qualitative realistic description of natural events, with a critical and distinct discontinuity in the wave front that establishes the Hh/Shh gradients. Modifications of the present model will be required to further refine it and to adapt it to different contexts. Such improvements will involve, for instance, developing a two-dimensional model taking into account later growth of the morphogenetic field, the modification of Gli isoforms, the effects of Sufu, Cdo/Ihog, Boc/Boi, the feedback regulatory loop between Gli1, Nanog and p53, or the integration of peptide growth factor-RTK-Ras/Akt signaling [23,78,82]. Also, the unidirectional ‘promotion’ observed in [83] refers to an optimal transportation problem that might best explained in the future by FLS rather than by linear diffusion since it is intrinsic to our model. FLS also predicts that Hh/Shh morphogenetic information is not freely available, oozing from its source, but rather that it is delivered in a restricted, possibly oriented manner. The recent description of nanotubes [72,74] and cytonemes [38,39,21,41], their implication in various kinds of signaling pathways [39,21,41] and our imaging of foreign cargo on a cell’s extensions raise the possibility that these structures could act as information channels or bridges, thus contributing to flux limitation. It is possible that extracellular or membrane binding partners mediate the localization of macroparticles to cell extensions. In this case, both mechanisms could contribute to flux limitation at different levels. Our detection of Shh/Hh complex macroparticles in such structures, rather that freely in the extracellular space or simply near to the cell’s membrane, provides evidence in this direction. These observations and the confirmation of the existence of limit velocities c detected in our Drosophila FRAP experiments (Fig. 1C–D), independent measurements

M. Verbeni et al. / Physics of Life Reviews 10 (2013) 457–475

469

[56,57], and on insect imaginal disk and human GBM cell extensions (Fig. 5I), raise the possibility that such extensions could have a more general role in morphogenetic information transfer between cells. The regulation of the velocity of macroparticle movement, c, modifies the speed of the concentration front, suggesting that a way to generate gradients of different sizes and shapes involves the modulation of macroparticle speeds. Modulation of the new parameter c that we highlight here may thus be central in the evolution of animal form and shape. Moreover, we suggest that FLS may be relevant not only to Shh/Hh signaling during development but also to other scenarios such as the control of morphogen action in stem cells and cancer. 4. Materials and methods Fly stocks. hh-Gal4 driver used in the ectopic expression experiments using the Gal4/UAS system. Transient expression of the UAS-constructs using the Gal4; tub-Gal80ts system was done by maintaining the fly crosses at 18 ◦ C and inactivating the Gal80ts for 24 to 30 hours at the restrictive temperature (29 ◦ C) before disks dissection. The pUAS-transgenes were: UAS-Hh-GFP (see[59]), UAS-Ihog-YFP [82]. Time lapse confocal image and FRAP analyses in the wing imaginal disks in tissue culture cells. A laser scanning confocal microscopes (Zeiss laser confocal and multifotón LSM710 coupled to an inverted AxioObserver (Zeiss) and Argon 458/488/514 nm laser and DPSS 561 nm laser) were used for confocal fluorescence imaging. Image J software (NIH USA) was used for image processing and determining fluorescence levels. All image time series were acquired with a 40x/1,30 EC “Plan-Neofluar” Oil DIC M27 (Zeiss) and using a Zeiss Zen 2008 program. Lentivectors. The fusion protein Ihog-YFP (described in [82]) was amplified using the forward primer 5 aata tctaga atg acg ctg ctc ac 3 and the reverse primer 5 ctga gtcgac tta cac gcc aac gct gtt 3 . Using the restriction sites for XbaI and SalI, the fragment was introduced into the lentivector Tween. The CDO sequence was amplified with the primers 5 act gac cgg tat gca tcc aga cct cgg 3 and 5 cag tac cgg tgt ctc ttg ggc ttg c 3 by PCR and introduced into the LL3.7, to build a CDO-GFP fusion protein. Numerical solutions of the FLS model. In order to numerically solve the eight coupled equations composing the FLS model, we wrote a program in FORTRAN77. The main problem we had to address was the interpolation of the flux limiter, introduced to reproduce propagation fronts, moving with finite velocity. The flux-limited equation, without reaction terms and with Neumann and weak Dirichlet boundary conditions is deeply analyzed in [49]. The discontinuity fronts appearing in the solutions are typically associated with nonlinear hyperbolic phenomena. Due to this behavior, we employed the conservative fifth-order finite difference WENO (Weighted Essentially Non-Oscillatory) scheme [70], to obtain piecewise polynomial reconstruction of the flux. To advance the solution in time, we used the optimal five-stage fourth-order strong stability preserving Runge–Kutta method of Spiteri and Ruuth [84]. To calculate numerical solutions we employed a grid with 500 spatial points and 5 × 107 time points. Cell culture and infection. The U 87 cell line was purchased from American Type Culture Collection (ATCC) and grown as specified. Cells were infected with 106 –107 viral stocks in the presence of 5 µg/ml polybrene. Cells were cultured on plastic dishes or on feeder layers of mouse embryonic fibroblasts inactivated with 10 µg/ml mitomycin C (Acros Organics) for 3 hours. 2-Photon laser-scanning microscopy and image analysis in tissue culture cells. Individual cell culture coverslips were transferred to a submerged chamber (Luigs and Neumann, Germany) and perfused with oxygenated artificial cerebro-spinal fluid (ACSF) warmed at 37 ◦ C (±0.5 ◦ C) at a rate of 1–2 ml/min. The ACSF contained (in mM): 124 NaCl, 3 KCl, 2 CaCl2, 1.3 MgSO4, 25 NaHCO3, 1.2 NaH2PO4, 10 D-glucose with pH = 7.34 and osmolarity of 300 mOsm when bubbled with 95% O2–5% CO2. Two-photon laser-scanning microscopy was performed using a trimscope scanning head (LaVision BioTec GmbH, Bielefeld, Germany) mounted on an upright BX51W I microscope (Olympus, Switzerland). The specimen was illuminated with 915 nm light from a tunable pulsed Ti:sapphire femtosecond Chameleon Ultra 2 laser (Coherent, Dieburg, Germany). All image time series were acquired with a 25X water immersion objective (Olympus, Switzerland). Emitted light was collected using photomultiplier tubes (Hamamatsu, Japan). Scanning and image acquisition were controlled using Imspector Pro software (LaVision BioTec

470

M. Verbeni et al. / Physics of Life Reviews 10 (2013) 457–475

GmbH, Bielefeld, Germany). For the movies tracking the vesicle movement, a single optical plane (200 × 200 µm, 1024 × 1024 pixels) was imaged at a frame rate ranging from 0.18 to 0.33 Hz. For image processing and quantification, image time series were first registered in 2D, using the stackreg package running in ImageJ (NIH, USA), and then loaded in Neuromantic (http://www.reading.ac.uk/neuromantic/) in order to label the movement trajectory of individual vesicle on cytonemes. The morphological files generated by Neuromantic were then quantified using the L-measure analysis software [85]. Acknowledgements We are grateful to Martí Ruiz-Altaba for comments on the manuscript and discussion and to Carmen Ibáñez from the Confocal Service of the Centro de Biología Molecular for their assistance in the in vivo Drosophila imagen. This work was partially supported by Projects MTM2008-05271 and MTM2011-23384 (to MV, OS, JS), BFU2008-03320/BMC, BFU2011-25987 and from Consolider Program CSD2007-00008 (to IG) from the Spanish MICINN; by Marie Curie FP6 (RTN 035528-2) and FP7 (ITN 238186) and by an institutional grant from the Fundación Areces to IG; by fellowships awarded by a Marie Curie FP6 (RTN 035528-2) contract to EM; by Junta de Andalucía Project P08-FQM-4267; by grants from the Swiss National Science Foundation (SNF professor grant number PP0033-119169 and Sinergia grant CRSII3-127456), the University of Geneva, the European Research Council (contract number ERC-2009-StG-243344-NEUROCHEMS), the Leenaards foundation and the European Molecular Biology Organization (young investigator program) to AC; by a long-term HFSP fellowship to ISC in ARA’s laboratory; and by grants from the Swiss National Science Foundation, Swiss Cancer League and funds from the Département d’Instruction Publique from the Canton et Republique de Genève to ARA. This project derives from sustained interactions over 5 years between J.S. and A.R.A. at the BIOMAT School at the University of Granada. Appendix A. Parameters used in numerical simulations In Table A.1, we report the values of the parameters used to numerically solve the system of partial derivative differential equations composing our FLS model. Most are taken from [7]. Notes. The values of the kinetic and binding parameters have been taken, when it was possible, from the literature, as indicated. Otherwise, we chose the values of the parameter in order to fulfill some basic requirements. The kinematic viscosity ν controls the concavity of the curves representing Shh concentration: an increase in the value of ν results in more concave curves. With the indicated value of ν, numerical solutions for [Shh] are concave until t ∼ 12 hours and then start to show a convex part. We adopted this value, considering the distribution profiles of Hh-GFP in [56] and of Shh-GFP in [57]. The entering flux, β, was chosen in order to obtain values for Shh concentration comparable to the ones obtained in [7]. The rates of Gli1Act and Ptc1cyt synthesis affect the quantity of receptors present on the cellular membrane. Actually, Gli1Act , due to its strong activator function, both exerts a positive autoregulation and induces ptc1. In turns, the concentration of Ptc1mem influences the spread of morphogens by sequestering Shh proteins. We aimed to have a gradient reaching ∼150 µm by 30 hours. For the same reason we choose the reported values of r and . The value of r is also based on the strong repressor activity of Gli3Rep [65]. As explained in the main text, in the region of the neural tube where the Shh pathway is active, the action of Gli1Act dominates, while in cells not yet reached by the signal, high levels of Gli3Rep assure that the target genes are kept silent. Our aim was to reproduce this opposite behavior of Gli1Act and Gli3Rep concentrations (see Figs. 3E and G). This feature of the pathway lead us to select the value of the rate of Gli3 synthesis reported in Table A.1. At t = 0 (see Table A.2), when there is no Shh reaching target cells, the concentrations of [Gli1Act ] is zero, as gene expression starts upon activation of the signal, but there must be preexisting Ptc1mem for it to act as the Shh/Hh receptor. The initial value of Ptcmem concentration was chosen in order to have a more than twofold increase upon signaling. Numerical solutions resulted to be insensitive to a change of the initial value of [Gli3Act ] between 0.5 and 10 nM. Initial value of [Gli3Rep ] was taken from [7].

M. Verbeni et al. / Physics of Life Reviews 10 (2013) 457–475

471

Table A.1 Parameter used in numerical simulations: definitions and values. Parameter

Description

Value

Source

c

Front propagation speed in FRAPS

mean value = 0.03 µm/s

Front propagation speed for FLS model Kinematic viscosity Value of the flux at the floor plate Dissociation constant for Shh-Ptc binding Backward rate constant for Shh-Ptc dissociation Forward rate constant for Shh-Ptc association Rate constant for internalization of the complex Shh·Ptc Rate constant for recycling of the complex Shh·Ptc to cellular surface Degradation rate constant for intracellular complex Shh·Ptc Rate constant for transport of newly produced Ptc to cellular membrane Rate of Gli1Act synthesis Degradation rate constant for all Gli proteins Rate of Ptc1 synthesis Rate of Gli3 synthesis Rate constant for the conversion of Gli3 to Gli3Rep Dissociation constant for Gli3 binding to DNA binding site Dissociation constant for Gli1 binding to DNA binding site Transcriptional repression Binding cooperativity

1.3 × 10−3 µm/s 5 × 10−9 cm2 /s 0.37 × 10−13 M cm s−1 1.5 × 10−9 M 0.1 min−1

Our measurements (Fig. 1D) [56,57] notes notes [69] [7]

6.67 × 107 M−1 min−1

koff /KShh

0.03 min−1 0.00402 min−1

0.03–0.3 min−1 for EGF [86] [5] for Dpp

0.00198 min−1

[7]

0.00036 min−1

[7]

2.74 × 10−11 M min−1 0.009 min−1

notes [7]

1.5 × 10−10 M min−1 1.5 × 10−10 M min−1 0.0117 min−1

notes notes [7]

8.3 × 10−10 M

[7]

KGli1 = KGli3

[7]

0.01 1.8

notes notes

ν β KShh koff kon kCin kCout kCdeg kcyt kG kdeg kP γg3 kg3r KGli3 KGli1 r 

Table A.2 Initial values of the concentrations. Variable

Value at t = 0

[Shh] [PtcShhmem ] [Ptcmem ] [PtcShhcyt ] [Ptccyt ] [Gli1Act ] [Gli3Act ] [Gli3Rep ]

0 nM 0 nM 0.21 nM 0 nM 0.01 nM 0 nM 1 nM 61.2 nM

Appendix B. The weight term Ptr We calculated the weight term Ptr appearing in the evolution equations for [Gli1Act ] and [Ptc1int ] employing the BEWARE method. This technique allows to evaluate the probability corresponding to transcriptional events occurring through the interactions of regulatory proteins (in our case Gli1Act , Gli3Act and Gli3Rep ) with specific DNA sequences (in our case, the regulatory region of the Shh-Gli target genes gli1 and ptc1). The first step consists in enumerating the possible configurations of transcription factors bound to DNA binding sites. Following Saha and Schaffer [7], we considered three binding sites, both for gli1 and ptc1 genes, resulting in 20 possible configurations.

472

M. Verbeni et al. / Physics of Life Reviews 10 (2013) 457–475

Each configurations is associated with a probability of transcription, which can be zero, if the configuration does not lead to gene transcription. This is the case, for example, of a configuration consisting only of repressors Gli3Rep or of a configuration with two repressors Gli3Rep and one activator Gli3Act , as the transcription factor Gli3Act has a weak activator function. The transcription probability corresponding to each configuration is evaluated as the relative concentration of regulatory proteins at stationary state. This fraction is calculated using statistical thermodynamics and depends on free energies of binding for each configurations, that is on dissociation constants KGli1 and KGli3 for binding of Gli transcription factors with DNA sites. In our model we consider KGli1 = KGli3 . Furthermore the relative concentration is multiplied by a repressor factor r < 1 for each repressor bound and by a factor  > 1 if two or more activator factor Gli1Act are bound, to take into account the cooperative [65] Gli1 function. In general, for a given configuration s, the expression for the corresponding transcription probability reads: Ps =

r k  (i−1+δi0 ) [Gli1Act ]i [Gli3Act ]j [Gli3Rep ]k , 20 l=1 Pl

where l ranges over the 20 possible configurations. The overall transcription probability Ptr is obtained summing over all the probability corresponding to configurations which actually contributes to gene activation. The explicit expression is given by: Ptr =

Prel Ztot

where 3 



 

Prel = 3  2 Gli1Act + Gli3Act + r Gli3Rep + KGli3 2

2





 

×  Gli1Act + Gli3Act + Gli1Act Gli3Act + r Gli3Rep + KGli3 and

3 3

2



2 Gli3Act Ztot = 3 2 Gli1Act + 3 Gli3Act + 3KGli3 Gli3Act + 3KGli3   2



3 + 3r Gli3Rep KGli3 + Gli3Act + KGli3

2  

+ 3r 2 Gli3Rep KGli3 + Gli3Act

3 2 





Gli3Act + r Gli3Rep + KGli3 + 3r 3 Gli3Rep + 3 Gli1Act 



 2

+ 3 Gli1Act Gli3Act + r Gli3Rep + KGli3 . In the equations defining the evolution of [Gli1Act ] and [Ptc1int ], concentrations [Gli3Act ] and [Gli3Rep ] are evaluated at time t, while [Gli1Act ] is evaluated at the delayed timed t − τ . Appendix C. Flux-limiters and optimal mass transportation At the same date (1992) that Rubin [87] questioned the Cattaneo law as an alternative to Fickian diffusion, proving that it violates the second principle of the thermodynamics, Ph. Rosenau [42] derived the flux-limited equation from the observation that the speed of sound is the highest admissible free velocity in a medium. As we have commented previously, there are different approaches to deduce the flux limited terms, in addition to the mentioned above [42]: From kinetic theory of multicellular growing system [47], from Hilbert methods in astrophysics [88,89], from optimal mass transportation [45,90], . . . We have chosen here an approach that might be less known: The optimal mass transportation. We are concerned with the optimal transport of masses from one location to the other (or the transport map from a mass density distribution to another one), where the optimality depends upon the context and then a cost function is applied to this transport process. The problems appear in several areas: Economics, probability, optimization, biology, meteorology and computer graphics, among others. The flux-limited equation can be derived [90] by means of an optimal mass transport theory, as a gradient flow of the Boltzmann entropy for the metrics corresponding to the cost function

M. Verbeni et al. / Physics of Life Reviews 10 (2013) 457–475

k(z) =

c2 (1 −

 1−

|z|2 ), c2

+∞,

473

if |z|  c, if |z| > c.

The idea of optimal mass transportation consists in defining the associated Wasserstein distance between two probability distributions ρ0 and ρ1 by      x −y h Wk (ρ0 , ρ1 ) := inf k dγ (x, y) , γ ∈Γ (ρ0 ,ρ1 ) h RN ×RN

where h > 0 and Γ (ρ0 , ρ1 ) is the set of probability measures in RN × RN whose marginals are ρ0 and ρ1 . Let F : [0, ∞) → [0, ∞) be a convex function and let P(RN ) be the set of probability density functions ρ : RN → [0, ∞). Starting from ρ0h := ρ0 ∈ P(RN ), we solve iteratively    h h inf hWk ρn−1 , ρ + F ρ(x) dx. ρ∈P (RN )

RN

This is a gradient descent with respect to the Wasserstein distance. Define ρ h (t) := ρnh for t ∈ [nh, (n + 1)h). The solution converges as h → 0+ to the solution of the nonlinear partial differential equation      ∂u |ξ |2 ∗ 2 1+ 2 −1 . = div u∇k ∗ ∇F (u) , k (ξ ) = c ∂t c Choosing F (r) = ln r, we deduce the relativistic heat equation   ∂u u∇u = ν div  . 2 ∂t |u|2 + νc2 |∇u|2 References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]

Turing AM. The chemical basis of morphogenesis. Philos Trans R Soc Lond B, Biol Sci 1952;237:37–72. Crick F. Diffusion in embryogenesis. Nature 1970;40:561–3. Meinhardt H. Space-dependent cell determination under the control of a morphogen gradient. J Theor Biol 1978;74:307–21. Wolpert L. Positional information and the spatial pattern of cellular differentiation. J Theor Biol 1969;25:1–47. Lander AD, Nie Q, Wan FY-M. Do morphogen gradients arise by diffusion? Dev Cell 2002;2:785–96. Lai K, Robertson MJ, Schaffer DV. The Sonic Hedgehog signaling system as a bistable genetic switch. Biophys J 2004;86:2748. Saha K, Schaffer DV. Signaling dynamics in Sonic Hedgehog tissue patterning. Development 2006;133:889–900. Kicheva K, Pantazis P, Bollenbach T, Kalaidzidis Y, Bittig T, Julicher F, et al. Science 2007;315:521. Neugebauer JM, Amack JD, Peterson AG, Bisgrove BW, Yost HJ. FGF signaling during embryo development regulates cilia length in diverse epithelia. Nature 2009;458:651–4. Schwank G, Restrepo S, Basler K. Growth regulation by Dpp: an essential role for Brinker and a non-essential role for graded signaling levels. Development 2008;135:4003–13. Müller P, Rogers KW, Jordan BM, Lee JS, Robson D, Ramanathan S, et al. Differential diffusivity of nodal and lefty underlies a reaction– diffusion patterning system. Science 2012;336:721–4. Harvey SA, Smith JC. Visualisation and quantification of morphogen gradient formation in the Zebrafish. PLoS Biol 2009;7(5):e1000101. Smith JC. Forming and interpreting gradients in the early xenopus embryo. Cold Spring Harb Perspect Biol 2009:1–12. Kondo S, Miura T. Reaction–Diffusion model as a framework for understanding biological pattern formation. Science 2010;329:1616–20. Roth S. Mathematics and biology: a Kantian view on the history of pattern formation theory. Dev Genes Evol 2011;221(5–6):255–79. John F. Partial differential equations. Springer; 1971. Einstein A. Zür Theorie der Brownschen Bewegung. Ann Phys 1906;19:371–81. Dessaud E, Yang LL, Hill K, Cox B, Ulloa F, Ribeiro A, et al. Interpretation of the Sonic Hedgehog morphogen gradient by a temporal adaptation mechanism. Nature 2007;450:717–20. Scherz PJ, McGlinn E, Nissim S, Tabin CJ. Extended exposure to Sonic Hedgehog is required for patterning the posterior digits of the vertebrate limb. Dev Biol 2007;308(2):343–54. Harfe BD, Scherz PJ, Nissim S, Tian H, McMahon AP, Tabin CJ. Evidence for an expansion-based temporal Shh gradient in specifying vertebrate digit identities. Cell 2004;118(4):517–28. Kornberg TB. The imperatives of context and contour for morphogen dispersion. Biophys J 2012;103:2252–6. Jessell TM. Neuronal specification in the spinal cord: inductive signals and transcriptional codes. Nat Rev Genet 2000;1:20–9.

474

M. Verbeni et al. / Physics of Life Reviews 10 (2013) 457–475

[23] Stecca B, Ruiz i Altaba A. Context-dependent regulation of the GLI code in cancer by Hedgehog and non-Hedgehog signals. J Mol Cell Biol 2010;2(2):84–95. [24] Gritli-Linde A, Lewis P, McMahon AP, Linde A. The whereabouts of a morphogen: direct evidence for short- and graded long-range activity of Hedgehog signaling peptides. Dev Biol 2001;236(2):364–86. [25] Lek M, Dias JM, Marklund U, Uhde CW, Kurdija S, Lei Q, et al. A homeodomain feedback circuit underlies step-function interpretation of a Shh morphogen gradient during ventral neural patterning. Development 2010;137:4051–60. [26] Parker DS, White MA, Ramos AI, Cohen BA, Barolo S. The cis-regulatory logic of Hedgehog gradient responses: key roles for gli binding affinity, competition, and cooperativity. Sci Signal 2011;4(176):ra38. [27] Litingtung Y, Dahn RD, Li Y, Fallon JF, Chiang C. Shh and Gli3 are dispensable for limb skeleton formation but regulate digit number and identity. Nature 2002;418:979–83. [28] te Welscher P, Fernandez-Teran M, Ros MA, Zeller R. Mutual genetic antagonism involving GLI3 and dHAND prepatterns the vertebrate limb bud mesenchyme prior to Shh signaling. Genes Dev 2002;16:421–6. [29] Guerrero I, Chiang C. A conserved mechanism of Hedgehog gradient formation by lipid modifications. Trends Cell Biol 2007;17:1–5. [30] Vyas N, Goswami D, Manonmani A, Sharma P, Ranganath H, VijayRaghavan K, et al. Nanoscale organization of Hedgehog is essential for long-range signaling. Cell 2008;133:1214–27. [31] Zeng X, Goetz JA, Suber LM, Scott Jr WJ, Schreiner CM, Robbins DJ. A freely diffusible form of Sonic Hedgehog mediates long-range signalling. Nature 2011;411:716–20. [32] Greco V, Hannus M, Eaton S. Argosomes: a potential vehicle for the spread of morphogens through epithelia. Cell 2001;106:633–45. [33] Panáková D, Sprong H, Marois E, Thiele C, Eaton S. Lipoprotein particles are required for Hedgehog and Wingless signalling. Nature 2005;435:58–65. [34] Callejo A, Culi J, Guerrero I. Patched, the receptor of Hedgehog, is a lipoprotein receptor. Proc Natl Acad Sci USA 2008;105:912–7. [35] Ruiz i Altaba A. Coexpression of HNF-3β and Isl-1/2 and mixed distribution of ventral cell types in the early neural tube. Int J Dev Biol 1996;40:1081–8. [36] Gurdon JB, Mitchell A, Mahony D. Direct and continuous assessment by cells of their position in a morphogen gradient. Nature 1995;376:520–1. [37] Meinhardt H, Prusinkiewicz P, Fowler DR. The algorithmic beauty of sea shells. Berlin: Springer-Verlag; 1998. [38] Ramirez-Weber FA, Kornberg TB. Cytonemes: cellular processes that project to the principal signaling center in Drosophila imaginal discs. Cell 1999;97(5):599–607. [39] Roy S, Hsiung F, Kornberg TB. Specificity of Drosophila cytonemes for distinct signaling pathways. Science 2011;15:354–8. [40] Gradilla AC, Guerrero I. Cytoneme-mediated cell–cell signaling during development. Cell Tissue Res 2013. [41] Sanders TA, Llagostera E, Barna M. Specialized filopodia direct long-range transport of Shh during vertebrate tissue patterning. Nature 2013;497:628–32. [42] Rosenau P. Tempered diffusion: a transport process with propagating front and inertial delay. Phys Rev A 1992;46:7371–4. [43] Chertock A, Kurganov A, Rosenau P. Formation of discontinuities in flux-saturated degenerate parabolic equations. Nonlinearity 2003;16:1875–98. [44] Kurganov A, Rosenau P. On reaction processes with saturating diffusion. Nonlinearity 2006;19:171–93. [45] Brenier Y. Extended Monge–Kantorovich theory. In: Caffarelli LA, Salsa S, editors. Optimal transportation and applications. Lectures given at the C.I.M.E. summer school help in Martina Franca. Lect Notes Math, vol. 1813. Springer-Verlag; 2003. p. 91–122. [46] McCann R, Puel M. Constructing a relativistic heat flow by transport time steps. Ann Inst Henri Poincaré, Anal Non Linéaire 2009;26:2539–80. [47] Bellomo N, Bellouquid A, Nieto J, Soler J. Multiscale biological tissue models and flux-limited chemotaxis for multicellular growing systems. Math Models Methods Appl Sci 2010;20(7):1179–207. [48] Andreu F, Caselles V, Mazón JM, Moll S. The Dirichlet problem associated to the relativistic heat equation. Math Ann 2010;347:135–99. [49] Calvo J, Mazón J, Soler J, Verbeni M. Qualitative properties of the solutions of a nonlinear flux-limited equation arising in the transport of morphogens. Math Models Methods Appl Sci 2011;21:893–937. [50] Caselles V. On the entropy conditions for some flux limited diffusion equations. J Differ Equ 2011;250:3311–48. [51] Campos J, Guerrero P, Sánchez O, Soler J. On the analysis of traveling waves to a nonlinear flux limited reaction–diffusion equation. Ann Inst Henri Poincaré, Anal Non Linéaire 2013;30(1):141–55. [52] Bellomo N, Soler J. On the mathematical theory of the dynamics of swarms viewed as complex systems. Math Models Methods Appl Sci 2012;22:1140006. [53] Bellomo N, Dogbè C. On the modelling of traffic and crowds – a survey of models, speculations, and perspectives. SIAM Rev 2011;53(3):409–63. [54] Schwank G, Dalessi S, Yang SF, Yagi R, de Lachapelle AM, Affolter M, et al. Formation of the long range Dpp morphogen gradient. PLoS Biol 2011;9(7):e1001111. [55] Kerszberg M, Wolpert L. Mechanisms for positional signalling by morphogen transport: a theoretical study. J Theor Biol 1998;191:103–14. [56] Su VF, Jones KA, Brodsky M, The I. Quantitative analysis of Hedgehog gradient formation using an inducible expression system. BMC Dev Biol 2007;7(43):1–15. [57] Chamberlain ChE, Jeong J, Guo Ch, Allen BL, McMahon AP. Notochord-derived Shh concentrates in close association with the apically positioned basal body in neural target cells and forms a dynamic gradient during neural patterning. Development 2008;135:1097–106. [58] Ingham PW, Fietz MJ. Quantitative effects of Hedgehog and decapentaplegic activity on the patterning of the Drosophila wing. Curr Biol 1995;5(4):432–40. [59] Torroja C, Gorfinkiel N, Guerrero I. Patched controls the Hedgehog gradient by endocytosis in a dynamin-dependent manner, but this internalization does not play a major role in signal transduction. Development 2004;131:2395–408.

M. Verbeni et al. / Physics of Life Reviews 10 (2013) 457–475

475

[60] Casali A. Self-Induced Patched receptor down-regulation modulates cell sensitivity to the Hedgehog morphogen gradient. Sci Signal 2010;3(136):1–8. [61] Chen Y, Struhl G. Dual roles for patched in sequestering and transducing Hedgehog. Cell 1996;87(3):553–63. [62] Aza-Blanc P, Ramírez-Webber FA, Laget MP, Schwartz C, Kornberg TB. Proteolysis that is inhibited by Hedgehog targets Cubitus interruptus protein to the nucleus and converts it to a repressor. Cell 1997;89:1043–53. [63] Alvarez-Medina R, Cayuso J, Okubo T, Takada S, Martí E. Wnt canonical pathway restricts graded Shh/Gli patterning activity through the regulation of Gli3 expression. Development 2008;135:237. [64] Lee J, Platt KA, Censullo P, Ruiz i Altaba A. Gli1 is a target of Sonic Hedgehog that induces ventral neural tube development. Development 1997;124:2537–52. [65] Ruiz i Altaba A. Combinatorial Gli gene function in floor plate and neuronal inductions by Sonic Hedgehog. Development 1998;125:2203–12. [66] Yao S, Lum L, Beachy PhA. The ihog cell-surface proteins bind Hedgehog and mediate pathway activation. Cell 2006;125:343–57. [67] Aza-Blanc P, Lin HY, Ruiz i Altaba A, Kornberg TB. Expression of the vertebrate Gli proteins in Drosophila reveals a distribution of activator and repressor activities. Development 2000;127(19):4293–301. [68] Casali A, Struhl G. Reading the Hedgehog morphogen gradient by measuring the ratio of bound to unbound Patched protein. Nature 2004;431:76–80. [69] Taipale J, Cooper MK, Maiti T, Beachy PA. Patched acts catalytically to suppress the activity of Smoothened. Nature 2002;418:892–7. [70] Jiang GS, Shu CW. Efficient implementation of weighted ENO schemes. J Comput Phys 1996;126:202–28. [71] Dessaud E, McMahon AP, Briscoe J. Pattern formation in the vertebrate neural tube: a Sonic Hedgehog morphogen-regulated transcriptional network. Development 2008;135:2489–503. [72] Sherer NM, Mothes W. Cytonemes and tunneling nanotubules in cell-cell communication and viral pathogenesis. Trends Cell Biol 2008;18(9):414–20. [73] Bilioni A, Sánchez-Hernández D, Callejo A, Gradilla AC, Ibáñez C, Mollica E, et al. Balancing Hedgehog, a retention and release equilibrium given by Dally, Ihog, Boi and shifted/DmWif. Dev Biol 2013;376(2):198–212. [74] Gerdes HH, Carvalho RN. Intercellular transfer mediated by tunneling nanotubes. Curr Opin Cell Biol 2008;20(4):470–5. [75] Rojas-Ríos P, Guerrero I, González-Reyes A. Cytoneme-mediated delivery of Hedgehog regulates the expression of bone morphogenetic proteins to maintain germline stem cells in Drosophila. PLoS Biol 2012;10(4):e1001298. [76] Beachy PhA, Hymowitz SG, Lazarus RA, Leahy DJ, Siebold C. Interactions between Hedgehog proteins and their binding partners come into view. Genes Dev 2010;24:2001–12. [77] Clement V, Sanchez P, De Tribolet N, Radovanovic I, Ruiz i Altaba A. Hedgehog-GLI signaling regulates human glioma growth, cancer stem cell self-renewal and tumorigenicity. Curr Biol 2007;17:165–72. [78] Zbinden M, Duquet A, Lorente-Trigos A, Ngwabyt SN, Borges I, Ruiz i Altaba A. NANOG regulates glioma stem cells and is essential in vivo acting in a cross-functional network with GLI1 and p53. EMBO J 2010;29:2659–74. [79] Ochs S. Fast transport of materials in mammalian nerve fibers. Science 1972;176(4032):252–60. [80] Meinhardt H. Models for the generation and interpretation of gradients. Cold Spring Harb Perspect Biol 2009;1:a001362. [81] Ruiz i Altaba A. Catching a Gli-mpse of Hedgehog. Cell 1997;90:193–6. [82] Callejo A, Bilioni A, Mollica E, Gorfinkiel N, Andrés G, Ibáñez C, et al. Dispatched mediates Hedgehog basolateral release to form the long-range morphogenetic gradient in the Drosophila wing disk epithelium. Proc Natl Acad Sci USA 2011;108(31):12591–8. [83] Udolph G, Luer K, Bossing T, Technau GM. Commitment of cns progenitors along the dorsoventral axis of Drosophila neuroectoderm. Science 1995;269:1278–81. [84] Spiteri RJ, Ruuth SJ. A new class of optimal high-order strong-stability-preserving time discretization methods. SIAM J Numer Anal 2002;40:469–91. [85] Scorcioni R, Polavaram S, Ascoli GA. L-Measure: a web-accessible tool for the analysis, comparison and search of digital reconstructions of neuronal morphologies. Nat Protoc 2008;3:866–76. [86] Lauffenburger DA, Linderman JJ. Receptors: models for binding, trafficking, and signaling. New York: Oxford University Press; 1993. [87] Rubin MB. Hyperbolic heat conduction and the second law. Int J Engng Sci 1992;30:1665–76. [88] Levermore CD, Pomraning GC. A flux-limited diffusion theory. Astrophys J 1981;248:321–34. [89] Coulombel J-F, Golse F, Goudon Th. Diffusion approximation and entropy based moment closure for kinetic equations. Asymptot Anal 2005;45(1–2):1–39. [90] Agueh M. Existence of solutions to degenerate parabolic equations via the Monge–Kantorovich theory. PhD thesis. Atlanta: Georgia Tech; 2001.