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

Bridging the gap between single-cell migration and collective dynamics

  1. Florian Thüroff
  2. Andriy Goychuk
  3. Matthias Reiter
  4. Erwin Frey  Is a corresponding author
  1. Ludwig-Maximilians-Universität München, Germany
Research Article
Cite this article as: eLife 2019;8:e46842 doi: 10.7554/eLife.46842
10 figures, 1 table and 2 additional files

Figures

Illustration of the computational model with the pertinent simulation steps.

(A) Illustration of a small cell cohort that adheres to a surface ((x,y)-plane). The polarization field, ϵ, is defined on the contact surface with the adhesion plane. The magnitude of the polarization field, which is indicated by the colorbar in Figure (C), encodes the local strength of cell-substrate adhesions and emulates the local mass of force-generating (pushing) cytoskeletal structures. Cell-cell adhesions are indicated in red. (B) Cytoskeletal structures respond to external mechanical stimuli through reaction networks involving different feedback loops. We greatly simplify these complex processes into two prototypic feedback loops, which break detailed balance and drive cell migration, as follows. The polarization field induces membrane protrusions and inhibits retractions. In turn, protrusions increase the polarization field (positive feedback) and therefore the likelihood of further protrusive activity, while retractions decrease the polarization field (negative feedback). In the absence of mechanochemical signals, the polarization field approaches its rest state. (C) Zoom-in to a common boundary shared between the substrate contact areas of three cells (bounded by the red lines), each represented by a contiguous set of occupied grid sites (hexagons). Top left: The upper right corner of the lower left cell (source cell) initiates a protrusion event against a neighboring element in the cell to its right (target cell), as indicated by the arrow, in an attempt to displace it. The success of each such attempted elementary event depends on the balance between contractile forces, cytoskeletal forces, and cell adhesion. Top right: If the protrusion event is successful, then the levels of regulatory factors are increased (decreased) in integer steps, at all lattice sites inside the source (target) cell that lie within a radius R of the accepted protrusion event (as indicated by the plus and minus signs). Bottom right: During the course of one MCS, different levels of regulatory factors accumulate locally within each cell, with positive levels of regulatory factors (green plus signs) promoting a build-up of cytoskeletal structures, negative levels of regulatory factors (red minus signs) causing degradation of cytoskeletal structures, and neutral levels of regulatory factors (white zero signs) causing relaxation towards a resting state, as indicated in the lower left image. The color code indicates local levels of cytoskeletal structures, ϵ.

Figure 2 with 2 supplements
Cell shape and persistence of migration as a function of cell polarizability.

(A) Mean-squared displacement (MSD) for single-cell movements at different maximum cell polarity Δϵ (stiffness parameters κP= 0.060, κA= 0.18; average polarization field ϵ0= 225; signaling radius R= 5; cell-substrate dissipation D= 0; cell-substrate adhesion penalty φ= 0; cytoskeletal update rate μ=0.1; 100 independent simulations for each set of parameters). Single cells perform a persistent random walk, i.e. they move ballistically (MSDτ2) for ττp, and diffusively (MSDτ) for ττp. Inset: Normalized velocity auto-correlation function for the same parameters as in the main figure. (B) Persistence time of directed cell migration plotted as a function of maximum cell polarity Δϵ, and perimeter stiffness κP (area stiffness κA=0.18; average polarization field ϵ0=225; signaling radius R=5; cell-substrate dissipation D= 0; cell-substrate adhesion penalty φ= 0; cytoskeletal update rate μ=0.1100 independent simulations for each set of parameters). The persistence time of the random walk increases with increasing cytoskeletal polarity and decreasing perimeter elasticity. (C) Cytoskeletal polarity also controls cell shapes, with crescent cell shapes (long persistence times) being observed at large cytoskeletal polartities, and round cell shapes (short persistence times) at small cytoskeletal polarities. Color code: cell polarization; cf. color bar in Figure 1C. (D) Single cell speed plotted as a function of maximum cell polarity Δϵ, and perimeter stiffness κP. (E) Single cell aspect ratio plotted as a function of maximum cell polarity Δϵ, and perimeter stiffness κP. (F) Speed and persistence time of single cells are correlated with the cell aspect ratio.

Figure 2—figure supplement 1
Role of substrate dissipation for cell shape and motility.

(A–D) Role of substrate dissipation for cells of varying maximum cell polarity Δϵ. The aspect ratio l+/l (A), the speed v (B), and the persistence time τp (C) as a function of substrate dissipation D for a series of values for Δϵ/κP indicated in the graphs. (D) Cell persistence and cell speed are correlated with the aspect ratio of the cell. Only cells with an aspect ratio larger than 2 are motile. In the simulations, the substrate dissipation was varied in the interval D[0, 50], and the maximum cell polarity Δϵ{20,30,40,50,60}. Fixed parameters: average polarization field ϵ0=225; area elasticity κA=0.18; membrane elasticity κP=0.060; cytoskeletal update rate μ=0.1; cell-substrate adhesion penalty φ=0100 independent simulations for each set of parameters. (E–F) Role of substrate dissipation for cells of varying membrane stiffness κP. The aspect ratio l+/l (E), the speed v (F), and the persistence time τp (G) as a function of substrate dissipation D for a series of values for Δϵ/κP indicated in the graphs. (H) Cell persistence and cell speed are correlated with the aspect ratio of the cell. Only cells with an aspect ratio larger than 2 are motile. In the simulations, the substrate dissipation was varied in the interval D[0, 50], and the membrane elasticity κP{0.054,0.057,0.060,0.063,0.066}. Fixed parameters: average polarization field ϵ0=225; area elasticity κA=0.18; maximum cell polarity Δϵ=50; cytoskeletal update rate μ=0.1; cell-substrate adhesion penalty φ=0100 independent simulations for each set of parameters.

Figure 2—video 1
Single cell motility and shape for different maximum cell polarities (κP=0.060, R=5).
Figure 3 with 1 supplement
Migratory behavior of single cells as a function of the cell’s signaling radius R at different values for the maximal cytoskeletal polarity Δϵ.

(Stiffness parameters κP=0.060κA=0.18; average polarization field ϵ0=225; cell-substrate dissipation D=0; cell-substrate adhesion penalty φ=0; cytoskeletal update rate μ=0.1100 independent simulations for each set of parameters). (A) The persistence times of directed migration of single cells exhibit a pronounced maximum at an optimal signaling radius, which depends on cell polarizability. (B) The shapes of single cells exhibit a pronounced maximal elongation at an optimal signaling radius, which depends on cell polarizability. (C) The signaling radius critically determines the synchronicity of internal cytoskeletal remodeling processes. Small signaling radii frequently lead to transient formation of mutually independent lamellipodia at different positions around the cell perimeter, thereby interrupting persistent motion (reducing persistence times of directed migration). Large signaling radii lead to structurally stable front-rear polarization profiles across the entire cell body (long persistence times of directed migration). Color code: cell polarization; cf. color bar in Figure 1C. (D) The speed of single cells does not drop to zero even when their persistence time of directed migration vanishes. This indicates single cell rotations. (E) The inverse curvature of the cell trajectories as a function of the signaling radius. (F) Depending on whether a cell migrates along its long axis (top) or short axis (bottom), it has to move a different projected contour length. If each protrusion takes roughly the same amount of time, then migration along the long axis (top; cell has to move a smaller projected contour length) allows for greater cell speeds than migration along the short axis (bottom; cell has to move a larger projected contour length).

Figure 3—video 1
Single cell motility and shape for different signaling radii (Δϵ=60, κP=0.060).
Figure 4 with 6 supplements
Phases of collective motion.

(4-cell systems; confinement radius r0=30.6; area stiffness κA=0.18; average polarization field ϵ0=225; signaling radius R=5; cytoskeletal update rate μ=0.1; cell-cell adhesion B=0; cell-cell dissipation ΔB=12; cell-substrate dissipation D=0; cell-substrate adhesion penalty φ=0 (r<r0), φ (r>r0); 100 independent simulations for each set of parameters). (A) Characteristic observables of collective cell rotation at different values of the cell perimeter stiffness parameter κP: mean (|ω|) and standard deviation (σω) of the magnitude of the cell cluster's angular velocity, and the standard deviation of the cell perimeter (σP). The statistics of collective cell motion depends only on the ratio of maximum cell polarity, Δϵ, to cell contractility, κP (specific polarizability). (B) Representative angular trajectories and (C) cell shapes (color code represents cell polarization; cf. Figure 1C) for the different parameter regimes as described in the main text. The cellular dynamics in the different parameter regimes are shown in Figure 4—video 1Figure 4—video 2 and Figure 4—video 3.

Figure 4—figure supplement 1
Collective motion for varying number of cells at low polarizability.

(N-cell systems; confinement radius r0=234N; stiffness parameters κP=0.060, κA=0.18; average polarization field ϵ0=225; maximum cell polarity Δϵ=28; signaling radius R=5; cytoskeletal update rate μ=0.1; cell-cell adhesion B=0; cell-cell dissipation ΔB=12; cell-substrate dissipation D=0; cell-substrate adhesion penalty φ=0 (r<r0), φ- (r>r0); 100 independent simulations for each set of parameters). For this choice of parameters 4-cell populations rotate in the 1-phase. We observe a similar behavior here: the cell clusters rotate slowly and reorient frequently. (A) Characteristic observables of collective cell rotation at different values of the cell perimeter stiffness parameter κP: mean (|ω|) and standard deviation (σω) of the angular velocity magnitude of cell motion, and the standard deviation of the cell perimeter (σP). The black line corresponds to a power-law fit of the form |ω|N-k/2r0-k with the fitted exponent k8/3. (B) Representative angular trajectories and (C) cell shapes (color code represents cell polarization; cf. Figure 1) for the different parameter regimes as described in the main text. .

Figure 4—figure supplement 2
Collective motion for varying number of cells at intermediate polarizability.

(N-cell systems; confinement radius r0=234N; stiffness parameters κP=0.060, κA=0.18; average polarization field ϵ0=225; maximum cell polarity Δϵ=50; signaling radius R=5; cytoskeletal update rate μ=0.1; cell-cell adhesion B=0; cell-cell dissipation ΔB=12; cell-substrate dissipation D=0; cell-substrate adhesion penalty φ=0 (r<r0), φ- (r>r0); 100 independent simulations for each set of parameters). For this choice of parameters 4-cell populations rotate in the 2-phase. We observe a similar behavior here: highly correlated rotations with no changes in rotational direction. (A) Characteristic observables of collective cell rotation at different values of the cell perimeter stiffness parameter κP: mean (|ω|) and standard deviation (σω) of the angular velocity magnitude of cell motion, and the standard deviation of the cell perimeter (σP). The black line corresponds to a power-law fit of the form |ω|N-1/2r0-1. (B) Representative angular trajectories and (C) cell shapes (color code represents cell polarization; cf. Figure 1) for the different parameter regimes as described in the main text. .

Figure 4—figure supplement 3
Collective motion for varying number of cells at high polarizability.

(N-cell systems; confinement radius r0=234N; stiffness parameters κP=0.060, κA=0.18; average polarization field ϵ0=225; maximum cell polarity Δϵ=70; signaling radius R=5; cytoskeletal update rate μ=0.1; cell-cell adhesion B=0; cell-cell dissipation ΔB=12; cell-substrate dissipation D=0; cell-substrate adhesion penalty φ=0 (r<r0), φ- (r>r0); 100 independent simulations for each set of parameters). For this choice of parameters 4-cell populations rotate in the 3-phase. We observe a similar behavior here: highly correlated rotations. (A) Characteristic observables of collective cell rotation at different values of the cell perimeter stiffness parameter κP: mean (|ω|) and standard deviation (σω) of the angular velocity magnitude of cell motion, and the standard deviation of the cell perimeter (σP). The black line corresponds to a power-law fit of the form |ω|N-1/2r0-1. (B) Representative angular trajectories and (C) cell shapes (color code represents cell polarization; cf. Figure 1) for the different parameter regimes as described in the main text.

Figure 4—video 1
Collective rotations of 4 cells in the 1-phase (Δϵ=28Δϵ/κP=467).
Figure 4—video 2
Collective rotations of 4 cells in the 2-phase (Δϵ=50Δϵ/κP=833).
Figure 4—video 3
Collective rotations of 4 cells in the 3-phase (Δϵ=70Δϵ/κP=1167).
Figure 5 with 3 supplements
Expansion of a confluent epithelial cell sheet after removal of boundaries positioned at x=±175 for two different parameter settings.

(Stiffness parameters κP=0.12κA=0.18; average polarization field ϵ0=35; signaling radius R=2; cytoskeletal update rate μ=0.1; cell-cell adhesion B=12; cell-cell dissipation ΔB=0; cell-substrate dissipation D=0; cell-substrate adhesion penalty φ=0100 independent simulations for each set of parameters). (A–D) Tissue expansion for a migration-dominated setup without explicit cell growth and mitosis. (3300-cell system; maximum cell polarity Δϵ=30). (E–H) Tissue expansion at low density and cell polarizability for a cell sheet comprised of dividing cells. (Initially a 2500-cell system; maximum cell polarity Δϵ=10; growth time Tg=180; division time Td=20; size threshold for cell growth AT=1Aref, where Aref is the size of a solitary cell in equilibrium). (A, E) Snapshots of the polarization field ϵ; cf. Figure 5—video 1 and Figure 5—video 2. (B, F) Kymographs showing the cell density averaged over the y-direction and (top) final snapshots of the cell density profiles. (C, G) Kymographs showing the component σxx of the stress tensor averaged over the y-direction and (top) final snapshots of the stress profiles. (D, H) Kymographs showing the component vx of the cell velocities averaged over the y-direction and (top) final snapshot of the velocity profiles.

Figure 5—figure supplement 1
Monolayer expansion depends on dissipation and cell polarizability.

Cell monolayer expansion depends on the cell-cell dissipation, cell-substrate dissipation, and maximum cell polarity. (Initially a 2500-cell system; stiffness parameters κP=0.12, κA=0.18; average polarization field ϵ0=35; maximum cell polarity Δϵ=10; signalling radius R=2; cytoskeletal update rate μ=0.1; cell-cell adhesion B=12; cell-substrate adhesion penalty φ=0; growth time Tg=180; division time Td=20; size threshold for cell growth AT=1Aref, where Aref is the size of a solitary cell in equilibrium; 100 independent simulations for each set of parameters). (A,B) Cell monolayer expansion depends on the cell-cell dissipation ΔB (maximum cell polarity Δϵ=10; cell-substrate dissipation D=0). (A) Maximal monolayer extension and roughness (we exclude an initial time interval of 200 MCS because it takes at least that long for first daughter cells to appear). Inset: Relative roughness of the spreading monolayer relative to its size. (B) Time traces for selected values of ΔB. (C,D) Cell monolayer expansion depends on the cell-substrate dissipation D (maximum cell polarity Δϵ=10; cell-cell dissipation ΔB=0). (C) Maximal monolayer extension and roughness (we exclude an initial time interval of 200 MCS because it takes at least that long for first daughter cells to appear). (D) Time traces for selected values of D. (E,F) Cell monolayer expansion depends on the maximum cell polarity Δϵ (cell-cell dissipation ΔB=0; cell-substrate dissipation D=0). (E) Maximal monolayer extension and roughness (we exclude an initial time interval of 200 MCS because it takes at least that long for first daughter cells to appear). (F) Time traces for selected values of Δϵ/κP. .

Figure 5—video 1
Motility-dominated tissue dynamics.
Figure 5—video 2
Proliferation-dominated tissue dynamics.
Figure 6 with 2 supplements
Expansion of a confluent epithelial cell sheet after removal of boundaries positioned at x=±175 for two different parameter settings that produce rough tissue fronts.

(Initially a 2500-cell system; stiffness parameters κP=0.10, κA=0.18; average polarization field ϵ0=35; maximum cell polarity Δϵ=20; signaling radius R=5; cytoskeletal update rate μ=0.1; cell-cell adhesion B=5; cell-cell dissipation ΔB=10; cell-substrate dissipation D=0; cell-substrate adhesion penalty φ=0; growth time Tg=180; division time Td=20; 100 independent simulations for each set of parameters). (A–D) Tissue expansion at low density and cell polarizability for a cell sheet comprised of quickly dividing cells. (Size threshold for cell growth AT=1.05Aref, where Aref is the size of a solitary cell in equilibrium). (E–H) Tissue expansion at low density and cell polarizability for a cell sheet comprised of slowly dividing cells. (Size threshold for cell growth AT=1.10Aref, where Aref is the size of a solitary cell in equilibrium). (A, E) Snapshots of the polarization field ϵ; cf. Figure 6—video 1 and Figure 6—video 2. (B, F) Kymographs showing the cell density averaged over the y-direction and (top) final snapshots of the cell density profiles. (C, G) Kymographs showing the component σxx of the stress tensor averaged over the y-direction and (top) final snapshots of the stress profiles. (D, H) Kymographs showing the component vx of the cell velocities averaged over the y-direction and (top) final snapshot of the velocity profiles.

Figure 6—video 1
Weak monolayer roughening (fingering) in motility-dominated tissue with quick proliferation.
Figure 6—video 2
Strong monolayer roughening in motility-dominated tissue with slow proliferation.
Appendix 1—figure 1
Illustration of the various sets defining a cell and its environment.

Grid sites occupied by cell α, i.e. its domain 𝒟(α), are indicated in red colors. The cell’s membrane sites, (α), are indicated by the lighter red color, the cell’s immediate neighborhood, 𝒩(α), is indicated in gray. Elementary events involving cell α always involve one grid site in (α) and one grid site in 𝒩(α). For the hexagonal lattices used in this work, each grid site 𝐱k is surrounded by 6 nearest neighbors which we collectively denote by 𝒩k.

Appendix 1—figure 2
Distribution of regulatory factors on the basis of accepted elementary events.

For ease of reference, grid rows have been numbered from 1 to 10. Left (A): Solid black lines indicate cells’ membrane positions after acceptance of the respective elementary event; colors indicate cellular domains before the respective elementary event has been accepted (gray: substrate; shades of yellow: cells). Blue and red circular arcs (of radius R) delineate areas of local increase or decrease in the level of regulatory factors, respectively. The following elementary events are depicted: (i) lower cell retracts (two grid sites in row 2); (ii) lower cell protrudes (row 5); (iii) upper cell protrudes (row 10). In addition, the following elementary events occur across the cell-cell boundary: (iv) retraction of upper cell leads to rupture of cell-cell contacts (row 6, right event); (v) either the lower cell protrudes and pushes the upper cell or the upper cell retracts and pulls on the lower cell (row 6, left event). Specifically, event (v) entails mechanical signaling between the upper and lower cell and, therefore, affects the distribution of regulatory factors in both cells. Right (B): Identical copy of the left image (A). Colors indicate local levels of regulatory factors F (blue: F is positive; white: F is zero; red: F is negative; gray: substrate site). Note, in particular, that a substrate grid site has been inserted where cell rupture occurred (row 6, right grid site). The following cases can be distinguished: (i) Grid site xk lies in the zone of influence of only positive (blue circles) or negative (red circles) chemical feedback, in which case the level of regulatory factors is positive or negative, respectively (e.g. red grid sites in row 2, or blue grid sites in row 5). (ii) Grid site xk lies outside of any zone of influence, in which case the level regulatory factors is zero (e.g. white grid sites in row 2). (iii) Grid site xk lies in the zone of influence of equally many positive and negative feedbacks, in which case the level of regulatory factors remains zero (e.g. fourth grid site in row 4). (iv) Grid site xk lies in a zone of predominantly positive or negative feedback, in which case the level of regulatory factors is positive or negative, respectively (e.g. third grid site in row 4). Recall that only the sign of F is of significance to update the cells’ polarization field; cf. Equation S13.

Appendix 1—figure 3
Cell-cell adhesion.

(A) Adhesive energy contribution in a cyclic process, where a protrusion of source cell α against target cell β is followed by the inverse retraction event. Both events involve a third party cell γ, leading to net energy dissipation after the cyclic process has been completed. Protrusion: (i) Three pre-existing cell-cell contacts between β and γ are torn apart (red dashed contacts); (ii) three new contacts between α and γ are formed; (iii) the contact length between source cell α and target cell β increases by one unit of length. This implies Δadh(𝒯pro)=(3Bβ,γ3Bα,γBα,β). Retraction: (i) Three pre-existing cell-cell contacts between α and γ are torn apart (red dashed contacts); (ii) three new contacts between β and γ are formed; (iii) the contact length between source cell α and target cell β decreases by one unit of length. This implies Δadh(𝒯ret)=(3Bα,γ3Bβ,γ+Bα,β). Altogether, this leads to Δadh(cycl)=Δadh(𝒯pro)+Δadh(𝒯 ret)=(3(ΔB)α,γ+3(ΔB)β,γ)0, i.e. a (non-negative) dissipative contribution, whose magnitude depends on the dissipation matrix (ΔB)α,β=Bα,βBα,β0. (B) Shear viscosity due to cell-cell adhesion. Consider two rows of adhesive cells sliding past each other as indicated in the figure (left row of cells moves up by one grid site; colors indicate different cells). The associated adhesion energy change (per cell) reads Δadh/nc=2(BB)0, where nc denotes the number of cells sliding past each other, and where we assumed cells of like type, i.e. Bα,βB and Bα,βB (αβ). The condition B>B, Equation S15e, thus implies positive friction associated with cellular shear flows, whose magnitude is proportional to the number of cells sliding past each other. Note that this shear viscosity vanishes for B=B, i.e. for zero dissipation matrix.

Appendix 2—figure 1
Role of area stiffness κA for cell size and motility.

(A) The cell area increases linearly with 1/κA. The aspect ratio (B), speed (C) and persistence (D) of the cell decrease with increasing cell size. In the simulations, the area elasticity was varied in the interval κA[0.09, 0.18], and the membrane elasticity was chosen from κP{0.054,0.057,0.060,0.063,0.066}. Fixed parameters: average cell polarization field ϵ0=225; maximum cell polarity Δϵ=50; signaling radius R=5; cytoskeletal update rate μ=0.1; cell-substrate dissipation D=0; cell-substrate adhesion penalty φ=0.

Tables

Table 1
Source and parameter files used for each figure.

All source and parameter files are found in Source data 1.

FigureSimulation codeProcessing codeParameters
Figure 2CPM_NoDivisionTrajectoryAnalysisSinglesingle_Q
Figure 2—figure supplement 1 (A-D)CPM_NoDivisionTrajectoryAnalysisSinglesingle_DQ
Figure 2—figure supplement 1 (E-H)CPM_NoDivisionTrajectoryAnalysisSinglesingle_DM
Figure 3CPM_NoDivisionTrajectoryAnalysisSinglesingle_R
Figure 4CPM_NoDivisionTrajectoryAnalysisCircularPatternrotation_Q
Figure 4—figure supplement 1CPM_NoDivisionTrajectoryAnalysisCircularPatternrotation_N_R1
Figure 4—figure supplement 2CPM_NoDivisionTrajectoryAnalysisCircularPatternrotation_N_R2
Figure 4—figure supplement 3CPM_NoDivisionTrajectoryAnalysisCircularPatternrotation_N_R3
Figure 5 (A-D)CPM_Divisionwound_nodiv
Figure 5 (E-H)CPM_Divisionwound_div
Figure 5—figure supplement 1 (A-B)CPM_Division_SupplementFrontAnalysiswound_div_A
Figure 5—figure supplement 1 (C-D)CPM_Division_SupplementFrontAnalysiswound_div_D
Figure 5—figure supplement 1 (E, F)CPM_Division_SupplementFrontAnalysiswound_div_Q
Figure 6 (A-D)CPM_Divisionwound_div_fing_1.0
Figure 6 (E-H)CPM_Divisionwound_div_fing_1.1
Appendix 2—figure 1CPM_NoDivisionTrajectoryAnalysisSinglesingle_A

Data availability

We have uploaded the source code used in the main part of our study as well as the one used in the appendix. Furthermore, we have provided the full list of parameters in the figure captions, as well as exemplary parameter files for all applicable figures.

Additional files

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)