A mechanochemical model recapitulates distinct vertebrate gastrulation modes

During vertebrate gastrulation, an embryo transforms from a layer of epithelial cells into a multilayered gastrula. This process requires the coordinated movements of hundreds to tens of thousands of cells, depending on the organism. In the chick embryo, patterns of actomyosin cables spanning several cells drive coordinated tissue flows. Here, we derive a minimal theoretical framework that couples actomyosin activity to global tissue flows. Our model predicts the onset and development of gastrulation flows in normal and experimentally perturbed chick embryos, mimicking different gastrulation modes as an active stress instability. Varying initial conditions and a parameter associated with active cell ingression, our model recapitulates distinct vertebrate gastrulation morphologies, consistent with recently published experiments in the chick embryo. Altogether, our results show how changes in the patterning of critical cell behaviors associated with different force-generating mechanisms contribute to distinct vertebrate gastrulation modes via a self-organizing mechanochemical process.


Mathematical Model
We derive our model in Cartesian coordinates and then provide an equivalent reformulation in polar coordinates.

Momentum Balance and Continuity Equation
We denote the tissue velocity field v(x, t) = [u(x, t), v(x, t)] and its vorticity scalar by ω(x, t).We recall the velocity gradient decomposition ∇v(x, t) = S(x, t) + W(x, t), (S1) with the rate-of-strain tensor S and the spin tensor W defined as We further decompose S into its isotropic and deviatoric parts as The viscous, active, total stresses and total stress trace are σ V = − pI + 2µS s σ A =m(x, t)[e(x, t) ⊗ e(x, t) − 1 2 I], e(x, t) = [cos φ(x, t), sin φ(x, t)] where µ denotes the shear viscosity, m denotes the active stress intensity arising from active myosin and φ the orientation of the myosin cables with respect to the x−axis.The friction exerted by the surrounding fluid on the epiblast is assumed to be very small, therefore we ignore substrate friction in our model.
Because viscosity is high (µ ≈ 9 × 10 3 P a • s), we ignore inertial terms and write the momentum balance as (S5) To close the system, we choose a simple continuity equation in which the flow divergence is proportional to the isotropic total stress (eq.(S4)) and active myosin concentration i.e.
where c has units [ P a • s] −1 and p 0 is a non-dimensional parameter modulating the contribution of active myosin to flow divergence.The intuition behind eq. ( S6) is that regions of high isotropic compressive stress and high active myosin, inducing isotropic apical cell contraction, exhibit negative flow divergence, i.e., cell ingression or folding.
The above relation then serves as a means of writing the local pressure in terms of the flow divergence and the active myosin concentration.Equation (S6) can also be derived considering the continuity equation of a two-dimensional confluent epithelium and using the above phenomenological arguments to express the source/sink terms as a function of ∇ • v and m.Substituting eq.(S6) into eq.(S5), we obtain the resulting force balance equations describing our active continuum without pressure terms (S7)

Active Stress Orientation
We assume the evolution of actomyosin cables orientation satisfies Jeffery's equation ( 42) which describes the direction evolution of axis-symmetric fibers under a low Reynold's number flow.In eq.(S8), φ t = ∂φ ∂t , λ is a parameter quantifying the sensitivity of actomyosin cables reorientation due to shear-induced rotations, S ij denotes the entries of S and the first two terms represent the advection term and the rigid-body rotation rate described by half of the flow vorticity.

Active Stress Magnitude
The dynamics of active myosin on actin filaments is characterized by catch bonds (23), i.e., the dissociation rate of actomyosin bonds depends on the tension along the actomyosin cable T c , and decreases exponentially for increasing tension.Laser ablation experiments demonstrate that supracellular actomyosin cables exhibit higher tension and higher levels of myosin than isolated cell junctions (52).We explain below our reasoning to estimate T c .

Junction-level modeling
Junction-level modeling of the feedback between total tension and active stress intensity supports active stress intensity as a good proxy for cable tension.Optical manipulation of apical cell-cell junctions in drosophila and chick embryo suggest that junctions are viscoelastic (53,54).Among simple viscoelastic models, measurements of junction deformation through time in chick embryos are best fit by a simple Maxwell model consisting of a spring and dashpot in series (54) (Figure S1A).Thus, we consider each junction as a Maxwell element in parallel with an active element generating additional active stresses (Figure S1B).Maxwellian viscoelasticity has been frequently used to model single-junction dynamics in drosophila (35,37,53).It is easy to show that the tension of a Maxwell element can be written as an elastic term with a continuously relaxing rest length (55).Consistent with the assumption of Maxwellian viscoelasticity and additional tension induced by myosin activity, the total junction tension is where K is a spring constant, r ij = r j − r i is the length between vertex positions, l ij is a dynamic rest length at which elastic tension vanishes, and m ij represents an active stress intensity due to myosin activity.The rest length relaxes to the junction length over time as where τ is the relaxation time.We model the dynamics of myosin activity (active stress intensity) in terms of mechanisms for myosin recruitment and tension-dependent dissociation where t r is the recruitment time scale of myosin from the cytoplasm, which has a given maximum concentration m n , and t −1 d e −koT ij represents the tension-dependent dissociation rate of actomyosin cables.χ has units [t −1 ] and accounts for the saturation of myosin activity on junctions, ensuring that active stress cannot accumulate without bound.Physically, available myosin may be depleted from basal pools or the number of active myosin motors that can operate on actin filaments may be limited.In this work, we remain agnostic as to which saturation mechanisms may operate and their relative contributions, so we capture saturation effects with a simple linear term.Later, we show that our results are not sensitive to the saturation mechanism.Note again that m, m n have units [P a] as they describe the stress magnitude arising from active myosin.From the time when active myosin is available to the time it generates active stresses, there is a time scale (t p ).We model this delay with the ratio (t p /t c ), where t c is the characteristic time scale of the system.
Equations (S9-S11) model the coupled dynamics of tension and myosin activity in a single junction (Figure S1B).To couple the dynamics of adjacent edges and model a supracellular actomyosin cable, we further assume relaxational vertex dynamics.Restricting the cable to one dimension, the dynamics of each vertex are dominated by its two neighboring cable edges as where ν represents effective friction resisting vertex displacement (37,35,53).The effective friction models the internal dissipative mechanisms such as the interaction of the cell junction with the cytosol and the viscous flow induced in the cytoskeleton due to the motion of the cell vertices (28).External dissipation by the surrounding fluid is assumed to be minimal.Thus, the length r ij between vertices experiences a strain rate where ∆ is the discrete Laplacian.In the continuum limit, we replace the discrete laplacian with ∆ l2 , where l is the characteristic junction length (37).This provides two coupled equations for the dynamics of tension and myosin activity: (S14) These equations can be combined to give a differential equation for the difference between T and m which can also be understood as a differential equation for the elastic (or passive) component of the tension.The first term causes elastic tension to propagate along cables.The second term causes elastic tension to diminish as the rest length relaxes.Where the propagation term dominates, we expect transient deviations between m and T .Otherwise, the dynamics tend to reach T = m.Simulating equations (S14) with an initial central peak of active stress intensity, we find that, when K l2 τ ν is small, T and m follow one another closely, as shown in Movie S13 and Figure S1C.
Supplementary Figure S1: Mechanical feedback along a 1D cable.A) Chick embryo experimental data from (54) shows that a junction deformed through optical manipulation exhibits a response that is best fit by Maxwellian viscoelasticity.B) Mechanical element for an actomyosin cable, consisting of a Maxwell element in parallel with an active element.C) Plots from the evolution of (S14) for sample parameters K l2 ν = 1, τ = 0.001, and χ

Tissue-level modeling and time scales
The tissue-scale behavior is fluid-like and hence well described by viscous and active stresses (eq.(S5)) at gastrulation time scales.Extending the above reasoning to the tissue scale, active myosin contributes to cable tension T c through an active stress contribution T A = e, σ A e = m 2 .At the junction level, T A corresponds to m ij , leaving K(r ij − l ij ) to be accounted for.One can also compute from the viscous stress T V = e, σ V e .However, σ V and T V reflect fluid-like tissue stresses arising from cell division and intercalation.If and how T V may affect junction-level cable tension is unclear.Regardless, we expect this effect to be of secondary importance because myosin-induced (active) tension (T A ) dominates T V (37).
We are not aware of estimates for K or ν, so we cannot precisely estimate K l2 ν .Moreover, the junction-level modeling presumes propagation over a 1D cable with no constraints on its extent.In the tissue, however, super-cellular cables have a limited extent spanning from 2 to 8 cells (5).Therefore, even if K and ν were known, propagation strength should be reduced by the effective length of super-cellular cables, which is small compared to the tissue scale.Last, the junction relaxation time has been estimated to be in the order of seconds in the chick epiblast (54).Altogether, the fast relaxation of the junction rest length and a small K l2 ν suggests that the elastic component of the tension is negligible at leading order.Following brief periods of rest length relaxation, mechanical feedback adjusts myosin activity to balance the tension of neighboring junctions (35).This is consistent with experimental observations that levels of myosin activity can serve as a readout for tension on cell junctions (2).Thus, in our coarsegrained continuum model, which neglects short-timescale elasticity and considers longer timescales over which the epithelium is effectively viscous, we expect and that m will exhibit some effective propagation of active stress intensity along cables in their predominant orientation φ.In two dimensions, the effects of this directional propagation can be captured by two directional derivatives along the cable, ∆ φ = D φ D φ , and with a coefficient ξ modulating the tissue-scale effect of tension propagation.The definition ∆ φ = D φ D φ is a simplification of the more general (e • ∇) 2 and assumes that the individual cables along which propagation occurs are approximately straight locally.
We reason that at the tissue scale, the dynamics of active stress intensity are governed by similar mechanisms of myosin recruitment and tension-dependent dissociation modeled at the junction scale.Using the above assumptions and eq.(S16), we model the evolution of active stress magnitude as The last two terms account for the advection of active stress intensity through material transport of actomyosin cables and for the directed propagation of active stress intensity along actomyosin cables with ∆ φ defined as the second directional derivative along the cable direction φ.

Model in Cartesian coordinates
The full model describing our active continuum can be summarized as

Non-dimensional model in Cartesian coordinates
Denoting by u c , x c , t c = x c /u c and m c = µu c /x c the characteristic velocity, characteristic length scale, characteristic time scale and the characteristic viscous shear stress, we rewrite eq.(S18) in non-dimensional compact form as p 0 describes the cell's propensity of ingressing from a myosin-induced apical isotropic contraction and can be thought of also as the ratio of isotropic to anisotropic active stresses.p 1 describes the ratio of the shear to bulk viscosity, p 2 is the alignment parameter that describes how the actomyosin cables tend to orient due to shear-induced rotations and is ≈ 1 for elongated fibers.p 3 represents the ratio of the myosin (or active stress intensity) recruitment time scale t r , from the ambient myosin pool of concentration m n /m c , to the characteristic t p time scale, multiplied by m n /m c .p 4 represents the ratio of the t p time scale to the dissociation time scale t d .p 5 is the product of the characteristic viscous shear stress and k o , which represents the inverse of the characteristic tension of the actomyosin cables (with units [P a] −1 ).p 6 is the strength of the simplest linear saturation term set by the ratio of t p to the characteristic dissociation time scale χ −1 controlling saturation.Practically, p 6 enables setting a desired saturation level for m. p 7 is the ratio between transport of the active stress magnitude via tension propagation (with effective diffusion constant ξ) and active stress magnitude transport via advection.
The first two terms of the force balance are a viscous force and a force arising from compressibility inhomogeneities.Instead, the last three terms are active forces induced by m and φ inhomogeneities.Finally, observing the first two terms in the last equation of (S19), we note that more myosin induces higher tension on the cables, which in turn decreases the dissociation rate: (m ↑) → T c ↑ (more positive) → diss.rate ↓ → myosin increases (m ↑) (Figure 2A).This observation suggests an instability, which will become clear in section 1.4.3.

Role of p 0
To gain insights on p 0 , we rewrite the force balance as From eq. (S6), we know that p 0 ≥ 0 reflects the expectation that high myosin induces isotropic apical cell contraction that causes active cell ingression.To elucidate the effect of p 0 on the relation between ∇[∇ • v] and ∇m, we neglect the shear force (−2p 1 ∆v), assume a simplified uniform ridge of active myosin in the y direction (Figure S2), and a spatially uniform cable distribution (i.e., ∇ • B = 0).Therefore, eq.(S21) simplifies to If cables are perpendicular to the myosin ridge (Figure S2A), as in the wild-type development (Figure 1E

One-dimensional model
To gain insights about the model (S19), we assume slow dynamics in the y−direction, ignore the φ−dynamics and set φ = 0, obtaining the following simplified one-dimensional model ) whose boundary conditions are inherited from the nonlinear problem, and correspond to zero Dirichlet boundary conditions for δu and zero flux boundary conditions for δm.We discretize the linearized PDE using a centered finite differencing scheme (56) and compute the generalized eigenvalue problem associated with eq.(S24).Supplementary Figure S3A top shows three equilibria m eq., while the bottom panel shows the corresponding 500 eigenvalues, with the largest real part, of the linearized operator.The spectrum shows that the lowest and highest equilibria are linearly stable while the intermediate equilibrium is linearly unstable for the chosen parameter values.We confirm these linearized results by solving the full nonlinear PDE from an initial Gaussian m distribution (Figure S3B) added to the unstable m eq.(Figure S3A top, red asterisk).From a biological perspective, a region of higher m 0 represents the mesendoderm precursor.At later times (Figure S3C), m undergoes a focusing instability, that in turn generates a 1D attractor corresponding to the PS and marked by a peak of negative divergence (red dashed curve).The complete time evolution is available as Movie S1.Finally, the velocity divergence u x shows a contracting region close to the PS as well as two symmetric expanding regions close to the boundary, consistent with the isotropic deformations observed in the Embryonic and Extra-Embryonic territories (2).Our simple 1D model reproduces several key features of amniote gastrulation.Yet, it is insufficient to reproduce the convergent extension and vortical patterns observed in experiments, which we investigate with our 2D model.
We note that eq.(S23) differs from the one-dimensional isotropic active gel instability in (57).Despite both associated with the dynamics of actomyosin on the cell cortex, eq.(S23) models myosin exchange with the cytoplasm as opposed to (57), which assume no exchange with the cytoplasm (eq.( 1) and [8] in (57)).Additionally, our model accounts for catch-bonds between actin and myosin in the form experimentally measured in (23), which is absent in (57).Finally, (57) accounts for friction with the substrate, absent in our case.For a more detailed biological explanation behind the myosin exchange with the cytoplasm and the catch-bond mechanism in our context, we refer the reader to Figure 3 of (5).

Model in polar coordinates
Because of the circular geometry of the Extra-embryonic tissue and the symmetric radial cellular motion at its outer boundary, we write model (S18) in polar coordinates.Supplementary Figure S4: Numerical mesh.Illustration of the discretization grid (first quadrant) used to solve eq.(S26).β denotes the cables orientation relative to α.
Denoting the velocity field by v = u(r, α, t)e r + v(r, α, t)e α and the cables orientation relative to α by β (Figure S4), we obtain

Numerical Scheme
We develop a finite-difference numerical scheme to solve eq.(S26) in MATLAB.First, we multiply the momentum balance equations by r 2 , leaving their solutions unaltered while avoiding numerical issues at grid points close to the origin.Then, we discretize spatial differential operators using centered finite-difference at second-order accuracy, except at the grid points (r 1 , α j ) (Figure S4), where we use forward finite-difference at second-order accuracy (56).
We solve the first two equations for u, v by inverting (once) the corresponding discretized differential operator, which we multiply at each time step to the discretized forcing vector arising from the updated β, m-dependent terms.We adopt a two-step Adams-Bashforth method to solve the last two PDEs for β and α.

Repeat step 1 and 2 for all
Outputs: u, v, β, m at all grid points {r i , α j , t k }.
To prevent the development of sharp derivatives in β at the origin due to the forward finite-difference scheme, we add a diffusion term in the β t equation with a non-dimensional prefactor of 10 −3 .In all our simulations, we use ∆r = 0.05, L = 1, ∆α = 4 • , ∆t = 5 × 10 −6 , T f ≈ 1.We monitor the accuracy of our numerical solver by substituting our solution at each time step into the momentum balance equations.We find deviations from zero of the order ≈ 10 −11 .

Induction of active stress intensity via tension propagation
To visualize the effect of directed propagation (∆ φ ) of active tension in polar coordinates, we solve the m t equation in (S26) activating only the diffusion term on the right-end side.Figure S5A shows the evolution of active stress intensity m starting with high values in the sickle-shaped region of mesendoderm precursors.At later times, active stress intensity propagates around the epiblast due to the tension-propagation mechanism captured by p 7 ∆ φ m (Figure S5B).In the full model, however, all terms in eq.(S17) contribute to the evolution of m.Because the tissue is fluidized, the propagation of active stress intensity is dominated by advection so that active stress intensity does not appear to spread to the full extent depicted in Figure S5.

Parameters, boundary and initial conditions: selection and sensitivity analysis
In all our simulations, we use a single set of parameter values summarized in table 1 We have selected parameter values using a combination of experiments, mechanistic arguments and numerical simulations; and have performed sensitivity analyses with respect to parameter changes as described below.
4.1 p 0 , p 1 , p 2 According to our arguments in Section 1.4.2,we set p 0 = 0 when active cell ingression is inhibited with drugs and p 0 = 2 otherwise.We use p 1 = 0.15 and note that increasing p 1 models a more compressible tissue.With higher p 1 , the typical counter-rotating vortices (polonaise movement) disappear due to increased cell ingression and less cell redistribution (see Fig. S11).From Jeffery's equation (S8), p 2 modulates the rotation rate due to shear of a material fiber, which is ≈ 1 for elongated fibers.Here we set p 2 = 0.9, and note that slight variations of this value does not affect our results, as shown in WT simulations Movie S14 (Movie S15) for p 2 = 0.6 (p 2 = 1).4.2 p 3 , p 4 , p 5 t p is the approximate time required to observe active stress at the tissue scale generated by active myosin cables.This phenomenon involves a complex cascade of operations from kinetic reactions to single-cell and multicellular mechanochemical processes.Here we estimate t p using the typical time for a junction to fully contract and induce an intercalation event driven by junctional myosin, a robust and experimentally accessible readout of tissue-scale active stress induced by cables.Experiments in the chick epiblast show that this process takes up to ≈ 20min (2).A typical time course is shown in fig 4F of (2) where it takes 60 minutes for three junctions to successively contract.See also supplementary videos 10 and 11 of (2) which show various examples of intercalating cells, where it takes 4-10 frames (1 frame is 3 minutes) for junctions to contract.This result is consistent with tissuescale mechanochemical models, where the actomyosin or active stress dynamics is the slowest (37,35,58), and a given configuration of active stress generates the corresponding tissue-scale velocity field by solving a static force balance (eq. 1 in S19).The recruitment and dissociation time scales are t r t d ≈ 10s in line with observations made in Drosophila (53).With these estimates, we set p 3 = 50 and p 4 = 100.p 5 is the coupling strength regulating myosin intensity and tension via the catch-bond feedback mechanisms.We are unaware of techniques to estimate p 5 , and set p 5 = 1.25.To test the sensitivity of our results with respect to p 5 , we show that a 20% reduction of the sensitivity to actomyosin dissociation rate to tension, induces a higher unstable myosin equilibrium concentration m eq.(Figure S6).The overall dynamic evolution of the system, however, remains qualitatively unchanged as shown in Movie S16, to be compared with Figure 2 and Movie S2.By contrast, a 25% increase of p 5 causes a bifurcation of m eq., in which case model evolution differs considerably from experiments, as shown in Movie S17.Similar conclusions hold if changes in p 3 and p 4 induce a bifurcation of m eq. .Our results, therefore, show that the active stress instability modulated by the catch-bond dynamics of actomyosin cables is essential for reproducing observed flow patterns.

p 6 , p 7
While myosin activity cannot increase without bounds, the precise mechanisms by which and degree to which it saturates on junctions are beyond the scope of our model.This is because we are interested in a finite-time developmental interval characterized by instability rather than an asymptotic equilibrium state.To this end, we add a linear saturation term −p 6 m as the simplest functional form to dampen positive feedback and create a higher stable equilibrium value for m (Figure S3A and Figure S7A).For our model parameters, the three m eq.persist for p 6 < 5.989, preserving the instability needed for active stress intensity to grow.To show the robustness of our results with respect to changes in p 6 , we show that even without (p 6 = 0) any saturation term (Figure S7B), the qualitative results of the 2D model do not change (Movie S18).
More generally, we find that slight variations in p 3 , p 4 , p 5 , p 6 do not alter the overall system behavior as long as the graph (Figure S3A) determining m eq does not change considerably (or undergoes a bifurcation).
Because the tissue is fluid-like and cells undergo embryo-scale motion, we argue that the advective transport of cells with active myosin dominates the induction of myosin via tension propagation on individual cables.Hence, we set p 7 = 0.01 throughout.Increasing the strength of active force propagation by 500% (p 7 = 0.05), the qualitative results of the 2D model (Movie S19) do not change, confirming the robustness of our results to parameters selection.

Boundary conditions
At the EP-EE boundary, the velocity is tangential to the boundary, ensuring that the EP area remains approximately constant over time, while in the EE region, the velocity is radial and roughly symmetric in α, explaining the expansion of the EE region over time (2).Because here we want to mechanistically predict what generates the observed flow patterns rather than enforcing it as an input, we do not select the boundary of our domain at the EP-EE interface, where the Dirichlet velocity boundary conditions would consist of a highly specific tangential flow profile.By contrast, we set our domain boundary ∂Ω (Figure 1A) slightly outside the EP-EE boundary, where, from the previous argument, the velocity is only radial and slightly higher than zero (we set |v b | = 0.2).Remarkably, our model recapitulates the correct vortical flows within the EP region and the associated tangential velocity behavior at the EP-EE interface.The most natural choice of boundary conditions for m, φ is no flux, as there is no experimental evidence of sources, sinks, or constant values of m, φ close to the EP-EE boundary.Mathematically, we set where n is the outer normal to ∂Ω.
To test the sensitivity of our results with respect to |v b |, we find that a 50% increase (decrease) of |v b | does not affect the qualitative behavior of our results, as shown in Movie S20 (Movie S21) to be compared with Movie S2.
Overall, increasing (decreasing) |v b | exerts a stronger (weaker) pulling of the PS towards the boundary.

Initial conditions
Consistent with experiments (Figure 1E at HH1), we set the initial orientation of actomyosin cables φ(x, t 0 ) along the tangential direction, and m(x, t 0 ) as a Gaussian distribution (Figure 2C left) added to the unstable m eq.equilibrium (Figure S3A), modeling the onset of high active myosin in the mesendoderm sickle territory.The Gaussian m(x, t 0 ) = m eq.+ Ae is centered at r c = 0.6, α c = 250 • , has amplitude A = 5 and standard deviations σ r = 0.06, σ α = 3σ r /r c .We select α c = 250 • to align the AP axis of the model (Figure 2C) with the AP axis from the experiment (Figure 2E).Changing α c only changes the AP axis orientation.We set σ α = 3σ r /r c so that the Gaussian width in the tangential direction is three times the width in the radial direction, reflecting the asymmetry of the mesendoderm sickle.
To test the sensitivity of our results to our Gaussian perturbation, we find that a 20% increase (decrease) in A Movie S22 (Movie S23) does not affect the qualitative behavior of our model.Similar conclusions hold for a 20% increase (decrease) in σ r Movie S24 (Movie S25) and for 20% increase (decrease) in r c Movie S26 (Movie S27).Increasing r c causes a stronger pulling of the PS towards the boundary.

Summary
In all our simulations of (S26) (i.e.eq. 1 in the main text) producing different gastrulation modes (Figures 2-5), we use a single set of parameter values and boundary conditions, summarized in table 2, except for the initial myosin distribution m(x, t 0 ) and the scalar, constant parameter p 0 , which controls the ratio of isotropic (ingression) to anisotropic (intercalation) active stresses, as described in detail in SI Sec.1.4.2.We change only p 0 and m(x, t 0 ).We set p 0 = 2 when there is active cell ingression and p 0 = 0 otherwise.Additionally, in the case of reduced mesendoderm precursors, we set m(x, t 0 ) as described in eq.(S27).By contrast, in the case of circular mesendoderm precursors (or ancestral condition), we use m(x, t 0 ) = m eq.+ Ae which models a Gaussian function of r, symmetric in α, added to the unstable m eq. and with the same parameters used in the wild type (Section 4.5).We summarize the combination of p 0 and m(x, t 0 ) to produce different gastrulation modes in Table 3.
5.1 Differences with the tensile-ring model from (3,46) An alternative modeling approach for amniotic gastrulation is proposed in (3,46) suggesting that a tensile ring at the EP-EE margin is sufficient to reproduce observation.We highlight below the main difference with our approach.
• Tensile ring model from (3) consists of a viscous flow fitted from data in space and time, F a acts on a ring located at EP-EE margin, and is fitted from data in space and time. (S29) Instead of fitting F a over space and time, (46) propose a mechanistic model to evolve F a at later times starting from an observed initial distribution.
By contrast, we model the fully coupled mechanochemical process and use information from experiments only at the initial time (eq.(S19)), hence providing a predictive mechanistic model.Specifically, we describe the process as a viscous flow      α∆v = ∇p + F a , first two eqs. in (S19), ∇ • v = g(F a , v, p), eq.(S6), ∂ t F a = f (F a , v), second two eqs. in (S19).

(S30)
Using only experimental actomyosin cables distribution at the initial time (HH1) and predicting experimental observations over the following twelve hours (up to HH3) is an important difference compared to the procedure proposed in (3) that requires data fitting over space and time.
• The active force F a in (3) acts only on a ring centered at the EP-EE margin.By contrast, we model F a in the whole embryo because actomyosin cables generating F a are abundantly distributed and evolve dynamically within the embryo consistent with experimental observations (Figure 1).At stages after HH1, experiments show that active myosin is dominant along the PS rather than on the EP-EE boundary (Figures 1D-E), consistent with our model prediction (Figure 2C).To support further this finding, in Figure 3 we show F a inferred from PIV velocity and eq.(S19) at HH3.The active forces are mainly perpendicular to the PS (Figures 3) rather than being organized on a ring-shaped structure.
• (3,46) do not address the forces that drive γ (eq.(S29)), which is fitted from data.By contrast, our model (eq.(S30)) explains and predicts them in time as a response to mechanical stress, recapitulating the flow divergence and areal changes that we see in experiments (Figures 2F-H), as well as the overall Dynamic Morpphoskeletons and other Lagrangian metrics (Figure 2 and Movies 2-3).
Our ability to predict four different perturbations in addition to the wild-type development (Figs.2,4 and Movies 2-12) supports the validity of our model and the underlying biological hypothesis.
Conditions and Initial Conditions, (S19) where B = e ⊗ e = cos 2 φ sin φ cos φ sin φ cos φ sin 2 φ , and m, φ, u, v, x, y, t are nondimensional.In eq.(S19), there are seven nondimensional parameters p 0 , p 1 = µc, p 2 = λ, p 3 = tp tr mn mc , p 4 = tp t d , p 5 = koµuc xc , p 6 = χt p , p 7 = ξ xcuc .(S20) at HH3 and Figure2Cright), eq.(S22) reveals that for all p 0 ≥ 0, ∇m and ∇[∇ • v] have opposite directions, reflecting the expectation that a ridge of high myosin induces a sink in the velocity field (Figure2H).By contrast, if cables are parallel to the myosin ridge (FigureS2B-C) as in the case of perturbation 3 (Figure4of (4) and Movie8), eq.(S22) reveals that for all p 0 > 1 (FigureS2B), ∇m and ∇[∇ • v] have opposite directions.On the contrary, if p 0 < 1, ∇m and ∇[∇ • v] would have the same direction (FigureS2C).Therefore, combining our model with experimental observations, we use p 0 > 1 when active cell ingression is present, and p 0 = 0 otherwise.Supplementary FigureS2: Influence of p 0 on the relationship between ∇m and ∇[∇ • v].We consider a simplified force balance (eq.(S22)) with no shear forces, a simple distribution m made of a ridge in the y direction and uniform actomyosin cables (i.e., ∇ • B = 0) perpendicular (A) or parallel (B, C) to the myosin ridge.Equation (S22) reveals that when cables are perpendicular (A) to the ridge, ∇m and ∇[∇ • v] have opposite directions for all values of p 0 > 0, reflecting the expectation that a ridge of high myosin induces a sink in the velocity field.When cables are perpendicular to the ridge, ∇m and ∇[∇ • v] have opposite directions for p 0 > 1 (B), and the same direction for p 0 < 1.

p 5 2
m m − p 6 m − m x u + p 7 m xx .(S23) Equation (S23) describes the dynamics along a straight line perpendicular to the primitive streak.We set Dirichlet boundary conditions on the velocity u(0, t) = −|u b |, u(L, t) = |u b | consistent with the symmetric cell migration in the extra-embryonic area, and no-flux boundary conditions for m: m x (0, t) = m x (L, t) = 0. We investigate the main features of eq.(S23) by seeking the constant fixed points for m which solve the implicit equation p 3 − (p 4 e − p 5 2 m + p 6 )m = 0. From the first equation, we compute the corresponding equilibrium velocity field u eq.= 2|u b | L x − |u b |.We then study the linear stability of these equilibrium solutions by computing the spectrum of the corresponding linearized PDE.Denoting by δh = [δu, δm] the perturbation from the equilibrium, we compute the linearized version of Eq. (S23) in operator form as Aδh = B(x)δh where Supplementary Figure S3: 1D model.(A) Top: Graph of p 3 − (p 4 e − p5 2 m + p 6 )m = 0 shows three m eq.marked with a blue circle, a red asterisk, and a green circle.Bottom: Truncated spectrum of the linearized PDE around corresponding m eq., u eq., consisting of the 500 eigenvalues with the highest real part.The spectrum at the intermediate equilibrium has positive real parts, implying that the intermediate equilibrium is linearly unstable.(B) The initial condition of the nonlinear 1D model using as initial myosin concentration a Gaussian added to the higher unstable m eq. .The blue curve shows the associated initial velocity u 0 , and the red dashed curve the initial velocity divergence u x0 .(C) Same as B at the final time.The full time evolution is available as Movie S1.Sample parameters: p 0 = 2, p 1 = 0.25, p 3 = 50, p 4 = 100, p 5 = 1.25, p 6 = 1, p 7 = 0.01, u b = 0.5, L = 4, dx = 0.008, dt = 0.0001.

Supplementary Figure S5 :
Visualization of directional propagation of active stress intensity.A) Initial and final distribution of active stress intensity m under evolution m t = p 7 ∆ φ m with no evolution of v or φ.B) Initial and final distribution of p 7 ∆ φ m.Active stress orientations are indicated by white (A) and black (B) line segments.Sample parameter: p 7 = 0.01.

Supplementary Figure S6 :
Sensitivity to p 5 .A) Exponential actomyosin dissociation rate p 4 e − p5 2 m for two values of p 5 .B) Reducing p 5 by 20%, almost doubles the unstable m eq .

Table 2 :
Parameter values used in all two-dimensional simulations.

Table 3 :
Summary of p 0 and m(x, t 0 ) used to generate different gastrulation modes.