Abstract
During embryogenesis tissue layers undergo morphogenetic flow rearranging and folding into specific shapes. While developmental biology has identified key genes and local cellular processes, global coordination of tissue remodeling at the organ scale remains unclear. Here, we combine in toto lightsheet microscopy of the Drosophila embryo with quantitative analysis and physical modeling to relate cellular flow with the patterns of force generation during the gastrulation process. We find that the complex spatiotemporal flow pattern can be predicted from the measured mesoscale myosin density and anisotropy using a simple, effective viscous model of the tissue, achieving close to 90% accuracy with one time dependent and two constant parameters. Our analysis uncovers the importance of a) spatial modulation of myosin distribution on the scale of the embryo and b) the nonlocality of its effect due to mechanical interaction of cells, demonstrating the need for the global perspective in the study of morphogenetic flow.
https://doi.org/10.7554/eLife.27454.001Introduction
Animal development is characterized by highly dynamic rearrangements of mechanically coupled cells. Such rearrangements must be tightly coordinated across the embryo to achieve normal morphology and organogenesis. During gastrulation of Drosophila, for example, the embryonic blastoderm – an epithelial monolayer of about 6000 cells on the surface of the embryo – undergoes a dramatic deformation that changes tissue topology and gives rise to the three germ layers. These processes involve a coherent flow of cells along the surface of the epithelial monolayer, which in turn drives folding and defines future shape of the embryo. The most prominent aspects of gastrulation are the formation of the ventral furrow which initiates the invagination and internalization of the mesoderm (Martin et al., 2009), and germband extension which involves convergent extension of the lateral ectoderm and the flow of the ventral germband onto the dorsal side of the embryo (Leptin, 1995). Both of these processes have been extensively studied, leading to the identification of developmental patterning genes specifically required for each process (Irvine and Wieschaus, 1994). Live imaging has also uncovered processspecific cell behaviors such as apical constriction of presumptive mesoderm cells during ventral furrow formation (Martin et al., 2009) and intercalation of neighboring cells in the lateral ectoderm during convergent extension (Zallen and Wieschaus, 2004; Bertet et al., 2004). These behaviors are associated with localized activity of the force generating nonmuscle myosin II (Martin et al., 2009; Irvine and Wieschaus, 1994; Zallen and Wieschaus, 2004; Bertet et al., 2004). However, despite considerable understanding of the local processes involved in such cellular rearrangements, a coherent picture of global morphogenetic flows has remained elusive (Butler et al., 2009; Lye et al., 2015; Blanchard et al., 2009).
Understanding how cell flows are coordinated across different cell populations requires distinguishing the roles of local cell behavior and longrange intercellular interactions. To what extent is the transformation of tissue driven locally by the processes associated with cells at that position? How important is the longrange interaction between different regions? In the context of the fly embryo, VF formation seems well explained locally by the apical area contraction of ventral mesoderm cells (Martin et al., 2009). On the other hand, the nonlocal interactions between the VF constriction (or posterior midgut invagination [Collinet et al., 2015]) and the convergent extension of lateral ectoderm remain a subject of active study (Collinet et al., 2015; Rauzi et al., 2015; Rickoll and Counce, 1980) which requires quantitative multiscale analysis.
There are two complementary approaches towards quantitative analysis of tissue flow. One approach focuses on cellscale behavior aiming to decompose tissue flow into specific cellular processes such as cellshape change and intercalation (Etournay et al., 2015; Bosveld et al., 2012) Alternatively one can ‘zoom out’, taking a continuum mechanics approach that aims to describe tissue flow on the whole organ scale (Landau et al., 2012; Prost et al., 2015; Marchetti et al., 2013). This coarsegrained mesoscopic perspective captures correlations in cell behavior which stem from intercellular interactions and the supracellular organization (Martin et al., 2009; Blankenship et al., 2006) of the cytoskeleton in epithelial tissues. In biophysics, the continuum mechanics approach has been developed to understand the behavior of active gels (Prost et al., 2015; Marchetti et al., 2013) during myosin driven viscoelastic flow (Mayer et al., 2010; Behrndt et al., 2012) and has been successfully used to model cortical flows in C. elegans zygotes at the firstcleavage state (Naganathan et al., 2014). Recent theoretical work (Noll et al., 2017) provides a bridge between cellbased and mesoscale continuum descriptions, focusing on the nontrivial consequences of stressdependent active cytoskeletal processes. Here, we shall use continuum mechanics approach to set up a framework for predicting global tissue flow at the whole organ level.
The main advantage of the continuum mechanics approach is its ability to capture key aspects of force balance associated with local deformation and flow. It allows to describe quantitatively, with only a few parameters, how the effect of local forcing spreads across a tissue. The tendency of cells to stick together and resist deformation results in a nonlocal relation between the myosin activity that drives the flow and actual flow velocities. The continuum mechanics approach therefore enables one to test different hypotheses helping to identify key contributing processes. For example, the question of local versus nonlocal response in the continuum mechanics approach translates into specific hypotheses regarded force balance: are the myosingenerated forces balanced locally by traction relative to a substrate or do they propagate within the epithelium layer through cell deformation and viscous coupling? Quantitative analysis can then be used to build, starting from the simplest model, a sequence of approximations that capture biological reality in increasing detail.
We shall describe below a novel image analysisbased approach that will use continuum mechanics to quantitatively relate different observables and will show that myosin distribution and anisotropy on mesoscopic scale is a fully adequate proxy of physical stress, thereby enabling a surprisingly predictive description of global flow. Specifically, we shall show that embryoscale tissue transformations during Drosophila gastrulation are represented by a temporal sequence of three topologically distinct flow field configurations. Each phase is accompanied by a characteristic spatial distribution of myosin molecular motors both on the basal as well as apical cell surface, which we quantify in terms of a coarsegrained ‘myosin tensor’ that captures both myosin concentration and anisotropy. To relate the observed global flow fields to myosin apical and basal distributions we assume that tissue flow is driven by stress proportional to the myosin tensor, and is effectively viscous with two parameters: effective shear and bulk viscosities, the latter controlling the compressible component of the flow. With a total of just three global parameters (only one of them time dependent), this simple model achieves remarkable agreement between predicted and measured spatiotemporal pattern of the flow. The analysis uncovers the importance of a) spatial modulation of myosin distribution and b) the longrange spreading of its effect due to mechanical interaction of cells. In particular, we find that transition to the germband extension phase of the flow is associated with the onset of effective areal incompressibility of the epithelium, which makes the relation of the flow and myosin forcing strongly nonlocal. Our quantitative analysis also reveals a new function for basal myosin in generating a dorsally directed flow and, combined with mutant analysis, points to an unconventional control mechanism of this function through twist dependent reduction of basal myosin levels on the ventral side. Finally, we shall argue that the ability to quantitatively describe the relation between the flow and its myosingenerated forcing provides a new approach to the study of the processes that control morphogenetic transformations, which can be used to disentangle novel control mechanisms such as mechanical feedback from the effects of gene expression patterning.
Results
To enable our study, we generated a pipeline that combines in toto light sheet microscopy (Krzic et al., 2012; Tomer et al., 2012) (Figure 1a), tissue cartography (Heemskerk and Streichan, 2015) (Figure 1b), and segmentationfree anisotropy detection to quantify global tissue flows, and myosin activation patterns (Figure 1—figure supplements 1–4). Using optical flow velocimetry applied to cylinder projections of the Surface of Interest (SOI) passing through cells below the apical cell surface (see SI for details Figure 1—figure supplement 5c), we find that tissue remodeling during Drosophila gastrulation is characterized by three simple flow field configurations (Figure 1c–e). The earliest flows start well before the ventral furrow (VF) forms, and are characterized by a dorsal sink and ventral source (Figure 1c). In contrast to the VF, no cells are internalized during this flow, but rather cells reduce cross section on the dorsal side (Figure 1—figure supplement 1c). As the VF forms, source and sink swap sides and a large group of cells internalize on the ventral side, as mesoderm precursors leave the surface of the blastoderm (Figure 1d). During germband extension (GBE), the flow pattern exhibits two saddles arranged on the dorsal and ventral sides as well as four vortices, two in the posterior and two in the anterior end (Figure 1e). Each of the three flow fields is accompanied by a typical spatial myosin configuration. The preVF flow associates with basal myosin that exhibits a pronounced Dorso Ventral (DV) asymmetry (Warn et al., 1980; Sokac and Wieschaus, 2008; Polyakov et al., 2014), with high levels of myosin on the dorsal and low levels on the ventral side (Figure 1f), while the apical pool appears uniform across the surface (Figure 1i, Figure 1—figure supplement 5d–g). The basal pool remains asymmetric during VF flow (Figure 1g), but the apical pool now also develops DV asymmetry in reversed orientation (Figure 1j). The asymmetry on the apical surface becomes further pronounced in the GBEphase (Figure 1h,k, Figure 1—figure supplements 6 and 7).
Global changes in myosin pools are a hallmark for transitions in flow field configuration (Figure 2a, Figure 1—figure supplements 1, 2, 3 and 4). Myosin is initially enriched in the basal pool, and as sink and source swap position, it begins to accumulate on the apical side. While the basal pool is isotropic (Figure 2—figure supplement 4a), cortical myosin on the apical cell surface is known to polarize during convergent extension (Zallen and Wieschaus, 2004; Bertet et al., 2004). To quantify this effect at the tissue level, we developed an automated segmentationfree anisotropy detection algorithm (Figure 2b, Figure 2—figure supplement 3a,b). Available methods for anisotropy detection mostly operate at the single cell level and construct a nematic tensor by integrating signal intensities along cell outlines (Aigouy et al., 2010). At the organismal scale membrane segmentation is costly, and often fails to define closed outlines of cells using only a polarized membrane marker. We overcome the need for fiduciary markers that increase experimental complexity by shifting the perspective to cell edges and using the Radon transform to implement a robust and rapid segmentationfree algorithm for detecting coursegrained anisotropy (Figure 2b) (Radon, 1917). Radon transforms integrate signal along lines of given orientation and normal distance from the origin. In this way, edges are mapped to peaks that reflect the total intensity along the length of an edge (Figure 2b) (see Appendix 1 for detail). Edge orientation and average myosin intensity are described by a 2 × 2 symmetric matrix (of rank one) defining the local ‘myosin tensor’ (Figure 2b). By averaging the resulting tensors in a given region, we obtain a quantitative description of local tissue anisotropy and overall levels that reflects the intensityweighted average of cell edges. The resulting coursegrained tensor has a nonzero trace, and thus can be separated into an isotropic and a traceless anisotropic part (Figure 3a, Figure 2—figure supplement 3b), from which we construct a measure for anisotropy (which is low in the basal pool Figure 2—figure supplement 4a). The anisotropic signal in the apical pool starts out low, but increases from about 8 min corresponding to late stage 6 (Figure 2a, Figure 2—figure supplement 4b,c). The anisotropy axis, readily computed by the eigenvectors of the myosin tensor, aligns well with local tangent to pair rule gene expression boundaries (Figure 2a,d, Figure 2—figure supplement 4d). This is the expected result given that anisotropies are thought to be driven by the patterned juxtaposition of pairrule gene expression (Zallen and Wieschaus, 2004).
We have examined and quantified tissue flow and myosin distribution in multiple (N = 22) wildtype embryos and found it highly reproducible (Figure 1—figure supplement 5). For the purpose of quantitative analysis presented below, we shall use the (suitably aligned) ‘ensemble’averaged flow and myosin distribution (see SI).
To relate myosin to stress, we assume signal intensity is proportional to myosin motor concentration and its local activity. The latter – pulling on cytoskeletal actin filaments – generates local force dipoles, which can be explicitly described in terms of local stress tensor (see Appendix for details) (Prost et al., 2015; Marchetti et al., 2013). On the coarse grained level, resulting stress would be defined by the activity weighted average over filament orientations and hence proportional to the myosin tensor as we define it. The resulting force per unit area of the epithelial layer is then proportional to the divergence of the myosin tensor (Landau et al., 2012). Note that the isotropic component of the myosin distribution (observed both in the apical and the basal pool) also generates a force that is proportional to the gradient of the measured concentration intensity profile (Figure 1f–k).
To relate myosin generated stress to morphogenetic flow, we assume that on the mesoscopic scale tissue flow is governed by effective viscoelasticity which arises from the mechanical properties of the underlying cytoskeletal network within the two dimensional epithelial layer of cells. This model assumes that on short time scales tissues respond elastically to mechanical perturbations (Bambardekar et al., 2015), yet on longer time scales elastic stress is relaxed through active rearrangement of the cytoskeleton as cells adapt to the imposed deformation. On the longer time scale tissue dynamics can be described by a twodimensional effective viscous flow equation with two effective viscosity parameters that (see Appendix 2) are directly related to the two elastic constants: shear modulus (controlling ‘sliding’ of cells relative to each other) and the planar bulk modulus (controlling areal compression or dilation) (Prost et al., 2015; Marchetti et al., 2013; Martin et al., 1972) (Figure 3a, Figure 3—figure supplement 1b). We note that effective viscosity spreads the impact of local forcing, generating a nonlocal response so the flow at any given point integrates the influence of forces acting all over the embryo. Inverting the equation using the finite element method, we obtain a quantitative prediction for the flow field generated by measured myosin localization patterns (see SI for details). Our model has only three global parameters: the ratio of effective viscosities, and the conversion factors relating normalized apical and basal myosin intensity to stress (Figure 3—figure supplement 1c). To keep the model as simple as possible, we do not allow spatial dependence of these parameters and keep conversion factors constant, leaving the ratio of viscosities as the only timedependent global fitting parameter (Figure 3—figure supplement 1b).
Even without spatial modulation of the parameters, the model achieves about 90% accurate description of the flow pattern before and after VF invagination (see Figure 3b). The main discrepancy of model predictions for preVF flow (see Figure 3c) is a displacement of sink and source positions along the AP axis by less than 30 µm. Prediction of GBE flow essentially agrees with measurements across the entire embryo, with the exception of a domain close to the vortices on the posterior end, due to a mismatch of fixedpoint location (Figure 3d). Remarkably, our model is even able to correctly predict subtle differences between anterior and posterior fixed points along the DV axis (Figure 3d). Measured flow is first dominated by sources and sinks that disappear later during GBE, suggesting that before and during VF invagination cells are less resistant to surface area compression than during GBE. Indeed, quantitative comparison with an independently measured flow field (Figure 1c–e) shows that the ν_{20} /ν_{27} ratio increases dramatically at the start of GBE phase (corresponding to the relative increase of the underlying 2D bulk modulus, see Appendix 2, Figure 3—figure supplement 1b) resulting in effective incompressibility of apical surface of cells. The temporal coincidence between completion of cellularization and increase of the bulk modulus provides an intriguing possible explanation of how the continuous transition in our time dependent variable might be realized. Poor agreement during VF invagination is due to a significant fraction of cells internalizing and thus leaving the surface. To account for this effect, we extend the model to allow a ‘cut’ in the lattice along ventral midline with an imposed inplane boundary force (perpendicular to the cut) representing the pulling effect of the VF (see SI for detail, Figure 3—figure supplement 1a). This relatively simple extension allows to recover ~90% accuracy (Figure 3be), illustrating how regional inhomogeneity associated with particular morphogenetic events could be quantitatively captured by suitable generalizations.
To evaluate the fit obtained in wild type embryos, we examined flows in mutant embryos in which the distribution of myosin is altered. Analysis based on tissue tectonics (Blanchard et al., 2009) has shown that strain rates in twist (twi) embryos, which lack the VF, exhibit slower kinetics compared to WT (Butler et al., 2009), however, the cause of this remains a subject of debate (Butler et al., 2009; Lye et al., 2015; Collinet et al., 2015). We have quantified the flow field and myosin activity patterns in twi mutants (Figure 4—figure supplement 1), and find that our model is able to accurately predict the flow profiles (Figure 4a). During early flow phases – corresponding to times of preVF flow in WT – DV asymmetry of the basal myosin pool is strongly reduced in comparison to WT, as is tissue movement towards the dorsal pole (Figure 4b, Figure 1—figure supplement 1, Figure 4—figure supplement 1a,d). Moreover, anisotropy of the apical myosin pool increases at a slower rate as compared to WT. As previously reported for strain rates (Butler et al., 2009), this is most pronounced for the first 20 min (Figure 4c, compare Figure 4 – fig. supplement = 0 with Figure 2—figure supplement 3e,f). In bcd nos tsl (bnt) embryos lacking all AP patterning, the early basal DV asymmetry is similar to WT, with only slightly reduced myosin asymmetries and dorsal movement (Figure 4b, Figure 2—figure supplement 1). At later stages, however, anisotropy of the apical myosin pool remains low and comparable to preVF WT levels. This result is expected given the uniform expression of pairrule genes in a bnt genetic background (Blankenship et al., 2006) (Figure 4c). Consistent with these myosin distributions, we see the early dorsal flow associated with basal myosin asymmetry but a failure to produce the more complex later flow patterns with their characteristic saddles and vortices. On a quantitative level our model’s predictive power for AP patterning deficient bnt mutant embryos is comparable to WT and twi mutants (Figure 4a).
Discussion
In summary, we have presented a simple biophysical model of morphogenetic flow that quantitatively describes complex tissue motion in terms of a hydrodynamic equation parameterized by two effective viscosities. The flow is driven by the stress defined by a linear superposition of two myosin tensors describing the apical and basal myosin pools. We propose that the basal myosin pool forms an isotropic and contiguous network (He et al., 2016), contracting in a similar fashion as purified actomyosin gels in vitro (Bendix et al., 2008; Alvarado et al., 2013). Imbalance within this network, caused by the twi dependent depletion of myosin on the ventral side, drives global dorsalward flow in the preVF phase, which continues to contribute until early GBE (Figure 4d). Interestingly, in silico perturbations indicate the local depletion (on the ventral side) has a global effect, most evidently manifested by a ‘sink’ on the dorsal side, which is lost in a simulation using the same parameters as WT, but no DV modulation of the basal myosin pool (Figure 4d, left panels). The apical pool decomposes into isotropic and anisotropic components. In addition to previously described accumulation of isotropic myosin in ventral regions (Martin et al., 2009), we observe a striking gradient of anisotropic apical myosin along the DV axis, reaching highest levels in lateral ectoderm and lowest levels in amnioserosa tissue at the dorsal pole. Because the force driving the flow arises from the nonuniformity of the stress, this modulation of myosin distribution is critical for the dynamics. While local rate of cell intercalation is often interpreted in terms of local myosin distribution on cellular and subcellular scales, our model shows that the local rate of strain is a result of the tissuewide distribution of forces generated by the spatial nonuniformity of myosin (mathematically described by the divergence of the myosin tensor) (Figure 4d, right panels). The importance of spatial modulation (Figure 4d) suggests a novel role of the dorsal signaling pathway in generation of GBE flow. Surprisingly, in twi mutants both the rate of increase as well as the peak myosin anisotropy are significantly reduced in the first 20 min of GBE flow (Figure 4c, Figure 4—figure supplement 1). The reduced intercalation and strain rates observed in these mutants has been previously reported (Butler et al., 2009), and interpreted in terms of possible generation of AP forces by the internalized VF (absent in twi mutants). Our model accounts for the reduced rate of strain in terms of the changes in spatial distribution and the reduced level of myosin anisotropy. This however brings up the question of how elimination of twi expression in the ventral mesoderm affects myosin anisotropy in the lateral ectoderm. We suggest that this effect may be due to mechanical feedback on myosin recruitment, which relates the later to local rate of strain. Through this ‘dynamic recruitment’ effect (Noll et al., 2017; FernandezGonzalez et al., 2009), changes in the ventral region that modify global flow patterns can affect myosin distribution and anisotropy in the lateral region. In this way, local modification of the myosin pattern can produce not only a nonlocal perturbation of the flow, but also a nonlocal perturbation of myosin distribution. The global nature of the flow is reinforced by the observed transition towards areal incompressibility at the onset of GBEflow, which together with reduced polarization kinetics and reduced strain rates observed in twi mutants, indicates that nonlocal consequences of stress generated largely in lateral ectoderm can account for the dorsal movement of the posterior midgut.
Taken together, our observations show that morphogenetic flow is a global response to local forcing which arises from the spatial modulation of myosin density and anisotropy. The latter is derived from the spatial patterns of developmental transcription factors, but we suggest may also involve mechanical feedback affecting recruitment of myosin. Our quantitative approach provides a framework for integrating the effect of local factors in the description of the global flow.
Materials and methods
Fly lines used
Request a detailed protocolHis2AvmCherry (Krzic et al., 2012), bcd^{e1}nos^{bn}tsl^{4}/TM3, halo twi^{ID96}/Cyo (twi^{ID96} is also known as twi [Martin et al., 2009]), sqhGFP klar (Martin et al., 2009), OregonR. Embryos where dechorionated following standard procedures, and mounted in agarose gels as previously described (Krzic et al., 2012).
Light sheet microscopy
Request a detailed protocolFluorescencebased live imaging was carried out on a MuVI SPIM (Krzic et al., 2012). Briefly, the optics consisted of two detection and illumination arms. Each detection arm forms a waterdipping epifluorescence microscope, consisting of an objective (Apo LWD 25x, NA 1.1, Nikon Instruments Inc.), a filter wheel (HS1032, Finger Lakes Instrumentation LLC), with emission filters (BLP01488R25, BLP02561R25, Semrock Inc.), tube lens (200 mm, Nikon Instruments Inc.), and an sCMOS camera (Zyla 4.2, Andor Technology plc.), with an effective pixel size of 0.26 µm. Each illumination arm consisted of a waterdipping objective (CFI Plan Fluor 10x, NA 0.3), a tube lens (200 mm, both Nikon Instruments Inc.), a scan lens (S4LFT0061/065, Sill optics GmbH and Co. KG), and a galvanometric scanner (6215 hr, Cambridge Technology Inc.), fed by lasers (06MLD 488 nm, Cobolt AB, and 561LS OBIS 561 nm, Coherent Inc.). Optical sectioning is achieved by translating the sample using a linear piezo stage (P629.1cd with E753 controller) sample rotation is performed with a rotational piezo stage (U628.03 with C867 controller) and a linear actuator (M231.17 with C863 controller, all Physik Instrumente GmbH and Co. KG).
Experiment control and data fusion
Request a detailed protocolStages and cameras are controlled using Micro Manager (Edelstein et al., 2014), to coordinate timelapse experiments, running on a Super Micro 7047GRTF Server, with 12 Core Intel Xeon 2.5 GHz, 64 GB PC3 RAM, and hardware Raid 0 with 7 2.0 TB SATA hard drives. Samples were recorded from two, by 90^{0} rotated views, at a typical optical sectioning of 1 µm, and temporal resolution of 75 s. As previously described (Krzic et al., 2012), MuVI SPIM optical stability allows a fusion strategy based on a diagnostic specimen. Recorded once per experiment, the diagnostic specimen is used to determine an initial guess for an affine transformation, which we feed into a rigid image registration algorithm (Klein et al., 2010), to fuse individual views, resulting in an isotropic resolution of. 26 µm in the registered image.
Surface of Interest extraction
Request a detailed protocolWe used tissue cartography to extract surfaces of interest (SOI) from embryos (Heemskerk and Streichan, 2015). Briefly, we identify the outline of the sample using the Ilastik detector to determine a point cloud for SOI construction. In a fitting step (implemented in the spherelike fitter), we create a smooth description of the SOI in terms of cylinder coordinates defined by AP axis and azimuth (Heemskerk and Streichan, 2015). Image intensity data are then projected onto a nested group of 5 layers (two normally evolved layers above and below the SOI), each three pixels apart, defining a 3–4 µm thick 'curved image stack'. For analysis the maximum intensity projection of nested layers was used. Although the shape of embryos of the same genotype is highly reproducible, small differences in the underlying point cloud can result in small differences of the SOI passing through the apical cell surface. To simplify comparison between embryos, we create a standard projection on a cylinder grid of fixed size, with the embryo surface oriented such that apical is left, posterior right, dorsal in the center, and ventral on top and bottom (Figure 1). Systematic distortions of measurements due to projecting the curved embryo surface to the plane are corrected using the metric tensor (Heemskerk and Streichan, 2015).
The apical surface is static, while the dynamic basal cell surface moves with the cellularization front. Projections of the latter could be created by reading signal on a surface obtained by evolving the apical SOI along its normal basal wards. However, small differences in cell height (<10% of a typical cell height, and ~1% of the embryo diameter), could result in small but systematic bias of the SOI around the cellularization front and impair projection quality. We avoid this problem by determining a new point cloud for each time point, for which we focus the ilastik detector on the interface between basal myosin and yolk. Our model approximates the embryo as a thin shell (see below), and hence as 2D surface. Therefore, we map the dynamic basal cell surface onto the cylinder grid of the static apical SOI.
Particle image velocimetry
Request a detailed protocolWe measure the flow field using the particle image velocimetry (PIV) method, that identifies local displacements between two time points (Adrian, 2005). Briefly, we implemented the phase correlation method that leverages favorable execution times of fast fourier transforms, to estimate local flow in a region of interest on the projections (Kuglin, 1975). To minimize effect from systematic distortions towards polar regions on cylinder projections, we adjust the size of the ROI according to the local metric strain, which we define as the deflation of the metric from flat space (Heemskerk and Streichan, 2015).
Reproducibility of the morphogenetic flow
Request a detailed protocolAlthough gastrulation in Drosophila is highly reproducible from embryo to embryo (Irvine and Wieschaus, 1994), in practical terms experiments are subject to a constant time shift, depending on the developmental stage of the sample at the start of imaging. Thus, we developed an automated routine that allowed us to identify a common time frame that the 36 (WT: N = 22, twist:7, bcd nos tsl:7) liveimaged embryos are registered to. Specifically, we introduced a constant time shift for a given flow field that minimizes the squared difference with respect to a reference flow field averaged over the embryo surface:
where $<{>}_{embryo}$denotes averaging across the embryo, ${\overrightarrow{v}}_{ref}$ is an arbitrary chosen reference from the ensemble, and ${\stackrel{\u20d1}{v}}_{i}$, ${t}_{off,i}$ denote the ith flow field and offset time respectively. In this way, we align samples to a chosen reference, in which we use the first occurrence of the cephalic furrow (CF) as a landmark indicating our choice for $t=0$ min. Within a given genotype, we could automatically determine the offset time. However, to align mutants to WT, we first aligned all mutant datasets, and then used landmarks such as the CF (twist), or the VF (bcd nos tsl), to define a common time frame as best as possible.
Time shifted accordingly, we created an ensemble average flow field for each genotype:
The magnitude of ensemble average is highly reproducible from embryo to embryo (note the small standard deviation Figure 1—figure supplement 5a). Flow trajectories during cellularization point towards the dorsal side (Figure 1—figure supplement 5b), showing persistent movement towards dorsal regions during preCF flow. This is accompanied by reduction of apical cell area in these regions (Figure 1—figure supplement 5c), as measured using confocal microscopy. While the length of preVF flow lines peaks on the anterior and posterior poles in WT, it is substantially reduced near poles in twist, and only mildly reduced in bcd nos tsl (Figure 1—figure supplement 5d–f). Together with the loss of basal DV asymmetry, this suggests that twist mediated reduction of basal myosin levels on the ventral side is responsible for dorsalward flow.
Myosin quantification
Intensity normalization
Request a detailed protocolUsing the imaging and preprocessing procedure as outlined above with samples expressing sqhGFP (Royou et al., 2002), we created projections of the apical and basal cell surfaces, with the goal of establishing a quantitative measurement of global myosin patterns (Figure 1—figure supplements 6 and 7). Ideally, quantification of signal intensities is carried out using identical conditions for each sample in the pool used for statistics, to minimize variability across samples. However, when performing in toto liveimaging, it is difficult to image more than one sample at a time and keep a high recording frequency. To minimize variability in a sequential recording scheme, we keep imaging conditions constant, but there are still possible variabilities in recorded signal intensity for biological but also technical reasons.
To account for such variability between experiments, we normalize recorded data (Figure 1—figure supplement 6b,b’,c). Signal intensity of all time points in a given experiment are summarized by normalizing the intensity distribution: upper and lower range are determined according to the $ll=0$ and $ul=95$percentile; normalization is done by subtracting the $ll$ and dividing by $(ulll)$, yielding a dimensionless normalized signal distribution (compare Figure 1—figure supplement 6b and b’). This strategy should not only allow for comparison on the same microscope, but also across microscopes, allowing for validation of in toto liveimaging from sequential experiments against fixed batches imaged e.g. on a confocal.
Basal myosin pool analysis via light sheet microscopy
Request a detailed protocolHere, we briefly outline the results for the DV asymmetry in the basal myosin pool that we reported in the main text. First, we time align intensity normalized basal projections as described above. Next we convolve each pullback with a Gaussian of width $\sigma ~3$ cell diameters to obtain basal myosin at the mesoscale (see discussion in model section below for definition). The results are then ensemble averaged to obtain ensemble myosin distribution as shown in Figure 1—figure supplement 6. To assess DV asymmetry, we focus on the region outlined by the black/white dashed line, where we first take an average along the AP axis and then compute average signal on dorsal side, and subtract from it the average signal on the ventral side. Repeating the outlined routine for all time points, we obtain the plot show in main text Figure 4b.
Basal myosin pool via confocal microscopy
Request a detailed protocolFigure 1—figure supplement 7c shows DV cross sections of fixed embryos stained for rb anti zipper and mouse anti dorsal, cut along the AP axis, and imaged on a confocal microscope. DV orientation of the samples is automatically determined based on the dorsal signal. To estimate the age of fixed embryos in relation to liveimaging data, we constructed a calibration curve for cell apicobasal height shown in Figure 1—figure supplement 7a. Using the known monotonic relation between cell height and age (Merrill et al., 1988), which we find lasts until about $8min$ after CF formation, we obtain estimate for the age of a fixed embryo based on measuring cell height. By segmenting the outline of basal myosin, we can then measure DV asymmetry in the same way as described for liveimaging data above.
A direct comparison between live imagingbased DV asymmetry measurement, and based on $N=345$ fixed DV cross sections from confocal shows that after applying normalization routine as described above, we find similar estimates for the DV asymmetry using both light sheet and confocal imaging (see Figure 1—figure supplement 7b). Note that uncertainties in the calibration curve propagate to exact age determination in the embryo, and thus increased fluctuations in DV asymmetry determined using confocal imaging.
Finite element implementation
Request a detailed protocolInversion of the continuum equation of state relating the coarsegrained myosin tensor and cellular flowfield was achieved using Finite Element Methods (FEM) in the weak formulation implemented within the FELICITY toolbox for MATLAB (Walker, 2017). Equations were inverted on a static triangular mesh representing the 'canonical' embryo surface produced via a point cloud (described above) subsequently turned into a smooth triangulation using MeshLab (Cignoni, 2008). As such, all objects within the equation of state had to be parameterized within the 3D embedding space of the mesh, which can be done by using the direction relation between projections and SOI (Heemskerk and Streichan, 2015). The only dynamic input to this inversion algorithm is the divergence of the myosin tensor. This was computed by interpolating the gradient of each Cartesian component of the tensor onto triangular faces of our mesh, producing a 3 × 3 × 3 object on each face. The partial trace of this object over directions within the tangent plane of the face result in the estimated divergence. This operation was repeated for all myosin pools used. All equations within the FEM software are projected onto the surface of our 3D mesh to manually ensure solutions only exist within the tangent plane.
To benchmark the quality of our FEM solver, we implemented the equation of motion on a sphere and tested results for known solutions, which confirmed our solver works within the expected numerical accuracy. To test dependence on discretization used, we compared our results using different meshes at varying mesh sizes, and found good agreement with all test cases.
In order to model internalization of the VF, we introduced the ability to add a 'cut' within the triangular mesh. Contraction of the tissue was modeled by manually introducing local force dipoles on edges within the mesh pointing along the bond. All vertices along the cut were given zero bulk modulus to allow for local tissue compression needed to simulate invagination. The location of the cut was estimated from the PIV flow fields. Specifically the ratio of the divergence of the velocity field to the velocity field's magnitude was used to estimate the spatial extent of the cut over times during ventral furrow formation.
Predictions for the flow field obtained via inversion are subject to an overall scale factor, that can't be determined by the model. To compare ensemble averaged flow field measurement $\stackrel{\u20d1}{v}\left(t\right)$ to model predictions $\stackrel{\u20d1}{u}\left(t\right)$in a quantitative fashion, we define a global measure for the spatial residual that is insensitive to such a scale factor. With the short hand notation $<\overrightarrow{u}>:=\sqrt{<\overrightarrow{u}{\left(x\right)}^{2}{>}_{embryo}}$ to define overall magnitude of the field u across the surface of the embryo ($<\overrightarrow{u}{\left(x\right)}^{2}{>}_{embryo}$ denotes averaging the space dependent field ${\stackrel{\u20d1}{u}\left(x\right)}^{2}$ across the embryo surface, so is not space dependent.), the residual is defined as
provides a spatial discrepancy map, indicating the prediction quality as a function of location on the embryo, that is insensitive to noise dominated fluctuations in domains of no flow (i.e. fixed points), as opposed to e.g. inner product.
Appendix 1
Segmentation free anisotropy detection using Radon transforms
When viewed in crosssection, the cell cortex in epithelia forms a polygon tiling, where the outline of each cell is well approximated by a closed sequence of edges. Anisotropic distribution of proteins is often characterized by homogeneously increased accumulation to cell edges of particular orientation, while it remains homogeneously low and comparable to background on other edges (Zallen and Wieschaus, 2004; Bertet et al., 2004; Blankenship et al., 2006). Note that typically the number of edges at low and high signal accumulation is of the same order of magnitude. Available methods to quantify cortical anisotropy mostly operate at the single cell level and construct a nematic tensor by integrating signal intensities along cell outlines (Aigouy et al., 2010). At the organismal scale membrane segmentation is costly, and for polarized markers low signal to noise on a significant number of edges often results in difficulties to close the cell circumference. We overcome the need for fiduciary markers that increase experiment complexity, by shifting perspective to cell edges and designing a robust and rapid segmentation free anisotropy detection algorithm.
Let's consider an image on a rectangular domain $\mathrm{\Omega}$, showing an edge of possibly nonuniform intensity, assigning each pixel at coordinates $\left(x,y\right)\in {\mathbb{R}}^{2}$ some intensity $I\left(x,y\right)\in \mathbb{R}$  as shown e.g. in Figure 2—figure supplement 1. Lets call this edge a linear signal or just signal. For simplicity, and without loss of generality, we assume the average intensity along the linear signal is $<{I}_{ls}>=1,$and average background intensity is $<{I}_{bg}>=0$. While under favorable conditions with high signal to noise ratio, it may be possible to identify the linear signal using conventional methods such as edge detection using a typical threshold, this task will be substantially more error prone at low signal to noise ratio.
Therefore, we decided to reformulate the problem using Radon transforms (Radon, 1917), that integrate signal along lines of a given angle and signed distance from the origin. Consider a line $L$ with normal ${\stackrel{\u20d1}{n}}_{L}=(\mathrm{cos}\alpha ,\mathrm{sin}\alpha )$, a signed normal distance $\delta \in \left\{\infty ,\infty \right\}$ away from the origin (e.g. shown in pink Figure 2—figure supplement 1, top left), which may be parameterized as follows
where $t\in ${0,1} is the parameter that marks position along the line $L$ (Figure 2—figure supplement 1). It is clear that all possible pairs of normal orientation and signed distance $(\alpha ,\delta )$, describes all possible lines covering the plane. The Radon transform of the image $I$  denoted $\mathcal{\mathcal{R}}I$  establishes a map from the rectangular domain $\mathrm{\Omega}$, on which the image is defined into the space spanned by line orientation and signed distance, where each line $L$ characterized by the pair $(\alpha ,\delta )$ is assigned the integrated(summed) intensity projection along the line  or in mathematical terms:
where the latter equation uses the explicit parametrization of the line $L$ in terms of $\alpha ,\delta ,t$ defined above. (The following is a verbal example for the equation above, and may be skipped): Figure 2—figure supplement 1 illustrates the algorithm that constructs the radon transform for the linear signal shown in white, one that has a normal orientation of ${45}^{0}$ and constant intensity ${I}_{ls}=1$ over a background of intensity ${I}_{bg}=1$. The magenta dashed line indicates a line $L$ with normal orientation $\alpha $ and a distance $\delta =0$ away from the origin shown in red. Integrating image intensity along such a line, for any orientation $\alpha \ne {45}^{0}$, we theoretically obtain only a contribution from the intersection with the linear signal, thus $\mathcal{\mathcal{R}}I(\alpha \ne {45}^{0},0)=1$. In contrast for $\alpha ={45}^{0}$, i.e. when the line is parallel and on top of the linear signal, the radon transform integrates the intensity along the entire line, returning the length $l$ of the linear signal: $\mathcal{\mathcal{R}}I(\alpha ={45}^{0},0)=l$. In our example, for $\delta \ne 0$ we will either have a single intersection of the line $L$ with the signal, again returning with the same outcome as above, or no intersection, where the radon transform returns 0, indicated by the dark blue regions in Figure 2—figure supplement 1.
Since the radon transform is linear, it follows that for any image that is the linear superposition of linear signals, the radon transform is the linear superposition of the radon transform of each linear signal. For example, if an image consisted of 2 linear signals (for example a four fold vertex), the resulting radon transform is the sum of the radon transform of each linear signal, with peak levels at orientation and signed distance of each linear signal. In this way, linear signals are mapped to peaks in the radon transform that reflect the total intensity along the length of the linear signal. Knowing the location of peaks in the radontransformed space allows us to reconstruct the position and orientation of each linear signal in the image. This simplifies the task of identifying the linear signal, as detection of peaks due to their compact structure is simpler than detection of edges, and the radon transformed signal will be less susceptible to fluctuations in the underlying data, and therefore enhance robustness.
To obtain average image intensity from the radon transform we use the fact that the height of the peak reflects the total intensity ${I}_{tot}$ along the linear signal to construct the average signal as $<I>={I}_{tot}/l$, where $l$ denotes the length of the linear signal. To determine the length of the linear signal, we need to determine the position where the linear signal begins and/or ends in the image. This is done by interpolating the signal along the orientation and position corresponding to the identified peak in the radon transform, and identifying discontinuities in the interpolated signal, which mark the boundaries of the linear signal. With beginning and end determined, we obtain the length as the magnitude of the vector connecting the two points. If beginning and end don't fall onto the boundary of the image, we perform an integration of a line restricted between these two endpoints, to obtain a more accurate estimate of the total line intensity.
Validation of segmentation free anisotropy detection
Pullbacks created with ImSAnE based on light sheet microscopy data of developing D. melanogaster embryos show the cortex of roughly 6000 cells, which would yield a complex radon transform pattern, effectively impeding peak detection. To simplify peak detection, we decided to carry out analysis in small blocks, where we focus on a local region of interest (ROI) that we sweep across the entire image. The size of the local region is chosen such that it can accommodate between 1–3 typical edge lengths ${l}_{edge}$. In this way, we typically obtain a small number (<10) of edges in the ROI, such that we could use the extendedmaxima transform (Soille, 2013) to identify the typically well separated peaks. To avoid double counting of edges  as they may appear from sweeping the ROI  we average detected lines of locally similar angles with a small angular difference $\Vert {\delta}_{\alpha}\Vert <\text{}\u03f5$ and distance $\delta <\text{}{l}_{edge}$.
Following the general prescription given above, we obtain local estimates for position, orientation, and average intensities of cell edges, demonstrated in Figure 2—figure supplement 2. A zoom on a local region indicated by the white box shows that this routine  fully automatically, and without need of fine tuning parameters  faithfully detects the orientation of cell edges, including dim edges with only minimally stronger signal compared to background. To further benchmark the quality of intensity estimates, we turned to known anisotropy of myosin during germband extension (Blankenship et al., 2006). Figure 2—figure supplement 1 shows on the right, that the presented algorithm finds the normalized intensity on AP edges (${90}^{0}$) is systematically higher than on DV edges (${0}^{0}$), reaching is maximum approximately 20 min after cephallic furrow formation, in good agreement with previous manual measurements (Blankenship et al., 2006).
Appendix 2
Mathematical model of tissue flow
Mechanics of epithelial tissue is largely defined by the properties of intracellular cytoskeletal cortexes, linked by cadherin mediated adherens junctions into a global transcellular mechanical network. This mechanical network is ‘active’ in the sense that it involves myosin motors that generate internal forces and can do work by contracting actinmyosin filaments under tension. A cytoskeletal network is also ‘adaptive’ in the sense that it can relax mechanical stress by reorganizing internally, on subcellular scale by recruiting (or releasing) myosin and other key molecular components, and on tissue scale, by allowing cells to rearrange locally (by T1 processes). The latter process allows cells to change neighbors and ‘flow’ relative to each other, while preserving the integrity of the epithelium and its cytoskeletal network. The detailed description of these complex cellular processes is not necessary for our present goals, which call merely for an approximate mesoscale description of tissue mechanics. By mesoscale we mean the scale of >3 cell diameters, still much smaller than the macroscopic scale of the embryo tissue flow. Microscopic complexity notwithstanding, on mesoscale we will think of tissue as a continuous medium and it will suffice for us to capture the facts that i) on short time scales tissue responds elastically to mechanical perturbations (Bambardekar et al., 2015) and ii) on the longer time scales elastic stress is relaxed through active rearrangement of the cytoskeleton as cells adapt to the imposed deformation. This is enough to define flow velocity in response to external and internal stress: on the timescale comparable to internal stress relaxation, tissue dynamics can be described by a generic viscoelasticity equation with two effective viscosity parameters which we shall now derive.
Shortterm elastic response means that incremental increase of strain causes an increase in stress, which can subsequently relax with the relaxation time ${\tau}_{R}$. To simplify the derivation we will present now, lets assume a flat twodimensional surface. Elastic stress is then governed by the Maxwell viscoelasticity and in the mesoscopic continuum approximation we have:
where $\stackrel{.}{x}=\frac{d}{dt}x$, so that $\stackrel{.}{u}}_{a$ is the rate of local displacement in the direction of spatial component a. Spatial derivatives of the $\stackrel{.}{u}}_{a$ vector define local rate of strain $\mathrm{\partial}}_{a}{\stackrel{.}{u}}_{b}+{\mathrm{\partial}}_{b}{\stackrel{.}{u}}_{a$. In addition, ${\delta}_{ab}$ is the Kronecker delta matrix and we have adopted the Einstein convention of summing over repeated indices (e.g. ${\partial}_{c}{u}_{c}={\sum}_{c}{\partial}_{c}{u}_{c}$). The first two terms on the right hand side describe generation of stress in proportion to the rate of strain ($\mu ,\lambda $ are the Lame coefficients parameterizing an elastic stressstrain relation) and the last term parameterizes relaxation of stress.
On the other hand, in the continuum approximation, tissue flow velocity ${{v}_{a}=\dot{u}}_{a}$ can be described by the (compressible) Stokes equation
where ${F}_{a}$ is external force (per unit area), $\rho $ is the density and ${\upsilon}_{0}$dynamic viscosity.
Now suppose that stress relaxation is sufficiently rapid to achieve quasiequilibrium
substituting ${\sigma}_{ab}$ into Equation 2 we have
Note that transient elasticity, parameterized by ${\tau}_{R}$, $\mu $ and $\lambda $, generates ‘effective viscosity’ ${\tau}_{R}\mu $, which can dominate the microscopic viscosity ${\upsilon}_{0}$ of the ‘fluid’ itself. This effective viscosity is in general anisotropic with ${\tau}_{R}\mu $ contributing to the shear viscosity which resists shear, acting to make flow more uniform, and the bulk viscosity ${\tau}_{R}\left(\mu +\lambda \right)$ which acts specifically on the on the compressible component of the flow.
If the effective viscosity is sufficiently high and flow changes sufficiently slowly, the inertial term can be neglected ($\rho {\stackrel{.}{v}}_{a}\ll$${{\upsilon \partial}_{b}}^{2}{v}_{a}$) and the quasistationary flow is defined by the balance of the external forcing and effective viscosity
where ${\upsilon}_{1}={(\upsilon}_{0}+{\tau}_{R}\mu )$ and ${\upsilon}_{2}={\tau}_{R}\left(\mu +\lambda \right)$. Since external force is related to external stress via ${F}_{a}={{{\partial}_{b}\sigma}_{ab}}^{ext}$ and we have argued that relevant ‘external stress’ that drives tissue flow in the embryo is proportional to the myosin tensor ${{\sigma}_{ab}}^{ext}~{m}_{ab}$ we have arrived at the equation we have used in the main text to relate tissue flow velocity with myosin distribution. (Note that the relation of ${{\sigma}_{ab}}^{ext}$ and ${m}_{ab}$ is due to the fact that myosin generates contractile forces, so that a local maximum of an isotropic myosin distribution ${m}_{ab}\left(r\right)={\delta}_{ab}m\left(r\right)$ would act just like a low pressure region, generating inward directed flow.)
Passive versus active mechanics
In the section above we put forward a simple effect model based on conventional Maxwell viscoelasticity. Epithelial tissues however are clearly not ordinary passive viscoelastic materials. Myosindriven rearrangement of the cytoskeleton is an active and adaptive process, which can be regulated on cellular and subcellular scales. We argue, however, that average flow on the mesoscopic scale can be usefully approximated by the passive viscoelastic model with suitable effective parameters. The ability of this simple model to capture highly nontrivial spatial flow patterns without the need for introducing spatial dependence of model parameters proves the validity of the approximation. Still, we do not expect this ‘passive viscoelasticity’ model to be complete. We anticipate that a more detailed description of myosin activity and myosin recruitment dynamics would be required in order to describe both the fast cellular scale fluctuations and the longterm myosin dynamics on the scale of the embryo. Other molecular/genetic factors will come into play as well: e.g. cadherin and other cytoskeletal components on small scale and transcription factors that guide morphogenesis on large scales. We plan to address these issues in the future.
References

1
Twenty years of particle image velocimetryExperiments in Fluids 39:159–169.https://doi.org/10.1007/s0034800509917
 2
 3
 4
 5

6
A quantitative analysis of contractility in active cytoskeletal protein networksBiophysical Journal 94:3126–3136.https://doi.org/10.1529/biophysj.107.117960
 7
 8
 9
 10
 11

12
Eurographics Italian Chapter ConferenceMeshLab: an Opensource mesh processing tool, Eurographics Italian Chapter Conference.

13
Local and tissuescale forces drive oriented junction growth during tissue extensionNature Cell Biology 17:1247–1258.https://doi.org/10.1038/ncb3226

14
Advanced methods of microscope control using μManager softwareJournal of Biological Methods 1:10.https://doi.org/10.14440/jbm.2014.36
 15

16
Myosin II dynamics are regulated by tension in intercalating cellsDevelopmental Cell 17:736–743.https://doi.org/10.1016/j.devcel.2009.09.003
 17

18
Tissue cartography: compressing bioimage data by dimensional reductionNature Methods 12:1139–1142.https://doi.org/10.1038/nmeth.3648

19
Cell intercalation during Drosophila germband extension and its regulation by pairrule segmentation genesDevelopment 120:827–841.

20
elastix: a toolbox for intensitybased medical image registrationIEEE Transactions on Medical Imaging 29:196–205.https://doi.org/10.1109/TMI.2009.2035616

21
Multiview lightsheet microscope for rapid in toto imagingNature Methods 9:730–733.https://doi.org/10.1038/nmeth.2064

22
Proceedings of the International Conference on Cybernetics and SocietyThe phase correlation image alignment method, Proceedings of the International Conference on Cybernetics and Society.
 23

24
Drosophila gastrulation: from pattern formation to morphogenesisAnnual Review of Cell and Developmental Biology 11:189–212.https://doi.org/10.1146/annurev.cb.11.110195.001201
 25

26
Hydrodynamics of soft active matterReviews of Modern Physics 85:1143–1189.https://doi.org/10.1103/RevModPhys.85.1143

27
Unified hydrodynamic theory for crystals, liquid crystals, and normal fluidsPhysical Review A 6:2401–2420.https://doi.org/10.1103/PhysRevA.6.2401
 28
 29

30
Requirements for autosomal gene activity during precellular stages of Drosophila melanogasterDevelopment 104:495–509.
 31
 32
 33
 34
 35

36
Uber die Bestimmung von Funktionen durch ihre Integralwerte langs gewissez MannigfaltigheitenBerSächsiche Akademie Der Wissenschaften 69:262–277.

37
Embryoscale tissue mechanics during Drosophila gastrulation movementsNature Communications 6:8677.https://doi.org/10.1038/ncomms9677

38
Morphogenesis in the embryo ofDrosophila melanogaster  Germ band extensionWilhelm Roux's Archives of Developmental Biology 188:163–177.https://doi.org/10.1007/BF00849045
 39
 40
 41
 42

43
Finite ELement Implementation and Computational Interface Tool for YouFELICITY.

44
Changes in the distribution of cortical myosin during the cellularization of the Drosophila embryoJournal of Embryology and Experimental Morphology 57:167–176.
 45
Decision letter

Frank JülicherReviewing Editor; MaxPlanckInstitute for the Physics of Complex Systems, Germany
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 "Global morphogenetic flow is accurately predicted by the spatial distribution of myosin motors" for consideration by eLife. Your article has been reviewed by three peer reviewers, including a member of our Board of Reviewing Editors, and the evaluation has been overseen by the Reviewing Editor and Anna Akhmanova as the Senior Editor.
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
Summary:
The authors perform a careful and thorough analysis of the collective cell movements during early development of the Drosophila embryo. Using light sheet microscopy and a quantitative image analysis the observed patterns of fluorescence intensity are mapped on a closed twodimensional surface representing the epithelial surface. In the case of fluorescently labelled myosin, patterns of myosin on the apical and on the basal side of the tissue have been quantified separately on the 2d surface. The work uses an elegant and original approach to quantify anisotropies of apical myosin using a Radon transform. In addition, flow fields are determined using PIV. The authors use the 2d distribution and anisotropies of apical and basal myosin to estimate patterns of active stress in the tissue. Using a continuum theory of an active gel, velocity fields are calculated and compared to the observed cell flow patterns. Strikingly, the authors show that the most important features of the cell flow can be accounted for from the knowledge of myosin distributions. Furthermore, the authors show that altered flow patterns in twi mutants can be explained by changed patterns of myosin distributions under mutant conditions. These are important results that provide a significant advance of the understanding of morphogenetic movements in the early fly embryo. The quantitative methods used are powerful and sophisticated. While this work is interesting and strong, there are several points that need to be addressed in a revision of the manuscript and several points need clarification.
Essential revisions:
1) A major concern as detailed below is that a lot of information is missing in the manuscript.
2) The authors do not show what the actual myosin distribution looks like, and there is little detail in how the 'coarsegrained' or 'smoothened' myosin distributions in Figure 1FK are obtained. The corresponding supplemental figure did not really help in this regard. If the data has been smoothened to produce these plots, the original data should also be shown. Scale bars should also be added to the figures. What is the 'contrast' in the myosin signals? What is the level of intensity between dorsal and ventral in relation to the average? What are the units in Figure 1F–K, Figure 2C, etc.?
3) The model used to calculate the flow fields is shown in Figure 3A. It is unclear why the isotropic basal and the anisotropic apical signals are used only. The apical myosin also has an isotropic signal. Is this not used and why not? Maybe the notation, which is not fully explained, leads to confusion here. Is the tensor m_{apical} traceless or does it contain an isotropic signal? This should be clarified. Presumably the isotropic and anisotropic contributions of apical myosin are contributing with two different factors. These factors should be defined and their values reported.
4) The quantification of anisotropies of the myosin distribution represents an interesting methodological advance. The presentation of the quantification of this anisotropy on the surface of the embryo is not very clear. Figure 2C and Figure 2D do not provide a full picture of the quantified anisotropy patterns. The authors should show a complete map of myosin anisotropy at different times.
5) Information about the anisotropy of the basal signal evaluated with the same method is missing. Why is the basal myosin anisotropy not contributing to the stress equation? Is it because basal myosin anisotropy is low?
6) In the appendix the anisotropy quantification using a Radon transform is explained. However, I did not see a clear definition of how the Radon transform information is used to determine the myosin tensor m_{apical} in Figure 3A.
7) If m(r) denotes a factor positively correlated to myosin concentration, the sign of the righthand side of Equation 6 may be wrong. With the current sign, a positive accumulation of myosin would drive an outward flow. It is essential that the authors clarify this.
8) In Figure 2D it is not clear what is plotted. The caption states "Eigenvectors of myosin tensor are plotted in cyan". It is not clear if the myosin tensor is a field defined on a regular grid, of only selected vectors are shown. Shown are bars that are hardly visible to the bare eye. They do not look like vectors and it is unclear to which points they are attached.
9) The authors should provide more information on the 2D shear and bulk viscosities ν and ν' and they ratio B. If one considers the height of the tissue layer as an explicit variable and treats the tissue as isotropic and incompressible in 3D, divergences in the 2D flow field couple to height changes and the ratio of the 2D bulk to the 2D shear viscosity would be 3 (see e.g. Batchelor, 2000). How does this compare to the values that are fitted here? Can the model be simplified to 2 parameters? Can B be constant in time (I assume not, Figure 1C–E)?
10) Related comment: Please provide the values of the parameters B, α, β in a table/plot in the main text together with confidence bounds. Are α and β constant in time? The 'true' number of model parameters remains unclear. Is this really just one parameter as stated in the abstract, or is it three as stated in the main text? It remains unclear whether the viscosity ratio B was different for each time point which would correspond to effectively many fit parameters. In the main text it is mentioned that B changes with time however this is unclear and needs to be clarified. The main text refers to Figure 1C–E about this but in the figure there is no information about B.
11) It seems the authors use a single value of B that is spaceindependent at each time point, but then allow B to instantaneously change everywhere and synchronously in space. The motivation for this remains unclear. How could that be achieved in the embryo?
12) In Figure 1C–E it is not stated in the caption (and not even in the supplement) which fluorescence signal (apical, basal myosin, combined?) is used to calculate the velocity field. This important information that should be found easily.
13) In Figure 3C, it is unclear whether the shown flow profiles are obtained with or without the "cut" in mesoderm? We would suggest adding a figure or schematic actually showing what the "cut in mesoderm" perturbation is doing.
14) One sentence states that the comparison between theoretical and experimental plots is done using normalised velocity fields rather than real velocity fields (caption of Figure 3B: "Fit residual, comparing predicted flow field with measured flow field normalized for magnitude"). The precise way the velocities were normalized remained unclear to the referees. Was the velocity normalized at each position separately or was only the overall magnitude of the velocity field normalized? The former case would be a problem because real velocity fields should be compared including velocity amplitude patters.
15) An important result is that "the model achieves about 90% accurate description of the flow pattern before and after VF invagination". It is quite unclear how this percentage of agreement is defined and determined. Also, in the caption of Figure 3 the residual is not defined. The definition of the residual given in the methods remains a bit unclear because the notation is not well explained (which quantities depend on position and which do not, how is <[…]>_{embryo} defined etc.).
16) The flow and myosin profiles in twist embryos are not reported. They should absolutely be shown.
17) It is not clear what information the theoretical analysis of mutants is bringing. Are the parameters used to fit the data the same as in WT?
18) It is misleading to refer to the model used in the calculations as viscoelastic. This is a description in the fluid limit and it should be referred to as such. Viscoelasticity does not enter; the corresponding viscoelastic relaxation time is not a parameter that is considered. Please correct corresponding references in the Abstract and main text, for example the language in the Results section (please add page and maybe line numbers in the revision) is confusing and the statement in the Abstract is misleading (surprisingly simple effective viscoelasticity model > viscous model). Of course, one can arrive at a fluidlike description by considering the longtime behavior of a viscoelastic material as the authors also show, and it is certainly ok to provide this derivation. But the theory that is used in the end is fluid and not being clear here can put the reader on a wrong track.
19) In Figure 3A a coordinate independent, covariant expression of the dynamic equations is presented. The equations given in the appendix are not covariant but depend on a coordinate system. It is therefore unclear whether the continuum theory applied to the problem is covariant or not. This should be clarified. Furthermore, the finite element method used is also not clear with respect to this point and not sufficiently explained. Do the results depend on the choice of the finite element discretization used or not? If not, has this been checked? If yes, what does this mean? We understand that the authors do not want to go into notational complexities here but the conceptual approach used should be clear.
https://doi.org/10.7554/eLife.27454.023Author response
Essential revisions:
1) A major concern as detailed below is that a lot of information is missing in the manuscript.
We thank the referees for the careful reading of the manuscript and suggested improvements. Below, we discuss the individual points of concern which have been addressed largely through modified, and additional Supplementary figures.
The main changes are:
Figure 1—figure supplement 1–4,
Old Figure 1—figure supplement 1–3 are now Figure 1—figure supplement 5–7,
2) The authors do not show what the actual myosin distribution looks like, and there is little detail in how the 'coarsegrained' or 'smoothened' myosin distributions in Figure 1F–K are obtained. The corresponding supplemental figure did not really help in this regard. If the data has been smoothened to produce these plots, the original data should also be shown. Scale bars should also be added to the figures. What is the 'contrast' in the myosin signals? What is the level of intensity between dorsal and ventral in relation to the average? What are the units in Figure 1F–K, Figure 2C, etc.?
We have added Figure 1—figure supplement 1, to demonstrate the initial steps of the imaging data analysis pipeline, consisting of defining embryo shape, apical, basal, and lateral surfaces. We explicitly define the “midplane” section, cutting across the lateral surfaces of cells, that was used as a common reference frame, and the apical and basal surfaces related by normal displacement. The individual surfaces are colorcoded and original signal on each surface shown. The midplane section surface was used to evaluate the flow field using the PIV method (as described in the Materials and methods section), which we indicated by adding a flow field next to it. To show the time course of myosin signal on apical and basal surfaces respectively, we added Figure 1—figure supplement 2, Figure 1—figure supplement 3 and Figure 1—figure supplement 4. We also provide measures for contrast, defined as root mean square contrast, shown in Figure 1—figure supplement 2A.
Original myosin signal is measured in arbitrary units. In the Material and methods section, we describe how we turn this into a dimensionless number by normalization. Therefore, all plots showing the myosin tensor have no units.
To better visualize the normalization we added panels B and B’ to Figure 1—figure supplement 2. Also, we added panel C, to report requested statistics on dorsal and lateral regions respectively.
We added scalebars to all projections of original embryo surfaces.”Pullbacks” (onto the cylinder) are locally distorted and thus a single scalebar would be misleading.
To explain the coarsegraining procedure, we added Figure 2—figure supplement 3, which in panel (A) shows outlines of cells, and the coarsegrained lattice. In panel (B) we show how the myosin tensor is defined on the mesoscale, which involves smoothing with a Gaussian kernel. All relevant quantities are indicated in the panel (see below comments).
3) The model used to calculate the flow fields is shown in Figure 3A. It is unclear why the isotropic basal and the anisotropic apical signals are used only. The apical myosin also has an isotropic signal. Is this not used and why not? Maybe the notation, which is not fully explained, leads to confusion here. Is the tensor m_{apical} traceless or does it contain an isotropic signal? This should be clarified. Presumably the isotropic and anisotropic contributions of apical myosin are contributing with two different factors. These factors should be defined and their values reported.
We thank the refees for pointing out the potentially confusing notation. Now the annotation for the model shown in Figure 3A shows that the isotropic basal pool is represented by the scalar field ‘m_{basal}’, (we explain below why we neglected anisotropic contributions). The apical pool is represented as the full tensor ‘m_{apical}’, which, as we explain in the text, is not traceless and has anisotropic and isotropic contributions, as shown in Figure 3A. To this end we explain in the annotation in Figure 3A how the tensor decomposes into an isotropic scalar (i.e. the trace), and an anisotropic traceless part. In this way there is no need for an additional factor distinguishing between isotropic and anisotropic contributions in the apical pool. The point is emphasized in the revised Results section:
“By averaging the resulting tensors in a given region, we obtain a quantitative description of local tissue anisotropy and overall levels that reflects the intensityweighted average of cell edges. The resulting tensor generally has a trace, and thus can be separated into an isotropic and a traceless anisotropic part
(Figure 3A).”
4) The quantification of anisotropies of the myosin distribution represents an interesting methodological advance. The presentation of the quantification of this anisotropy on the surface of the embryo is not very clear. Figure 2C and Figure 2D do not provide a full picture of the quantified anisotropy patterns. The authors should show a complete map of myosin anisotropy at different times.
We thank the referees for these suggestions and added Figure 2—figure supplement 3, where we explain how the coarsegrained myosin tensor on the mesoscale is obtained. Panel (B) emphasizes the distinction between isotropic and anisotropic parts, to provide a clear connection to the quantitative measure of anisotropy in basal and apical pools shown in panels (C and D). We further provide a sequence of representative timepoints for the apical surface, showing that myosin anisotropy is strongest in the germband and peaks around 17 minutes. For better intuition, we also added plots on the apical surface of the embryo in panel (E). For a complete picture of Figure 2D, we demonstrate a time series of the principal axes of the anisotropic myosin tensor.
5) Information about the anisotropy of the basal signal evaluated with the same method is missing. Why is the basal myosin anisotropy not contributing to the stress equation? Is it because basal myosin anisotropy is low?
As correctly guessed by the referees, we neglect basal pool anisotropy in our study, because it is small. We added Figure 2—figure supplement 3C to demonstrate this.
6) In the appendix the anisotropy quantification using a Radon transform is explained. However, I did not see a clear definition of how the Radon transform information is used to determine the myosin tensor m_apical in Figure 3A.
Figure 2B outlines the concept of how we describe edge orientation r(i), based on the detected angle α using the radon transform. Below it we show how the myosin tensor m_{ab} on each edge is computed. To make the connection with the model clearer, we changed the notation, to follow the convention in Figure 3A, and indicated that we compute the tensor for each edge (i) but only on the apical surface by specifically referring to m_{apical(i)}.
To clarify how we compute the myosin tensor on the mesoscale, we added Figure 2—figure supplement 3, which demonstrates the conversion from nematic m(i) obtained by Radon transform on cell edges to a regular lattice covering the embryo surface. Panel (A) shows a sketch of a cell array in green, and points in the regular lattice covering the embryo in black. From this we compute distance between a given edge i, and lattice site x, which is used to construct a Gaussian weight in the definition of the mesoscale myosin tensor, explained in panel (B).
7) If m(r) denotes a factor positively correlated to myosin concentration, the sign of the righthand side of Equation 6 may be wrong. With the current sign, a positive accumulation of myosin would drive an outward flow. It is essential that the authors clarify this.
Of course, the referee is correct, in that a positive accumulation of myosin should create an inward flow. We have corrected the misprint in the equation.
8) In Figure 2D it is not clear what is plotted. The caption states "Eigenvectors of myosin tensor are plotted in cyan.". It is not clear if the myosin tensor is a field defined on a regular grid, of only selected vectors are shown. Shown are bars that are hardly visible to the bare eye. They do not look like vectors and it is unclear to which points they are attached.
We added Figure 2—figure supplement 3 that explains how the myosin tensor field is defined. In the caption to Figure 2D we now explain that we show the axis parallel to the leading eigenvectors of the tensor field, which we only plot along even skipped stripes for simplicity of comparison.
9) The authors should provide more information on the 2D shear and bulk viscosities ν and ν' and they ratio B. If one considers the height of the tissue layer as an explicit variable and treats the tissue as isotropic and incompressible in 3D, divergences in the 2D flow field couple to height changes and the ratio of the 2D bulk to the 2D shear viscosity would be 3 (see e.g. Batchelor, 2000). How does this compare to the values that are fitted here? Can the model be simplified to 2 parameters? Can B be constant in time (I assume not, Figure 1C–E)?
We have reworded the Discussion section in the text to clarify that viscosities in question are “effective” quantities that relate to the shear and bulk moduli and the rate of stress relaxation (Figure 3A, more fully explained in Appendix 2) rather than derived from the isotropic and incompressible 3D fluid as suggested by the referee. As explained in the response to comment 10 and 11 below, B(t) (we updated the notation, to reflect which variables depend on space or time) is time dependent. However, it is low through much of cellularization and only increases significantly from the onset of cephalic furrow formation (T= 0 min) (see Figure 3—figure supplement 1B). Cells already completed much of their elongation before (compare Figure 1—figure supplement 3A). As outlined in the response to comment 11 below, we instead speculate the steady increase of incompressibility relates to completion of cellularization, which occurs concomitant with the apparent increase of the bulk modulus.
10) Related comment: Please provide the values of the parameters B, α, β in a table/plot in the main text together with confidence bounds. Are α and β constant in time? The 'true' number of model parameters remains unclear. Is this really just one parameter as stated in the abstract, or is it three as stated in the main text? It remains unclear whether the viscosity ratio B was different for each time point which would correspond to effectively many fit parameters. In the main text it is mentioned that B changes with time however this is unclear and needs to be clarified. The main text refers to Figure 1C–E about this but in the figure there is no information about B.
We thank the referees for pointing out the inconsistency between our original Abstract and the main text. We corrected the statement that there is “only one parameter”, to now say there is” only one time dependent and two constant parameters”. Indeed, B is a number that changes with time. In the added Figure 3—figure supplement 1C, we elaborate on the number of model parameters, and explain that they are not spatially dependent. We stress that the complex spatial structure of the flow at any given time is accurately fitted with only one parameter, as the timeindependent constants are determined by other timeslices. We also show how B depends on time, and provide values including confidence intervals of the constant conversion factors α and β.
11) It seems the authors use a single value of B that is spaceindependent at each time point, but then allow B to instantaneously change everywhere and synchronously in space. The motivation for this remains unclear. How could that be achieved in the embryo?
Updated Figure 3A annotation now indicates that apical and basal myosin is a field, while B(t) is a spaceindependent but time dependent scalar.
To improve presentation of this aspect in our model, we added Figure 3—figure supplement 1. Panel (B) shows the smooth temporal dependence of B, starting out low, and increasing slowly with cellularization. After gastrulation is completed, B is constant high. In the text, we also added the following statement:
“The temporal coincidence between completion of cellularization and increase of the bulk modulus provides an intriguing possible explanation of how the continuous transition in our time dependent variable might be realized.”
12) In Figure 1C–E it is not stated in the caption (and not even in the supplement) which fluorescence signal (apical, basal myosin, combined?) is used to calculate the velocity field. This important information that should be found easily.
We have added a figure, Figure 1—figure supplement 1, explaining the imaging pipeline used to extract the flow field, shown in panel (I) and define apical and basal myosin distributions
13) In Figure 3C, it is unclear whether the shown flow profiles are obtained with or without the "cut" in mesoderm? We would suggest adding a figure or schematic actually showing what the "cut in mesoderm" perturbation is doing.
Following this suggestion and added a schematic, showing the cut in the mesoderm in the new Figure 3—figure supplement 1A. We also explain in the caption that flow profile shown in Figure 3E was produced using the ‘cut’.
14) One sentence states that the comparison between theoretical and experimental plots is done using normalised velocity fields rather than real velocity fields (caption of Figure 3B: "Fit residual, comparing predicted flow field with measured flow field normalized for magnitude"). The precise way the velocities were normalized remained unclear to the referees. Was the velocity normalized at each position separately or was only the overall magnitude of the velocity field normalized? The former case would be a problem because real velocity fields should be compared including velocity amplitude patters.
The caption of Figure 3B now reads:
“Fit residual, comparing predicted flow field with measured flow field (see subsection “Finite element implementation” for a detailed definition of the residual) as a function of time. Both fields are normalized for average magnitude. The average magnitude of predicted velocity field defines one of our fitting parameters.”
We also clarified our notation in the residual definition, to avoid potential
confusion (see response to comment 15).
15) An important result is that "the model achieves about 90% accurate description of the flow pattern before and after VF invagination". It is quite unclear how this percentage of agreement is defined and determined. Also, in the caption of Figure 3 the residual is not defined. The definition of the residual given in the methods remains a bit unclear because the notation is not well explained (which quantities depend on position and which do not, how is <[...]>_{embryo} defined etc.).
The caption of Figure 3 points towards the Materials and methods section for more detail. To make this link clearer, we now emphasize the residual is defined in subsection “Finite element implementation”.
In subsection “Finite element implementation”, we improve notation interpretation, explain the averaging procedure across the embryo, and emphasize which quantities in the residual depend on space. The SI text now defines:
With the short hand notation $<\overrightarrow{u}>:=\sqrt{<\overrightarrow{u}{\left(x\right)}^{2}{>}_{embryo}}$ to define overall magnitude of the field u across the surface of the embryo $<\overrightarrow{u}{\left(x\right)}^{2}{>}_{embryo}$ denotes averaging the space dependent field ${\stackrel{\u20d1}{u}\left(x\right)}^{2}$ across the embryo surface, so is not space dependent.), the residual is defined as
16) The flow and myosin profiles in twist embryos are not reported. They should absolutely be shown.
We have added Figure 4—figure supplement 1 to demonstrate representative time points of the flow field, basal, apical isotropic, and apical anisotropic myosin in twist mutants.
17) It is not clear what information the theoretical analysis of mutants is bringing. Are the parameters used to fit the data the same as in WT?
We thank the referees for pointing out the missing description. The caption of Figure 4D now states that model parameters are as previously determined for WT, and we emphasize in the Discussion section:
“Interestingly insilico perturbations indicate the local depletion (on the ventral side) has a global effect, most evidently manifested by a “sink” on the dorsal side, which is lost in a simulation using the same parameters as WT, but no DV modulation of the basal myosin pool (Figure 4D, left panels).”
18) It is misleading to refer to the model used in the calculations as viscoelastic. This is a description in the fluid limit and it should be referred to as such. Viscoelasticity does not enter; the corresponding viscoelastic relaxation time is not a parameter that is considered. Please correct corresponding references in the Abstract and main text, for example the language in the Results section (please add page and maybe line numbers in the revision) is confusing and the statement in the Abstract is misleading (surprisingly simple effective viscoelasticity model > viscous model). Of course, one can arrive at a fluidlike description by considering the longtime behavior of a viscoelastic material as the authors also show, and it is certainly ok to provide this derivation. But the theory that is used in the end is fluid and not being clear here can put the reader on a wrong track.
The reviewer is correct: while microscopic description of the tissue may be more complex, on the time scales of interest, cellular flow is effectively viscous. Yet, we feel that it is important to note that “effectively viscous behavior” does not mean that the tissue is a simple viscous fluid, which it is not. (Note point #9 illustrates possible misunderstandings along this line of thought.) Our reference to the underlying viscoelastic model was intended to communicate that point and our derivation in Appendix 2: Mathematical model of tissue flow”relates the longitudinal and transverse effective viscosities to the shear and bulk elastic moduli and the stress relaxation rate. We have however made corrections suggested by the reviewer.
19) In Figure 3A a coordinate independent, covariant expression of the dynamic equations is presented. The equations given in the appendix are not covariant but depend on a coordinate system. It is therefore unclear whether the continuum theory applied to the problem is covariant or not. This should be clarified. Furthermore, the finite element method used is also not clear with respect to this point and not sufficiently explained. Do the results depend on the choice of the finite element discretization used or not? If not, has this been checked? If yes, what does this mean? We understand that the authors do not want to go into notational complexities here but the conceptual approach used should be clear.
The equation presented in Figure 3A shows a general expression of the continuum model in terms of Laplace operator, divergence and gradient, which is valid on flat but also curved surfaces. Of course, these operators need to be correctly implemented (on a curved surface), which in our case is accomplished computationally via finite element method on a triangulated surface. For the derivation of our equation in the Appendix 2, however we assumed for simplicity a flat configuration, which greatly simplifies notation, and helps the reader to develop a better intuition. In the Materials and methods section we now also explain how we tested accuracy of our finite element solver and dependence on the discretization used. In subsection “Finite element implementation”, we added:
“To benchmark the quality of our FEM solver, we implemented the equation of motion on a sphere and tested results for known solutions, which confirmed our solver works within the expected numerical accuracy. To test dependence on discretization used, we tested our results using different meshes at varying mesh sizes and found good agreement with all test cases.”
https://doi.org/10.7554/eLife.27454.024Article and author information
Author details
Funding
National Science Foundation (PHY1220616)
 Boris I Shraiman
Howard Hughes Medical Institute
 Eric F Wieschaus
National Institutes of Health (NICHD 1K99HD088708)
 Sebastian J Streichan
Gordon and Betty Moore Foundation (GBMF #2919)
 Boris I Shraiman
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Trudi Schüpbach as well as members of Shraiman and Wieschaus labs for helpful discussions; Lars Hufnagel for initial support with the MuVi SPIM setup; and Reba Samanta for preparing histological samples used to create Figure 1—figure supplement 7. This work was funded by the Gordon and Betty Moore Foundation grant #2919 (SJS), NSF PHY1220616 (BIS), grant #1K99HD08870801 from National Institute of Child Health and Human Development (SJS). EFW is an investigator with the Howard Hughes Medical Institute.
Reviewing Editor
 Frank Jülicher, MaxPlanckInstitute for the Physics of Complex Systems, Germany
Publication history
 Received: April 5, 2017
 Accepted: January 20, 2018
 Accepted Manuscript published: February 9, 2018 (version 1)
 Version of Record published: March 8, 2018 (version 2)
Copyright
© 2018, Streichan 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

 3,768
 Page views

 695
 Downloads

 35
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, Scopus, PubMed Central.