1. Developmental Biology
  2. Physics of Living Systems
Download icon

Size-dependent patterns of cell proliferation and migration in freely-expanding epithelia

  1. Matthew A Heinrich
  2. Ricard Alert
  3. Julienne M LaChance
  4. Tom J Zajdel
  5. Andrej Košmrlj
  6. Daniel J Cohen  Is a corresponding author
  1. Department of Mechanical and Aerospace Engineering, Princeton University, United States
  2. Lewis-Sigler Institute for Integrative Genomics, Princeton University, United States
  3. Princeton Center for Theoretical Science, Princeton University, United States
  4. Princeton Institute for the Science and Technology of Materials (PRISM), Princeton University, United States
Research Article
  • Cited 0
  • Views 1,266
  • Annotations
Cite this article as: eLife 2020;9:e58945 doi: 10.7554/eLife.58945

Abstract

The coordination of cell proliferation and migration in growing tissues is crucial in development and regeneration but remains poorly understood. Here, we find that, while expanding with an edge speed independent of initial conditions, millimeter-scale epithelial monolayers exhibit internal patterns of proliferation and migration that depend not on the current but on the initial tissue size, indicating memory effects. Specifically, the core of large tissues becomes very dense, almost quiescent, and ceases cell-cycle progression. In contrast, initially-smaller tissues develop a local minimum of cell density and a tissue-spanning vortex. To explain vortex formation, we propose an active polar fluid model with a feedback between cell polarization and tissue flow. Taken together, our findings suggest that expanding epithelia decouple their internal and edge regions, which enables robust expansion dynamics despite the presence of size- and history-dependent patterns in the tissue interior.

eLife digest

Cells do not exist in isolation. Instead, they form tissues, where individual cells make contact with their neighbors and form microscopic ‘architectures’. Epithelia are a type of tissue where cells are arranged in flat sheets, and are found in organs such as the lining of the kidney or the skin.

Tissues need to grow, especially early in life. If tissues are damaged – for example, if the skin is cut or grazed – cells also need to divide (to create new healthy cells) and move as a group (to close the wound). Such coordinated motions result in cells exhibiting distinct group behaviors, similar to those observed within crowds of people or schools of fish. If coordination breaks down, problems can happen such as uncoordinated tissue growth seen in cancer.

However, how cell movements are coordinated is still not fully understand. For example, researchers know that cells’ positions within a group can determine how they behave, meaning that even the same type of cell could behave differently at the edge or center of a tissue. This suggests that the initial size and shape of a tissue should influence its subsequent growth and behavior; however, the nature of this influence is still largely unknown. Heinrich et al. therefore wanted to determine the differences in the way larger and smaller tissues grow.

Microscope imaging was used to track the growth of circular, artificial tissues made from single-layered sheets of dog kidney cells grown in the laboratory. Comparing how quickly the tissues expanded revealed that the area of tissue circles that started out smaller increased at a much faster rate than that of tissue circles that were larger to begin with. This turned out to be because the edges of the tissues grew at a constant speed, independent of their initial size or shape, but circles with a smaller area have a larger proportion of cells on their edges. The motions of the cells at the center of the tissues had no effect on how the edges of the tissue grew. A final observation was that the way tissues of a given size behaved depended on whether they had grown to be that size, or they started off that big.

These results shed light on how groups of cells interact in growing tissues. In the future, this information could be used to predict how different tissues grow over time, potentially helping scientists engineer better artificial tissues or organs for transplantation.

Introduction

Writing in 1859, physiologist Rudolf Virchow presented the concept of the ‘Zellenstaat’ or ‘Cell State,’ describing tissues as ‘a society of cells, a tiny well-ordered state’ (Virchow, 1855). This social framework motivated Abercrombie and Heaysman, 1954 work on cellular behavior that elucidated how encounters between cells can regulate locomotion and proliferation via contact inhibition. Since then, concerted interdisciplinary effort has been brought to bear on understanding how cell-cell interactions give rise to the complex collective behaviors driving so many crucial biological processes. One of the most foundational collective behaviors is collective cell migration—the directed, coordinated motion of cellular ensembles that enables phenomena such as gastrulation, wound healing, and tumor invasion (Friedl and Gilmour, 2009). Given this importance, considerable effort spanning biology, engineering, and physics has been directed towards understanding how local cellular interactions can give rise to globally coordinated motions (Alert and Trepat, 2020; Hakim and Silberzan, 2017).

Studies of collective cell migration are most often performed using epithelial tissues due to their fundamental role in multicellular organisms and strong cell-cell adhesion, which in turn gives rise to elegant, cohesive motion. Moreover, given that epithelia naturally form surfaces in vivo, studying epithelial layers in vitro has a physiological basis that can inform our understanding of processes such as healing (Poujade et al., 2007), envelopment (Steinberg, 2007), and boundary formation (Dahmann et al., 2011). These features have made epithelia both the gold standard in collective cell migration studies, and one of the most well-studied models for biological collective behaviors.

Due to the complexity of collective behaviors, much effort has gone towards reductionist assays that restrict degrees of freedom and ensemble size to simplify analysis and interpretation. One such approach is to confine a tissue within predefined boundaries using micropatterning to create adhesive and non-adhesive regions (Doxzen et al., 2013; Deforet et al., 2014; Notbohm et al., 2016; Pérez-González et al., 2019; Peyret et al., 2019; Petrolli et al., 2019). Such confinement mimics certain in vivo contexts such as constrained tumors as well as aspects of compartmentalization during morphogenesis (Lecuit and Lenne, 2007). Alternately, many studies have explored the expansion of tissues that initially grow into confluence within confinement but are later allowed to migrate into free space upon removal of a barrier. A popular assay of this type relies on rectangular strips of tissue that are allowed to expand in one or both directions (Poujade et al., 2007; Trepat et al., 2009; Petitjean et al., 2010; Reffay et al., 2011; Nnetu et al., 2012; Serra-Picamal et al., 2012; Zhang et al., 2017; Uroz et al., 2018; Tlili et al., 2018), where averaging along the length of the strip can reveal coordinated population-level behaviors such as complex migration patterns, non-uniform traction force fields, and traveling mechanical waves. Other studies have focused on the isotropic expansion of micro-scale (< 500 μm diameter) circular tissues using the barrier stencil technique (Jang et al., 2017) as well as photoswitchable substrates (Rolli et al., 2012). Still more work has explored approaches to induce directional migration, from geometric cues to applied electric fields (Vedula et al., 2012; Cohen et al., 2014).

In contrast to micro-scale confinement assays, other work has focused on large, freely-expanding tissues of uncontrolled initial size and shape, which grow from either single cells (Puliafito et al., 2012; Huergo et al., 2011) or cell-containing droplets (Lee et al., 2013; Beaune et al., 2014). Related experiments track long-term growth of cell colonies via images taken once per day over several days, but this low temporal resolution cannot access timescales over which migration is important (Huergo et al., 2011; Simpson et al., 2013). Thus, there is still a lack of assays to study long-term expansion and growth of large-scale tissues with precisely-controlled initial conditions, especially initial tissue size, shape, and density.

To address this gap, we leveraged bench-top tissue patterning (Poujade et al., 2007; Cohen et al., 2016) to precisely pattern macro-scale circular epithelia of two sizes (>1 mm in diameter) and performed long-term, high frequency, time-lapse imaging after release of a barrier. To elucidate the consequences of size effects on the tissue, we tracked every cell, relating the overall expansion kinetics to cell migration speed, cell density, and cell-cycle dynamics. We find that, whereas the tissue edge dynamics is independent of the initial conditions, the tissue bulk exhibits size-dependent patterns of cell proliferation and migration, including large-scale vortices accompanied by dynamic density profiles. Together, these data comprise the first comprehensive study of macro-scale, long-term epithelial expansion, and our findings demonstrate the importance of exploring collective cell migration across a wider range of contexts, scales, and constraints.

Results

Expansion of millimeter-scale epithelia of different sizes and shapes

We began by characterizing the overall expansion and growth of tissues with the same cell density but different initial diameters of 1.7 mm and 3.4 mm (a 4X difference in area, with tissues hereafter referred to as either ‘small’ or ‘large’), using an MDCK cell line stably expressing the 2-color FUCCI cell-cycle marker (Sakaue-Sawano et al., 2008; Streichan et al., 2014; Uroz et al., 2018; Beaune et al., 2014; Benham-Pyle et al., 2016). We patterned the tissues by culturing cells in small and large circular silicone stencils for ∼18 hr (Cohen et al., 2016; Poujade et al., 2007), whereupon stencils were removed and tissues were allowed to freely expand for 46 hr (Figure 1A, Figure 1—video 1), while images were collected at 20 min intervals using automated microscopy (see Materials and Methods). Our cell seeding conditions and incubation period were deliberately tuned to ensure that the stencils did not induce contact inhibition of proliferation prior to stencil removal (checking FUCCI to ensure the tissue was not arrested in G1). Upon stencil removal, tissues expanded while maintaining their overall circular shape throughout the 2 day experiment. Unless otherwise noted, cell density at stencil removal was ∼2700 cells/mm2, a value consistent with active and growing confluent MDCK epithelia (Streichan et al., 2014; Uroz et al., 2018).

Figure 1 with 5 supplements see all
Expansion dynamics of millimeter-size cell monolayers.

(A) Footprint throughout 46 hr growth period of representative small (left) and large (right) circular tissues, with the tissue outlines drawn at 4 h increments. Initial diameters were 1.7 mm and 3.4 mm. (B) Small circles exhibit faster relative area, A(t)/A0, increase than large circles, where A0 and A(t) are the areas of tissues at the beginning of the experiment and at time t, respectively. Purple points show the relative area increase, A(t+t0)/A(t0), of small tissues from the time t0=30 h when they reached the size of the large circles. (C) Average tissue density ρ(t)=N(t)/A(t) has non-monotonic evolution in small tissues but monotonically increases in large tissues, where N(t) is the number of cells in a tissue at time t. (D) Edge radial velocity vr is largely independent of initial tissue size and cell density. We grouped initial cell densities as ρ1=[2350,3050] cells/mm2, ρ2=[1650,2350] cells/mm2, and ρ3=[1300,1650] cells/mm2. (E) Experimental data on tissue shape and model fits. Assuming a constant migration speed vn in direction normal to the edge, we can predict the area expansion dynamics of elliptical tissues with different aspect ratios. The model fits our data for all tissues with vn29.5 µm/hr, yielding normalized χ2 values of 0.79, 0.13, and 0.06 for aspect ratios of 8, 4, and 1 respectively (χ2< 1 indicates a good fit; see Materials and methods). In B, data are from n = 16 tissues across five independent experiments (small and large circles). In C, n = 11 across four experiments for small circles, and n = 9 across three experiments for large circles. In D, n = 16 across five independent experiments for small and large circles, ρ=ρ1; n = 13 across three experiments for small circles, ρ=ρ2; and n = 11 across three experiments for small circles, ρ=ρ3. In E, n = 4 across two experiments for a/b = 1 and a/b = 4, and n = 5 across two experiments for a/b = 8. Shaded regions correspond to standard deviations.

First, we measured relative areal increase (Figure 1B) and relative cell number increase (Figure 1—figure supplement 1) of small and large tissues. By 46 hr, small and large tissues had increased in area by 6.4X and 3.3X, respectively, while cell number increased by 9.2X and 5.5X, respectively. Since proliferation outpaces area expansion in long-term growth, average tissue density increased by the end of the experiment. The evolution of average tissue density was more complex, however, as small tissues experienced a density decrease from 4 to 12 hr while large tissues exhibited a monotonic increase in cell density (Figure 1C). Accordingly, at any given time after stencil removal, large tissues had a higher density than small tissues. Non-monotonic density evolution has been observed in thin epithelial strips (Poujade et al., 2007) and likely arises from competition between migration and proliferation dynamics, which we discuss later.

We then related area expansion to the kinematics of the tissue edge. To quantify edge motion, we calculated the average radial velocity of the tissue boundary, vr(t), at 1 hr intervals over 46 hr (Materials and methods). We found that vr is independent of both tissue size and a wide range of initial cell densities, in all cases reaching ∼30 μm/h after ∼16 hr (Figure 1D). Before reaching this constant edge velocity, vr ramps up during the first 8 hr after stencil removal, and, notably, overshoots its long-time value by almost 30%. We hypothesize that the overshoot is due to the formation of fast multicellular finger-like protrusions that emerge at the tissue edge in the early stages of expansion and then diminish (Figure 1—video 2). This hypothesis is supported by a recent model showing that edge acceleration (as observed during the first 8 hr in Figure 1D) leads to finger formation (Alert et al., 2019). It is remarkable that the edge radial velocity vr(t) is independent of the initial tissue size and density, especially considering that cell density evolution shows opposite trends at early stages of expansion for small and large tissues (Figure 1C). This observation suggests that the early stages of epithelial expansion are primarily driven by cell migration rather than proliferation or density-dependent decompression and cell spreading.

The observation that vr is independent of tissue size ought to explain why small tissues have faster relative area expansions than large tissues. We hypothesized that the relation between tissue size and areal increase could be attributed primarily to the perimeter-to-area ratio. Assuming a constant edge velocity vn normal to the tissue boundary, the tissue area increases as dA=Pvndt, where P is the perimeter of tissue and dt is a small time interval. Thus, the relative area increase dA/A=(P/A)vndt scales as the perimeter-to-area ratio, which is inversely proportional to the radius for circular tissues, so the relative area increases faster for smaller tissues (Figure 1B).

To verify that the perimeter-to-area ratio is proportional to the relative area increase, we analyzed elliptical tissues with the same area and cell density but different perimeters (Figure 1—video 3). Increasing the perimeter-to-area ratio of a tissue by increasing its aspect ratio indeed leads to faster relative area expansion (Figure 1E). A simple, edge-driven expansion model with linear increase of the tissue major and minor axes predicts A(t)/A(0)=(a+vnt)(b+vnt)/(ab), where a and b are the initial major and minor axes of the tissue. This model fits our data well assuming the same edge speed vn29.5 μm/h for all tissues (Figure 1E). This observation suggests that edge speed is mostly independent of edge curvature. However, we measure a smaller edge speed at the major axes of ellipses, which are high-curvature points with radius of curvature rc0.75 mm (Figure 1—figure supplement 2). Such high curvatures are concentrated around the major axes of our elliptical tissues. However, most of the tissue edge has a smaller curvature, and therefore advances at a curvature-independent speed. Further, even high curvature regions blunt due to expansion over time (see Figure 1—video 3). As a result, our model with a single edge speed vn29.5 μm/h is sufficient to capture the area expansion of both circular and elliptical tissues (Figure 1E).

Together, our findings demonstrate that epithelial shape and size determine area expansion dynamics via the perimeter-to-area ratio. This relationship results from the fact that tissues exhibit a constant, size-independent, migration-driven edge speed normal to tissue boundary. Since initial tissue size does not affect boundary dynamics, but does impact the relative growth and expansion of the tissue, we hypothesize that cells in the tissue bulk exhibit tissue size-dependent behaviors.

Spatiotemporal dynamics of migration speed and radial velocity

Having demonstrated the role of the boundary in the expansion of large-scale epithelia, we sought to relate tissue areal expansion rate to internal collective cell migration dynamics. We used Particle-Image-Velocimetry (PIV, Materials and methods) to obtain flow fields describing cell migration within freely expanding epithelia (Poujade et al., 2007; Petitjean et al., 2010; Angelini et al., 2010; Cohen et al., 2014; Aoki et al., 2017). We constructed kymographs (Materials and Methods) to display the full spatiotemporal flow patterns of the tissue (Figure 2A,B; Serra-Picamal et al., 2012; Zhang et al., 2017), averaging over the angular direction and over 16 tissues (for representative kymographs, see Figure 2—figure supplement 1). We also separately show time evolution (Figure 2C) and spatial profiles (Figure 2D) of speed and radial velocity to compare small and large tissues.

Figure 2 with 1 supplement see all
Speed and radial velocity in inner and outer tissue zones.

(A,B) Average kymographs of (A) speed and (B) radial velocity vr throughout expansion for small (left) and large (right) tissues. (C) Evolution of the average speed of boundary (top) and center (bottom) zones, defined as regions extending ∼200 μm from the tissue center and tissue edge, respectively. This width of the zones corresponds approximately to the velocity-velocity correlation length for MDCK cells (Petitjean et al., 2010). While the speed in the edge zone remains high in both small and large tissues, the speed in the center zone begins to decrease ∼24 hr sooner in large tissues than in small tissues, as the central zone of the small tissues has particularly high speed from 18 to 36 hr. (D) Profiles of speed (top) and radial velocity (bottom) at 36 hr, from the edge of the tissue inwards. Arrows indicate that the tissues are indexed from the edge of the tissue inwards. All data are from n = 16 tissues across five independent experiments (small and large circles). Speed and radial velocity profiles of large and small tissues match closely for the first 500 μm from the tissue edge. The average difference between the profiles in this zone is 0.39 μm/h (speed) and 0.27 μm/h (radial velocity), respectively, while the smallest standard deviation for any point in either profile is 0.56 μm/h.

Kymographs of speed and radial velocity reveal the existence of an edge region of fast, outward, radial cell motion (Figure 2A,B), with speeds similar to the radial edge velocity reported in Figure 1D. Up to ∼500 μm from the tissue edge, the speed and radial velocity profiles are practically identical for small and large tissues (Figure 2D), showing that cell motion near the tissue edge is independent of tissue size.

The tissue centers, in contrast, exhibit size-dependent behaviors. For both small and large tissues, a wave front of cell speed and radial velocity propagates toward the tissue centers at ∼90 μm/h (Figure 2A and B, dashed lines). This is approximately 3X faster than the tissue edge speed, consistent with previously described waves of strain rate in cell monolayers (Serra-Picamal et al., 2012). Soon after the wave of radial velocity reaches the center, it retreats, leaving a region of low radial velocity that increases in extent in the center of both small and large tissues (Figure 2B). This decrease of radial velocity is accompanied by a reduction in cell speed in the center of large tissues but not in small tissues, in which cell speed remains high until 36 hr (Figure 2A and C Bottom). We examine the behavior of this high-speed but low-radial-velocity central region of small tissues in the next section.

Emergence of large-scale vortices

The propagation of low radial velocity out from the center of small tissues coincides with the formation and expansion of a millimeter-scale, persistent vortex (see Figure 3A, Figure 3—video 1 for representative vortex). These large vortices are observed in both small and large tissues (Figure 3—video 2), but they only reach tissue-spanning sizes in small tissues.

Figure 3 with 5 supplements see all
Vortex formation in expanding tissues.

(A) Vortical flows seen from 10 hr traces of cell trajectories in small (left) and large (right) tissues. We color each trajectory according to its local orientation. (B) Growth rate of perturbations of wave vector modulus q around the unpolarized state of the tissue bulk, Equation 3. Perturbations with wavelength longer than 2π/qc grow (Ω>0), leading to large-scale spontaneous flows in the tissue bulk. We show curves for the following values of the polarity-velocity coupling parameter: νs=0,1,2,3,4 mm−1. For the remaining parameters, we took Ta=100 Pa/μm, ξ=100 Pa⋅s/μm2, η=25 MPa⋅s, γ=10 kPa⋅s, a=20 Pa, K=10 nN, as estimated in Pérez-González et al., 2019. (C) Average kymographs of vorticity show that the vortex in small tissues appears in the center and expands to >1 mm (n = 16), while vorticity in large tissues is present only during the early stages of tissue expansion (n = 16). The black bars indicate a characteristic vortex size. (D) Characteristic vortex size (marker size), time (horizontal axis), and intensity (vertical axis) of each tissue’s maximal vortex intensity. Small tissue vortices are generally more intense, with p<0.0001. (E) For small tissues, the time of maximal vortex intensity decreases with the initial cell density.

To visualize the form and scale of these vortices, we tracked individual cell motion and colored cell trajectories according to their orientation (Püspöki et al., 2016) for a representative small and large tissue vortex (see Figure 3A and Materials and Methods). We plotted trajectories for the time periods that the vortex was most apparent, which was 20–40 hr in the small tissue (Figure 3A, left) and 10–30 hr in the large tissue (Figure 3A, right). During the vortex period in small tissues, cell trajectories are primarily radial in the boundary zone, but mainly tangential in the entire central zone (Figure 3A left, see Figure 3—figure supplement 1 for vortex trajectory quantification).

To understand the emergence of the vortices, we build on a continuum physical model of tissue spreading that describes the cell monolayer as a two-dimensional compressible active polar fluid (Blanch-Mercader et al., 2017; Pérez-González et al., 2019; Alert et al., 2019). Consistent with our velocity measurements (Figure 2C), we assume that cells at the edge zone are radially polarized and motile, whereas cells in the bulk of the tissue are unpolarized and non-motile. We describe cell polarization at a coarse-grained level via a polarity field 𝐩 that obeys the following dynamics (Alert and Trepat, 2020):

(1) t𝐩=𝐡γ+νs𝐯.

Here, γ is the rotational viscosity that damps polarity changes. Respectively, 𝐡=-a𝐩+K2𝐩 is the so-called molecular field that governs polarity relaxation: the first term drives the polarity to zero, and the second term opposes spatial variation of the polarity field. As a result of these terms, the radial polarity at the tissue edge decays over a length scale Lc=K/a into the tissue bulk.

With respect to previous models of tissue spreading, we add the last term in Equation 1, which couples the polarity to the tissue velocity field 𝐯. This coupling is a generic property of active polar fluids interacting with a substrate (Brotto et al., 2013; Kumar et al., 2014; Oriola et al., 2017; Maitra et al., 2020). Previous works in agent-based models showed that similar polarity-velocity alignment interactions (Alert and Trepat, 2020) can lead to waves (Petrolli et al., 2019), flocking transitions (Szabó et al., 2006; Henkes et al., 2011; Basan et al., 2013; Malinverno et al., 2017; Giavazzi et al., 2018), and vortical flows (Rappel et al., 1999; Camley et al., 2014; Li and Sun, 2014; Segerer et al., 2015; Barton et al., 2017; Lin et al., 2018) in small, confined, and polarized tissues. Here, using a continuum model, we propose that cell polarity not only aligns with but is also generated by tissue flow, and we ask whether this polarity-velocity coupling can lead to large-scale spontaneous flows in the unpolarized bulk of unconfined tissues.

To determine the flow field 𝐯, we impose a balance between internal viscous stresses in the tissue, with viscosity η, and external cell-substrate forces, including viscous friction with coefficient ξ, active traction forces with coefficient Ta, and the cell-substrate forces associated with the polarity-velocity coupling νs:

(2) η2𝐯=ξ𝐯-Ta𝐩-νs𝐡.

This force balance predicts that even if cell polarity, and hence active traction forces, are localized to a narrow boundary layer of width Lc50 μm (Blanch-Mercader et al., 2017; Pérez-González et al., 2019), cell flow can penetrate a length λ=η/ξ into the tissue. Based on our measurements (Figure 2D), we estimate λ0.5-1 mm, which is larger than the velocity correlation length of 200 μm in the tissue bulk (Petitjean et al., 2010).

A linear stability analysis of Equations 1 and 2 shows that perturbations of wave number q around the quiescent (𝐯=0) and unpolarized (𝐩=0) state grow with a rate

(3) Ω(q)=-aγ(1+Lc2q2)+Taνs-aνs2(1+Lc2q2)ξ(1+λ2q2).

This result shows that, if Taνs>a(ξ/γ+νs2), the unpolarized state of an active polar fluid described by Equations 1 and 2 is unstable (Ω>0) to perturbations of wavelength longer than a critical value 2π/qc given by Ω(qc)=0 (Figure 3B). This analysis suggests that, for tissues larger than this critical value 2π/qc, the quiescent tissue bulk becomes unstable and starts to flow spontaneously at large scales, consistent with the emergence of large-scale vortices. The mechanism of this instability is the positive feedback between flow-induced cell polarization and the flows due to migration of polarized cells. The fact that a critical size of the order of millimeters is required for this long-wavelength instability might explain why large-scale vortices have not been observed in previous studies, which considered smaller tissues.

Vortex kinematics

To quantify the kinematics of the large-scale vortical flows, we obtained the vorticity field ω(𝐫,𝐭)=×𝐯(𝐫,𝐭). Before averaging over tissues, we took the dominant direction of rotation of each tissue to correspond to positive vorticity. This direction was counterclockwise in 51.5% of tissues and clockwise in 49.5% of tissues, with a sample size of 68. With this convention, the vortex core always has positive vorticity. Accordingly, the outer region of the vortex exhibits negative vorticity (Figure 3C, see Figure 3—figure supplement 2 for kymographs and heatmaps of vorticity representative tissues), which corresponds to the counter-rotation that occurs when the central vortical flow transitions to the outer radial flow (Figure 3A, left). We define a characteristic vortex radius as the radial position of the center of the negative-vorticity region, which is ∼1 mm at 36 hr in small tissues (Figure 3C, black bars).

To analyze vortex dynamics across different tissues with varying vortex positioning, and to quantitatively capture the onset and strength of vortices, we calculated the enstrophy spectrum (q,t)=|ω~(𝐪,t)|2, where ω~(𝐪,t)=(d𝐫/A)ω(𝐫,t)ei𝐪𝐫 are the spatial Fourier components of the vorticity field ω(𝐫,t) (Alert et al., 2020). The enstrophy spectrum is the power spectral density of the vorticity field as a function of the wave-vector modulus q, and therefore provides a measure of the vortex intensity at a length scale 2π/q. The kymographs of the enstrophy spectrum show that most of the vortex’s intensity is found at a characteristic length scale of ∼1 mm (Figure 3—figure supplement 3).

For each tissue we characterized the maximal vortex strength by the maximum value of (q,t) as well as its associated wavelength 2π/q and time of occurrence. We represented these three quantities on a scatter plot, which shows that vortices in small tissues have generally higher intensity than those in large tissues (Figure 3D). Vortices in small tissues are also larger relative to tissue size, since the absolute size of vortices in small and large tissues is similar (Figure 3D). Furthermore, vortex strength peaks several hours later in small tissues than in large tissues (Figure 3D). We hypothesized that this difference is due to large tissues featuring a faster density increase than small tissues (Figure 1C). To test this hypothesis, we varied the initial cell density of small tissues and observed that the time of maximum vortex intensity decreases with increasing density (Figure 3E, Figure 3—figure supplement 3). These results prompted us to examine spatiotemporal cell density evolution.

Spatiotemporal dynamics of cell density

Given that cell density appears to affect vortex formation and is known to control contact inhibition of locomotion and proliferation (Schnyder et al., 2020), we explored the spatiotemporal evolution of cell density. Constructing average kymographs in the same way as for speed, radial velocity, and vorticity, we observe that the vortex region in the center of small tissues is accompanied by an unexpected local density minimum (Figure 4A). Strikingly, snapshots of small and large tissues reveal that large-scale vortices occur in low-density regions, regardless of location within the tissue (Figure 4—figure supplement 1). However, given that vortices in large tissues are often off-centered, the low-density region does not appear in their average kymograph of cell density (Figure 4A).

Figure 4 with 1 supplement see all
Spatiotemporal dynamics of cell density during epithelial expansion.

(A) Averaged kymographs of cell density for small (left, n = 11) and large (right, n = 9) tissues. Small tissues develop a central low-density region that persist more than 20 hr. (B) Cell density, ρ, at the center of large tissues increases gradually, while cell density at the center of small tissues has non-monotonic evolution. (C) For different initial tissue sizes and densities, the evolution of the cell density, ρ, at the boundary zone converges to similar values at about 12 hr, which coincides with the end of the overshoot of edge radial velocity in Figure 1D. Center and boundary zones are defined as in Figure 2B. (D) Simulated evolution of cell densities obtained from the numerical solution of the continuity equation using the average radial velocity measurements vr(r,t) (Figure 2B) and a uniform and constant cell proliferation rate corresponding to a 16 h cell doubling time. In (B,C) the initial cell density ranges and number of replicates ρ1, ρ2, and ρ3 are the same as in Figure 1D.

To investigate the effects of initial conditions, we tracked the density evolution of the center and boundary zones across tissues with different starting densities and sizes, grouping initial densities into three ranges as before (Figure 4B and C). As with the average density in Figure 1C, the density monotonically increases in large tissues centers but is non-monotonic in small tissues. Notably, the cell density at the center of small tissues of different initial cell densities reach a common minimum during the 16–32 hr time period (Figure 4B), which includes the vortex onset time. At the boundary zone, the long-time evolution of the cell density is independent of initial tissue size and density (Figure 4C). This common long-time evolution is reached at about 12 hr (Figure 4C), which coincides with the time at which the edge radial velocity stabilizes upon the overshoot (Figure 1D).

To understand the unexpected transient density decrease at the center of small tissues, we sought to explain it as the result of combined advective transport based on the measured radial flow fields 𝐯r(𝐫,t) and homogeneous cell proliferation at a rate k(𝐫,t)=k0 throughout the tissue. To test this hypothesis, we solved the continuity equation for the cell density field ρ(𝐫,t),

(4) ρt=-(ρ𝐯)+k0ρ,

using the average radial velocity profiles vr(r,t) measured by PIV (Figure 2D), and a proliferation rate k0=1.04 h−1, which corresponds to a cell doubling time of 16 hr (Materials and methods). This minimal model recapitulates the major features of the evolving density profiles for both small and large tissues (compare Figure 4D with Figure 4A). Therefore, the unexpected formation of a central low-density region results from the combination of outward tissue flow and proliferation within the colony. However, further research is required to determine the biophysical origin of the non-monotonic density evolution. Moreover, having assumed a density-independent proliferation rate, our model predicts a cell density in the center of large tissues higher than the one measured at the end of the experiment, and it does not quantitatively reproduce the cell density profiles at the edge regions. These discrepancies suggest that more complex cell proliferation behavior is required to fully recapitulate the density dynamics in expanding cell monolayers.

Spatiotemporal dynamics of cell cycle

To better understand how tissue expansion affects cell proliferation, we analyzed the spatiotemporal dynamics of cell-cycle state. Our cells stably express the FUCCI markers, meaning that cells in the G0-G1-S phase of the cell cycle (referred to here as G1) fluoresce in red (shown as magenta), and cells in the S-G2-M phase of the cell cycle (referred to here as G2) fluoresce in green (Sakaue-Sawano et al., 2008). Additionally, immediately-post-mitotic cells do not fluoresce and appear dark. Small and large tissues are initially well mixed with green and magenta cells, confirming that cells are actively cycling throughout the tissue at the time of stencil removal (Figure 5—figure supplement 1). During tissue expansion, spatiotemporal patterns of cell-cycling behavior emerge (Figure 5A, Figure 5—video 1).

Figure 5 with 2 supplements see all
Coordinated spatiotemporal cell-cycle dynamics.

Transition from the G1 (magenta) to the G2 (green) phase of the cell cycle corresponds to DNA replication (during S phase). Subsequently, a cell proceeds to mitosis (M phase, dark), and eventually back to the G1 phase upon cell division. (A) Fluorescence images of the Fucci marker of cell-cycle state at the end of the experiment (46 hr) of representative small and large tissues overlaid with nuclei positions (gray). The boundary zone of both tissues has more cells in the G2 than in the G1 phase, along with a substantial proportion of dark cells (inset). Scale bars 1 mm. (B) Average kymographs (small, n = 5; large, n = 11) of cell-cycle-state fraction. In small tissues, a G1-dominated transition zone, which appears as a vertical magenta streak from 16 hr onward, is interposed between G2-dominated center and edge zones. While the size of small tissues from 30 to 46 hr matches that of large tissues from 0 to 16 hr (dashed boxes), cell-cycle states between these times are clearly distinct. (C) Fraction of cell-cycle states in the boundary zone. (D) Fraction of cell-cycle states in the center zone. Center and boundary zones are defined as in Figure 2. For C and D, n = 5 for small tissues and n = 11 for large tissues. (E) Scatter plot of density and speed, with color indicating the fraction of cells at G1 and G2, corresponding to each PIV pixel of the final timepoint of a representative small (left) and large (right) tissue.

To quantitatively investigate these cell-cycle patterns, we obtained the local fractions of G1, G2, and post-mitotic cells by evaluating cell cycle state for each cell nucleus (see Materials and Methods). We then overlaid kymographs of the G1 and G2 cell-cycle-state fractions (Figure 5B) and plotted the time evolution of G1, G2, and post-mitotic fractions together (Figure 5C,D). Immediately after stencil removal, we observe a cell division pulse in all tissues, which manifests in a decrease in G2 and increase in post-mitotic fraction (Figure 5C,D). After about 12 hr of tissue expansion, the boundary region becomes primarily populated by rapidly-cycling cells (Figure 5B,C), which results in a predominance of cells in this region that either have recently divided (post-mitotic, black) or are likely to divide soon (G2, green). The high numbers of post-mitotic cells indicate that cells in G1 rapidly proceed to mitosis. Given that the edge radial speed overshoots during the first 12 hr of tissue expansion (Figure 1D), future work is necessary to characterize the effect of cell cycling on edge motion at early stages of expansion.

In the central region of small tissues (Figure 5B left, D left), we observe cell-cycling dynamics similar to the boundary region. Thus, in the tissue-spanning vortex of small tissues, cells are also rapidly cycling. The fraction of cells in G1 only starts to increase at ∼40 hr (Figure 5D left), coinciding with the weakening of the vortex (Figure 3C left). In contrast, the center zone of large tissues undergoes strong cell-cycle arrest at the G1-G2 transition at about 30 hr, also coinciding with the weakening of the vortex in large tissues (Figure 5B right, D right). Cells already past G1 at this time continue to division and re-enter G1, evidenced by the steady increase in local fraction of G1 accompanied by a steady decrease in G2 after 30 hr. Similar cell-cycle arrests were previously reported both in growing epithelia (Streichan et al., 2014) and in spreading 3D cell aggregates (Beaune et al., 2014). Before the onset of cell-cycle arrest, the center of large tissues exhibits large-scale coordinated cell-cycling dynamics in the form of anti-phase oscillations, with peaks in G2 fraction accompanied by troughs in G1 fraction (Figure 5B right, D right).

Finally, we sought to link cell-cycle dynamics to the kinematics of tissue expansion by studying correlations between local measurements of cell cycle, cell speed, and cell density (Figure 5E). Here, each point represents one PIV window, with color indicating its average cell-cycle state. As expected, cell speed is negatively correlated with cell density. Further, in large tissues, the cell-cycle state transitions from G1-dominated to G2-dominated when cell density increases above ∼5000 cells/mm2 and cell speed falls below ∼12 μm/h (Figure 5E right). In this regime, the decrease of cell speed with increasing cell density bears similarities to previously-reported glass transitions and contact inhibition of locomotion (Angelini et al., 2011; Zimmermann et al., 2016; Garcia et al., 2015). Small tissues, by contrast, lack the G1-dominated, slow, high-density cell population (Figure 5E, left) found in the center of large tissues. Taken together, our findings emphasize that cell cycling, cell flow, and cell density patterns are inextricably linked and depend on the initial size of an expanding tissue.

Discussion

We began this study by asking how changes in initial size affect the long-term expansion and growth of millimeter-scale epithelia. By means of high spatiotemporal resolution imaging and precisely controlled initial conditions, our assays systematically dissected tissue expansion and growth from the overall boundary kinematics (Figure 1) to the internal flow patterns (Figures 2, 3 and 4) and cell-cycle dynamics (Figure 5). While we demonstrated that ‘small’ tissues increase in area relatively much faster than do ‘large’ tissues, our data suggest a surprising and stark decoupling of the outer and inner regions of an expanding epithelium. Notably, the behaviors of the edge zones are largely independent of tissue size, cell density, and history, while interior dynamics depend strongly on these factors.

Unexpectedly, the overall tissue growth and expansion dynamics (Figure 1) could be attributed to one dominant feature: these epithelia expanded at the same edge speed regardless of initial tissue size, shape, and cell density. The only exception is the major axes of ellipses, where the normal edge speed is smaller when the radius of curvature is rc<0.75 mm. This observation, combined with the fact that the velocity penetration length is 500 mm (Figure 2D), suggests that a tissue must be 1 mm in diameter for the tissue edge to move independently of bulk flows. As a result of this robust edge motion, the areal expansion rate of the tissue is dictated by its perimeter-to-area ratio. To further emphasize the decoupling of the boundary and internal dynamics of epithelia, consider that the key findings in Figure 1 neither predict nor depend upon the radically different internal dynamics we observed within ‘small’ and ‘large’ tissues. For instance, despite the roiling vortices occupying large portions of ‘small’ tissues and the pronounced, large-scale contact inhibition of ‘large’ tissues–two antithetical phenomena–no hints of these behaviors can be detected in the motion of the boundary.

Critically, the type and timing of internal dynamics are dictated not by the current size but by the expansion history of a given tissue. While a small tissue eventually expands to reach the initial size of a large tissue, it exhibits different internal dynamics from the large tissue at this size (Figure 6). This difference in internal dynamics is perhaps easiest to observe in spatiotemporal evolution of cell cycle (Figure 6D, Figure 5B dashed boxes); the small-tissue footprint from 30 to 46 hr closely matches the large-tissue footprint from 0 to 16 hr, but the cell cycle distribution during these time periods bears almost no similarities. This applies as well to other important bulk properties of the tissue (Figure 6A–C), as cell cycle is tightly linked to cell speeds and density (Figure 5E). For example, at equal current sizes, the center of initially-small tissues features high vorticity with decreasing cell speed whereas initially-large tissues exhibit low vorticity and increasing cell speed (Figure 6A,B). Respectively, at equal current sizes, while absolute cell densities in the tissue centers share some overlap, it is notable that the rate of density change at the tissue center is increasing faster in initially-small tissues than in initially-large tissues (Figure 6C). However, the most striking differences in cell density evolution occur not at equal current sizes but during the early stages of tissue expansion: whereas the cell density at the center of large tissues increases at all times, the center of small tissues features a marked density decrease between ∼8 and ∼24 hr (Figure 4A,B). Overall, while edge dynamics are stereotyped and conserved across different sizes, our findings suggest that initial tissue size impacts the bulk dynamics by altering the constraints under which the tissue grows. We expect that tissues with sizes between our two choices would exhibit similar edge dynamics and internal patterns that cross over between our small and large tissues.

Initial tissue size, rather than current tissue size, determines the internal dynamics of expanding epithelia.

Here, we quantify the internal state of the tissue in terms of the cell speed in tissue center (A), maximal vortex power (B), cell density in tissue center (C), and fraction of cells in the G1 phase of the cell cycle in tissue center (D). At late times, initially-small tissues reach radii that initially-large tissues had at early times. When they have the same current size (overlap region in between dashed lines), initially-small and initially-large tissues have distinct internal dynamics of cell migration and cell proliferation. The tissue center zone in A, C, and D was defined as in Figure 2.

The vortices are a particularly striking example of such size- and history-dependent internal patterns (Figure 3, Figure 6B). Our active fluid model suggests that the vortices emerge from a dynamical instability of the tissue bulk, which occurs when the tissue reaches a critical size. Thus, whereas the instability itself is a bulk phenomenon independent of the tissue edge, edge-driven expansion allows small tissues to reach the critical size that triggers the instability. In addition, our data suggest a strong correlation between vortex formation and the development of non-monotonic density profiles. Not only did small tissues exhibit co-occurrence of vortices with density decreases in the tissue center, but also off-center vortices in large tissues always co-localized with a local density decrease (Figure 4—figure supplement 1). Our model does not currently describe cell density, and hence cannot explain the relationship between vortex formation and local density decreases. Thus, our experimental findings call for the development of more detailed models that couple cell density to both the velocity and the polarity fields, accounting for how density gradients influence cell polarization (Alert and Trepat, 2020).

The pronounced decoupling between boundary and internal dynamics in epithelia confers stability to the overall expansion of the tissue, making it robust to a wide range of internal perturbations. From the perspective of collective behavior, we speculate that such robust boundary dynamics may be beneficial in a tissue such as an epithelium whose teleology is to continuously expand from its free edges to sheath organ surfaces. Further, the ability to accurately predict epithelial expansion with a single parameter, the edge speed, will have practical uses in experimental design and tissue-engineering applications. Finally, given that many of the phenomena presented here only occurred due to the millimetric scale of our unconfined tissues and the long duration of the experiments, our results showcase the value of pushing the boundaries of large-scale, long-term studies on freely-expanding tissues.

Materials and methods

Key resources table
Reagent type (species)
or resource
DesignationSource or referenceIdentifiersAdditional information
Cell line (canine)MDCK-FUCCIStreichan et al., 2014.
Wildtype: ECACC-00062107
N/ARFP signal is G1, GFP is G2.
OtherDMEM - low glucoseSigma-Aldrich, IncCat.D5523
Otherfetal bovine serumAtlanta BiologicalsCat.S11550
OtherSilicone stencil material,
250 µm thick.
Stockwell ElastomericsCat.HT6240-40DTissue patterning material
Software, algorithmFIJINIH ImageJ Project
Software, algorithmMATLABMathworks, Inc2019AAll code compatible
back to at least 2015B
Software, algorithmMachine Learning ToolsLaboratory code LaChance and Cohen, 2020—see ReferencesFull code on GitHub

Cell culture

Request a detailed protocol

All experiments were performed with MDCK-II cells expressing the FUCCI cell-cycle marker system as received from: Streichan et al., 2014. After treatment with Mycoplasma Removal Agent (MPI Biological), cells tested negative for mycoplasma (MycoProbe, R and D Systems). We cultured cells in MDCK media consisting of low-glucose (1 g/L) DMEM with phenol red (Gibco, USA), 1 g/L sodium bicarbonate, 1% streptomycin/penicillin, and 10% FBS (Atlanta Biological, USA). Cells were maintained at 37°C and 5% CO2 in humidified air.

Tissue patterning

Request a detailed protocol

We coated tissue-culture plastic dishes (BD Falcon, USA) with type-IV collagen (MilliporeSigma, USA) by incubating 150 μL of 50 μg/mL collagen on the dish under a glass coverslip for 30 min at 37°C, washing three times with deionized distilled water (DI), and allowing the dish to air-dry. We then fabricated silicone stencils with cutouts of desired shape and size and transferred the stencils to the collagen coated surface of the dishes. Stencils were cut from 250 μm thick silicone (Bisco HT-6240, Stockwell Elastomers) using a Silhouette Cameo vinyl cutter (Silhouette, USA). We then seeded the individual stencils with cells suspended in media at 1000 cells/mL. Suspended cells were concentrated at 2.25×106 cells/mL and pipetted cells into the stencils at the appropriate volume. Care was taken not to disturb the collagen coating with the pipette tip. To allow attachment of cells to the collagen matrix, we incubated the cells in the stencils for 30 min in a humidified chamber before flooding the dish with media. We then incubated the cells for an additional 18 hr to allow the cells to form monolayers in the stencils, after which the stencils were removed with tweezers. Imaging began 30 min after stencil removal. Media without phenol red was used throughout seeding and imaging to reduce background signal during fluorescence imaging.

Live-cell time-lapse imaging

Request a detailed protocol

All imaging was performed with a 4X phase contrast objective on an automated, inverted Nikon Ti2 with environmental control (37°C and humidified 5% CO2) using NIS Elements software and a Nikon Qi2 CMOS camera. Phase contrast images were captured every 20 min, while RFP/GFP channels were captured every 60 min at 25% lamp power (Sola SE, Lumencor, USA) and 500 ms exposure time. No phototoxicity was observed under these conditions for up to 48 hr. Final images were composited from 4 × 4 montages of each dish using NIS Elements.

Tissue edge radial velocity

Request a detailed protocol

Tissues were segmented to make binary masks using a custom MATLAB (Mathworks) script. Tissue edge radial velocity was measured from the binary masks within more than 200 discrete sectors of the tissue; the edge radial velocity of all sectors were averaged to arrive at the tissue average edge radial velocity. Radial velocity at each sector was calculated for each timepoint as the rate of change of the average extent of the boundary pixels of the sector, utilized a rolling average of 3 timepoints (1 hr) to account for capture phase offsets resulting from capturing phase and fluorescence images at different frequencies. Sectors originated from the center of each tissue at the initial timepoint and were ∼20 µm wide at the edge of the tissue at the starting point.

Radius of curvature for the major and minor axes of elliptical tissues

Request a detailed protocol

Curvature at the major and minor axes of growing tissues was approximated at each time-point by fitting an ellipse to the tissue footprint and taking the radius of curvature at the minor and major axes as b2/a and a2/b, respectively, where a is the major semi-axis length and b is the minor semi-axis length.

Statistical tests and goodness of fit

Request a detailed protocol

Normalized χ2 values in Figure 1E were calculated as 1Ni=1N(ui-μi)2σi2, where N is the number of time-points in the curve, ui are the model predictions, and μi and σi are the mean and standard deviation of the measured values, respectively. With these definitions, a fit with χ2<1 is good.

The P-value in Figure 3D was calculated using a Mann-Whitney U test, and the two-tailed p-value of p<10-4 indicates that the large and small vortex power data indeed come from different populations.

Cell counts

Request a detailed protocol

The FUCCI system contains a period after M-phase where cells go dark, making FUCCI unreliable for cell counting. Instead, we developed and trained a convolutional neural network to reproduce nuclei from 4X phase contrast images using our in-house Fluorescence Reconstruction Microscopy tool (LaChance and Cohen, 2020) . The output of this neural network was then segmented in ImageJ to determine nuclei footprints and centroids.

Tissue PIV and density measurements

Request a detailed protocol

Tissue velocity vector fields were calculated from 2 × 2 resized phase contrast image sequences using the free MATLAB package PIVLab (Thielicke and Stamhuis, 2014) with the FFT window deformation algorithm. We used a 1 st pass window size of 64 × 64 pixels and second pass of 32 × 32 pixels, with 50% pixel overlaps. This resulted in a 115 × 115 μm window. The window size was chosen to be smaller than the velocity-velocity correlation length but large enough to enable fast computation of PIV fields for many tissues. As seen in Figure 2—figure supplement 1, using a window size of 57 × 57 μm, which contains only a few cells, yields higher resolution velocity fields but does not qualitatively affect the measured speed and radial velocity. We focus on large-scale features of the velocity field, which are not affected by choosing a smaller PIV window size.

Local density was also calculated for each PIV window by counting the number of approximate nucleus centroids in that window. Data from PIV were smoothed in time with a moving average of 3 time points centered at each timepoint as before.

Average kymographs

Request a detailed protocol

First, we constructed kymographs for individual tissues using distance from the tissue center as the spatial index for each measurement window corresponding to a kymograph pixel. We did not plot kymograph pixels for which more than 95% of the measurements at that distance were beyond the tissue footprint. We then averaged the individual tissue kymographs, aligning by the centers.

Trajectory colorization

Request a detailed protocol

We first generated a plot of all relevant trajectories (Tinevez et al., 2017) colorized randomly in grayscale using a custom MATLAB (Mathworks) script. We then used the Fiji plugin OrientationJ on this plot to colorize the resulting image according to orientation (Püspöki et al., 2016).

Cell density simulation

Request a detailed protocol

To test whether the observed spatiotemporal evolution of density ρ(r,t) could be explained by flow of material (rather than divisions, extrusions, and cell death), we solved the continuity equation for a homogenous tissue in a circular geometry with spatiotemporal evolution of average radial velocity vr(r,t) as measured from PIV in experiments ( Figure 2B). The continuity equation is

(5) ρt=-𝐣+k0ρ,

where a homogeneous cell proliferation rate k0=1.04h-1 is assumed throughout the tissue, which corresponds to the cell doubling time of 16 hr. The current density is 𝐣=ρ𝐯𝐫-Dρ, where we included a diffusion term with a small diffusion constant D=0.22mm2/h for numerical stability.

The continuity Equation (5) was discretized using the finite volume method (Eymard et al., 2000), which is briefly summarized below. The tissue domain was divided into an inner circle Ω0 of radius r1/2=12Δr and circular annuli Ωi with inner radii ri-1/2=(i-12)Δr and outer radii ri+1/2=(i+12)Δr, respectively, where i=1,2,3, and Δr=115μm corresponds to the width of 1 window in the PIV analysis. The continuity Equation (5) was then integrated over the inner circle Ω0 and circular annuli Ωi as

(6a) 1A00r1/2(2πrdr)ρt=1A00r1/2(2πrdr)[j+k0ρ],
(6b) 1Airi1/2ri+1/2(2πrdr)ρt=1Airi1/2ri+1/2(2πrdr)[j+k0ρ],

where A0=πr1/22 is the area of the inner circle Ω0 and Ai=πri+1/22-πri-1/22 is the area of the circular annulus Ωi. The integrals in Equation (6a, b) can be approximated as

(7a)ρ(0,t)t=2πA0r1/2j(r1/2,t)+k0ρ(0,t),(7b)ρ(ri,t)t=2πAi[ri+1/2j(ri+1/2,t)ri1/2j(ri1/2,t)]+k0ρ(ri,t).

Here, density profiles ρ(ri,t) are evaluated at ri=iΔr for all i=0,1,2,. Current densities are evaluated as j(ri+1/2,t)=ρ(ri+1/2,t)vr(ri+1/2,t)-D[ρ(ri+1,t)-ρ(ri,t)]/Δr for all i=0,1,2,, where ρ(ri+1/2,t)=[ρ(ri,t)+ρ(ri+1,t)]/2 and vr(ri+1/2,t)=[vr(ri,t)+vr(ri+1,t)]/2. Density profiles ρ(ri,t) were then obtained by integrating Equation (7) with the forward Euler method using a time step Δt=20 min to align with experimental data collection of radial velocity profiles vr(ri,t) from Figure 2B. The initial conditions were ρ(ri,0)=2700 cells/mm2 for ri<rtissue and ρ(ri,0)=0cells/mm2 for ri>rtissue, where rtissue is the radius of tissue at the beginning of experiment. For comparison with experimental data (see Figure 4), we thresholded the kymographs of simulated density at 100cells/mm2, which corresponds to much lower density than a confluent tissue.

Relating local cell density to vortex centers

Request a detailed protocol

For panels (E) and (F) in Figure 4—figure supplement 1, we applied a Fourier low-pass filter on vorticity fields, retaining only large-scale vorticity fluctuations (with wavelengths longer than 1 mm). We excluded the tissue edge region (500 μm from the boundary) that is outward polarized and does not exhibit vortical flows. Each point in panels (E) and (F) corresponds to a point in the filtered vorticity field, plotted against the cell density in that point.

Cell cycle analysis

Request a detailed protocol

The Fucci system consists of an RFP and GFP fused to proteins Cdt1 and Geminin, respectively (Sakaue-Sawano et al., 2008). Cdt1 levels are high during G1 and low during the rest of the cell cycle, while Geminin levels are high during the S, G2, and M phases (Sakaue-Sawano et al., 2008; Streichan et al., 2014). After capturing the appropriate fluorescence images, preprocessing was implemented identically for GFP and RFP channels to normalize channel histograms. To determine local cell cycle fraction, we determined the median value of RFP and GFP signal for each cell nucleus and manually selected thresholds for RFP and GFP signals separately to classify cell cycle for each cell as G0-G1-S (RFP above threshold), S-G2-M (RFP below threshold and GFP above threshold), or postmitotic (RFP and GFP below threshold). Local cell cycle fraction of each state could then be easily computed for each PIV pixel. Note that S phase (both RFP and GFP signals above threshold) did not prove to be a reliable feature for segmentation.

Code and data availability

Request a detailed protocol

Data for representative small, large, and ellipse tissues (Heinrich et al., 2020) and analysis Matlab scripts (Heinrich, 2020) have been made available (copy archived at https://github.com/elifesciences-publications/FreelyExpandingTissues).

References

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
    Finite Volume Methods
    1. R Eymard
    2. T Gallouet
    3. R Herbin
    (2000)
    In: J. L Lions, P Ciarlet, editors. Handbook of Numerical Analysis, 7. Universite Pierre et Marie Curie. pp. 713–1020.
    https://doi.org/10.1016/S1570-8659(00)07005-8
  21. 21
  22. 22
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
  48. 48
  49. 49
  50. 50
  51. 51
  52. 52
  53. 53
  54. 54
  55. 55
  56. 56
  57. 57
  58. 58
  59. 59
  60. 60
  61. 61
  62. 62
  63. 63
  64. 64
  65. 65
    Cellular-Pathologie
    1. R Virchow
    (1855)
    Archiv Für Pathologische Anatomie Und Physiologie Und Für Klinische Medicin 8:3–39.
    https://doi.org/10.1007/BF01935312
  66. 66
  67. 67

Decision letter

  1. Jody Rosenblatt
    Reviewing Editor; King's College London, United Kingdom
  2. Didier YR Stainier
    Senior Editor; Max Planck Institute for Heart and Lung Research, Germany
  3. Alexandre Kabla
    Reviewer; University of Cambridge, United Kingdom

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

This study investigates how a cellular colony's initial size affects its expansion through cell migration and proliferation, which become decoupled from its original size. This analysis highlights the importance of emerging vortices within the colonies and how anisotropy can affect tissue boundaries within cellular monolayers. While these studies were conducted in flat tissue culture monolayers, they may impact how cellular patterns in tissues may arise.

Decision letter after peer review:

Thank you for submitting your article "Size-dependent patterns of cell proliferation and migration in freely-expanding epithelia" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Didier Stainier as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Alexandre Kabla (Reviewer #3).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

We would like to draw your attention to changes in our revision policy that we have made in response to COVID-19 (https://elifesciences.org/articles/57162). Specifically, we are asking editors to accept without delay manuscripts, like yours, that they judge can stand as eLife papers without additional data, even if they feel that they would make the manuscript stronger. Thus the revisions requested below only address clarity and presentation.

Summary:

The manuscript by Heinrich et al. investigates how a tissue's initial size affects its expansion through cell migration and proliferation, which become decoupled from its original size. This analysis highlights the importance of emerging vortices within the colonies and how anisotropy can affect tissue boundaries within cellular monolayers. Overall, all three reviewers found your paper to be well done and of potential interest to a variety of scientists and felt your approach was insightful and quite thorough.

Revisions:

Together, the reviewers converged on a few points that we felt would not require necessarily any extra experimentation but could be obtained from your current data. These points focus on:

1) Study more systematically the relationship between curvature at the boundary and radial velocity, comparing round and oval shapes. Furthermore, how large does the tissue need to be for the radial velocity to be constant?

2) Directly compare round tissues that have the same current size but originate from different original sizes to determine if differences derive from original size or current size.

3) Please add more statistical analysis, instead of generalised statements.

4) Please add more methods regarding PIV.

5) Please add more discussion regarding the limitations of the model.

While these are the points, we would like you to address, I've not deleted the individual responses, as they may be more precise in addressing how you will do this, given I did not serve as a reviewer, just an editor here.

Reviewer #1:

This manuscript by Heinrich et al. explores how the initial size of a tissue affects its expansion through cell migration and proliferation. In their system, cell dynamics and proliferation at the outer edge of the tissue are decoupled from its original size. However, the initial size of the culture has an effect on how the inner core of the tissue behaves. The manuscript is well written, and the work is interesting and original with potential implications for development. However, I believe that some additional work is needed to make the paper suitable for eLife:

1) Although for the most part the data are nicely quantified, for a journal of this category statistical tests should be used to draw conclusions. For example,: “This model fits our data well”. This is qualitative and the authors should provide some objective metric of the goodness of fit. Same concern in subsection “Spatiotemporal dynamics of migration speed and radial velocity” were the authors use qualitative words like “similar” or “practically identical” rather than statistical tests to compare their data. This is a general comment that applies to many other comparisons established in the text, for example those related to Figure 3D in subsection “Vortex kinematics” or the correlation in Figure 3E.

2) In Figure 1—video 3, it seems that the oval tissues grow anisotropically with areas of greater curvature growing slower. I do appreciate that small and large tissues have different curvatures yet the same radial velocity, but what is the radial velocity for the experiments shown in 1E? Is it the same across tissues with different shapes? Is it homogeneous within tissues with variable curvature (higher aspect ratio)? Given that video, could changes in curvature (or other factors) affect the speed of migration and consequently the expansion of those tissues?

3) I think this is briefly mentioned in the Discussion but it is important to directly compare the behaviour of large and small tissues not only over time but also when their sizes are equivalent to demonstrate that any effect arises from the original size of the tissues and not from their current size. There should be a thorough analysis of this in the main Results section.

4) It would be good to back up the statement that vortices nucleate in low-density regions regardless of their location, with some quantifications rather than a couple of representative examples.

Reviewer #2:

In this study, Heinrinch et al. analyzed population dynamics of epithelial cells. They focused on the role of colony size, cell cycle and tissue expansion. Overall, they provide an elegant and interesting experimental study to describe the expansion of cellular monolayers. I could recommend the publication of the paper since it would be of interest for a large community. In particular, the emergence of vortices inside the colonies is very interesting and their results could shed new light on the formation of tissue boundaries within cellular monolayers. The analysis of the vortices is very complete and rigorous, and the cell cycle study is well addressed with clear quantification. I previously reviewed this paper for another journal, and I think that the authors made some efforts to improve their manuscript. In particular, the modeling part provides an interesting new piece of work. I still have a few comments to clarify and strengthen some of the points presented here before publication.

a) The role of cell cycle and expansion should be more discussed. In particular, the role of colony size and cell cycle remains largely unexplained. Experiments using synchronized cells or blocking proliferation could be helpful.

b) Additional experiments using drugs blocking cell division, acto-myosin contractility, actin polymerisation and perturbing cell-cell junctions may be useful to provide a deeper understanding.

c) It would be interesting to provide data on single cell shapes (anisotropy, spreading…) for small and large colonies over space and analyze the potential changes of cellular shapes as a function of their position in the monolayers. Do the authors observe changes in the cell shape anisotropy on small and large tissues? Along this line, the analysis of height of the colonies in the different configurations would be helpful. Do cells exhibit any changes in their apico-basal polarity depending on the size of the tissue and their position that could explain the different behaviors?

d) The sudden decrease in radial velocity could be due to cell division since the duration of cell cycle is about 16 hours in MDCK. It is surprising that this time of inflection is constant for both small and large colonies with varying density. Any explanation? Have the authors tried this for large colonies of varying density as they did for small ones?

e) Could the authors check if the propagation length is constant in both small and large tissues? How does it compare with the correlation length of MDCK cells?

Reviewer #3:

The paper from Heinrich et al. presents a rich analysis of the expansion and proliferation of epithelia. The authors use a commonly studied model system for wound healing and collective migration, MDCK epithelia on a flat substrate, and carefully analyse the behaviour of large (millimetre size) unconfined patches of cells. They map the displacement field of the cells across the tissue as well as monitor the cell cycle. This is to my knowledge the first paper to combine such analysis in a systematic manner at a large scale. They report novel behaviours, such as the emergence of different patterns of behaviours in the rim (radial movements) and in the bulk (rotational movements). More importantly, they clearly demonstrate that size matters in the way such partition of behaviour occurs, and in how this patterning evolves over time. A model provides a physical interpretation of some of these observations, although not all aspects of the work are well explained at this point, as mentioned by the authors themselves. This is a very nice paper which will prompt further discussions (and work) in the community. I recommend its publication.

I have very few concerns about the work. The paper is generally well written and a pleasure to read. Clearly a lot of work has been put into it, and although I would want to know more about certain aspects out of curiosity, it would be unfair to expect this to be included here.

1) The paper is specifically about large "millimetre-scale" epithelia, and I was wondering how large they need to be to match this criterion? When is it too small and what sets this length scale? I would expect that small patches would not be able to expand in size with the same radial velocity as the large populations. Do you have any experimental evidence about this?

2) Another point I would recommend clarifying concerns the calculation of the speed (e.g. on Figure 2). The PIV method section indicates what the spatial resolution is (~100µm by 100µm) but the temporal resolution could be made more explicit. I would expect the results to depend on these length-scales and timescales. How were these windows selected? Could data be included to show that the conclusions are robust with respect to this choice?

3) The modelling presented in the paper does a great job summarising part of the results. However, more could be said about what it fails to reproduce and why.

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "Size-dependent patterns of cell proliferation and migration in freely-expanding epithelia" for further consideration by eLife. Your revised article has been evaluated by Didier Stainier (Senior Editor) and a Reviewing Editor.

The manuscript has been improved but there are some remaining issues that need to be addressed before acceptance, as outlined below:

1) Regarding the relationship between curvature and edge speed, the authors find that the normal velocity is not constant along the ellipse edge, but they use a constant velocity in their model. They should address this briefly in the Results section and explain why they are able to make the assumption of constant velocity in the model. For example, by measuring if the radius of curvature along the ellipses is mainly greater than 1mm (which I imagine it is) and therefore a constant velocity can be assumed as shown in Figure 1—figure supplement 2B.

2) Since a major claim is that the behaviour of the tissues was all due to their original size rather than their current size, it is important to directly compare the behaviour of tissues that were originally large and small at the time point. Is the cell density significantly different for large and small tissues when they are the same size? It is hard to see from the figure. Also, Figure 6 is entirely addressed in the Discussion section rather than in the results, despite its importance. Why is it not in the Results section?

https://doi.org/10.7554/eLife.58945.sa1

Author response

Revisions:

Together, the reviewers converged on a few points that we felt would not require necessarily any extra experimentation but could be obtained from your current data. These points focus on:

1) Study more systematically the relationship between curvature at the boundary and radial velocity, comparing round and oval shapes. Furthermore, how large does the tissue need to be for the radial velocity to be constant?

We performed new analyses specifically relating curvature to normal expansion velocity including our elliptical tissue data. These data are seen in Figure 1—figure supplement 2. They show normal velocity as a function of curvature and demonstrate that the normal velocity is independent of curvature except in the extreme curvature cases on the highest aspect ratio ellipses. A fuller description can be found in paragraph four of the Results section, and we also provide further context in paragraph two of the Discussion section. With regards to the minimum tissue size needed for constant radial velocity, our analyses suggest tissue diameters of at least 1 mm given that we calculate the stable boundary zone is ~500 µm. This distinction has been added to the Discussion section.

2) Directly compare round tissues that have the same current size but originate from different original sizes to determine if differences derive from original size or current size.

This is an important consideration and insightful question that frequently comes up during presentations of this work, so we devoted an additional figure (Figure 6) to more quantitatively explore this. Our new analyses specifically explore the question of tissues of similar current size but different initial sizes with respect to migration speed, vortex power, cell density, and cell cycle (Figure 6). These new analyses further highlight the importance of history and context in understanding tissue behaviors and make it clear that initial tissue size and mechanical history, rather than current size, drives tissue behavior. A fuller discussion can be found in paragraph three of the Discussion section.

3) Please add more statistical analysis, instead of generalised statements.

Proper statistics are critical, so we have conducted new statistical analyses for all comparative studies. No prior claims have changed, but we now have more quantitative standards to support them. Specific analyses are as follows. Model goodness-of-fit is now evaluated using Chi-squared analysis, which can be seen in Figure 1E, and interpreted as Chi-squared values < 1 representing reasonable fits. Comparison of large and small tissue speeds and radial velocities were assessed by comparing the difference between the datasets compared to their respective standard deviations, as seen in Figure 2D. T-testing via the Mann-Whitney test was performed to compare distributions for vorticity (Figure 3D) and highlight the strongly statistically significant difference between them (p < 0.0001). Further analysis of the spatial correlation between density and vorticity resulted in a new panel (Figure 4—figure supplement 1E, F) displaying the relationships more quantitatively and across the whole dataset. All analyses are discussed in the Materials and methods section.

4) Please add more methods regarding PIV.

We appreciated this common as selection and assessment of PIV parameters is non-trivial and deserves more detailed discussion to aid in reproducibility. To address this, we carefully evaluated the effects of PIV parameters on our results and verified that changing the interrogation box did not appreciably alter the large-scale features or structures of the flow fields. That the larger PIV window is sufficient is especially useful because it can dramatically reduce computation time for large times relative to a smaller window size. These data are now presented in Figure 2—figure supplement 1, and the approach is described in the Materials and methods section.

5) Please add more discussion regarding the limitations of the model.

We discussed the limitations of the model in paragraph four of the Discussion section. We now explain that our model does not account for the cell density field, whereas our data show that vortex formation co-occurs with low-density regions. Again, our model is not intended to be exhaustive, but to capture key mechanistic details with as few parameters as possible. Our data and simulation results therefore call for the development of more detailed models to capture the relationship between vortex formation and cell density.

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

The manuscript has been improved but there are some remaining issues that need to be addressed before acceptance, as outlined below:

1) Regarding the relationship between curvature and edge speed, the authors find that the normal velocity is not constant along the ellipse edge, but they use a constant velocity in their model. They should address this briefly in the Results section and explain why they are able to make the assumption of constant velocity in the model. For example, by measuring if the radius of curvature along the ellipses is mainly greater than 1mm (which I imagine it is) and therefore a constant velocity can be assumed as shown in Figure 1—figure supplement 2B.

This is important to clarify, so we have the following to address it. As noted, high curvature only applies to a small region of the 1:8 ellipses, which will blunt over time as it grows, so our constant velocity model still fits the overall area expansion data well.

“Such high curvatures are concentrated around the major axes of our elliptical tissues. However, most of the tissue edge has a smaller curvature, and therefore advances at a curvature-independent speed. Further, even high curvature regions blunt due to expansion over time (see Figure 1—video 3). As a result, our model with a single edge speed 𝑣𝑛 ≃ 29.5 𝜇m/h is sufficient to capture the area expansion of both circular and elliptical tissues (Figure 1E).”

2) Since a major claim is that the behaviour of the tissues was all due to their original size rather than their current size, it is important to directly compare the behaviour of tissues that were originally large and small at the time point. Is the cell density significantly different for large and small tissues when they are the same size? It is hard to see from the figure.

It is true that the difference in density is less apparent in Figure 6C, and we have added additional discussion specifically exploring this point and emphasizing both the rate of change of density and the fact that the primary density phenotypes occur in the early stages of growth where smaller tissues experience a drop in density while larger tissues experience a monotonic increase. These are summarized below.

“…at equal current sizes, while absolute cell densities in the tissue centers share some overlap, it is notable that the rate of density change at the tissue center is increasing faster in initially-small tissues than in initially-large tissues (Figure 6C). However, the most striking differences in cell density evolution occur not at equal current sizes but during the early stages of tissue expansion: whereas the cell density at the center of large tissues increases at all times, the center of small tissues features a marked density decrease between ∼ 8 and ∼ 24 h (Figure 4A,B).”

Also, Figure 6 is entirely addressed in the Discussion section rather than in the results, despite its importance. Why is it not in the Results section?

We felt strongly that Figure 6 belongs in the Discussion section because it represents a compilation and comparison of data already presented in the prior figures in different forms. In other words, Figure 6 does not show new measurements but rather reanalyzes the data in all the previous figures in a way that showcases history effects. Thus, it is a figure that very much engages with, and emerges from, the points raised in the Discussion section relating to history effects, putting our results into a broader context. Therefore, we feel that moving Figure 6 to the Results section would weaken its impact and the Discussion section overall. That said, we can move it to the Results if necessary.

https://doi.org/10.7554/eLife.58945.sa2

Article and author information

Author details

  1. Matthew A Heinrich

    Department of Mechanical and Aerospace Engineering, Princeton University, Princeton, United States
    Contribution
    Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-9041-5554
  2. Ricard Alert

    1. Lewis-Sigler Institute for Integrative Genomics, Princeton University, Princeton, United States
    2. Princeton Center for Theoretical Science, Princeton University, Princeton, United States
    Contribution
    Formal analysis, Validation, Methodology, Writing - original draft, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1885-9177
  3. Julienne M LaChance

    Department of Mechanical and Aerospace Engineering, Princeton University, Princeton, United States
    Contribution
    Data curation, Software, Formal analysis, Investigation, Visualization, Methodology
    Competing interests
    No competing interests declared
  4. Tom J Zajdel

    Department of Mechanical and Aerospace Engineering, Princeton University, Princeton, United States
    Contribution
    Formal analysis, Investigation, Visualization, Methodology, Writing - original draft
    Competing interests
    No competing interests declared
  5. Andrej Košmrlj

    1. Department of Mechanical and Aerospace Engineering, Princeton University, Princeton, United States
    2. Princeton Institute for the Science and Technology of Materials (PRISM), Princeton University, Princeton, United States
    Contribution
    Formal analysis, Supervision, Funding acquisition, Methodology, Writing - original draft, Project administration, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-6137-9200
  6. Daniel J Cohen

    Department of Mechanical and Aerospace Engineering, Princeton University, Princeton, United States
    Contribution
    Conceptualization, Resources, Data curation, Software, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing - original draft, Project administration, Writing - review and editing
    For correspondence
    danielcohen@princeton.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-5819-1135

Funding

National Institutes of Health (1 R35GM133574-01)

  • Daniel J Cohen

Human Frontier Science Program (LT000475/2018-C)

  • Ricard Alert

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

Acknowledgements

DJC acknowledges the National Institutes of Health R35 GM133574-01. RA acknowledges support from the Human Frontiers of Science Program (LT000475/2018-C).

Senior Editor

  1. Didier YR Stainier, Max Planck Institute for Heart and Lung Research, Germany

Reviewing Editor

  1. Jody Rosenblatt, King's College London, United Kingdom

Reviewer

  1. Alexandre Kabla, University of Cambridge, United Kingdom

Publication history

  1. Received: May 15, 2020
  2. Accepted: August 18, 2020
  3. Accepted Manuscript published: August 19, 2020 (version 1)
  4. Version of Record published: September 17, 2020 (version 2)

Copyright

© 2020, Heinrich et al.

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

Metrics

  • 1,266
    Page views
  • 240
    Downloads
  • 0
    Citations

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

Download links

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

Downloads (link to download the article as PDF)

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

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

Further reading

    1. Developmental Biology
    2. Neuroscience
    Sarah McCallum et al.
    Research Article Updated

    The presence and identity of neural progenitors in the enteric nervous system (ENS) of vertebrates is a matter of intense debate. Here, we demonstrate that the non-neuronal ENS cell compartment of teleosts shares molecular and morphological characteristics with mammalian enteric glia but cannot be identified by the expression of canonical glial markers. However, unlike their mammalian counterparts, which are generally quiescent and do not undergo neuronal differentiation during homeostasis, we show that a relatively high proportion of zebrafish enteric glia proliferate under physiological conditions giving rise to progeny that differentiate into enteric neurons. We also provide evidence that, similar to brain neural stem cells, the activation and neuronal differentiation of enteric glia are regulated by Notch signalling. Our experiments reveal remarkable similarities between enteric glia and brain neural stem cells in teleosts and open new possibilities for use of mammalian enteric glia as a potential source of neurons to restore the activity of intestinal neural circuits compromised by injury or disease.

    1. Cancer Biology
    2. Developmental Biology
    Erik Oliemuller et al.
    Research Article Updated

    SOX11 is an embryonic mammary epithelial marker that is normally silenced prior to birth. High SOX11 levels in breast tumours are significantly associated with distant metastasis and poor outcome in breast cancer patients. Here, we show that SOX11 confers distinct features to ER-negative DCIS.com breast cancer cells, leading to populations enriched with highly plastic hybrid epithelial/mesenchymal cells, which display invasive features and alterations in metastatic tropism when xenografted into mice. We found that SOX11+DCIS tumour cells metastasize to brain and bone at greater frequency and to lungs at lower frequency compared to cells with lower SOX11 levels. High levels of SOX11 leads to the expression of markers associated with mesenchymal state and embryonic cellular phenotypes. Our results suggest that SOX11 may be a potential biomarker for breast tumours with elevated risk of developing metastases and may require more aggressive therapies.