The proneural wave in the Drosophila optic lobe is driven by an excitable reaction-diffusion mechanism

Abstract

In living organisms, self-organised waves of signalling activity propagate spatiotemporal information within tissues. During the development of the largest component of the visual processing centre of the Drosophila brain, a travelling wave of proneural gene expression initiates neurogenesis in the larval optic lobe primordium and drives the sequential transition of neuroepithelial cells into neuroblasts. Here, we propose that this ‘proneural wave’ is driven by an excitable reaction-diffusion system involving epidermal growth factor receptor (EGFR) signalling interacting with the proneural gene l’sc. Within this framework, a propagating transition zone emerges from molecular feedback and diffusion. Ectopic activation of EGFR signalling in clones within the neuroepithelium demonstrates that a transition wave can be excited anywhere in the tissue by inducing signalling activity, consistent with a key prediction of the model. Our model illuminates the physical and molecular underpinnings of proneural wave progression and suggests a generic mechanism for regulating the sequential differentiation of tissues.

https://doi.org/10.7554/eLife.40919.001

Introduction

The development of multicellular organisms relies on a multitude of transient coordination processes that provide the spatiotemporal cues for cell fate decision-making and thereby ensure that tissues are specified with the correct size, pattern and composition (Perrimon et al., 2012; Oates et al., 2012; Sato et al., 2013). In one strategy, large-scale patterning is engineered by self-organised concentration waves of biomolecular fate determinants that travel across tissues through intercellular exchange and the regulation of gene expression. Such travelling waves, which are viable carriers of spatiotemporal information, are a ubiquitous feature of developmental pattern formation, where they arise through different underlying mechanisms, from coordinated intracellular oscillations (Oates et al., 2012; Jörg et al., 2015; Hubaud et al., 2017; Verd et al., 2018) to self-organised reaction-diffusion processes (Lubensky et al., 2011; Formosa-Jordan et al., 2012; Fried et al., 2016; Gavish et al., 2016; Corson et al., 2017).

During the development of the fruitfly Drosophila melanogaster, a propagating wave of gene expression orchestrates the patterning of the largest component of the visual processing centre: neuroepithelial cells in the optic lobe of the larval brain divide symmetrically, expanding the progenitor pool, and then undergo a sequential transition into asymmetrically dividing neuroblasts, which generate the neurons of the medulla (Figure 1a,b) (Egger et al., 2007; Yasugi et al., 2008; Egger et al., 2010; Yasugi et al., 2010). The transition from neuroepithelial cells to neuroblasts occurs at a ‘transition zone’ that sweeps from one side of the optic lobe to the other and is marked by the expression of the proneural gene lethal of scute (l’sc) (Yasugi et al., 2008). This ‘proneural wave’ is controlled by the coordinated action of different signalling pathways: epidermal growth factor receptor (EGFR) signalling and Delta-Notch signalling (Figure 1c–f) (Yasugi et al., 2010). The sequential nature of the transition is crucial to generate populations of cells of different developmental ages that give rise to a diverse array of terminally differentiated medulla neurons (Li et al., 2013; Sato et al., 2013; Suzuki et al., 2013; Erclik et al., 2017). The transition zone exhibits localised EGFR signalling as well as expression of l’sc (Figure 1b,c). Absence of EGFR signalling leads to loss of the differentiation wave, indicating that EGFR signalling is a key component for proneural wave progression (Yasugi et al., 2010). The neuroepithelium exhibits low levels of Notch signalling activity (Egger et al., 2011; Weng et al., 2012). However, Notch activity peaks directly before the transition from neuroepithelial cell to neuroblast, drops during the transition and then is restored upon neuroblast transformation (Figure 1c) (Contreras et al., 2018). In addition, coordinating roles are played by the JAK/STAT and Fat-Hippo pathways, which are broadly expressed in the neuroepithelium and prevent premature and ectopic transition of the neuroepithelium (Yasugi et al., 2008; Yasugi et al., 2010; Wang et al., 2011a; Reddy et al., 2010; Kawamori et al., 2011; Weng et al., 2012; Tanaka et al., 2018).

Molecular basis for the proneural wave in the Drosophila optic lobe.

(a) Schematic depiction of the Drosophila larva at the late 3rd instar stage when the proneural wave is transforming the neuroepithelium into medulla neuroblasts. (b) Optic lobe in a lateral view showing the neuroepithelium (labelled with Notch intracellular domain (NICD), white), the transition zone (L’sc, cyan) and the neuroblasts (Dpn, red). (c) L’sc expression and Notch signalling activity around the transition zone. Top: Magnification of the region outlined in (b), showing neuroblasts (Dpn, red), L’sc expression (cyan) and the neuroepithelium (NICD, white). Middle: Confocal image showing that Notch signalling activity (HLH-mgamma, purple) increases just before the transition zone (marked by L’sc, cyan), drops during the transition and then increases again in neuroblasts. Bottom: The proneural wave, characterised by expression of L’sc as well as EGF receptor (EGFR) and Notch signalling activity, sequentially converts the neuroepithelium into neuroblasts. (d) EGFR signalling in the transition zone activates expression of the transmembrane protein Rhomboid, which in turn cleaves the membrane-tethered form of the EGFR ligand Spitz (mSpi) to generate its active secreted form (sSpi). (The shaded region depicts an individual cell in the neuroepithelium.) sSpi can bind to the EGFR on the same cell and neighbouring cells. (e) Delta-Notch signalling is a contact-dependent signalling pathway active in both the neuroepithelium and the neuroblasts. The Delta ligand binds to Notch receptors on adjacent cells upon which their intracellular domain (NICD) is cleaved. The NICD regulates target genes, which, in turn, affects expression of Delta. (f) Active EGFR signalling promotes the expression of L’sc within the same cell, which is sufficient for the neuroepithelium to neuroblast transition and which in turn downregulates EGFR signalling.

https://doi.org/10.7554/eLife.40919.002

The question of how the specific functional feedbacks of EGFR signalling and proneural gene expression generate a localised propagating transition zone requires a mechanistic explanation of wave progression based on molecular feedbacks and signalling cascades. Such a description should explain (i) the dynamic nature of the wave, (ii) the emergence of a localised transition zone with spatially confined expression of the proneural gene l’sc and (iii) the specific profiles of gene expression and signalling activity around the transition zone. Moreover, the nature and function of the interaction of these components with Delta-Notch signalling, more commonly associated with lateral inhibition of neighbouring cells, is poorly understood, see Appendix 3. While a recent effort of a phenomenological description of the proneural wave (Sato et al., 2016) has started to model the coarse-grained aspects of proneural wave progression, the emergence of some major characteristics of the wave (such as spatially confined proneural gene expression in a localised transition zone) has not been addressed. Here we propose a model of signalling activity and proneural gene expression that describes the emergence of the proneural wave. Within this framework, the neuroepithelium behaves as an excitable medium in which changes in gene expression at the tissue boundary initiate a spontaneous wave of signalling activity that effects the transition from neuroepithelium to neuroblasts.

Results

Travelling front model of EGFR signalling activity

To develop the model, we first considered interactions between L’sc expression and associated signalling pathways within the transition zone. Previously, it was proposed that sequential induction of EGFR signalling is responsible for the progression of the proneural wave (Yasugi et al., 2008). EGFR signalling activates the expression of L’sc (Yasugi et al., 2010), which is sufficient to drive the neuroepithelium to neuroblast (NE to NB) transition (Yasugi et al., 2008; Contreras et al., 2018). The EGFR is activated by binding the secreted form of its ligand, Spitz. Secreted Spitz is generated by cleavage of a membrane-bound precursor by the transmembrane protease, Rhomboid (Klämbt, 2000). EGFR, Rhomboid and secreted Spitz together form an autocrine positive feedback loop (Figure 1d) (Wiley et al., 2003; Sato et al., 2013). In a first step, we noted that the dynamics of EGFR signalling alone has features that are sufficient to enable such a sequential induction and produce a travelling front of EGFR signalling activity; a feature notably absent in recent attempts to model the proneural wave, which also require further components to stabilise the propagating EGFR signalling front (Sato et al., 2016). In a minimal model based on the EGFR/Rhomboid/Spitz positive feedback loop, EGFR signalling activity is represented by a single component ‘E’ (e.g., the local cellular concentration of the active form of Spitz) that is diffusible between cells and involves the aforementioned positive feedback (Figure 1d and Figure 2a; Appendix 1),

(1) ϕEt=η2ϕE+μh(ϕE)-kϕE.
Dynamics of the proneural wave as an excitable reaction-diffusion system.

(a) Minimal model of EGFR signalling. The dynamics of the EGFR/Rhomboid/Spitz feedback loop is condensed in a single component ‘E’ (green), which is diffusible between cells and able to self-activate (for details, see Appendix 1). This single component represents a proxy for the activity of the feedback loop shown in Figure 1c, for example the local concentration of the active form of Spitz. The corresponding reaction-diffusion system, Equation 1, can give rise to a propagating front that leaves behind an elevated signalling state. Plots show the numerical solution of Equation 1 in a one-dimensional representation of tissue (for simplicity) at two time points with initially elevated levels of E at the left-hand side of the domain. Specifying position in units of the diffusion length η/k and time in units of the decay time k-1, the remaining chosen parameters are μ=4 and n=3. (b) Model of EGFR signalling interacting with the proneural gene l’sc. EGFR signalling activates L’sc expression (component ‘L’, blue), which effectively inhibits EGFR signalling by driving the NE to NB transition (for details, see Appendix 2). The corresponding reaction-diffusion system, Equation 2, can give rise to a propagating localised pulse of signalling activity and proneural gene expression corresponding to the transition zone. Parameters for E are the same as in panel a; parameters for L are μL=0.4, kL=0.2. (c) Schematic depiction of the mechanism giving rise to a localised transition zone, shown in panel b. Diffusion of signalling components (E, green) into the neuroepithelium leads to activation of the positive feedback loop, which locally excites signalling and proneural gene expression (L, blue) (Materials and methods). The excitation ceases as downregulation of signalling occurs, a consequence of the transition triggered by L’sc expression. (d) Regulatory network of the refined model including Delta-Notch signalling (D and N) and a local variable Ω indicating the cell state (Ω=0 indicates neuroepithelial cells and Ω=1 indicates neuroblasts; for details, see Appendix 3). Each shaded cell indicates one lattice site corresponding to one cell of the tissue. (e) Simulation of the integrated model of EGFR signalling, L’sc expression, Delta-Notch signalling and the NE to NB transition in a one-dimensional array of cells. The emerging spatial signalling and gene expression profile is characterised by a pulse of EGFR signalling, L’sc and Delta, and a drop in Notch signalling activity within the transition zone. The drop in Notch is preceded by a pulse of Notch signalling activity (pink arrowheads), which is due to a local lateral inhibition effect mediated by Delta-Notch signalling. Parameters are given in Appendix 3—table 1 except for η=0.03.

https://doi.org/10.7554/eLife.40919.003

Here ϕE denotes the local strength of E activity, η is the effective diffusion constant, μ is the gain rate in signalling activity due to positive feedback, k is the decay rate, and h(ϕ)=ϕn/(1+ϕn) is a Hill function parameterising the nonlinear positive feedback. Generically, reaction-diffusion systems involving diffusion and self-activation are known to support travelling bistable fronts that leave behind an elevated signalling state (Figure 2a; Video 1; Appendix 1) (Muratov and Shvartsman, 2004; Keener and Sneyd, 2009; Graham et al., 2010; Tayar et al., 2015). EGFR signalling is a natural candidate for being a key driver of the proneural wave and experimental evidence has shown that EGFR signalling is both necessary and sufficient for wave progression (Yasugi et al., 2010).

Video 1
Travelling EGFR signalling front.

The movie shows the simulation of a one-dimensional version of the EGFR signalling model Equation 1 corresponding to the snapshots shown in Figure 2a. All simulation parameters as in Figure 2a.

https://doi.org/10.7554/eLife.40919.004

Travelling pulse model of EGFR signalling and proneural gene expression

However, notably, the EGFR/Rhomboid/Spitz feedback loop does not remain active in the wake of the travelling wave, but remains spatially confined as Rhomboid is expressed only transiently in the travelling transition zone (Yasugi et al., 2010; Sato et al., 2013). Therefore, in a second step, we considered the influence of the proneural gene l’sc, represented by a second component ‘L’ in our model. Elevated EGFR signalling activates L’sc expression (Yasugi et al., 2010), which is sufficient to drive the NE to NB transition (Yasugi et al., 2008; Contreras et al., 2018). In this minimal model, L’sc downregulates EGFR signalling as a consequence of the transition, leading to an indirect negative feedback (Figure 1f). The corresponding reaction-diffusion system for the local strength of E and L activity, ϕE and ϕL, is given by

(2) ϕEt=η2ϕE+μEh(ϕE)[1h(ϕL)]kEϕE ,ϕLt=μLh(ϕE)kLϕL ,

where μi (with i=E,L) indicate production rates and ki denote decay rates. Simulations of Equation 2 demonstrate that this type of feedback is sufficient to describe a travelling localised pulse of signalling activity and L’sc expression through the tissue (Figure 2b,c; Video 2; Appendix 2). Notably, the dynamics of L’sc (L) alter the bistable signalling behaviour of EGFR signalling (E) into an excitable one: once sufficiently perturbed by diffusion from an adjacent cell with an elevated signalling state, the intracellular reaction dynamics produces a transient expression pulse that downregulates itself as a result of the NE to NB cell fate transition (Appendix 2).

Video 2
Travelling EGFR signalling pulse and proneural gene expression.

The movie shows the simulation of a one-dimensional version of the model of EGFR signalling interacting with the proneural gene l’sc, Equation 2, corresponding to the snapshots shown in main text Figure 2b. All simulation parameters as in Figure 2b.

https://doi.org/10.7554/eLife.40919.005

Integrated model of the proneural wave to include EGFR-L’sc-Notch interactions

We next aimed to develop a more refined model that could be challenged by experiment and compared with previously published data. Such a model necessarily includes Delta-Notch signalling, which has been shown to influence how long cells remain in the L’sc expressing state (Yasugi et al., 2010; Wang et al., 2011b; Weng et al., 2012). As a mediator of lateral inhibition, Delta-Notch signalling is often associated with the emergence of ‘salt-and-pepper’-like patterns of cell fate (Bray, 2006; Shaya and Sprinzak, 2011). However, this pattern is not seen during proneural wave progression and the reasons for this are not clear (Egger et al., 2010; Pérez-Gómez et al., 2013; Sato et al., 2016). To address this question, we extended our minimal model to include canonical Delta-Notch interactions (Figure 2d) (Collier et al., 1996; Bray, 2006; Simakov and Pismen, 2013): (i) trans-activation of Notch by Delta, (ii) downregulation of Delta by Notch within the same cell and (iii) cis-inhibition (downregulation of Notch by Delta in the same cell). The model incorporates interactions between the Delta-Notch signalling pathway, EGFR signalling and L’sc expression, namely, upregulation of Delta through EGFR signalling (Yasugi et al., 2010), upregulation of EGFR signalling through Notch signalling (Yasugi et al., 2010), downregulation of L’sc through Notch signalling (Reddy et al., 2010) and downregulation of Notch expression through L’sc (Egger et al., 2010). Despite these complex interactions, the functional ‘module’ comprising EGFR signalling and L’sc expression still remains the driver of the wave (green box in Figure 2d), while Delta-Notch signalling acts to provide further timing cues for the transition and to prevent premature differentiation (pink box in Figure 2d) (Egger et al., 2010; Reddy et al., 2010; Yasugi et al., 2010). The integrated model also includes an explicit representation of the cell state dynamics during the transition between neuroepithelial cell to neuroblast. The NE to NB transition is promoted by L’sc (Yasugi et al., 2008) and downregulation of Notch in the presence of EGFR signalling (Yasugi et al., 2010; Weng et al., 2012). The mathematical details of the refined model are given in Appendix 3.

Congruence with experimental data

In addition to the emergence of a propagating transition zone, the integrated model also yielded predictions on the spatial profiles of signalling activity and gene expression (Figure 2e; Figure 4a; Video 3; Video 4; Appendix 4), which were in striking agreement with features observed in prior experiments. First, EGFR signalling, as well as L’sc and Delta expression, was found to be elevated only in the transition zone (Figure 1c) (Egger et al., 2010; Yasugi et al., 2010). Second, a peak of Notch activity is observed slightly in advance of the transition zone (pink arrowheads in Figure 2e) followed by a sharp drop in Notch activity (Figure 1c) (Egger et al., 2011; Orihara-Ono et al., 2011; Weng et al., 2012; Contreras et al., 2018). According to the model, lateral inhibition promotes high-Delta/low-Notch and low-Delta/high-Notch states in adjacent cells, and leads to a travelling ‘laterally inhibited’ cell state as the wave progresses. By contrast, the drop in Notch levels in transitioning cells arises due to cis-inhibition in our model as a consequence of Delta binding to Notch within the same cell, as has been shown experimentally (Reddy et al., 2010; Weng et al., 2012; Contreras et al., 2018).

Video 3
Travelling proneural wave in the integrated model on a 1D array.

The movie shows a simulation of the proneural wave model Equations 17–19 simulated on a one-dimensional array of cells. All simulation parameters as in Figure 2e.

https://doi.org/10.7554/eLife.40919.006
Video 4
Travelling proneural wave in the integrated model on a 2D hexagonal lattice.

The movie shows a simulation of the proneural wave model Equations 17–19 simulated on a hexagonal lattice with circular geometry with a radius of 15 lattice sites. All simulation parameters are given in Appendix 3—table 1.

https://doi.org/10.7554/eLife.40919.007

To challenge the model further, we checked whether the documented effects of misregulating EGFR signalling, Notch signalling or L’sc expression (Yasugi et al., 2010) could be reproduced. For example, clones in which EGFR signalling has been constitutively activated tend to advance the transition zone within the clone, while the absence of EGFR signalling leads to loss of the proneural wave (Figure 3a,b) (Yasugi et al., 2010). To this end, we simulated the model dynamics on a two-dimensional hexagonal lattice that mimics the topology of the neuroepithelium. Consistent with experiment, the model captured the acceleration, delay or loss of the proneural wavefront within a clone depending upon its genetic makeup (Figure 3a–f; Appendix 6).

Simulations of the proneural wavefront as well as clones (outlined cells) capturing different mutant and transgenic conditions.

(a) Knockout of EGFR signalling (μE=0 within the clone), (b) EGFR signalling constitutively active (E signalling always active within the clone), (c) L’sc knockout (μL=0 within the clone), (d) L’sc constitutively active (L synthesis always active within the clone), (e) Notch downregulation (β=0 within the clone), (f) Notch upregulation (additional N synthesis with rate β/2 within the clone). In all panels, white arrowheads indicate advancements and retardations of the wavefront as compared to wildtype tissue due to the respective genetic alterations of the clones. The system given by Equations 17–19 was numerically simulated on a 20×20 hexagonal lattice with initially localised levels of E in the first three columns at the left boundary of the system so that the wave travels to the right. All other parameters are given in Appendix 3—table 1. All shown simulation snapshots are taken at time t=25, except for panel B, which is taken at t=17.5. The column ‘Experiment’ shows sketches of experiments with mutant and transgenic clones and animals and refer to the corresponding original literature.

https://doi.org/10.7554/eLife.40919.008

Low levels of Notch signalling activity are observed both in the neuroepithelium and in neuroblasts but not at the transition zone (Figure 1c). It has been shown experimentally that Notch activity is required to maintain neuroepithelial cell fate and the loss of Notch results in premature transformation into neuroblasts (Egger et al., 2010; Ngo et al., 2010; Reddy et al., 2010; Yasugi et al., 2010; Orihara-Ono et al., 2011; Wang et al., 2011b; Pérez-Gómez et al., 2013). Intriguingly, despite the observation of active Notch signalling, there is no evidence of lateral inhibition in the neuroepithelium. Lateral inhibition causes neighbouring cells to acquire complementary cell fates and so results in the emergence of a ‘salt-and-pepper pattern’ of Notch signalling activity (Bray, 2006; Shaya and Sprinzak, 2011). The reason for the absence of salt-and-pepper Notch signalling in the neuroepithelium is not clear. Notably, our model predicts that the basal level of Notch activity observed in the the neuroepithelium could be the reason for the suppression of lateral inhibition patterns outside the transition zone. In our model, basal levels of Notch activity in the neuroepithelium lead to a spatially homogeneous ‘oversaturation’ that prevents the Delta levels from rising before being activated by EGFR signalling. This is the case even in the presence of biochemical fluctuations (Figure 4a). However, if basal Notch levels are lowered to small values compared to the threshold levels for activation and inhibition in our model, we indeed recapitulate the salt-and-pepper patterns that are a consequence of lateral inhibition (Figure 4b). An analytical argument for the suppression of lateral inhibitions through basal Notch activity is given in Appendix 5.

Basal Notch activity suppresses lateral inhibition patterns.

The panels show snapshots of the proneural wave model Equation 17Equation 20 simulated on a hexagonal lattice with circular geometry with a radius of 15 lattice sites. White arrows indicate the direction of wave progression. Transient lateral inhibition patterns can occur if basal Notch levels are low compared to the thresholds for activation and inhibition of the Delta–Notch interactions: (a) In the scenario with basal Notch activity, lateral inhibition patterns are suppressed (basal Notch gain rate β=10). (b) In the scenario with downregulated basal Notch activity, lateral inhibition patterns appear (basal Notch gain rate β=1). Other parameters are given in Appendix 3—table 1; both panels are simulated with biochemical noise strength γ/μE=0.5 (see Equation 20). Initial conditions were localised elevated levels of E in those outer boundary cells that have angles between π/3 and 5π/3 as measured from the center of the circular lattice. (c) Downregulation of Notch levels by expressing Notch RNAi in clones does not result in the emergence of a salt-and-pepper expression pattern of Delta (pink). Clone outlines are marked by white dotted lines. (d) Expressing Notch RNAi in clones results in the complete loss of detectable Notch (N(intracellular domain, ICD), purple) within the clones. Clone outlines are marked by white dotted lines. (e) Phase diagram for the occurrence of lateral inhibition in the two-cell system (for details, see Appendix 5). Here, β denotes the basal production rate and λ denotes the gain rate. (c) and (d) are single section confocal images, scale bars represent 20 μm.

https://doi.org/10.7554/eLife.40919.009

We tested this prediction of the model by lowering Notch levels in the neuroepithelium but we did not observe ‘salt-and-pepper’ patterns of Delta/Notch expression within clones expressing an RNAi against Notch (Figure 4c). However, the absence of the emergence of lateral inhibition is likely due to the complete loss of detectable Notch in cells expressing the RNAi (Figure 4d), while the reduction of Notch levels in the model prediction is more subtle (Figure 4b). Referring to the ‘phase diagram’ in Figure 4e, it can be seen that both basal and Delta-regulated Notch activity need to be in the appropriate range for lateral inhibition patterns to occur, which is difficult to achieve experimentally. Furthermore, our model entails that Notch downregulation is a necessary (but not generally sufficient) condition for inducing salt-and-pepper patterns.

Dependence of wave speed and transition zone width on kinetic rate parameters

Next, we asked which aspects of the signalling and gene expression changes in our model have the largest effect on two important features of the system: the speed of the proneural wave and the width of the transition zone. To this end, we performed a sensitivity analysis on the kinetic rate parameters, as detailed in Appendix 7.

This analysis, based on the so-called ‘Morris method’ (Morris, 1991; Campolongo et al., 2007; Wu et al., 2013), entails a resampling of the parameter space of the model and yields three indices for each probed parameter. These indices indicate the impact of each parameter on the assessed output. The Morris indices m and m* describe the impact of the respective parameter on the output, with m including positive and negative effects (which may cancel each other as the parameter is varied) and m* the overall absolute effect (Campolongo et al., 2007). The third index σ measures the non-linearity of the parameter/output relation and/or interactions with other parameters (Wu et al., 2013). Detailed definitions of the respective indices are given in Appendix 7.

Sensitivity analysis of the model.

Plots show the Morris indices m, m* and σ as described in the main text and Appendix 7, indicating the effect of a parameter on (a) the wave speed and (b) the width of the transition zone. The indices m and m* indicate the influence of a parameter on the respective output with m comprising positive and negative effects and m* measuring the absolute effect, whereas non-zero values of σ indicate a nonlinear influence and/or interactions with other parameters. The μi and ki denote the gain and decay rates for the respective components i=E,L,D,N, η denotes the diffusion constant of the component ‘E’ and β denotes the basal Notch gain rate (see Equation 17).

https://doi.org/10.7554/eLife.40919.010

Here we probed the effects of the kinetic rate parameters of EGFR signalling, Delta-Notch signalling and L’sc expression on the propagation speed of the proneural wave and the width of the transition zone. As expected, this analysis showed that the diffusion, gain and decay rate of EGFR signalling (η, μE and kE) are the key regulators of wave speed, since EGFR signalling constitutes the driver of the wave in our model (Figure 5a). In contrast, L’sc gene expression (parametrised by μL and kL) had almost no effect on wave speed since its dynamics is ‘pulled’ by the EGFR signalling front. Interestingly, the basal gain rate of Notch signalling (β) as well as its decay rate (kN) play another prominent role in setting the wave speed. This is consistent with experimental data showing that Notch signalling promotes EGFR signalling at the transition zone, leading to a reinforcement of activation as the wave arrives at undifferentiated cells (Yasugi et al., 2010). Considering the width of the transition zone (Figure 5b), we found that while all parameters had some effect, L’sc expression clearly had a tightening effect as it promoted differentiation and therefore leads to a faster termination of differentiation. In contrast, basal Notch signalling in the epithelium (described by β and kN) tended to enlarge the width of the transition zone in our sensitivity analysis. Indeed, a recent study showed that overexpression of Notch at the transition zone extends the width of the L’sc stripe and delays the transformation into neuroblasts (Contreras et al., 2018), providing support for this prediction of the model. Subsequently, it follows that a complementary prediction of the model would be that a reduction of Notch at the transition zone would decrease the width of the L’sc stripe. We tested this prediction by knocking down Notch at the transition zone in clones (Figure 6). Within clones expressing Notch RNAi, the proneural wave was not only accelerated (Figure 3), as observed previously in Notch mutant clones (Egger et al., 2010; Reddy et al., 2010; Yasugi et al., 2010), but the width of the transition zone also appeared smaller (yellow arrowheads in Figure 6). In summary, this sensitivity analysis suggests that EGFR and Notch signalling are the key regulators of wave speed and width of the transition zone while L’sc expression provides an additional acceleration of the transition of the wave and therefore has a negative influence on the width of the transition zone.

As in Notch mutant clones (Egger et al., 2010; Reddy et al., 2010; Yasugi et al., 2010), the proneural wave is accelerated in clones expressing Notch RNAi, see also Figure 3e.

(a) Expression of Notch RNAi resulted in the downregulation of Notch (N(ICD), purple) and accelerated the transformation (which requires L’sc, cyan) of neuroepithelial cells to neuroblasts (Dpn, red) within clones. (b) The downregulation of Notch appears to decrease the width of the transition zone, as assessed by L’sc (cyan) and Delta (pink) expression within Notch RNAi clones. Dotted white lines mark clone boundaries and yellow arrows indicate the position of the transition zone within Notch RNAi clones. Images are single section confocal slices, scale bars represent 20 μm.

https://doi.org/10.7554/eLife.40919.011

Ectopic excitation of the transition in vivo

A signature feature of the model dynamics is that, through interactions between L’sc and EGFR signalling, the neuroepithelium functions as an excitable medium. As such, the model predicts that local induction of EGFR signalling would initiate a circular (target-like) transition wave (Figure 7a; Video 5). To test whether a transition wave in the neuroepithelium could be excited at a position remote from the proneural wave, we induced clones expressing a downstream effector of the EGFR signalling pathway, Pointed P1 (PntP1; Figure 7b,c; Materials and methods). We found upregulation of L’sc at clonal boundaries and expression of Dpn within the clone, suggesting the ectopic generation of neuroblasts within the epithelium, that is, a NE to NB transition (Figure 7b,c). Our results agree with previous experiments showing the induction of neuroblasts within the neuroepithelium in response to ectopic EGFR signalling (Yasugi et al., 2010) and are in striking agreement with model simulations based on the same perturbation (Figure 7a; Video 5; Supplementary Text).

Video 5
Travelling proneural wave with ectopic activation of EGFR signalling within clones.

The movie shows a simulation of the proneural wave model Equations 17–19 on the same lattice as in Video 3 but with four single-cell clones with constitutively active EGFR signalling. The movie corresponds to the snapshots shown in Figure 7a.

https://doi.org/10.7554/eLife.40919.012
Integrated model of EGFR signalling, L’sc expression, Delta-Notch signalling and the NE to NB transition predicts key features of proneural wave progression in wildtype tissue and following perturbation.

The effect of constitutively active EGFR signalling outside the transition zone: Comparison between the model prediction and experiment. (a) Snapshots of a model simulation on a two-dimensional hexagonal lattice representing the neuroepithelium with randomly distributed clones derived as target waves centred on the site (cell) in which EGFR signalling has been activated (see Materials and methods and Video 5). Cyan indicates levels of L’sc, red indicates neuroblasts (NB) and the grey grid the neuroepithelium (NE). The third column shows a merged image. (b,c) Ectopic expression of PntP1 within the neuroepithelium induces L’sc expression and a NE to NB transition. Clones expressing PntP1, a downstream effector of the EGFR signalling pathway, are indicated by white outlines or green arrowheads. Clones that merge with the transition zone are marked with white arrowheads. Clones within the neuroepithelium that are clearly separated from the transition zone, are marked by green arrowheads; clones that merge with the transition zone are marked with white arrowheads. L’sc is labelled in cyan; neuroblasts are labelled in red by expression of the Hes family transcription factor Deadpan (Dpn).

https://doi.org/10.7554/eLife.40919.013

Discussion

Our findings suggest that the proneural wave involves the activation of an excitatory pulse of signalling activity and gene expression, giving rise to a tightly-regulated propagating transition zone. In contrast with Turing-based activator-inhibitor mechanisms (Turing, 1952), which typically comprise fast diffusible inhibitors, the reaction-diffusion system described here is based on strictly local inhibition (Figure 2b,c). The role of sequential patterning by the proneural wave is to ensure the correct timing and composition of the neuroblast population (Bertet et al., 2017). A similar process of sequential patterning occurs during the progression of the morphogenetic furrow in the Drosophila eye (Roignant and Treisman, 2009; Lubensky et al., 2011; Formosa-Jordan et al., 2012; Wartlick et al., 2014; Fried et al., 2016; Gavish et al., 2016). However, the progression of the morphogenetic furrow also entails transient growth as well as subsequent photoreceptor patterning and differentiation to generate ommatidia (Sato et al., 2013).

In our model, which focuses on the driving mechanism behind the proneural wave, we have refrained from considering additional signalling pathways that neither play key roles in driving the proneural wave nor exhibit strong signatures of bidirectional feedbacks (in contrast to EGFR and Delta-Notch signalling). These include the JAK/STAT and Hippo pathways, which serve important roles in modulating proneural wave progression but are not actively involved in propagating the transition zone through a reaction-diffusion-like mechanism (Yasugi et al., 2008; Reddy et al., 2010; Yasugi et al., 2010; Kawamori et al., 2011; Wang et al., 2011a; Weng et al., 2012).

On a mechanistic level, the excitable propagation behaviour illuminated here provides a mechanism to capture the transient and localised activity of the proneural gene l’sc and EGFR signalling, as well as robustness against fluctuating signalling activity and gene expression. In contrast to a recent model (Sato et al., 2016), our model also implies that neither differentiation nor proneural gene expression is required for the transient stabilisation of EGFR signalling activity that is required to advance the wave front. In the context of vertebrate somitogenesis, intracellular excitability has recently been suggested to underlie the emergence of genetic oscillations (Hubaud et al., 2017). The appearance here of excitability in the context of a propagating front of gene expression suggests that such a mechanism may serve more widely as a generic and robust strategy to achieve sequential transition waves in developing tissues.

Materials and methods

Fly strains

Request a detailed protocol

Flies were raised on standard cornmeal medium at 25°C. Strains used were:

yw, hsFLP; FRT40A, tub-GAL80/CyO, ActGFP; tubP-GAL4,UAS-mCD8-GFP/TM6B

w; FRT40A; UAS-Pnt-P1/TM6B

w; FRT40A/CyO; N RNAi/TM6B (N RNAi lines used were BL33611 and BL33616)

HLH-mgamma-GFP (Almeida and Bray, 2005) was used to report active Notch signalling.

UAS-Pnt-P1 clones and N RNAi clones

Request a detailed protocol

To generate clones in the developing optic lobe, larvae were collected 48–50 hr after egg laying and were heat shocked for 20 min at 37°C. Larvae were then dissected 50 or 60 hr after clone induction.

Immunohistochemistry

Request a detailed protocol

Larval brains were dissected in PBS and fixed for 20 min at room temperature in 4% formaldehyde and fixation buffer (PBS, 5 mM MgCl2, 0.5 mM EGTA). After fixation, brains were rinsed and washed in 0.3% PBS Triton X100 (PBT). Samples were blocked with 10% normal goat serum (NGS) in 0.3% PBT at room temperature and incubated with the primary antibodies overnight at 4°C. Brains were then washed in 0.3% PBT and incubated with the secondary antibodies overnight at 4°C. Brains were washed in 0.3% PBT and mounted in Vectashield (Vector Laboratories, Burlingame, CA, USA). The following primary antibodies and dilutions were used: guinea pig anti-Dpn (1:10,000) and rat anti-L’sc (1:5,000) (Caygill and Brand, 2017), chicken anti-GFP (1:2,000) from Abcam, mouse anti-Delta (1:50, C594.9B) from DSHB and mouse anti-Notch (intracellular domain, ICD) (1:50, C17.9C6) from DSHB. Fluorescently conjugated secondary antibodies Alexa405, Alexa488, Alexa546 and Alexa633 (all 1:200) from Life Technologies.

Images were acquired with a Leica TCS SP8 confocal microscope (Leica Microsystems, Wetzlar, Germany) and analysed with Fiji (Schindelin et al., 2012). Figures and illustrations were assembled using Adobe Photoshop CS3 and Adobe Illustrator CS3 (Adobe Systems, San Jose, CA, USA).

Appendix 1

Bistable EGFR signalling fronts

We describe the details of our theoretical model by building up an integrated model of the proneural wave starting from a simple system of EGFR signalling activity. Reaction-diffusion models with different types of genetic interactions and correspondingly different phenomenologies have been used to describe the dynamics of EGFR signalling in other contexts such as Drosophila oogenesis (Shvartsman et al., 2002; Zartman et al., 2009; Simakov et al., 2012; Fauré et al., 2014) and the Drosophila retina (Pennington and Lubensky, 2010).

Model formulation

To illustrate the mechanism of wave propagation, we first show how EGFR signalling activity alone may give rise to a propagating signalling front. As outlined in the main text, we consider the dynamics of one signalling component E, which encapsulates the collective dynamics of the EGFR/Rhomboid/Spitz network (see below) and, hence, involves positive feedback and effective diffusion on the tissue level (see Figure 1c and Figure 2a). We describe the component E by a continuous field that defines the signalling activity, ϕ=ϕ(𝐱,t), where 𝐱 denotes the position coordinate within the tissue and t the time. (Throughout this text, we consider one-dimensional systems for illustration purposes and two-dimensional systems to describe the proneural wave in vivo. Hence, the position variable 𝐱 and the nabla operator refer to 𝐱=(x1,x2) and =(/x1,/x2) in the case of two dimensions and 𝐱=x and =/x in the case of one dimension.) The corresponding reaction-diffusion system is given by

(3) ϕt=η2ϕ+r(ϕ) ,

where η is the effective diffusion constant and r(ϕ) is the reaction term describing the intracellular feedback,

(4) r(ϕ)=μh(ϕΦ)-kϕ.

Here, μ denotes the gain rate in signalling activity due to positive feedback, k the degradation rate and h is a monotonically increasing function describing the nonlinear positive feedback with Φ being the threshold activity. Here, we choose a function of the Hill type (Novák and Tyson, 2008),

(5) h(ϕ)=ϕn1+ϕn,

where n is the ‘Hill exponent’ characterising the nonlinearity of the feedback. Systems of the type given by Equations 3–5 are known to exhibit solutions in the form of a propagating front for appropriate parameters (Keener and Sneyd, 2009). For illustrative purposes, we initially consider the dynamics of signalling activity in one spatial dimension on a domain of length with ‘no-flux’ boundary conditions, (ϕ/x)|x=0=0=(ϕ/x)|x=. A numerical example of the system specified by Equations 3–5 is shown in Figure 2a. A localised concentration beyond a certain threshold concentration at the boundary x=0 initiates a travelling front of E, which propagates at constant velocity, leaving behind elevated levels of sustained signalling activity. (For the example shown in Figure 2a, we choose an initial condition of the form ϕ(x,0)=e-x2/x02 with x02=5η/k.)

This behaviour is familiar in reaction-diffusion systems with reaction terms of the form given by Equation 4 and can be understood by studying the local reaction dynamics in a phase portrait: Appendix 1—figure 1a shows the reaction term r, which describes the local net loss/gain rate of signalling activity E. For sufficiently large gain rates, the function has three equilibrium points ϕi* for which the net loss/gain rate vanishes (dots in Appendix 1—figure 1a). Therefore, once the system has reached such an equilibrium, it remains there until perturbed. The two equilibrium points ϕ0* and ϕ2* are stable (black dots), that is small perturbations of the signalling activity will drive the system back to these equilibrium points—higher levels lead to a net loss whereas smaller levels lead to a net gain (black arrows in Appendix 1—figure 1a). The two equilibria are separated by an unstable equilibrium ϕ1* (white dot). If the system is only slightly perturbed around ϕ=ϕ0* (black dot), the reaction dynamics will drive the system back towards this equilibrium, because degradation is stronger than self-activation in this region, r<0. However, if the signalling activity rises above ϕ1* (e.g., due to diffusion from the neighboring cell), self-activation will outcompete degradation, r>0, and the signalling activity will be driven towards the other stable equilibrium ϕ2*. This reaction-diffusion mechanism is known as a bistable front (Keener and Sneyd, 2009).

Appendix 1—figure 1
Key features of reaction dynamics leading to bistable front propagation.

(a) Reaction term ρ as given by Equation 6. Dots indicate the fixed points ϕi* with i=0,1,2 for which ρ(ϕi*)=0. Filled dots indicate stable fixed points, the open dot indicates the unstable fixed point. (b) Potential U associated with the reaction term ρ shown in panel A and defined by ρ=-U/ϕ. Parameters in both panels are λ=4 and n=3. (c) Reaction term ρ as given by Equation 6 for Hill exponents n=2 (light green), n=4 (green), n=10 (dark green), and the limiting case n (dashed black), given by Equation 9, for λ=4. (d) Front velocity c as a function of λ. Numerical solutions obtained from simulations of Equation 6 with Hill exponent n=2 and analytical approximation Equation 11 for the n limit. (e) Example of the front profile ϕ given by Equation 10 with λ=4.

https://doi.org/10.7554/eLife.40919.016

Further insight into the dynamics of the front propagation can be gained using analytical techniques (Keener and Sneyd, 2009). Non-dimensionalising the model by rescaling ϕϕ/Φ, tkt, and xk/ηx, the corresponding reaction-diffusion system in one dimension is given by

(6) ϕt=2ϕx2+ρ(ϕ),ρ(ϕ)=λh(ϕ)-ϕ.

Here, λ=μ/(kΦ) is the only (dimensionless) parameter and ρ is the dimensionless reaction term.

Anticipating that the system gives rise to a travelling front with velocity c and a stationary front profile ϕ, we make the ansatz ϕ(x,t)=ϕ(x-ct) in Equation 6, which yields the ordinary differential equation

(7) ϕ′′+cϕ+ρ(ϕ)=0,

where the prime denotes the derivative with respect to the argument y=x-ct of ϕ.

Multiplying this equation by ϕ and integrating from - to yields an implicit expression for the velocity of the front,

(8) c=ΔU(ϕ)2dy,

where ΔU=limy[U(ϕ(y))-U(ϕ(-y))] with the potential U(ϕ) defined by ρ=-U/ϕ (see Appendix 1—figure 1b). Analytical solutions to Equation 8 only exist for special functional forms of the reaction term ρ. To obtain analytical insights into the wave speed, we consider the limiting case of large Hill exponents, n in Equation 5. This corresponds to a switch-like response of the gain rate once the activation threshold is reached. In this case, the reaction term acquires a piecewise linear functional form,

(9) ρ(ϕ)=λΘ(ϕ-1)-ϕ,

where Θ is the Heaviside step function with the convention Θ(0)=1/2 (see Appendix 1—figure 1c). This function has three equilibria for λ>1. For this reaction term, the front profile and front velocity can be calculated analytically via Equation 7 and Equation 8. The front profile is given by

(10) ϕ(y)={e-λ-1yy0λ+(1-λ)ey/λ-1y<0,

where we have fixed the arbitrary position of the co-moving reference frame by imposing the condition ϕ(0)=1. Appendix 1—figure 1e shows an example of the front solution Equation 10. The velocity of the front is given by

(11) c=λ-2λ-1.

Expanding in the limit of large λ, dropping all orders higher than λ-1/2 and restoring the original parameter dependence by multiplying with the velocity scale ηk, we find

(12) cημΦ-32ηΦμk.

Equation 12 implies that the front velocity increases with increasing diffusion constant η and synthesis rate μ as well as decreasing degradation rate k and decreasing activation threshold ΦAppendix 1—figure 1d shows that Equation 11 is a good approximation for the velocity of the front even when compared to numerical simulations with a finite Hill exponent n.

A one-component model reproduces the key features of more detailed models

Previously, we sought to capture the spatio-temporal dynamics of EGFR signalling by a single component, described by the signalling activity ϕ, even though EGFR signalling comprises three components: the EGF receptor, the transmembrane factor Rhomboid and the EGFR ligand Spitz (see Figure 1d). To show that the system above entails the essential features of more detailed descriptions, we now consider the kinetics of these three components and show that a corresponding model gives rise to the same qualitative dynamics. We consider a dimensionally reduced reaction-diffusion system where ψE denotes the concentration of bound EGFRs, ψR represents the concentration of Rhomboid, and ψS is the concentration of the secreted active form of Spitz (sSpi),

(13) ψEt=λEHE(ψS)ψE ,ψRt=λRHR(ψE)ψR ,ψSt=2ψSx2+λSψRψS .

Here, λE is the binding rate of EGFRs, λR is the synthesis rate of Rhomboid, and λS is the secretion rate of Spitz; the functions HE and HR describe saturation of EGFR binding and saturation of the Rhomboid synthesis rate, respectively, and are assumed to be of the same qualitative form as the function h, given by Equation 5. Only the equation for the component S contains a diffusion term as only Spitz is secreted. The negative terms account for receptor unbinding and degradation of gene products; for simplicity, we have chosen identical unit decay rates and concentration thresholds. Numerical examples of Equation 13 show that the positive feedback of E, R, and S can generate a bistable front, see Appendix 1—figure 2.

Appendix 1—figure 2
Numerical example of the three-component system Equation 13.

Different curves show ψE (black), ψR (green), and ψS (blue). The two panels show the time points t=30 (top) and t=60 (bottom). Functions and parameters are hR(ϕ)=hE(ϕ)=h(ϕ), given by Equation 5 with n=3, and λi=4 for i=E,R,S. Boundary conditions are (ψi/x)|x=0=0=(ψi/x)|x=.

https://doi.org/10.7554/eLife.40919.017

Since each component is required to activate the next one in the loop, all three components show the same behaviour in terms of their localisation. Moreover, since the activation of one component is sufficient to activate the entire loop, diffusibility of S confers a diffusion-like effect to the entire E–R–S system. Therefore, the parameter μ of the one-component model Equation 3 and Equation 4 is to be interpreted as a composite parameter that characterises the strength of the positive feedback of the entire cycle. Likewise, the parameter k represents the combined effects of degradation and receptor unbinding, and η describes the effective diffusibility, mediated by diffusion of component S. Formally, a connection to the one-component model Equation 6 can be established when the dynamics of E and R is sufficiently fast such that they are always close to their stationary values determined by ψE/t=0 and ψR/t=0 while the front is travelling. In this case, we can solve for ψE and ψR and obtain a closed equation for ϕψS,

(14) ϕt=2ϕx2+λh(ϕ)-ϕ,

where λλS and h(ϕ)=λRHR(λEHE(ϕ)). Note that Equation 14 is formally equivalent to Equation 6; if HE and HR are sigmoidal functions with the qualitive shape of Equation 5, then so is h.

https://doi.org/10.7554/eLife.40919.015

Appendix 2

Excitable dynamics from EGFR–proneural interactions

While the minimal model Equations 3–5 leads to the emergence of a front that leaves behind an elevated signalling state, the transition zone of the proneural wave is characterised by localised EGFR signalling and proneural gene expression. We now show how such a travelling localised pulse arises if we take into account interactions between EGFR signalling and the proneural gene l’sc. The interactions between EGFR signalling and L’sc are schematically represented in Figure 1f and Figure 2b: the component E activates the expression of L while the component L effectively downregulates E as a consequence of the transition that is induced. The corresponding reaction-diffusion equations for the two fields ϕE and ϕL are given by

(15) ϕEt=2ϕE+ρE(ϕE,ϕL) ,ϕLt=ρL(ϕE,ϕL) ,

with the reaction terms

(16) ρE=λEh(ϕE)h¯(ϕL)ϕE ,ρL=λLh(ϕE)ϕL ,

where we have introduced the Hill function h¯(ϕ)1-h(ϕ) which describes an inhibitory effect. For simplicity, we have also considered identical unit concentration thresholds for activation and inhibition. Here, λE and λL are the gain rates of the components E and L, respectively. Note that in the absence of L (ϕL=0), the reaction term ρE reduces to the one of the one-component model given in Equation 6, ρE(ϕE,0)=ρ(ϕE) (Appendix 2—figure 1). Again, we consider a finite one-dimensional domain of length and no-flux boundary conditions, (ϕE/x)|x=0=0=(ϕE/x)|x= and (ϕL/x)|x=0=0=(ϕL/x)|x=.

Appendix 2—figure 1
Reaction dynamics of the two-component model Equation 15 and Equation 16.

(a) Reaction term of the component E in the absence of L, given by ρE(ϕE,0) in Equation 16. Dots indicate fixed points for which ρE(ϕE,0)=0. (b) Full local reaction dynamics for the two-component model. Vector field F=(ρE,ρL) as given by Equation 16. Dots indicate points with ϕL=0 and ρE=0. Parameters in both plots are λE=λL=4 and n=3. Coloured curves show the nullclines for E (green) and L (blue).

https://doi.org/10.7554/eLife.40919.019

Figure 2b in the main text displays a numerical example of the system Equation 15 and Equation 16. Again, starting from a localised concentration of the component E at x=0, a localised pulse of E and L travels through the system at a constant speed. (For the example shown in Figure 2b, we choose an initial condition of the form ϕ(x,0)=e-x2/x02 with x02=10η/k.) Again, this behaviour can be elucidated by studying the local reaction dynamics, which is now given by the vector field 𝐅=(ρE,ρL), defined by Equation 16) (Appendix 2—figure 1). Note that now only the point (ϕE,ϕL)=(0,0) is a stable equilibrium. The black and white dots indicate points with ϕL=0 that satisfy ρE=0 as in the one-component system (compare to Appendix 1—figure 1a). However, note that only the point (0,0) also satisfies ρL=0. Again, diffusion from a neighboring cell will increase the levels of E. When the concentration level marked by the white point is exceeded, the reaction dynamics will elevate the levels of E and L until a turning point is reached when downregulation of E is sufficiently strong to suppress its positive self-feedback. Thus, the system finally returns to the fixed point (ϕE,ϕL)=(0,0) which marks the end of the pulse. Hence, the model Equation 15 and Equation 16 leads to propagation of a pulse of E and L at a defined speed.

https://doi.org/10.7554/eLife.40919.018

Appendix 3

Integrated model of the proneural wave

The system described in Appendix 2 invokes a core mechanism giving rise to a propagating transition zone. As motivated in the main text, we now extend our model to include Delta–Notch interactions, based on the classical description advanced by (Collier et al., 1996) and also include cis-inhibition (Sprinzak et al., 2010; Shaya and Sprinzak, 2011).

An earlier attempt to describe the proneural wave advanced by (Sato et al., 2016) has focused on a phenomenological description of the proneural wave, which, already in its general approach, differs from our model. The corresponding model is, in major parts, built around observed patterns of gene expression and patterning phenomena, rather than starting from the molecular underpinnings: there, proneural gene expression is not an independent dynamic ingredient; rather, the cell state (NE or NB) and the state of proneural gene expression is combined in a single abstract variable, making it impossible to address the effects of alterations in proneural gene expression independently of the NE to NB transition. Hence, in (Sato et al., 2016), the general driving mechanism of the wave is fundamentally different from our model: it relies on EGFR signalling inducing a change in cell state and, conversely, a change in cell state transiently driving EGFR signalling activity through a phenomenological prescription. In contrast to our model, this renders autocrine EGFR signalling without intrinsic bistability and therefore without the ability to autonomously drive the wavefront. Moreover, the resulting model is unstable with respect to additive fluctuations in gene expression and signalling activity; slight misexpression or perturbations in signalling activity will result in an immediate premature differentiation, as discussed in (Sato et al., 2016). Delta-Notch interactions are incorporated as a subsystem without intrinsic multistability usually found necessary in attempts to describe the emergence of lateral inhibition phenomena (Collier et al., 1996).

Since Delta–Notch interactions can give rise to lateral inhibition, that is stable low-Delta/high-Notch and high-Delta/low-Notch states in adjacent cells, instead of a continuous description of the tissue, we now consider a lattice where each lattice site represents a cell. We consider the signalling and gene activities ϕ𝐱i(t) with i=E,L,D,N describing EGFR signalling activity, L’sc expression, and Delta and Notch, respectively. The index x indicates the lattice site. Furthermore, we introduce the cell state Ω𝐱(t) which takes values from 0 to 1, where Ω𝐱=0 indicates that cell x is a neuroepithelial cell and Ω𝐱=1 indicates that cell x is a neuroblast. Figure 2d shows the regulatory network of the model. In contrast to the effective inhibition of signalling by the proneural gene l’sc that we considered before, we now include the more realistic shutdown of signalling as a consequence of differentiation. Moreover, motivated by the presence of low levels of Notch signalling in the neuroepithelium and the neuroblasts (Egger et al., 2010; Orihara-Ono et al., 2011), we include a basal source of Notch that is independent of trans-activation by Delta in adjacent cells. As in the previous sections, we here consider identical reference concentrations for activation and inhibition and rescale all concentrations by this reference concentration, so that the fields ϕ𝐱i with i=E,L,D,N are dimensionless. Moreover, we consider identical degradation constants for all four components and rescale time by the degradation time. The corresponding dynamical equations are given by

(17) dϕxEdt=η[Δ^ϕE]x+μE(h(ϕxE)+h(ϕxN/Φ1))h¯(2Ωx)kEϕxE,dϕxLdt=μLh(ϕxE)h¯(ϕxN/Φ2)kLϕxL,dϕxDdt=μD[h(ϕxE)+h¯(ϕxN)]h¯(2Ωx)kDϕxD,dϕxNdt=[β+μNh([Σ^ϕD]x)]h¯(ϕxD)h¯(ϕxL)kNϕxN.

Here, the parameters μi denote gain rates, ki denote decay rates, β denotes the basal gain rate of the component N, and the Φi denote threshold levels for positive and negative feedbacks. The operators Δ^ and Σ^ are the lattice Laplacian and the sum over concentrations of neighboring lattice sites, respectively, and defined by [Δ^ϕ]𝐱=yU𝐱(ϕ𝐲-ϕ𝐱) and [Σ^ϕ]𝐱=yU𝐱ϕ𝐲, where U𝐱 is the set of neighbours of site x. The dynamics of the cell state Ω𝐱 is given by

(18) dΩ𝐱dt=-dVdΩ(Ω𝐱)+fint(ϕ𝐱E,ϕ𝐱L,ϕ𝐱D,ϕ𝐱N),

where the function V is a ‘potential’ for the cell state which ensures that in the absence of signalling and proneural gene expression, Ω𝐱 has two stable equilibria: Ω𝐱=0 (neuroepithelium) and Ω𝐱=1 (neuroblast) (see Appendix 3—figure 1). The qualitative features of our model do not depend on the exact choice of V. The term fint acts as a ‘force’ that triggers the transition from one state to the other depending on the local signalling activity and proneural gene expression. The functional form of fint is based on the observations that (i) proneural gene expression seems to be sufficient but not necessary for the transition while (ii) Notch expression seems to keep cells in their proneural state (Yasugi et al., 2010). Therefore, we here choose

(19) V(Ω)=Ω2(1Ω)2/4,fint(ϕE,ϕL,ϕD,ϕN)=h(ϕL/Φint)+h¯(ϕN/Φint)h(ϕE/Φint).

The cell state potential V has wells at Ω=0 (neuroepithelial state) and Ω=1 (neuroblast state) (see Appendix 3—figure 1), so that both states are stable and a transition occurs when the hill in between both states can be overcome. This choice of fint leads to a transition of the cell state from neuroepithelial (Ω=0) to neuroblast (Ω=1) when (i) L exceeds the threshold level Φint and/or (ii) E exceeds and N drops below the threshold level.

Appendix 3—figure 1
The cell state potential V, given by Equation 19 has two wells, corresponding to the neuroepithelial state (NE, Ω=0) and the neuroblast state (NB, Ω=1).
https://doi.org/10.7554/eLife.40919.021
Appendix 3—table 1
Reference parameter set used for the model Equations 17–19.
https://doi.org/10.7554/eLife.40919.022
Parameter(s)ValueDescriptionAffected components
η0.02diffusion constantE
μE, μL, β10gain ratesE, L, N
μD, μN5gain ratesD, N
kE, kL, kD, kN1decay ratesE, L, D, N
Φ1100Notch thresholdE
Φ20.5Notch thresholdL
Φint10threshold for differentiationΩ
n3Hill exponentE, L, D, N
γ(as indicated)biochemical noise strengthE, L, D, N
https://doi.org/10.7554/eLife.40919.020

Appendix 4

Robustness of the model against fluctuations and disorder

Robustness against biochemical noise

Gene expression and biochemical reactions typically suffer from fluctuations due to small numbers of molecules involved (Tsimring, 2014). To achieve reliable morphological results, any biochemical mechanism governing morphogenetic processes must be robust against such types of noise. From an analytical point of view, stability of the zero signalling fixed point in our model (see Appendix 2—figure 1) ensures that a cell does not differentiate prematurely due to a certain degree of noise in proneural gene expression or fluctuating signalling activity. To numerically demonstrate that our model is robust against fluctuations in molecule concentrations, we performed simulations of the system in the presence of biochemical fluctuations in all four components. The noisy system is given by Equations 17–19 with each dynamical equation replaced according to

(20) dϕidtdϕidt+γξi(t),(i=E,L,D,N)

where γ denotes the noise strength and ξi denotes Gaussian white noise characterised by the expectation values ξi(t)=0 and ξi(t)ξj(t)=δijδ(t-t). Furthermore, δij denotes the Kronecker delta and δ(t) the Dirac delta distribution.

Appendix 4—figure 1 shows numerical examples of the system for different noise strengths γ. In these examples, the wave robustly travels from the left to the right for small and intermediate noise levels (as compared to the gain rate μE) (Appendix 4—figure 1a,b), whereas premature differentiation is only observed for large noise levels that introduce fluctuations comparable to physiological concentrations (Appendix 4—figure 1c). Finally, random differentiation throughout the tissue only occurs if the system is dominated by fluctuations (Appendix 4—figure 1d).

Appendix 4—figure 1
Response of the dynamics to biochemical fluctuations.

Panels show snapshots of the system for different noise strengths: (a) γ/μE=0, (b) γ/μE=0.75, (c) γ/μE=1, (d) γ/μE=1.5. All other parameters are given in Appendix 3—table 1. The system given by Equations 17–20 was simulated on a hexagonal lattice with circular geometry with a radius of 15 lattice sites. Initial conditions were localised elevated levels of E in those outer boundary cells that have angles between π/3 and 5π/3 as measured from the center of the circular lattice. The respective simulation panels show snapshots of the activity of EGFR signalling (green), L’sc expression (blue), Delta activity (magenta) and Notch activity (pink), as well as the cell state Ω, for which black indicates neuroepithelium (Ω=0) and red indicates neuroblasts (Ω=1). Colour intensity indicates the local gene expression levels or signalling activities, respectively. The snapshots show the time t=12.5.

https://doi.org/10.7554/eLife.40919.024

Robustness against lattice defects

To test whether the proposed mechanism is robust with respect to disordered lattice structures, we considered a site-diluted version of our model, in which randomly chosen sites in the hexagonal lattice were ‘deactivated’, that is removed from the dynamics and the coupling topology. Simulations showed that while the presence of disorder leads to a local distortion of the propagating front profile, the overall mechanism remains intact and leads to robust and sequential differentiation of the neuroepithelial tissue (Appendix 4—figure 2a). Even regions which are partially ‘shielded’ by defects eventually become differentiated as the diffusion-mediated propagation has no intrinsic directionality and is able to reach such regions when the wave has surrounded the corresponding region (Appendix 4—figure 2). Moreover, the overall phenomenology of mutant and transgenic clones remains intact on such imperfect lattices, as illustrated using the ‘Notch upregulation’ clone (cf. Figure 3f and Appendix 4—figure 2b).

Appendix 4—figure 2
Simulated proneural wave on a lattice with random ‘defects’ (black sites).

(a) Proneural wave propagation on a hexagonal lattice with circular geometry. Parameters and initial conditions as in Appendix 4—figure 1 with γ=0. The snapshots show the time t=10 (top row), t=25 (middle row) and t=40 (bottom row). (b)Notch upregulation clone as shown in Figure 3f, but on a heavily site-diluted lattice. The white arrowhead indicates retarded differentation, cf. Figure 3f.

https://doi.org/10.7554/eLife.40919.025
https://doi.org/10.7554/eLife.40919.023

Appendix 5

Suppression of lateral inhibition patterns

An analytical argument demonstrating how basal Notch levels suppress lateral inhibition can be made in a simple picture involving Delta–Notch interactions in two cells x=1,2,

(21) dϕxDdt=h¯(ϕxN)ϕxD ,dϕxNdt=β+λϕx¯DϕxN ,

where x¯ refers to the respective other cell and h¯(ϕ)1-h(ϕ) with h being the Hill function Equation 5, as before. Here, λ is the gain rate for Notch and β indicates the basal production. For simplicity, we consider a linear positive feedback in the dynamics of Notch and the Hill exponent n=2 for h¯. In steady state, where dϕ𝐱D/dt=0=dϕ𝐱N/dt, we can eliminate ϕ1N and ϕ2N to obtain

(22) ϕ𝐱D=Γ(ϕ𝐱¯D),

where Γ(ϕ)=h¯(β+λϕ). From this, it follows that both ϕ1D and ϕ2D satisfy ϕ𝐱D=Γ(Γ(ϕ𝐱D)). Among other solutions, this equation has two solutions of the form ϕ±=p±q with

(23) p=1211+β2-βλ,q=p2-1+β2λ2,

which, in the case that they are real, correspond to the low-Delta/high-Notch and high-Delta/low-Notch states in adjacent cells as they satisfy ϕ±=Γ(ϕ).

However, bistability only exists if both solutions are real, which is the case for q>0. From Equation 23, we find that this corresponds to

(24) λ>2(β2+1)(β+β2+1).

Figure 4e shows the corresponding phase diagram for the occurrence of lateral inhibition in the two-cell system. Therefore, a basal production term, if strong enough, can prevent lateral inhibition patterns.

https://doi.org/10.7554/eLife.40919.026

Appendix 6

Description of mutant and transgenic clones

Simulation of clones

To capture experimental scenarios in which mutant or transgenic clones were induced in the neuroepithelium, we modified the model Equation 17 accordingly. In our model, a clone is defined by a cell or a group of cells within the simulated tissue which has altered kinetic rate parameters or a different initial condition, depending on the type of experimental perturbation (Figure 3). For the case of downregulation of proneural gene or signalling factors, the synthesis or binding rate of the respective gene or receptor in the clone cells is decreased or set to zero, as indicated in the caption of Figure 3. For the case of upregulation, which in all considered cases corresponds to a constitutively active gene, we added a source term to the clone that leads to constant synthesis and furthermore set the initial condition of the clone to the elevated steady-state concentration of the respective gene or signalling activity.

Simulation of Figure 7a

To simulate the effects of clones in which EGFR signalling is constitutively activated within the neuroepithelium (Figure 7), we simulated the model on a hexagonal lattice with circular boundaries (Video 4). The initial condition was set to ϕE|t=0=μE/2 in a one-dimensional array of cells in the outermost cell layer that have angles between π/3 and 5π/3 as measured from the center of the circular lattice. Moreover, we arbitrarily selected four lattice sites and endowed them with a constant production of E by adding the term μEh¯(2Ω𝐱) to the reaction dynamics of the component E in Equation 17; this mimics the constitutively active EGFR signalling. Figure 7a shows a snapshot at time t=14. All parameter values are given in Appendix 3—table 1.

https://doi.org/10.7554/eLife.40919.027

Appendix 7

Sensitivity analysis of the model

Morris method for global sensitivity analyses

To test the sensitivity of key observables on model parameters, we here employ the so-called Morris method, a widely used method for global sensitivity analyses (Morris, 1991; Campolongo et al., 2007; Wu et al., 2013). To briefly summarise, for a fixed model observable 𝒪 and a given set of parameters θ1,,θn, the Morris method consists in repeatedly sampling a discretised parameter space (or subspace of interest) of the model and, for each parameter i, calculating the so-called ‘elementary effects’

(25) ei=𝒪(θ1,,θi+Δ,,θn)-𝒪(θ1,,θn)Δ,

That is the finite-difference quotient of the output with respect to the parameter, given a finite step size Δ that is chosen adequately; for standard choices and further details on the method, see, for example (Wu et al., 2013). This sampling procedure yields a distribution Pi(e) of elementary effects for each parameter i, from which the following sensitivity indices are computed,

(26) mi=ei,mi*=|e|i,σ=e2i-ei2,

where i denotes the expectation value under the distribution Pi. The interpretation of these indices is given in the main text.

Probed observables and parameters

As output observables we here choose the linear propagation velocity v of the proneural wave and the width w of the transition zone. We formally define these quantities as follows for the lattice-based full proneural wave model Equations 17–19. By x, we denote the extension of the system in the direction of the travelling wave. We define Ω¯x=-1yΩxy as the average of the cell state Ω in the direction perpendicular to the wave, where is the extension of the lattice in the perpendicular direction. A wave with constant velocity leads to a proportionally linear increase in number of neuroblasts, so that the wave velocity v (in lattice sites per unit time) is given by

(27) v=ddtxΩ¯x.

Practically, we determine v as the slope obtained from a linear fit of xΩ¯x in the linear regime.

We define the width w of the transition zone via the spatial spread of transitioning cells, that is those with 0<Ω<1. Formally, this width can be defined as

(28) w=[2πx(Ω¯x-Ω¯x-1)2]-1.

The discrete derivative, Ω¯x-Ω¯x-1, which measures the steepness of the profile, is non-zero only in the transition zone. For example, for a Gaussian profile of the discrete derivative, Ω¯x-Ω¯x-1e-x/2σ2 with a variance σ21, Equation 28 yields w=σ. To avoid confounding effects by initial and boundary conditions in model simulations, we use the temporal median of w as a proxy for the width of the transition zone.

We compute the Morris indices Equation 26 for the dependence of v and w on the kinetic and diffusion properties of the integrated model, that is on the parameters η, β and μi, ki for i=E,L,D,N and allow them to vary between the 0.2-fold and 5-fold reference value given in Appendix 3—table 1, while keeping all other parameters fixed to their values given in Appendix 3—table 1. 55000 parameter samples were used to compute expectation values.

https://doi.org/10.7554/eLife.40919.028

Data availability

All data generated or analysed during this study are included in the manuscript.

References

  1. Book
    1. Bertet C
    2. Çelik A
    3. Wernet MF
    (2017)
    The Developmental Origin of Cell Type Diversity in the Drosophila Visual System
    Cham, Switzerland: Springer International Publishing.
    1. Shvartsman SY
    2. Muratov CB
    3. Lauffenburger DA
    (2002)
    Modeling and computational analysis of EGF receptor-mediated cell communication in Drosophila oogenesis
    Development 129:2577–2589.
    1. Turing AM
    (1952) The chemical basis of morphogenesis
    Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 237:37–72.
    https://doi.org/10.1098/rstb.1952.0012

Article and author information

Author details

  1. David J Jörg

    1. Cavendish Laboratory, Department of Physics, University of Cambridge, Cambridge, United Kingdom
    2. The Wellcome Trust/Cancer Research UK Gurdon Institute, University of Cambridge, Cambridge, United Kingdom
    Contribution
    Conceptualization, Software, Formal analysis, Investigation, Visualization, Methodology, Writing—original draft, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-5960-0260
  2. Elizabeth E Caygill

    1. The Wellcome Trust/Cancer Research UK Gurdon Institute, University of Cambridge, Cambridge, United Kingdom
    2. Department of Physiology, Development and Neuroscience, University of Cambridge, Cambridge, United Kingdom
    Contribution
    Investigation, Visualization, Methodology, Writing—original draft, Writing—review and editing
    Competing interests
    No competing interests declared
  3. Anna E Hakes

    1. The Wellcome Trust/Cancer Research UK Gurdon Institute, University of Cambridge, Cambridge, United Kingdom
    2. Department of Physiology, Development and Neuroscience, University of Cambridge, Cambridge, United Kingdom
    Contribution
    Investigation, Visualization, Methodology, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8664-1014
  4. Esteban G Contreras

    The Wellcome Trust/Cancer Research UK Gurdon Institute, University of Cambridge, Cambridge, United Kingdom
    Contribution
    Investigation, Visualization, Methodology, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-0934-741X
  5. Andrea H Brand

    The Wellcome Trust/Cancer Research UK Gurdon Institute, University of Cambridge, Cambridge, United Kingdom
    Contribution
    Conceptualization, Supervision, Funding acquisition, Investigation, Methodology, Writing—original draft, Project administration, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2089-6954
  6. Benjamin D Simons

    1. Cavendish Laboratory, Department of Physics, University of Cambridge, Cambridge, United Kingdom
    2. The Wellcome Trust/Cancer Research UK Gurdon Institute, University of Cambridge, Cambridge, United Kingdom
    3. The Wellcome Trust/Medical Research Council Stem Cell Institute, University of Cambridge, Cambridge, United Kingdom
    Contribution
    Conceptualization, Supervision, Funding acquisition, Investigation, Methodology, Writing—original draft, Project administration, Writing—review and editing
    For correspondence
    bds10@cam.ac.uk
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-3875-7071

Funding

Wellcome (092096)

  • David J Jörg
  • Elizabeth E Caygill
  • Anna E Hakes

Cancer Research UK (C6946/A14492)

  • David J Jörg
  • Elizabeth E Caygill
  • Anna E Hakes
  • Esteban G Contreras
  • Andrea H Brand
  • Benjamin D Simons

Biotechnology and Biological Sciences Research Council (Project Grant BB/L007800/1)

  • Elizabeth E Caygill
  • Andrea H Brand

Wellcome (Senior Investigator Award 103792)

  • Anna E Hakes
  • Esteban G Contreras
  • Andrea H Brand

Royal Society (Darwin Trust Research Professorship)

  • Andrea H Brand

Wellcome (Senior Investigator Award 098357)

  • Benjamin D Simons

Royal Society (EP Abraham Research Professorship RP\R1\180165)

  • Benjamin D Simons

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We are grateful to Pau Formosa-Jordan, Claude Desplan, François Schweisguth and members of the Simons and Brand labs for useful discussions.

Version history

  1. Received: August 9, 2018
  2. Accepted: February 8, 2019
  3. Version of Record published: February 22, 2019 (version 1)

Copyright

© 2019, Jörg et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,986
    views
  • 299
    downloads
  • 13
    citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. David J Jörg
  2. Elizabeth E Caygill
  3. Anna E Hakes
  4. Esteban G Contreras
  5. Andrea H Brand
  6. Benjamin D Simons
(2019)
The proneural wave in the Drosophila optic lobe is driven by an excitable reaction-diffusion mechanism
eLife 8:e40919.
https://doi.org/10.7554/eLife.40919

Share this article

https://doi.org/10.7554/eLife.40919

Further reading

    1. Computational and Systems Biology
    2. Developmental Biology
    Gang Xue, Xiaoyi Zhang ... Zhiyuan Li
    Research Article

    Organisms utilize gene regulatory networks (GRN) to make fate decisions, but the regulatory mechanisms of transcription factors (TF) in GRNs are exceedingly intricate. A longstanding question in this field is how these tangled interactions synergistically contribute to decision-making procedures. To comprehensively understand the role of regulatory logic in cell fate decisions, we constructed a logic-incorporated GRN model and examined its behavior under two distinct driving forces (noise-driven and signal-driven). Under the noise-driven mode, we distilled the relationship among fate bias, regulatory logic, and noise profile. Under the signal-driven mode, we bridged regulatory logic and progression-accuracy trade-off, and uncovered distinctive trajectories of reprogramming influenced by logic motifs. In differentiation, we characterized a special logic-dependent priming stage by the solution landscape. Finally, we applied our findings to decipher three biological instances: hematopoiesis, embryogenesis, and trans-differentiation. Orthogonal to the classical analysis of expression profile, we harnessed noise patterns to construct the GRN corresponding to fate transition. Our work presents a generalizable framework for top-down fate-decision studies and a practical approach to the taxonomy of cell fate decisions.

    1. Developmental Biology
    2. Evolutionary Biology
    Zhuqing Wang, Yue Wang ... Wei Yan
    Research Article

    Despite rapid evolution across eutherian mammals, the X-linked MIR-506 family miRNAs are located in a region flanked by two highly conserved protein-coding genes (SLITRK2 and FMR1) on the X chromosome. Intriguingly, these miRNAs are predominantly expressed in the testis, suggesting a potential role in spermatogenesis and male fertility. Here, we report that the X-linked MIR-506 family miRNAs were derived from the MER91C DNA transposons. Selective inactivation of individual miRNAs or clusters caused no discernible defects, but simultaneous ablation of five clusters containing 19 members of the MIR-506 family led to reduced male fertility in mice. Despite normal sperm counts, motility, and morphology, the KO sperm were less competitive than wild-type sperm when subjected to a polyandrous mating scheme. Transcriptomic and bioinformatic analyses revealed that these X-linked MIR-506 family miRNAs, in addition to targeting a set of conserved genes, have more targets that are critical for spermatogenesis and embryonic development during evolution. Our data suggest that the MIR-506 family miRNAs function to enhance sperm competitiveness and reproductive fitness of the male by finetuning gene expression during spermatogenesis.