Rectified random cell motility as a mechanism for embryo elongation

Most vertebrate species exhibit an elongated body axis. This characteristic pattern is established during embryogenesis as the tissues progressively form in an anterior-to-posterior direction. Microsurgical ablation of the posterior presomitic mesoderm (PSM), which contains the precursors of skeletal muscles and the axial skeleton, severely reduces posterior elongation movements, indicating that this tissue plays a major role in the control of posterior elongation of the embryonic axis. Analysis of cell motility in the chicken embryo PSM (Bénazéraf et al., 2010) shows that there is an anterior-posterior gradient in the activity of cells. However, locally, the motility inside the PSM of chicken and zebrafish embryos is manifested by random, undirected, Brownian-like cellular motion (Bénazéraf et al., 2010; Lawton et al., 2013). These random, diffusive movements contrast with the oriented cell intercalation movements controlling elongation of the anterior parts of the embryo. Within the PSM, this motility gradient is downstream of a chemical gradient of the secreted fibroblast growth factor FGF, as shown schematically in Fig. 1A. FGF is known to play an important role in cell motility (Delfini et al., 2005); indeed, increasing FGF concentration in the PSM causes cell motility to increase without any orientation preference (Bénazéraf et al., 2010), reducing the speed of elongation. Earlier qualitative models focused on the role of the motility gradient as a driving force for elongation (Bénazéraf et al., 2010), but how this motility gradient is related to elongation is still unclear.

The observation that at a local level cell motility is random and undirected needs to be reconciled with the emergence of global body elongation. As new cells enter the PSM, they are exposed to a high concentration of FGF and become highly motile, but do not move in an oriented manner. Owing to the confinement of motility driven by the presence of relatively immobile and stiff lateral tissues (Bénazéraf et al., 2017), the cells generate an effective mechanical pressure. Thus, a potential mechanism for the rectification that leads to the observed directional cell velocity and body elongation could be this mechanical pressure. Because the expression of FGF is highest at the posterior PSM and decreases away from it, there is a natural reduction in motility anteriorly, consistent with observations. Eventually, the effects of anterior adhesion dominate over motility and this causes cellular condensation into somites.

To quantify these processes, here we use experimental observations to measure the effective diffusivity of cells as well as their advection speed as a function of their location relative to the last formed somite. Our observations suggest a minimal microscopic cellular description of a zone of proliferating cells with high motility, which we use to develop both a quantitative cellular model and an equivalent macroscopic continuum theory, thus providing a theoretical and computational framework for body elongation. These complementary approaches yield simple expressions for the speed and scale of body elongation that are consistent with our experimental measurements, with potential implications for our understanding of outgrowth-driven morphogenesis in other settings in which FGF-driven gradients in motility and contractility drive the extension of the limb bud (Gros et al., 2010) and the gut (Nerurkar et al., 2019).

Our experiments were carried out with chicken embryos at Hamburger-Hamilton (HH) stages 10-11 (Hamburger and Hamilton, 1992), corresponding to the period when the elongation of the embryo is most substantial (Denans et al., 2015). A time series of an elongating PSM is shown in Fig. 1B (see also Movie 1). To measure the elongation rate, we registered the movement of the embryo with respect to the last-formed somite at the beginning of the experiment (t=0) as depicted by the black dashed line in Fig. 1B, and tracked the advancement of the Hensen's node, L_e(t) as a function of time. In Fig. 1C, we see that L_e(t) increases linearly with time, with a mean elongation rate V=(2.8±0.3)×10⁻² μm/s (averaged from five embryos).

To evaluate the role of cell motility on overall body elongation, we then examined the movement of cells by electroporating the PSM cells with fluorescent reporters that label cell nuclei specifically (Hatakeyama and Shimamura, 2008), as shown in Fig. 1D (see also Movie 1). As a first approximation, the motion of the cells can be considered two-dimensional in the anterior-posterior and medio-lateral directions because the relative dorso-ventral depth of the PSM is small. In the reference frame of the last-formed somite at t=0, during a fixed acquisition time of 4 h, posterior cell trajectories showed a larger net displacement than the anterior ones, consistent with prior experiments (Bénazéraf et al., 2010). To quantify the variations of cell motility along the body axis, we divided the PSM into three regions (anterior, middle and posterior), as shown schematically in Fig. 1E, and obtained the mean square displacement (MSD) of the cells given by < Δr²(t)>, where defines the distance that the cell travels in a lag time t, as shown in the inset of Fig. 1F (Qian et al., 1991; Wirtz, 2009).

Our observations quantify how the cells move in space and time along the PSM. When new cells are added close to the tail bud (TB), they are highly motile but move in random directions; as they move further away from the zone of proliferation in the TB, on average they gradually slow down and stop moving.

Next, we provide microscopic and macroscopic chemo-mechanical mathematical and computational models based on the observed gradient in motility and the confinement due to the more rigid lateral tissues, and show that they explain the observed rectified cell motility and the resulting posterior body elongation.

The TB boundary is modeled as a wall of immobile cells attached to their neighbors by elastic springs with spring constants, k_chain=2k, and allowed to move posteriorly as a result of the mechanical pressure exerted by the motile cells anterior to it. Our simulations show that after an initial transient state, the motion of the wall reaches a steady state whereby cells added at a constant rate cause the wall to move at a constant velocity (Movie 2). The quasi-one-dimensional dynamics manifest in the simulations arises because of the impermeable boundary conditions in the direction perpendicular to the direction of elongation.

The fraction of motile cells is much larger near the moving wall, and as one moves anteriorly, this density falls off quickly because of the decrease in FGF activity. Changing FGF activity by changing τ, and thus varying the probability with which each individual cells stops moving (and therefore the total fraction of moving cells), changes the velocity and motility profiles (see Fig. 2C,D).

The first equation describes the variations in the density of actively motile cells that also gradually decay at a rate τ, whereas the second equation characterizes how the stresses associated with the active pressure generated by the random motility of the cells causes the tissue to respond and push the TB. Here, η is the viscosity of the PSM and ξ is the viscous friction associated with motion of the elongating PSM relative to the surrounding tissues (endoderm, ectoderm, neural tube and lateral plate). We note that these friction coefficients define a characteristic viscous screening length l_ξ∼(η/ξ)^1/2 that we assume is of the order of a few cell sizes, i.e. l_ξ∼a.

The equations need a closure relation for the pressure; a minimal relation assumes that the active pressure is proportional to the density of motile cells p∼αρ, consistent with both our microscopic model, with ⁠, and also consistent with earlier simulations of active Brownian particles (Mallory et al., 2014). Although other relations of the form p∼αρ^q are plausible, we choose to stay with a minimal model here.

To understand the dependence of the solution of Eqns 11 and 12 on the problem parameters, we can rewrite the equations in a form that depends on five dimensionless variables:

Here Π₁ is the ratio of the cell size a and the elongation length scale ⁠, Π₂ is the ratio of the maximum active stress αρ₀ and the external viscous friction ξVL∼ξV²τ, and Π₃ is the ratio of the internal viscous stress ηV/L and the active stress αρ₀. We note that the parameter ⁠, where l_ξ is the viscous screening length defined ealier (after Eqn 12). Finally, there are two dimensionless parameters associated with the boundary conditions Π₄, which is the ratio of the external stress at the PSM boundary F and the active stress αρ₀, and Π₅, which is the scaled ratio of the internal cell diffusivity and the rate of addition of cells at the PSM boundary ⁠.

Using the following parameter values Π₁=0.01, Π₂=2, Π₃=1, Π₄≈0.005, Π₅=0.003, we solve the initial boundary value problem (Eqns 17-21) using a finite difference method. In Fig. 3B,C, we see that our calculated profiles for the velocity v(x) and diffusivity compare well with the experimentally observed profiles, and are suggestive of a simple exponential law (see also Fig. S5). In Fig. 3B, we show that the system evolves into a steady state following a short transient (insets), with an active zone near the TB that moves at constant speed. Both the cellular motility (Fig. 3B) and diffusivity (Fig. 3C) decay towards the anterior where the FGF degrades. Comparing our numerical results for the moving average of cell velocity 〈v〉 and the rescaled diffusivity with the experimentally determined values as a function of position relative to Hensen's node (see Materials and Methods) we see that the profiles compare well. Furthermore, on converting the results to dimensional values, we can estimate the active stress αρ₀≈5−50 Pa, qualitatively similar to measurements on amniote embryos (Zhou et al., 2009; Serwane et al., 2017). From the simulations, we can also estimate the typical length associated with motility decay L∼1200 μm and the length scale l_ξ∼1800 μm associated with friction, which are of the same order as the experimental estimates, L=Vτ∼260−520 μm and ⁠, 3-15 times the cell size.

A sensitivity analysis of our continuum model (see section 3 in the supplementary information) shows that the parameters Π₁, the ratio of the cell size to the elongation length scale that arises from activity, and Π₃, the ratio of viscous stresses to the active stress, affect the velocity and motility profile significantly. The other parameters, Π₂, Π₄ and Π₅, have a negligible effect on our results over a range relevant for our problem (for the range of parameter values used, see Table S3; for the associated results, see Figs S6-S8). The dependence of the length scales corresponding to the typical motility and velocity profiles as a function of the dimensionless parameters Π₁ and Π₃ is shown in Fig. S10.

In vertebrate embryos, posterior structures are formed sequentially by a combination of cell proliferation and cell motility that together leads to body elongation. Although it has long been observed that the elongation process involves the posterior displacement of the TB, the physical mechanisms underlying this have not been clearly elucidated. Here, we have quantified this process experimentally and theoretically and shown that it occurs as a result of two effects: the addition of motile cells at a boundary (the TB) and confinement both laterally and anteriorly. These two effects lead to the rectification of random cell diffusivity, which generates the forces underlying elongation and the emergence of characteristic velocity and activity/diffusivity length scales. Simple cellular models and a one-dimensional continuum framework capture the essence of the behavior in terms of a pair of fundamental dimensionless parameters that characterize the ratio of the active driving stress to the internal and external viscous stresses, and are consistent with a narrow boundary layer of activity near the TB, in agreement with observations. Natural extensions of the models include adding the effects of dimensionality, such as the lateral growth in the absence of FGF decay (Bénazéraf et al., 2010), two-dimensional compressive stresses that generate vortical flows near the TB (Xiong et al., 2020), and the gradients in viscosity from the anterior to the posterior region (Mongera et al., 2018). Such multidimensional models could also be modified to describe other embryonic outgrowths processes that exhibit graded diffusive behavior of cells downstream of FGF signaling, such as the budding of vertebrate limbs (Gros et al., 2010) and the onset of gut elongation (Nerurkar et al., 2019).

Fertilized chicken eggs were obtained from a commercial provider (Les Couvoirs de l'Est, Willgottheim, France) and incubated at 37°C in a humidified incubator. After 24 h, HH stage 4-5 embryos were mounted on filter paper and transferred ventral side up to 35 mm agar/albumen petri dishes for injection (Chapman et al., 2001). Electroporation of the PSM was performed using H2B-mCherry or H2B-Venus nuclear markers as described previously (Bénazéraf et al., 2010). Electroporated embryos were returned to incubator and left to grow to HH stage 10-11 before imaging.

The imaging procedure used here was similar to a previously described procedure (Czirók et al., 2002; Rupp et al., 2003). Briefly, the embryos were transferred to custom-made six-well observation chambers containing agar/albumin gel, and positioned ventral side up. The time-lapse imaging was performed at 37°C using a motorized upright microscope (Leica DMR, Leica Microsystems) with a 10× objective (N.A. 0.3) and a CCD digital camera (QImaging Retiga 1300i) at a rate of 10 frames/h. At each time point, brightfield and fluorescent images of the embryo were taken at three various fields to cover the total length of the PSM. Post-acquisition processing was performed on images as described previously (Czirók et al., 2002; Rupp et al., 2003) to obtain a 2D time series. Cell tracking and trajectory analysis was performed on fluorescent images using custom-made MATLAB (MathWorks) routines. Trajectories were analyzed on 4-h long movies. To account for cell position modification due to tissue drift, we divided our movies to 1 h segments. The location of each cell was defined with respect to the node at the beginning of each 1 h trajectory. For each cell trajectory, the MSD was calculated and adjusted with Eqn 1 to obtain D and v. Further, the moving average filter, Smooth (MATLAB), was used for the plots of D and v as a function of position to obtain Fig. 3.

The discrete cell simulations were performed by discretizing Eqn 2 using the Euler–Maruyama method and integrating in time. In the initial state, all the particles are arranged in a square lattice (particles are slightly perturbed from the lattice) in a region that has dimensions of 80a×160a. At a given time, the particles can be at x∈[0,X_max], y∈[0,Y_max], where X_max is the position of the rigid boundary and Y_max=160a is the lateral width. We use reduced, dimensionless units in which all energies are given in terms of a typical energy w: U*=U/w, all lengths in terms of the typical cell size a, r*=r/a and masses in terms of a mass unit M, m*=m/M. Using these fundamental units, we can also reduce the effective temperature, i.e. the active random motility ⁠, time ⁠, fraction ρ*=Na³/V (V is the volume) and any other physical quantity of interest, where and ρ* are all dimensionless.