The proneural wave in the Drosophila optic lobe is driven by an excitable reactiondiffusion mechanism
Abstract
In living organisms, selforganised 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 reactiondiffusion 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.001Introduction
The development of multicellular organisms relies on a multitude of transient coordination processes that provide the spatiotemporal cues for cell fate decisionmaking 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, largescale patterning is engineered by selforganised 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 selforganised reactiondiffusion processes (Lubensky et al., 2011; FormosaJordan 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 DeltaNotch 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 FatHippo 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).
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 DeltaNotch 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 coarsegrained 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 membranebound 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),
Here ${\varphi}^{\mathrm{E}}$ denotes the local strength of E activity, $\eta $ is the effective diffusion constant, $\mu $ is the gain rate in signalling activity due to positive feedback, $k$ is the decay rate, and $h(\varphi )={\varphi}^{n}/(1+{\varphi}^{n})$ is a Hill function parameterising the nonlinear positive feedback. Generically, reactiondiffusion systems involving diffusion and selfactivation 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).
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 reactiondiffusion system for the local strength of E and L activity, ${\varphi}^{\mathrm{E}}$ and ${\varphi}^{\mathrm{L}}$, is given by
where ${\mu}_{i}$ (with $i=\mathrm{E},\mathrm{L}$) indicate production rates and ${k}_{i}$ 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).
Integrated model of the proneural wave to include EGFRL’scNotch 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 DeltaNotch 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, DeltaNotch signalling is often associated with the emergence of ‘saltandpepper’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érezGómez et al., 2013; Sato et al., 2016). To address this question, we extended our minimal model to include canonical DeltaNotch interactions (Figure 2d) (Collier et al., 1996; Bray, 2006; Simakov and Pismen, 2013): (i) transactivation of Notch by Delta, (ii) downregulation of Delta by Notch within the same cell and (iii) cisinhibition (downregulation of Notch by Delta in the same cell). The model incorporates interactions between the DeltaNotch 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 DeltaNotch 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; OriharaOno et al., 2011; Weng et al., 2012; Contreras et al., 2018). According to the model, lateral inhibition promotes highDelta/lowNotch and lowDelta/highNotch 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 cisinhibition 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).
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 twodimensional 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).
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; OriharaOno et al., 2011; Wang et al., 2011b; PérezGó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 ‘saltandpepper pattern’ of Notch signalling activity (Bray, 2006; Shaya and Sprinzak, 2011). The reason for the absence of saltandpepper 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 saltandpepper 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.
We tested this prediction of the model by lowering Notch levels in the neuroepithelium but we did not observe ‘saltandpepper’ 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 Deltaregulated 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 saltandpepper 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 socalled ‘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 $\sigma $ measures the nonlinearity 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.
Here we probed the effects of the kinetic rate parameters of EGFR signalling, DeltaNotch 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 ($\eta $, ${\mu}_{\mathrm{E}}$ and ${k}_{\mathrm{E}}$) 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 ${\mu}_{\mathrm{L}}$ and ${k}_{\mathrm{L}}$) 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 ($\beta $) as well as its decay rate (${k}_{\mathrm{N}}$) 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 $\beta $ and ${k}_{\mathrm{N}}$) 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.
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 (targetlike) 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).
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 tightlyregulated propagating transition zone. In contrast with Turingbased activatorinhibitor mechanisms (Turing, 1952), which typically comprise fast diffusible inhibitors, the reactiondiffusion 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; FormosaJordan 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 DeltaNotch 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 reactiondiffusionlike 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 protocolFlies were raised on standard cornmeal medium at 25°C. Strains used were:
yw, hsFLP; FRT40A, tubGAL80/CyO, ActGFP; tubPGAL4,UASmCD8GFP/TM6B
w; FRT40A; UASPntP1/TM6B
w; FRT40A/CyO; N RNAi/TM6B (N RNAi lines used were BL33611 and BL33616)
HLHmgammaGFP (Almeida and Bray, 2005) was used to report active Notch signalling.
UASPntP1 clones and N RNAi clones
Request a detailed protocolTo 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 protocolLarval 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 antiDpn (1:10,000) and rat antiL’sc (1:5,000) (Caygill and Brand, 2017), chicken antiGFP (1:2,000) from Abcam, mouse antiDelta (1:50, C594.9B) from DSHB and mouse antiNotch (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. Reactiondiffusion 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, $\varphi =\varphi (\mathbf{\mathbf{x}},t)$, where $\mathbf{\mathbf{x}}$ denotes the position coordinate within the tissue and $t$ the time. (Throughout this text, we consider onedimensional systems for illustration purposes and twodimensional systems to describe the proneural wave in vivo. Hence, the position variable $\mathbf{\mathbf{x}}$ and the nabla operator $\nabla $ refer to $\mathbf{\mathbf{x}}=({x}_{1},{x}_{2})$ and $\nabla =(\partial /\partial {x}_{1},\partial /\partial {x}_{2})$ in the case of two dimensions and $\mathbf{\mathbf{x}}=x$ and $\nabla =\partial /\partial x$ in the case of one dimension.) The corresponding reactiondiffusion system is given by
where $\eta $ is the effective diffusion constant and $r(\varphi )$ is the reaction term describing the intracellular feedback,
Here, $\mu $ 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 $\mathrm{\Phi}$ being the threshold activity. Here, we choose a function of the Hill type (Novák and Tyson, 2008),
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 $\mathrm{\ell}$ with ‘noflux’ boundary conditions, ${(\partial \varphi /\partial x)}_{x=0}=0={(\partial \varphi /\partial x)}_{x=\mathrm{\ell}}$. 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 $\varphi (x,0)={\mathrm{e}}^{{x}^{2}/{x}_{0}^{2}}$ with ${x}_{0}^{2}=5\eta /k$.)
This behaviour is familiar in reactiondiffusion 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 ${\varphi}_{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 ${\varphi}_{0}^{*}$ and ${\varphi}_{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 ${\varphi}_{1}^{*}$ (white dot). If the system is only slightly perturbed around $\varphi ={\varphi}_{0}^{*}$ (black dot), the reaction dynamics will drive the system back towards this equilibrium, because degradation is stronger than selfactivation in this region, $r<0$. However, if the signalling activity rises above ${\varphi}_{1}^{*}$ (e.g., due to diffusion from the neighboring cell), selfactivation will outcompete degradation, $r>0$, and the signalling activity will be driven towards the other stable equilibrium ${\varphi}_{2}^{*}$. This reactiondiffusion mechanism is known as a bistable front (Keener and Sneyd, 2009).
Further insight into the dynamics of the front propagation can be gained using analytical techniques (Keener and Sneyd, 2009). Nondimensionalising the model by rescaling $\varphi \to \varphi /\mathrm{\Phi}$, $t\to kt$, and $x\to \sqrt{k/\eta}x$, the corresponding reactiondiffusion system in one dimension is given by
Here, $\lambda =\mu /(k\mathrm{\Phi})$ is the only (dimensionless) parameter and $\rho $ is the dimensionless reaction term.
Anticipating that the system gives rise to a travelling front with velocity $c$ and a stationary front profile $\varphi $, we make the ansatz $\varphi (x,t)=\varphi (xct)$ in Equation 6, which yields the ordinary differential equation
where the prime denotes the derivative with respect to the argument $y=xct$ of $\varphi $.
Multiplying this equation by ${\varphi}^{\prime}$ and integrating from $\mathrm{\infty}$ to $\mathrm{\infty}$ yields an implicit expression for the velocity of the front,
where $\mathrm{\Delta}U={\mathrm{lim}}_{y\to \mathrm{\infty}}[U(\varphi (y))U(\varphi (y))]$ with the potential $U(\varphi )$ defined by $\rho =\partial U/\partial \varphi $ (see Appendix 1—figure 1b). Analytical solutions to Equation 8 only exist for special functional forms of the reaction term $\rho $. To obtain analytical insights into the wave speed, we consider the limiting case of large Hill exponents, $n\to \mathrm{\infty}$ in Equation 5. This corresponds to a switchlike response of the gain rate once the activation threshold is reached. In this case, the reaction term acquires a piecewise linear functional form,
where $\mathrm{\Theta}$ is the Heaviside step function with the convention $\mathrm{\Theta}(0)=1/2$ (see Appendix 1—figure 1c). This function has three equilibria for $\lambda >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
where we have fixed the arbitrary position of the comoving reference frame by imposing the condition $\varphi (0)=1$. Appendix 1—figure 1e shows an example of the front solution Equation 10. The velocity of the front is given by
Expanding in the limit of large $\lambda $, dropping all orders higher than ${\lambda}^{1/2}$ and restoring the original parameter dependence by multiplying with the velocity scale $\sqrt{\eta k}$, we find
Equation 12 implies that the front velocity increases with increasing diffusion constant $\eta $ and synthesis rate $\mu $ as well as decreasing degradation rate $k$ and decreasing activation threshold $\mathrm{\Phi}$. 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 onecomponent model reproduces the key features of more detailed models
Previously, we sought to capture the spatiotemporal dynamics of EGFR signalling by a single component, described by the signalling activity $\varphi $, 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 reactiondiffusion system where ${\psi}^{\mathrm{E}}$ denotes the concentration of bound EGFRs, ${\psi}^{\mathrm{R}}$ represents the concentration of Rhomboid, and ${\psi}^{\mathrm{S}}$ is the concentration of the secreted active form of Spitz (sSpi),
Here, ${\lambda}_{\mathrm{E}}$ is the binding rate of EGFRs, ${\lambda}_{\mathrm{R}}$ is the synthesis rate of Rhomboid, and ${\lambda}_{\mathrm{S}}$ is the secretion rate of Spitz; the functions ${H}_{\mathrm{E}}$ and ${H}_{\mathrm{R}}$ 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.
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 diffusionlike effect to the entire E–R–S system. Therefore, the parameter $\mu $ of the onecomponent 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 $\eta $ describes the effective diffusibility, mediated by diffusion of component S. Formally, a connection to the onecomponent 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 $\partial {\psi}^{\mathrm{E}}/\partial t=0$ and $\partial {\psi}^{\mathrm{R}}/\partial t=0$ while the front is travelling. In this case, we can solve for ${\psi}^{\mathrm{E}}$ and ${\psi}^{\mathrm{R}}$ and obtain a closed equation for $\varphi \equiv {\psi}^{\mathrm{S}}$,
where $\lambda \equiv {\lambda}_{\mathrm{S}}$ and $h(\varphi )={\lambda}_{\mathrm{R}}{H}_{\mathrm{R}}({\lambda}_{\mathrm{E}}{H}_{\mathrm{E}}(\varphi ))$. Note that Equation 14 is formally equivalent to Equation 6; if ${H}_{\mathrm{E}}$ and ${H}_{\mathrm{R}}$ are sigmoidal functions with the qualitive shape of Equation 5, then so is $h$.
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 reactiondiffusion equations for the two fields ${\varphi}^{\mathrm{E}}$ and ${\varphi}^{\mathrm{L}}$ are given by
with the reaction terms
where we have introduced the Hill function $\overline{h}(\varphi )\equiv 1h(\varphi )$ which describes an inhibitory effect. For simplicity, we have also considered identical unit concentration thresholds for activation and inhibition. Here, ${\lambda}_{\mathrm{E}}$ and ${\lambda}_{\mathrm{L}}$ are the gain rates of the components E and L, respectively. Note that in the absence of L (${\varphi}^{\mathrm{L}}=0$), the reaction term ${\rho}_{\mathrm{E}}$ reduces to the one of the onecomponent model given in Equation 6, ${\rho}_{\mathrm{E}}({\varphi}^{\mathrm{E}},0)=\rho ({\varphi}^{\mathrm{E}})$ (Appendix 2—figure 1). Again, we consider a finite onedimensional domain of length $\mathrm{\ell}$ and noflux boundary conditions, ${(\partial {\varphi}^{\mathrm{E}}/\partial x)}_{x=0}=0={(\partial {\varphi}^{\mathrm{E}}/\partial x)}_{x=\mathrm{\ell}}$ and ${(\partial {\varphi}^{\mathrm{L}}/\partial x)}_{x=0}=0={(\partial {\varphi}^{\mathrm{L}}/\partial x)}_{x=\mathrm{\ell}}$.
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 $\varphi (x,0)={\mathrm{e}}^{{x}^{2}/{x}_{0}^{2}}$ with ${x}_{0}^{2}=10\eta /k$.) Again, this behaviour can be elucidated by studying the local reaction dynamics, which is now given by the vector field $\mathbf{\mathbf{F}}=({\rho}_{\mathrm{E}},{\rho}_{\mathrm{L}})$, defined by Equation 16) (Appendix 2—figure 1). Note that now only the point $({\varphi}^{\mathrm{E}},{\varphi}^{\mathrm{L}})=(0,0)$ is a stable equilibrium. The black and white dots indicate points with ${\varphi}^{\mathrm{L}}=0$ that satisfy ${\rho}_{\mathrm{E}}=0$ as in the onecomponent system (compare to Appendix 1—figure 1a). However, note that only the point $(0,0)$ also satisfies ${\rho}_{\mathrm{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 selffeedback. Thus, the system finally returns to the fixed point $({\varphi}^{\mathrm{E}},{\varphi}^{\mathrm{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.
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 cisinhibition (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). DeltaNotch 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 lowDelta/highNotch and highDelta/lowNotch 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 ${\varphi}_{\mathbf{\mathbf{x}}}^{i}(t)$ with $i=\mathrm{E},\mathrm{L},\mathrm{D},\mathrm{N}$ describing EGFR signalling activity, L’sc expression, and Delta and Notch, respectively. The index $\mathbf{x}$ indicates the lattice site. Furthermore, we introduce the cell state ${\mathrm{\Omega}}_{\mathbf{\mathbf{x}}}(t)$ which takes values from 0 to 1, where ${\mathrm{\Omega}}_{\mathbf{\mathbf{x}}}=0$ indicates that cell $\mathbf{x}$ is a neuroepithelial cell and ${\mathrm{\Omega}}_{\mathbf{\mathbf{x}}}=1$ indicates that cell $\mathbf{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; OriharaOno et al., 2011), we include a basal source of Notch that is independent of transactivation 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 ${\varphi}_{\mathbf{\mathbf{x}}}^{i}$ with $i=\mathrm{E},\mathrm{L},\mathrm{D},\mathrm{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
Here, the parameters ${\mu}_{i}$ denote gain rates, ${k}_{i}$ denote decay rates, $\beta $ denotes the basal gain rate of the component N, and the ${\mathrm{\Phi}}_{i}$ denote threshold levels for positive and negative feedbacks. The operators $\widehat{\mathrm{\Delta}}$ and $\widehat{\mathrm{\Sigma}}$ are the lattice Laplacian and the sum over concentrations of neighboring lattice sites, respectively, and defined by ${[\widehat{\mathrm{\Delta}}\varphi ]}_{\mathbf{\mathbf{x}}}={\sum}_{y\in {U}_{\mathbf{\mathbf{x}}}}({\varphi}_{\mathbf{\mathbf{y}}}{\varphi}_{\mathbf{\mathbf{x}}})$ and ${[\widehat{\mathrm{\Sigma}}\varphi ]}_{\mathbf{\mathbf{x}}}={\sum}_{y\in {U}_{\mathbf{\mathbf{x}}}}{\varphi}_{\mathbf{\mathbf{y}}}$, where ${U}_{\mathbf{\mathbf{x}}}$ is the set of neighbours of site $x$. The dynamics of the cell state ${\mathrm{\Omega}}_{\mathbf{\mathbf{x}}}$ is given by
where the function $V$ is a ‘potential’ for the cell state which ensures that in the absence of signalling and proneural gene expression, ${\mathrm{\Omega}}_{\mathbf{\mathbf{x}}}$ has two stable equilibria: ${\mathrm{\Omega}}_{\mathbf{\mathbf{x}}}=0$ (neuroepithelium) and ${\mathrm{\Omega}}_{\mathbf{\mathbf{x}}}=1$ (neuroblast) (see Appendix 3—figure 1). The qualitative features of our model do not depend on the exact choice of $V$. The term ${f}_{\mathrm{int}}$ 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 ${f}_{\mathrm{int}}$ 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
The cell state potential $V$ has wells at $\mathrm{\Omega}=0$ (neuroepithelial state) and $\mathrm{\Omega}=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 ${f}_{\mathrm{int}}$ leads to a transition of the cell state from neuroepithelial ($\mathrm{\Omega}=0$) to neuroblast ($\mathrm{\Omega}=1$) when (i) L exceeds the threshold level ${\mathrm{\Phi}}_{\mathrm{int}}$ and/or (ii) E exceeds and N drops below the threshold level.
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
where $\gamma $ denotes the noise strength and ${\xi}_{i}$ denotes Gaussian white noise characterised by the expectation values $\u27e8{\xi}_{i}(t)\u27e9=0$ and $\u27e8{\xi}_{i}(t){\xi}_{j}({t}^{\prime})\u27e9={\delta}_{ij}\delta (t{t}^{\prime})$. Furthermore, ${\delta}_{ij}$ denotes the Kronecker delta and $\delta (t)$ the Dirac delta distribution.
Appendix 4—figure 1 shows numerical examples of the system for different noise strengths $\gamma $. 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 ${\mu}_{\mathrm{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).
Robustness against lattice defects
To test whether the proposed mechanism is robust with respect to disordered lattice structures, we considered a sitediluted 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 diffusionmediated 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 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 $\mathbf{x}=1,2$,
where $\overline{\mathbf{x}}$ refers to the respective other cell and $\overline{h}(\varphi )\equiv 1h(\varphi )$ with $h$ being the Hill function Equation 5, as before. Here, $\lambda $ is the gain rate for Notch and $\beta $ indicates the basal production. For simplicity, we consider a linear positive feedback in the dynamics of Notch and the Hill exponent $n=2$ for $\overline{h}$. In steady state, where $\mathrm{d}{\varphi}_{\mathbf{\mathbf{x}}}^{\mathrm{D}}/\mathrm{d}t=0=\mathrm{d}{\varphi}_{\mathbf{\mathbf{x}}}^{\mathrm{N}}/\mathrm{d}t$, we can eliminate ${\varphi}_{1}^{\mathrm{N}}$ and ${\varphi}_{2}^{\mathrm{N}}$ to obtain
where $\mathrm{\Gamma}(\varphi )=\overline{h}(\beta +\lambda \varphi )$. From this, it follows that both ${\varphi}_{1}^{\mathrm{D}}$ and ${\varphi}_{2}^{\mathrm{D}}$ satisfy ${\varphi}_{\mathbf{\mathbf{x}}}^{\mathrm{D}}=\mathrm{\Gamma}(\mathrm{\Gamma}({\varphi}_{\mathbf{\mathbf{x}}}^{\mathrm{D}}))$. Among other solutions, this equation has two solutions of the form ${\varphi}_{\pm}=p\pm \sqrt{q}$ with
which, in the case that they are real, correspond to the lowDelta/highNotch and highDelta/lowNotch states in adjacent cells as they satisfy ${\varphi}_{\pm}=\mathrm{\Gamma}({\varphi}_{\mp})$.
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
Figure 4e shows the corresponding phase diagram for the occurrence of lateral inhibition in the twocell system. Therefore, a basal production term, if strong enough, can prevent lateral inhibition patterns.
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 steadystate 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 ${{\varphi}^{\mathrm{E}}}_{t=0}={\mu}_{\mathrm{E}}/2$ in a onedimensional array of cells in the outermost cell layer that have angles between $\pi /3$ and $5\pi /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 ${\mu}_{\mathrm{E}}\overline{h}(2{\mathrm{\Omega}}_{\mathbf{\mathbf{x}}})$ 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.
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 socalled 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 $\mathcal{O}$ and a given set of parameters ${\theta}_{1},\mathrm{\dots},{\theta}_{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 socalled ‘elementary effects’
That is the finitedifference quotient of the output with respect to the parameter, given a finite step size $\mathrm{\Delta}$ 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 ${P}_{i}(e)$ of elementary effects for each parameter $i$, from which the following sensitivity indices are computed,
where ${\u27e8\cdot \u27e9}_{i}$ denotes the expectation value under the distribution ${P}_{i}$. 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 latticebased 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 ${\overline{\mathrm{\Omega}}}_{x}={\mathrm{\ell}}_{\u27c2}^{1}{\sum}_{y}{\mathrm{\Omega}}_{xy}$ as the average of the cell state $\mathrm{\Omega}$ in the direction perpendicular to the wave, where ${\mathrm{\ell}}_{\u27c2}$ 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
Practically, we determine $v$ as the slope obtained from a linear fit of ${\sum}_{x}{\overline{\mathrm{\Omega}}}_{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<\mathrm{\Omega}<1$. Formally, this width can be defined as
The discrete derivative, ${\overline{\mathrm{\Omega}}}_{x}{\overline{\mathrm{\Omega}}}_{x1}$, which measures the steepness of the profile, is nonzero only in the transition zone. For example, for a Gaussian profile of the discrete derivative, ${\overline{\mathrm{\Omega}}}_{x}{\overline{\mathrm{\Omega}}}_{x1}\propto {\mathrm{e}}^{x/2{\sigma}^{2}}$ with a variance ${\sigma}^{2}\gg 1$, Equation 28 yields $w=\sigma $. 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 $\eta $, $\beta $ and ${\mu}_{i}$, ${k}_{i}$ for $i=\mathrm{E},\mathrm{L},\mathrm{D},\mathrm{N}$ and allow them to vary between the 0.2fold and 5fold 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.
Data availability
All data generated or analysed during this study are included in the manuscript.
References

Regulation of postembryonic neuroblasts by Drosophila grainyheadMechanisms of Development 122:1282–1293.https://doi.org/10.1016/j.mod.2005.08.004

BookThe Developmental Origin of Cell Type Diversity in the Drosophila Visual SystemCham, Switzerland: Springer International Publishing.

Notch signalling: a simple pathway becomes complexNature Reviews Molecular Cell Biology 7:678–689.https://doi.org/10.1038/nrm2009

An effective screening design for sensitivity analysis of large modelsEnvironmental Modelling & Software 22:1509–1518.https://doi.org/10.1016/j.envsoft.2006.10.004

Pattern formation by lateral inhibition with feedback: a mathematical model of deltanotch intercellular signallingJournal of Theoretical Biology 183:429–446.https://doi.org/10.1006/jtbi.1996.0233

A discrete model of Drosophila eggshell patterning reveals cellautonomous and juxtacrine effectsPLOS Computational Biology 10:e1003527.https://doi.org/10.1371/journal.pcbi.1003527

Regulation of neuronal differentiation at the neurogenic wavefrontDevelopment 139:2321–2329.https://doi.org/10.1242/dev.076406

A model of the Spatiotemporal dynamics of Drosophila eye disc developmentPLOS Computational Biology 12:e1005052.https://doi.org/10.1371/journal.pcbi.1005052

Continuum theory of gene expression waves during vertebrate segmentationNew Journal of Physics 17:093042.https://doi.org/10.1088/13672630/17/9/093042

Fat/Hippo pathway regulates the progress of neural differentiation signaling in the Drosophila optic lobeDevelopment, Growth & Differentiation 53:653–667.https://doi.org/10.1111/j.1440169X.2011.01279.x

BookMathematical Physiology (2nd ed)Springer Science+Business Media.https://doi.org/10.1007/9780387758473

EGF receptor signalling: the importance of presentationCurrent Biology 10:R388–R391.https://doi.org/10.1016/S09609822(00)004851

Signal propagation and failure in discrete autocrine relaysPhysical Review Letters 93:500–504.https://doi.org/10.1103/PhysRevLett.93.118101

Design principles of biochemical oscillatorsNature Reviews Molecular Cell Biology 9:981–991.https://doi.org/10.1038/nrm2530

Switch and template pattern formation in a discrete reactiondiffusion system inspired by the Drosophila eyeThe European Physical Journal E 33:129–148.https://doi.org/10.1140/epje/i2010106476

Signaling mechanisms controlling cell fate and embryonic patterningCold Spring Harbor Perspectives in Biology 4:a005975.https://doi.org/10.1101/cshperspect.a005975

Pattern formation in the Drosophila eye discThe International Journal of Developmental Biology 53:795–804.https://doi.org/10.1387/ijdb.072483jr

Waves of differentiation in the fly visual systemDevelopmental Biology 380:1–11.https://doi.org/10.1016/j.ydbio.2013.04.007

Fiji: an opensource platform for biologicalimage analysisNature Methods 9:676–682.https://doi.org/10.1038/nmeth.2019

From Notch signaling to finegrained patterning: Modeling meets experimentsCurrent Opinion in Genetics & Development 21:732–739.https://doi.org/10.1016/j.gde.2011.07.007

Modeling and computational analysis of EGF receptormediated cell communication in Drosophila oogenesisDevelopment 129:2577–2589.

Noise in biologyReports on Progress in Physics 77:026601–026630.https://doi.org/10.1088/00344885/77/2/026601

The chemical basis of morphogenesisPhilosophical Transactions of the Royal Society of London. Series B, Biological Sciences 237:37–72.https://doi.org/10.1098/rstb.1952.0012

Role of JAK/STAT signaling in neuroepithelial stem cell maintenance and proliferation in the Drosophila optic lobeBiochemical and Biophysical Research Communications 410:714–720.https://doi.org/10.1016/j.bbrc.2011.05.119

Changes in Notch signaling coordinates maintenance and differentiation of the Drosophila larval optic lobe neuroepitheliaDevelopmental Neurobiology 72:1376–1390.https://doi.org/10.1002/dneu.20995

Sensitivity analysis of infectious disease models: methods, advances and their applicationJournal of The Royal Society Interface 10:20121018.https://doi.org/10.1098/rsif.2012.1018
Decision letter

Naama BarkaiSenior and Reviewing Editor; Weizmann Institute of Science, Israel
In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.
Thank you for submitting your article "The proneural wave in the Drosophila optic lobe is driven by an excitable reactiondiffusion mechanism" for consideration by eLife. Your article has been reviewed by two reviewers and the evaluation has been overseen by Naama Barkai as the Reviewing and Senior Editor. The reviewers have opted to remain anonymous.
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
It is uniformly agreed that the system is interesting and the theoretical work valid. There are several issues, however, that will need to be addressed.
1) Motivating the study: given that the system was previously studied (Sato paper), it may be necessary to better define why a new model is 'needed', and what new biological insights it can bring to the table. For example, what are the mysteries that are still not understood and can be better explained by a more indepth quantitative study? This is discussed in Appendix 3 but should be explained also in the main text.
2) Related to point 1, not all statements in Appendix 3 are supported – e.g. the statement that the Sato model is not robust. Also, the new model may not be fundamentally different from the Sato model, but rather extends it: models clearly make useful simplifications. The current work may 'open' some 'blackbox' and by this provide deeper insights. This does not reduce the significance of the study. Please address that.
3) The analysis of Appendix 5 is interesting and should be part of the main text. The lack of lateral inhibition in this case presents a key biological insight that distinguishes the present study. The experimental test proposed by reviewer #1 below is therefore important.
4) Please explain all assumptions, interactions and parameters, as requested below.
Finally, please address all comments and suggestions in the detailed reviews below.
Reviewer #1:
This paper describes a model to explain the progression of a wave of differentiation that crosses a neuroepithelium to generate a brain structure in the optic lobes of Drosophila. This type of wave has been described (and modeled) quite extensively, in particular the morphogenetic furrow that patterns the Drosophila eye. In this case, we are dealing with the neurogenic wave that transforms epithelial cells into neural stem cells (neuroblasts) and produces one of the optic lobe structures. In fact, a recent paper had also modeled the same wave of differentiation and this paper is discussed here.
This type of patterning mechanisms is likely to be fairly common in biology and it is thus quite important to understand its fundamental aspects, which means that a comparison with the morphogenetic furrow is important to emphasize the common points and points of divergence. This will therefore be an important paper for physicists and for more theoretically oriented biologists who want to understand the rules of patterning.
This being said, as a biologist, I do not feel completely qualified to comment on the mathematics of the model. Yet, I think that the authors should have made more efforts to make more accessible the aspects of the model in which a biologist might be interested and might want to use for his/her own research. For example, it is quite difficult to find a clear description in one place of all the interactions (positive and negative) between pathway components that are actually used in the model. Figure 2D comes closest to doing so but, as with other parts of the paper, it refers to an appendix that is disconnected from the main text. The authors should therefore explain in much more precise terms the exact interactions and must justify why each arrow exists, and or why some of them (which have also been described as important for the process) have been excluded.
In particular:
 What is the evidence to support Notch promotion of EGFR signaling?
 Why isn't proneural promotion of Dl expression included?
 Is Notch repression of cell fate conversion necessary on top of the Notch inhibition of proneural factors?
 The authors do not even mention the data on JAK/STAT and Hippo affecting proneural wave progression. Why do they ignore this?
Some of these justifications have been relegated to appendices and it would make the paper much more important to the biologists if much of Appendix 3 (and similarly for Appendix 5) were to be part of the main text. The authors must do this.
There are several comparisons with the recent paper from the Sato lab, with assertions that the current model makes more accurate predictions. However, both papers use relatively crude assessments of phenotype: acceleration or retardation of the proneural wave. From the model simulations in Figure 3, the authors show that there are many 'flavors' to proneural acceleration (or retardation), i.e. different parameters can produce similar outcomes. If the authors really want to say their model is more accurate or better, then they should test model predictions appropriately (i.e. look at EGFR activity and Notch and Dl levels, in addition to Dpn and/or L'sc for the experimental conditions detailed in Figure 3. All of this should be quite easy and would add significant validation to the model.
Furthermore, the authors should also explain why their model of proneural wave progression adds to previous similar models of wave progression, in particular the morphogenetic furrow in the eye disc, also in Drosophila. One key feature that is distinct with the proneural wave is the lack of Notchmediated lateral inhibition. As this was a puzzling question, the explanation that this is due to basal low levels of Notch in the neuroepithelium is very attractive. A prediction of the model (Appendix 5, Figure 1B) is that reducing Notch levels in the neuroepithelium should result in a saltandpepper pattern typical of lateral inhibition. This would be a key result that can very easily be tested experimentally by using a neuroepithelial driver to knock down to different extents Notch using RNAi. This would provide a qualitatively different prediction from just acceleration or retardation of the wave and would convince the biologist that the model is useful to address this point.
Reviewer #2:
The presented article proposes a mathematical model for a proneural wave in the Drosophila optical lobe as driven by an excitable reaction diffusion mechanism, that explains observed data and predicts a new experiment. The biochemical nature of the wave is due to a complex process, involving interaction of EGFR and DeltaNotch signaling pathways. Their interplay establishes a transition zone, that travels over the tissue and triggers a differentiation wave.
Although a phenomenological model of this process exists, the authors propose a new approach. This is justified as it helps to analyze at greater depth how this biochemical process, with interaction between multiple components, results in the wave. As such the article is an example for how mainly model driven approaches may help refining the view on a developmentally relevant, but poorly understood process.
The structure of the article is well developed. The appendix, in particular for the explanation of traveling front and pulse models, is didactically useful for a readership that has a background in life sciences.
Possible improvements:
In general, the theory is compared to data at a qualitative level. If possible, it would help to develop a more quantitative interface to experiments. For example, Figure 2E discusses a transient activation of notch signaling that is mentioned to also occur in the optic lobe. The model predicts a striking spatial profile, that may be contrasted to a similar profile of data for a few of the genes involved.
Further, Figure 3 shows a comparison of model simulations on a 2D lattice to experiment. One of the striking differences is that the simulation is based on an ordered hexagonal lattice, while data looks like a strongly disordered array. Order may be less important in case of strong effects such as shown in panels A,B,D. It remains unclear however, how the underlying lattice order will affect the more subtle effects shown in the cell state variable for panels C,E,F.
In Appendix 3, the authors explain the concrete realization that H_{0} is not important to capture the qualitative features of the model. However, in the spirit of the welldeveloped didactic tone in the preceding appendices, the model might become more transparent to a broader audience, if key features of H_{0} are further developed. The mathematics is clear, but additional graphics and a brief discussion would help to guide the intuition of this aspect.
The following Appendix 4 presents a valuable indepth analysis about robustness against fluctuations. However, it is difficult to find a more detailed discussion about choice of parameters presented in Appendix 3—Table 1. Are these parameters chosen in a physiologically relevant regime? What's the effect of changes to the model? For future experimentalists aiming to test this model, it may be good to provide a range of possible parameter values in which the presented features will be valid. For example, is there a microscopic motivation for the choice of the Hill coefficient n = 3?
This is an exciting model, and it would be very helpful for future development if the model could provide additional predictions that are potentially difficult to test using current technology. For example, could the model be used to predict what sets the speed of the wave? This point is not essential to support the major conclusions of the article. But in the interest of theory motivating new experiments, such an addition might be helpful.
https://doi.org/10.7554/eLife.40919.030Author response
It is uniformly agreed that the system is interesting and the theoretical work valid. There are several issues, however, that will need to be addressed,
1) Motivating the study: given that the system was previously studied (Sato paper), it may be necessary to better define why a new model is 'needed', and what new biological insights it can bring to the table. For example, what are the mysteries that are still not understood and can be better explained by a more indepth quantitative study? This is discussed in Appendix 3 but should be explained also in the main text.
We have expanded the Introduction and Discussion section of our manuscript as described below.
First, the paper by Sato et al., has clearly pioneered a modelling approach on proneural wave progression, however, it remains vague on the molecular basis of the wave. It states that EGFR signalling triggers differentiation, while differentiation triggers EGFR signalling activity, resulting in essentially a twostate system that propagates by construction. It does not address the emergence of major characteristics of the wave such as spatially confined proneural gene expression in a localised transition zone. Moreover, a concrete molecular mechanism and/or cascade that could be pinpointed (and therefore falsified) or reenacted, does not exist. To highlight this more clearly, we have expanded the corresponding section in the Introduction, mentioning the lack of understanding of the following four major points:
i) the dynamical nature of the wave,
ii) the emergence of a localised transition zone with spatially confined expression of the proneural gene, l’sc
iii) the specific profiles of gene expression and signalling activity around the transition zone,
iv) the effects of molecular kinetics on wave propagation,
In the Discussion section, we have added new sections to highlight the specific explanatory and mechanistic novelties.
In addition, the Sato model exhibits some rather puzzling aspects in its dynamical implementation, such as the linearity of DeltaNotch interactions and its lack of resilience to noise. Rather than provide a critical discussion of these more detailed aspects in the main text, this comparison is discussed in Appendix 3.
2) Related to point 1, not all statements in Appendix 3 are supported – e.g. the statement that the Sato model is not robust.
This point refers to the following sentence in our manuscript: "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."
Indeed, this statement is already made, in a more formal way, in the original paper by Sato et al., 2016: “Because the equation for A is unstable at A = 0, the simple addition of noise induces spontaneous differentiation apart from the wave front.”
(Here, A is the variable encoding both the cell state and proneural gene expression; in the words of the authors: “A is an abstract value specifying the state of differentiation […]. We assume that A is intimately related to the expression levels of ASC proteins.” The instability arises due to the fact that any perturbation of A away from 0 drives A towards the differentiated state if EGF expression is sufficiently large.)
Since a repetition of the formal argument already given in the original reference would be redundant, we have added "as discussed in Sato et al., (2016)." to the corresponding sentence in Appendix 3.
Also, the new model may not be fundamentally different from the Sato model, but rather extends it: models clearly make useful simplifications. The current work may 'open' some 'blackbox' and by this provide deeper insights. This does not reduce the significance of the study. Please address that.
Our model is formulated in a framework of reactiondiffusion systems similar to that used in Sato et al. and we do not dispute that this is a sensible and natural framework to describe the process at hand. However, we do not believe that our model is an extension or unfolding of the Sato model. The Sato model is, in vital parts, based on the phenomenology of the wave. To exemplify this, let us look at the formal mechanism that leads to stabilisation of the wave in Sato et al., (2016). EGF is diffusible but the wave can only propagate if EGF activity is directly or indirectly activated:
“However, EGF ligand and Dl are only produced at the interface between undifferentiated and differentiated cells but not produced in differentiated NBs. Thus, we include the EGF and Dl production terms a_{e} A(A_{0} – A) and a_{d} A(A_{0} – A), respectively.” (Sato et al., 2016).
It is not straightforward to imagine the biochemical basis of such an interaction, which activates EGF signalling for intermediate values of the abstract differentiation variable A. Importantly, our model neither builds upon nor expands on such or similar notions.
Instead, our model is constructed by gradually building up the phenomenology of the wave from known molecular interactions. To our knowledge, the idea that autocrine EGF signalling drives the proneural wavefront and that differentiation shuts down this driver mechanism and renders it excitable has not been presented before. It is certainly not implicit in the model of Sato et al., where in fact the molecular mechanism of wave propagation remains rather cryptic for the above reasons. We would therefore argue that our model provides conceptually different insights rather than a refinement of previous work.
3) The analysis of Appendix 5 is interesting and should be part of the main text. The lack of lateral inhibition in this case presents a key biological insight that distinguishes the present study. The experimental test proposed by reviewer #1 below is therefore important.
We have moved the corresponding sections about the suppression of lateral inhibition patterns and the role of cisinhibition to the main text, together with the corresponding simulation figures, which constitute the new Figure 4 in the main text.
4) Please explain all assumptions, interactions and parameters, as requested below.
In the revised manuscript, we have partitioned the original section on the full model of the proneural wave into two separate sections, one of which explains the structure of the full model (including a reference for each interaction) shown in Figure 2D.
Finally, please address all comments and suggestions in the detailed review below.
Reviewer #1:
[…] This will therefore be an important paper for physicists and for more theoretically oriented biologists who want to understand the rules of patterning.
This being said, as a biologist, I do not feel completely qualified to comment on the mathematics of the model. Yet, I think that the authors should have made more efforts to make more accessible the aspects of the model in which a biologist might be interested and might want to use for his/her own research. For example, it is quite difficult to find a clear description in one place of all the interactions (positive and negative) between pathway components that are actually used in the model. Figure 2D comes closest to doing so but, as with other parts of the paper, it refers to an appendix that is disconnected from the main text. The authors should therefore explain in much more precise terms the exact interactions and must justify why each arrow exists, and or why some of them (which have also been described as important for the process) have been excluded.
Please see our response to point 4 from the reviewing editor.
In particular:
 What is the evidence to support Notch promotion of EGFR signaling?
It has been shown previously that expression of N^{act} results in the activation of PntP1 (Yasugi et al., 2010, see Figure 6D,D').
 Why isn't proneural promotion of Dl expression included?
As a starting point for representing DeltaNotch signalling in our model, we had originally chosen the wellestablished model of (Collier et al., 1996), which includes the effective feedback between Notch and Delta (via proneural genes) as a direct interaction between effective variables. In fact, we have checked in simulations that an additional explicit promotion of Dl expression through proneural genes does not alter the phenomenology of the model. This is due to the fact that such an additional interaction would only reinforce an already existing feedback: Notch would suppress l'sc, which, in turn, would lead to decreased activation of Delta, i.e., would provide a negative feedback that is already present. A high degree of functional redundancy with respect to additional interactions is common in fairly complex genetic regulatory networks, which generically fall into the “sloppiness universality class” (see, e.g., the excellent papers by Gutenkunst et al., 2007; Daniels et al., 2008; Waterfall et al., 2009). To maintain the clarity of what is already a complex model, we have therefore opted to avoid such manifest redundancies.
 Is Notch repression of cell fate conversion necessary on top of the Notch inhibition of proneural factors?
It is known from the literature that expression of the proneural gene l'sc is sufficient but not necessary for the NEtoNB transition (Yasugi et al., 2008; Egger et al., 2010; Sato et al., 2013). Furthermore, Notch is thought to regulate the width of the transition zone, making the spatially confined region within the transition zone (in which Notch activity is downregulated) a prime candidate for a region "permissive" for transition. In this spirit, Notch acts as an inhibitor of the transition in our model. Otherwise, EGFR signalling alone would be able to initiate the transition, for which we have no reason to assume is the case.
 The authors do not even mention the data on JAK/STAT and Hippo affecting proneural wave progression. Why do they ignore this?
While we are aware that the JAK/STAT and Hippo pathways modulate proneural wave progression, neither seems to be a crucial driver of the wave itself, nor do their expression patterns suggest a concomitant motion with the wave. Rather, the broad expression of their ligands and/or receptors in the neuroepithelium suggests that they act as modulators of proneural wave progression that prevent premature transition (in the case of JAK/STAT) and to regulate the mitotic activity of neuroepithelial cells (in the case of Hippo) (Yasugi et al., 2008; Reddy et al., 2010; Sato et al., 2013). This is why we refrained from including them in our manuscript. In the revised manuscript, we have added to the Discussion section a paragraph on these pathways, elaborating on their role in modulating proneural wave progression, while not being drivers of the wave, as suggested by experiments.
Some of these justifications have been relegated to appendices and it would make the paper much more important to the biologists if much of Appendix 3 (and similarly for Appendix 5) were to be part of the main text. The authors must do this.
As explained above, we have partitioned the "full model" section in the manuscript into two sections, the first of which ("Integrated model of the proneural wave") is devoted to explaining all model interactions including references. The second section ("Congruence with experimental data") contains the considerations about the suppression of lateral inhibition patterns that were originally part of Appendix 5, as suggested by the referee. We believe that this increases the clarity of the manuscript and highlights the significance of the findings. Nevertheless, we have opted to keep the formal, more mathematical aspects of the full model formulation and the suppression of lateral inhibition in the respective appendices, where they are still available to the interested reader.
There are several comparisons with the recent paper from the Sato lab, with assertions that the current model makes more accurate predictions. However, both papers use relatively crude assessments of phenotype: acceleration or retardation of the proneural wave. From the model simulations in Figure 3, the authors show that there are many 'flavors' to proneural acceleration (or retardation), i.e. different parameters can produce similar outcomes. If the authors really want to say their model is more accurate or better, then they should test model predictions appropriately (i.e. look at EGFR activity and Notch and Dl levels, in addition to Dpn and/or L'sc for the experimental conditions detailed in Figure 3. All of this should be quite easy and would add significant validation to the model.
We have performed a general sensitivity analysis of the model, which illuminates the influence of a parameter on key observables of the system such as the proneural wave speed and the width of the transition zone. Many of these parameters refer to interactions that are not present in the Sato model, such as the ones entailed by independent variables for proneural gene expression and differentiation status.
Furthermore, the authors should also explain why their model of proneural wave progression adds to previous similar models of wave progression, in particular the morphogenetic furrow in the eye disc, also in Drosophila.
In the revised manuscript, we have expanded the paragraph in the Discussion mentioning progression of the morphogenetic furrow as a similar process, but also pointing out that the eye disc exhibits quite different features with regard to growth and deformation of the underlying tissue as well as the network structure that leads to propagation of the furrow.
One key feature that is distinct with the proneural wave is the lack of Notchmediated lateral inhibition. As this was a puzzling question, the explanation that this is due to basal low levels of Notch in the neuroepithelium is very attractive. A prediction of the model (Appendix 5, Figure 1B) is that reducing Notch levels in the neuroepithelium should result in a saltandpepper pattern typical of lateral inhibition. This would be a key result that can very easily be tested experimentally by using a neuroepithelial driver to knock down to different extents Notch using RNAi. This would provide a qualitatively different prediction from just acceleration or retardation of the wave and would convince the biologist that the model is useful to address this point.
We tested this prediction by lowering Notch levels in the neuroepithelium but did not observe 'saltandpepper' patterns of Delta/Notch expression within clones expressing Notch RNAi (see the new Figure 4C and D). However, the absence of the emergence of lateral inhibition is likely due to the complete loss of detectable Notch in cells expressing Notch RNAi, while the reduction of Notch levels in the model prediction is more subtle. Referring to the 'phase diagram' in the new Figure 4E, it can be seen that both basal and Deltaregulated 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 saltandpepper patterns.
Reviewer #2:
[…] The structure of the article is well developed. The appendix, in particular for the explanation of traveling front and pulse models, is didactically useful for a readership that has a background in life sciences.
Possible improvements:
In general, the theory is compared to data at a qualitative level. If possible, it would help to develop a more quantitative interface to experiments. For example, Figure 2E discusses a transient activation of notch signaling that is mentioned to also occur in the optic lobe. The model predicts a striking spatial profile, that may be contrasted to a similar profile of data for a few of the genes involved.
We have included a new image to show this more clearly. Please see Figure 1C.
Further, Figure3 shows a comparison of model simulations on a 2D lattice to experiment. One of the striking differences is that the simulation is based on an ordered hexagonal lattice, while data looks like a strongly disordered array. Order may be less important in case of strong effects such as shown in panels A,B,D. It remains unclear however, how the underlying lattice order will affect the more subtle effects shown in the cell state variable for panels C,E,F.
While the reviewer is right that, in general, the underlying lattice order can have strong effects on patterning phenomena, especially in systems with repulsive interactions where defects and frustration effects can play an important role (e.g., in those cases where lateral inhibition does occur!), the lattice structure is less important in diffusive systems with suppressed repulsion such as the one reported here. Hence, we have no reason to assume that the effects of clones illustrated here would be systematically different on different lattice structures. If this were the case, even the results of clonal experiments would depend on the specific local topology of the clone for which we have no indication. To illustrate this point, we have performed simulations on sitediluted lattices that introduce random defects in the lattice structure. These simulations have been included in the robustness analysis in Appendix 4, showing that both the general mechanism of wave propagation as well as the phenomenology of the more subtle clonal effects remain intact, as illustrated using the clonal scenario in which Notch has been upregulated.
In Appendix 3, the authors explain the concrete realization that H_{0} is not important to capture the qualitative features of the model. However, in the spirit of the welldeveloped didactic tone in the preceding appendices, the model might become more transparent to a broader audience, if key features of H_{0} are further developed. The mathematics is clear, but additional graphics and a brief discussion would help to guide the intuition of this aspect.
To improve the motivation for the choice of this function that encapsulates the two possible cell states (neuroepithelial and neuroblast), we have replaced the function H_{0} by a 'potential' function whose derivative is H_{0}. We have added an additional figure in Appendix 3, where it is referred to, explaining the interpretation of this potential function in the spirit of the preceding appendices, cf. Appendix 1—Figure 1B. We believe this has improved the clarity of the model description.
The following Appendix 4 presents a valuable indepth analysis about robustness against fluctuations. However, it is difficult to find a more detailed discussion about choice of parameters presented in Appendix 3—Table 1. Are these parameters chosen in a physiologically relevant regime? What's the effect of changes to the model? For future experimentalists aiming to test this model, it may be good to provide a range of possible parameter values in which the presented features will be valid. For example, is there a microscopic motivation for the choice of the Hill coefficient n = 3?
This is an exciting model, and it would be very helpful for future development if the model could provide additional predictions that are potentially difficult to test using current technology. For example, could the model be used to predict what sets the speed of the wave? This point is not essential to support the major conclusions of the article. But in the interest of theory motivating new experiments, such an addition might be helpful.
First, we would like to mention that the core phenomenology of the model is very robust and does not depend on the precise choice of the model parameters unless some clear separation of time scales (e.g., in the kinetic parameters of the different components) is introduced, which can prevent proneural wave progression. This also holds for the Hill coefficient, which sets the steepness of the feedback.
Having said that, we found it interesting to probe the influence of the model parameters in a systematic way, especially with respect to key observables such as the speed of the wave (as mentioned by the reviewer) but also the width of the transition zone. To this end, we have performed a comprehensive sensitivity analysis of the model using the socalled 'Morris method' (Morris, 1991). This method evaluates the impact of the kinetic rate parameters and diffusion constants of the different model components on the output quantities of interest. The results of this new analysis are presented in the main text and Figure 5 of the revised manuscript with technical details deferred to the new Appendix 7. We believe that this additional analysis strengthens the predictive value of the model and can guide the direction of further experiments. We thank the reviewer for this suggestion.
https://doi.org/10.7554/eLife.40919.031Article and author information
Author details
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 FormosaJordan, Claude Desplan, François Schweisguth and members of the Simons and Brand labs for useful discussions.
Senior and Reviewing Editor
 Naama Barkai, Weizmann Institute of Science, Israel
Publication history
 Received: August 9, 2018
 Accepted: February 8, 2019
 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,767
 Page views

 280
 Downloads

 10
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Further reading

 Cancer Biology
 Developmental Biology
Gainoffunction mutations in the proteintyrosine phosphatase SHP2 are the most frequently occurring mutations in sporadic juvenile myelomonocytic leukemia (JMML) and JMMLlike myeloproliferative neoplasm (MPN) associated with Noonan syndrome (NS). Hematopoietic stem and progenitor cells (HSPCs) are the disease propagating cells of JMML. Here, we explored transcriptomes of HSPCs with SHP2 mutations derived from JMML patients and a novel NS zebrafish model. In addition to major NS traits, CRISPR/Cas9 knockin Shp2^{D61G} mutant zebrafish recapitulated a JMMLlike MPN phenotype, including myeloid lineage hyperproliferation, ex vivo growth of myeloid colonies, and in vivo transplantability of HSPCs. Singlecell mRNA sequencing of HSPCs from Shp2^{D61G} zebrafish embryos and bulk sequencing of HSPCs from JMML patients revealed an overlapping inflammatory gene expression pattern. Strikingly, an antiinflammatory agent rescued JMMLlike MPN in Shp2^{D61G} zebrafish embryos. Our results indicate that a common inflammatory response was triggered in the HSPCs from sporadic JMML patients and syndromic NS zebrafish, which potentiated MPN and may represent a future target for JMML therapies.

 Developmental Biology
Zebrafish are an established research organism that has made many contributions to our understanding of vertebrate tissue and organ development, yet there are still significant gaps in our understanding of the genes that regulate gonad development, sex, and reproduction. Unlike the development of many organs, such as the brain and heart that form during the first few days of development, zebrafish gonads do not begin to form until the larval stage (≥5 dpf). Thus, forward genetic screens have identified very few genes required for gonad development. In addition, bulk RNA sequencing studies which identify genes expressed in the gonads do not have the resolution necessary to define minor cell populations that may play significant roles in development and function of these organs. To overcome these limitations, we have used singlecell RNA sequencing to determine the transcriptomes of cells isolated from juvenile zebrafish ovaries. This resulted in the profiles of 10,658 germ cells and 14,431 somatic cells. Our germ cell data represents all developmental stages from germline stem cells to early meiotic oocytes. Our somatic cell data represents all known somatic cell types, including follicle cells, theca cells and ovarian stromal cells. Further analysis revealed an unexpected number of cell subpopulations within these broadly defined cell types. To further define their functional significance, we determined the location of these cell subpopulations within the ovary. Finally, we used gene knockout experiments to determine the roles of foxl2l and wnt9b for oocyte development and sex determination and/or differentiation, respectively. Our results reveal novel insights into zebrafish ovarian development and function and the transcriptome profiles will provide a valuable resource for future studies.