Bridging scales.

Actin dynamics is studied at many levels: from association of actin monomers to form filaments [33], to biophysical force production by F-actin [34], the assembly and branching of F-actin [35], and the resultant shape and motion of cells [36, 37].

An excellent review of actin dynamics models in single cell migration, [38] provides a solid bridge between molecular scale and cell scale phenomena. The authors explain the main model formalisms and show how models can be used to understand experimental data. A review of the link from actin and its properties to cell mechanosensing and behavior is given in [39]. Reviews of the literature on actin-based cell migration reveal a mature field, where quantitative methods and theory have ripened in tandem. This synergy has benefited greatly from biochemical pioneers, such as Tom Pollard (Yale University), who helped to nurture an appreciation for mathematical and physical modeling in the community.

Levels of detail.

The pair of papers [33, 40] aptly illustrates the dichotomy between highly detailed computational models [40] and more conceptual minimal models [33]. On one hand, the highly detailed [40] synthesizes a large body of experimental data for actin assembly and branching, including many actin-related component. At the opposite extreme, the models of [33], as well as later papers [36, 37] emphasize physical principles, universal properties, and general insights. In some sense, the overall predictions of both model types can be compared. The minimal models are more tractable for analysis, parameter sweeps, and overall insights, but are harder to connect to detailed molecular biology experiments. The complexity of highly detailed computational models (matching molecular experiments) makes it harder to navigate their results, which are more like a “1-to-1” map. Such maps may be important for those who “live in the neighborhood,” but are essentially baffling for newcomers. While distinct, one could argue that these approaches are complementary, so that overall insights can be obtained from one, and specific details from the other. The dichotomy between the detailed and the simplified models will reappear throughout this review, serving as 1 hallmark of distinct views.

Links with experiments.

A review of the experiments and theory for actin dynamics in the motility of keratocytes is provided in [41]. A retrospective that emphasizes the importance of theory, and case-study where theory has led the way is described in [42]. Notable contributions linking experiments to theory for cell shape include [36, 37], where shapes of cells are classified and related to actin branching and treadmilling.

Bridging scales.

Several papers link the dynamics of GTPases to cytoskeletal assembly and cell shape. Among these are experimental [50, 51] and computational [5, 52] works. A review of signaling that organizes the cell front and rear in Dictyostelium discoideum is given in [53], who also survey models of these as excitable systems. Some of these papers are described in fuller detail below.

Levels of detail.

GTPases act like molecular switches that cycle on and off the cell membrane and interact via crosstalk through effector proteins. A large collection of proteins participate in GTPase signaling: Guanine nucleotide exchange factors (GEFs) activate and GTPase-activating proteins (GAPs) inactivate the GTPases with varying degrees of specificity. GTPases are also sequestered in the cytosol by binding to guanine nucleotide dissociation inhibitors (GDIs). In some models, notably [54, 55], the details of the binding and mechanisms of activation–inactivation are lumped into phenomenological terms such as Michaelis–Menten or Hill functions. In [55], emphasis on cross-talk of Cdc42, Rac, and Rho and on spatial polarization comes at the expense of the molecular steps themselves (some of which remain unknown). In [56], more detail on such steps is modeled and eventually reduced to 3 ordinary differential equations (ODEs) using a quasi-steady-state approach. The authors focus on the role of a Cdc42-GEF in oscillatory phenotype in neuronal growth cone motility. In contrast, [57] concentrate on GDI binding of GTPases in a highly detailed computational model (originally crafted in BioNetGen and Virtual Cell).

Crosstalk between GTPases is modeled in great detail by [58], who constructed a Boolean model for epidermal growth factor (EGF) signaling to Rac and Rho. Their model has 38 intermediates and multiple reactions that activate or inhibit components. (See also [59], who combines PhysiCell with MaBoss to allow the modeling of intracellular signaling networks as Boolean networks). At lower level of detail are Rac-Rho mutually antagonistic models that leave out the forest of interacting nodes [60].

Simpler models for single GTPase include [61] for cell polarity and [62, 63] for cell shape. There, assumptions are made to condense underlying complex interactions into simpler “rules” of behavior. For example, “high GTPase activity leads to outwards protrusion force at the cell edge” (in the case of Rac) or “contraction” (in the case of Rho) [64, 65]. An advantage of this simplification is that simple models consisting of few partial differential equations (PDEs) allow for well-honed methods of applied mathematics, including nonlinear dynamics, stability theory, asymptotic analysis, and/or bifurcation theory to explore model predictions, parameter dependence, and regimes of behavior.

Serendipitously, simplified signaling models, like [61], occasionally also expose interesting mathematical structures to study. A case in point is “wave-pinning” [61, 66]: a wave of GTPase activity initiated at 1 end of a cell stalls before reaching the opposite end, resulting in robust “cell polarization.” Simplified models also permit numerical implementation in more complex geometries. For example, in [67], a wave-pinning model for a single GTPase in a single 3D cell is implemented by bulk diffusion methods. The cell does not deform, but the spatial localization of the GTPase is quantified in a fully 3D geometry. (See also [68], who considered neutrophil fluidization, and [69], who considered polarization in a 3D sphere).

At a gradually increasing scale are models that include not just GTPases (Rac, Rho, and Cdc42) but also other layers (phosphoinositides, actin, Arp2/3, and myosin) that interact to regulate cell polarity in response to stimuli [70], or dictate cell shape and directed motility [5, 52]. In the latter 2 papers, a single moving “keratocyte” is represented using a Cellular Potts Model (CPM), an energy minimization agent-based platform for computing cell shapes. It is shown that GTPase signaling can account for cell polarization, reorientation [52], and resolution of conflicting cues or obstacles [5]. (Compare with [71], who employ phase-field computational methods toward similar goals).

By way of comparison, in another more detailed approach, in [72], a cell is modeled as having a solid boundary that is moved using a Stefan condition. For the chemical signaling, Rho, Rac, 2 species of GEFs, F-actin, and G-actin are included. The model accounts for the basic repertoire of neutrophil motility.

All in all, models of GTPase signaling have taught us several valuable lessons, some with universal ramifications. First, modeling has provided clues to the functional significance of the seemly strange cycling of GTPases between membrane and cytosol: namely, the separation of diffusion rates resulting from these distinct compartments could be playing a role in pattern formation processes, an ingredient enabling GTPase localization or patterning in cells. Second, the stripped-down models have shown that chemical polarization in a cell need not depend on crosstalk between multiple types of GTPases—it can be set up by a single member of the family, given sufficient positive feedback and some depletion of its cellular pool [61]. Groups are still occasionally rediscovering on their own, the link between cell size and cell polarization that was implicit in [66], suggesting the need for more expository reviews of mathematical results. An important concept introduced in [5], but not yet fully recognized in the community, is the synergy between GTPase dynamics and its effect of cell shape and boundary curvature of the cell edge. Simply put, the “motion by curvature” of the chemical system interacts with boundary conditions to accelerate the dynamics. Altogether, the appreciation of GTPase signaling has benefited greatly from a host of distinct modeling and mathematical approaches.

Links with experiment.

A review of the links between models of cell migration and experimental data (image processing, cell tracking, and feature extraction from 1 cell to many) is given in [73]. Here we focus more specifically on experiments that highlight intracellular signaling.

The work by [70] provides data for the reorientation and polarity responses of HeLa cells from various starting cell states (polar, anti-polar, or nonpolar). The authors showed that an internal circuit of signaling (Cdc42, Rac, Rho, and phosphoinositides) could account qualitatively for the observed responses of these cells in microfluidic channels, with an externally controllable response to Rac. In both this and the follow-up [74], the 1D geometry of the channels helps to reduce the geometric complexity of cell shape, allowing for a better match to 1D spatial model representations.

How are models for GTPase crosstalk experimentally linked to cell morphology and motility dynamics? In experimental results of [51], the authors showed that signaling circuit of mutually antagonistic Rac and Rho could affect not only the dynamics of F-actin, but also shapes and migration of mesenchymal breast cancer cells. Together with a mathematical model for the Rac-Rho interactions, they were able to manipulate the Rac-Rho competition, demonstrate bistable states, and show that manipulating the system by inhibiting PAK (a kinase that mediates inhibition of Rho by Rac) displayed hysteresis characteristic of bistable systems. The significance of this paper is that it demonstrates a direct link between a simple hypothetical model for the way that Rac and Rho GTPases operate in a cell and the next level up, that of overall cell morphology.

Sometimes, individual papers provide only part of the story, but taken as a whole, a series of papers gives a broader view. The sequence of work in [50, 60, 75, 76] explore how Rac-Rho mutual antagonism is linked to cell morphology. In [75], the authors experimentally manipulated Rac and Rho activities to show spread or contracted cells. Interestingly, they found that combining constitutively active (CA) Rac and Rho simultaneously results in mixed morphologies in human glioma cells. The fact that high Rho and high Rac activity produces a mixture of possible coexisting stable steady state cell shapes was independently predicted in a purely modeling study by [60]. Related experiments by [50] on cell shapes were also later modeled and explained in a follow-up paper by [76]. These papers demonstrate that relatively simple stripped-down depictions of cellular signaling “modules” (such as Rac-Rho) can account for important and unexpected observations at the level of the cell as a whole.

The paper [77] explores how 3D collagen microtracks and confinement affect cell migration. The authors find that the degree of cell-extracellular matrix (ECM) interactions are key determinants of speed, morphology, and cell-generate substrate strains during motility.

In recent times, a clear link has been established between GTPase activity and mechanical tension experienced by cells. A pioneering experimental paper that showed this 2-way feedback is [78]. The effects of forces on the Rho family proteins, including the involvement of GEFs that activate Rac1, RhoA, or Cdc42 are reviewed in [79]. Some GEFs respond to cyclic stretch, and others to tensile force or shear stress and substrate stiffness. Rac1 and Cdc42 are activated by stretching of adhesion bonds [80]. Rho is mainly used in maintenance of focal adhesions, but it appears to play a prominent role in cell–matrix interactions.

Strain and strain gradients affect cell orientation. Single cell experiments with cyclic stretch are described in [81]. The authors suggest that the Rho-ROCK pathway that regulates myosin light-chain activity is responsible for sensing and responding to strain gradients.

Overall, it appears that the link between mechanical stimuli and GTPase signaling is still young, providing ample opportunities for creative modeling. So far, there is no consensus on what are the “takeaways” from the models so far. Further, it would appear that this gap should be filled if we are to successfully bridge between 1 cell and many, since the mechanics of cell collisions, as well as cell collective migration, entail sensing of both chemical and mechanical stimuli in the interacting cells.

Bridging scales.

Phase-field methods have been used to model cell shape and motility in 2 spatial dimensions (2D) [71]. The outline of a cell is represented by a level curve of some function ϕ(x, y, t). For example, in [88], a minimal model for Rac signaling was used to describe the shape and motility dynamics of a single phase-field cell. The model is then extended to pairwise and cell–train interactions in a narrow 2D strip. Interestingly, the authors had to assume an intracellular signaling “inhibitor” that is activated at cell–cell contacts to obtain appropriate behavior. This is one of the earliest models spanning the three levels of organization, from subcellular, to single and multiple cells. The model exploits a reasonably simple level of detail at each stage to reproduce several behaviors, such as reversals, cells walking-past one another, cell sticking, and forming a “cell–train.” We feel that this paper provides good prototype to be emulated by community members: a clear and well-studied intracellular system, linked to evolving cell shape, with feedback from cell–cell interactions. The paper leads to natural follow-up questions, amenable to experimental investigation: what cell components play the role of the putative “inhibitor”?

Models for neural crest cell (NCC) swarms in vivo and in vitro were proposed in [87, 89, 90], focusing on distinct behaviors of leader and follower cells. (See “Links with experiments” for a summary). For the same NCC system, the paper [7] developed a multiscale computational model linking signaling to cell group migration. A Rac-Rho GTPase circuit was represented by differential equations (ODEs) at nodes on the perimeter of a deforming polygonal cell. Co-attraction was represented as a non local tendency of cells to cluster, short range contact inhibition by up-regulating Rho at nodes in contact with other cells. Stochastic noise in Rac activity generated characteristic run lengths and reorientation seen in typical NCCs.

In [7], it was shown that these underlying biochemical systems could account for CIL. This study can be compared with the GTPase and inhibitor system described above, [91]. Coherent migration of a group of n cells, each of diameter d, was possible along a confinement corridor of width on the order of . This paper, [7], has a wider corridor and larger numbers of cells than [88], but confinement serves a similar function of enhancing directed cell movement. In real embryos, such “corridors” could be defined by permissive and inhibitory regions for NCC migration, as shown experimentally in developing embryos [92]. As a comparison to the models of [7, 88], in [93] a “rule-based” approach was applied to the same problem of NCC group migration, without bridging to the biochemical circuit.

In somewhat related work, [8], 2D cells are described by viscoelastic triangular finite elements. Cells can adhere to one another using linear springs. At each cell vertex, intracellular signalling controls cell protrusion and cell–matrix adhesion. The signaling networks consists of Rac and several other components, including the downstream kinase PAK, the focal adhesion protein Paxillin, the cell adhesion molecule cadherin, and Merlin, a protein implicated in contact inhibition of cells in a feedback loop with Rac1 and PAK. Key findings of this study were (1) cells elongate on stiffer substrates; and (2) stiffer substrates allow for better directional guidance and collective cell migration. Interestingly, the authors state that the principal role of intracellular signaling is to coordinate cell movement direction over distances of ≈ 500μm.

Levels of detail.

There are many “simple” SPP models in the literature on collective cell behavior, reviewed, for example, in [94, 95]. It is relatively straightforward to create rules that result in collective behaviors, but harder to do so when more details are included.

In 1 SPP cell model in a 1D spatial domain, [25], CIL and force-induced cell repolarization result in properties of collective cell behavior. Rules governing pairwise interactions and probabilities for CIL are derived from experiments by [96] described in the next section. It turns out that there is an optimal number of cells in a group for migration persistence. This paper joins [70, 74, 96] and others in proposing and using simplified 1D geometry to examine the basics of cell–cell interactions. The synergy between papers [25, 96] is very helpful, and points to the need for more paired experiment-model studies.

SPP models are easily extended to higher dimensions, and/or more realistic pairwise forces. For example, [29] considers a cell swarming model in 2D that is closely related to spherical cell models in [27, 97]. While [29] neglects cell–cell friction, the paper couples the cells’ compressibility, which plays the role of a local force of repulsion, to adhesion due to long protrusions such as filopodia that leads to nonlocality.

The highlight of [29] is the inclusion of these non local cell–cell attraction forces. In this way, the paper makes a seminal link between extended cell–cell interactions and a body of theoretical work on non local interactions in swarms or physical particles. The same paper also features the concept of H-stability (preservation of bounded density as the number of particles increases). This idea is used to ascertain correct cell spacing in swarming models, facilitating links to experimental observations. This paper hence imports a panoply of techniques, methods of analysis, and results, from theoretical work on non local models dating back to the 1990s [98], and inspires new ways of interpreting experiments.

Parallels with macroscopic swarms and flocks are useful, so long as we remember that these operate in vastly different friction regimes. (Cells experience an overdamped regime, where inertial forces are negligible). In yet another example of a cell swarming model, the paper [99] accounts for the development of stripes in zebrafish pigmentation patterns. The model comprises 2 cell types, with cell differentiation and death, and successfully reproduces a wide array of observed wild-type and mutant patterns.

In a different approach at a lower level of detail, [100] describes a group of border cells in Drosophila as a single sphere to determine the role of cluster size in chemotactic migration speed. The authors assume a spherical cluster moving by Stokes’ law in a viscous environment. The migration force is proportional to cell surface receptor occupancy in an exponential gradient of attractant. The results show good fits to both normal cells and cells that have a deficiency of receptors. The same biological problem is treated in [101], using a force-based approach for individual spherical cells (rather than cluster as a whole as in [100]) in 2D, followed by 3D (for small clusters of up to 8 cells). The computation incorporates cellular adhesion and repulsion, with some stochastic effects.

At the next level of detail, the shape of interacting cells is explicitly represented. The papers [91, 102] apply phase-field methods to model cell shapes in 2D with intracellular signaling. Interestingly, the signaling system is itself “minimal:” The authors assume basic “wave-pinning” GTPase dynamics, as in [61], with additional stochastic elements. The levels of internal GTPase give rise to either protruding or retracting forces on the cell edge. Then, [91] models the rotational cell motion when 2 cells interact (compare with [31] who showed a similar behavior in CPM cells). For coordinated movement, the directions of these gradients need to be linked. A key finding of [91] was that the emergence of rotational motion is strongly dependent on the type of cell–cell coordination mechanism (which includes the previous “inhibitor” field and other phenomenological cell polarization mechanisms).

We see gradual ascent in the numbers of cells considered and level of detail. We now find groups addressing the question “What can cells do well together that they do poorly on their own?” For example, collective chemotaxis was shown to be better in groups than in single cells [18, 103, 104]. (See also the related work by [22] for larger cell groups, reviewed in the next section).

Overall, our sampling of the literature reveals a paucity of models at the small-cell group level, relative to subcellular and tissue-scale modeling. We believe that papers such as [29] open opportunities for swarm-centered modelers to make a contribution to the cellular biology realm, by translating numerous results of physics and mathematics, suitably reinterpreted, to cells. We also contend that understanding the ways that cell interact in small groups is an important and relatively tractable step to formulating more informed tissue-scale models.

Links with experiments.

Pairwise interactions of MTLn3 rat breast adenocarcinoma cells inside microfluidic channels with chemotactic cues were quantified in [74]. These experiments provide important information for the spatio temporal evolution of GTPase activity during cell collisions. When cells collided in [74], they either stalled, adhered and moved together, or exhibited CIL. More recently, [105] studied how mechanics, and cell polarity drive collective cell migration in larger groups of cells. More experiments of this type would be very important in laying the groundwork for understanding pairwise cell interactions and in directly informing the development of models. Such models could then be used to interpret and deepen experimental observations.

In a similar vein is [96], a study of CIL experiments on micropatterned surfaces, where data were collected for the probability of cell repolarization after contact with another cell. The authors also explored persistence of cell trains, bridging from pairs to larger groups. An agent-based model was then used to predict cell–train outcomes based on the observed repolarization probabilities.

A model that replicates the experimental setup of [106] appears in [107]. Flat cylindrical model cells are assumed to align ECM fibers. Primarily a single-cell model, it is also extended to 2 interacting cells where, interestingly, leader-follower behavior is seen: The leader creates a “track” of aligned ECM fibers that the second cell can follow.

A number of works shed light on the molecular basis of cell–cell interactions. For instance, Wnt signaling is a ubiquitous pathway that regulates polarity, migration, and many other cellular processes. Wnt signaling at cell–cell interfaces that locally upregulates Rho is studied in [108]. Ephrins form cell–cell ligand–receptor bonds that funnel signals to GTPases to control cell repulsion or attraction. Involvement of Ephrin receptors in CIL [109], effect of tension on CIL [110, 111], and on the activation of RhoA [112] have been similarly studied.

Experiments on migration of primordial progenitor cells in the developing zebrafish [113] established the importance of guidance cues such as repulsion due to ephrins. Physical barriers also guide cell migration and organ placement. Actin was visualized for cells colliding with barriers expressing ephrin receptors. A particle-based model with an attraction–repulsion Lennard-Jones potential was then used to test the effects of reflecting and nonreflecting boundaries. (The Lennard-Jones potential is composed of power functions of the form r⁻¹² and r⁻⁶ for local repulsion and long-range attraction).

NCCs chemotax toward some cell types in a chase-and-run behavior [114]. For example, NCCs chase the placode cells, namely cells that are destined to form feathers or teeth or hair follicles. Before contact, each cell has high Rac at its front. Upon contact, the cell-adhesion protein N-cadherin, at the cell–cell interface, inhibits Rac in both cells leading to separation and escape of the placodes. Once the cells separate, chemotaxis is reestablished.

The group of Danuser [115] used high-resolution traction-force microscopy (TFM) to measure traction forces for a single cell and small groups of 2 to 6 epithelial cells. They then used a thin-plate finite element model (assuming homogeneous elasticity) to reconstruct forces along the interfaces between cells in a cluster. It was shown that the forces were correlated with E-cadherin localization. The authors reduce the cell group to a network representation, with vectors connecting cell centroids.

A recent paper explores the relationship between initial tissue size and the dynamics of tissue growth [116].

Experiments on the migration of chick neural crest cells in vivo and in vitro are described by several groups, including that of Kulesa. His experiments have been directly linked to models previously mentioned [87, 89, 90]. The group has investigated how cells reconcile mixed chemotactic signals, as in breast cancer cells in [74], and how information is shared between cells. The differences between leaders and followers is of recent interest in the group.

While papers such as [25, 117] argue eloquently for the need to carry out small-group cell experiments, we find this, too, to be a missed opportunity. Possibly, the in vitro work that this entails is viewed as less compelling than in vivo experiments, and so, less persuasive from the grantsmanship perspective. Hence, there are still relatively few experiments to map out the links between cell signaling and cell–cell collisions, between events inside the cell and the small-scale swarms that form. We recommend this area as a promising one for experiment-model synergy.

Bridging scales.

Here we concentrate on bridging between models for biochemical signaling and multicellular behavior. In the next section, we will focus on bridges between 1 or a few cells to larger populations and from simple to more detailed geometry (points, spheres, and cell shapes) and dimensionality (1D, 2D, and 3D).

The modeling paper by [13] is a vertex-based computation of epithelial dynamics, where force balance of cell vertices includes an active contractile force on the cell perimeter due to myosin. ODEs track Rho activation by cell perimeter stretching and downstream myosin activation (by phosphorylation) at a given vertex. Vertex motion is a result of force balance. Friction and passive forces are derived from energy penalties for area constraint and from the energy of adhesion to neighbor edges. This aspect of the paper resembles the energy-based approach in CPM computations, e.g., in [16]. The authors quantified the effects of cell density on correlation lengths, and regimes of behavior such as streaming, contractile waves, group movement, rotation in a circular domain, presence of vortices, and non uniformity of myosin distribution.

Some similarities are shared by [13, 16]. Both are tissue-based computations, with ODEs for GTPase biochemistry, an assumption that GTPase leads to cell contraction or expansion, and feedback from cell tension back to the activation of the GTPase. Both papers find fluctuations in cell shape or size from the full feedback between signaling and contractility (Fig 2C in [13], and the comparable Fig 3B in [16]). The signaling biochemistry is assigned to vertices in [13] versus the cell interior in [16]. In [16], a single cell is first linked to a 1D chain of cells, and then to 2D CPM cells, where waves of cell contraction are observed. The levels of detail in these 2 papers are comparable, though [13] also includes predictions for myosin, unlike [16]. See also [147] for a similar approach restricted to 1D cell chains.

In the same general class, the work of [10] links cell contraction (in a vertex-based hexagonal cell with 6 spokes) to biochemical details for a circuit of proteins known to regulate myosin contractility in Drosophila dorsal closure. This model is more detailed on the specific known interactions of protease-activated receptor (PAR) proteins that form a negative feedback loop with actomyosin (Baz, Par-6, aPKC), including a set of 9 ODEs for components linked to myosin dynamics along cell edges and spokes. The authors account for several phases in the developmental process, including oscillations observed in the tissue (using time delays) and the eventual contraction of the tissue. By way of comparison, [148] has a greater level of detail of the cell shape, including more viscoelastic spokes and edges, but essentially no biochemical detail. An example of 3D vertex-based simulations with cell rearrangement and out-of-plane bending can be found in [149].

Levels of detail.

The gene regulatory network involved in Drosophila ventral furrow formation is modeled at a fine level of molecular detail in [150] using a Boolean modeling approach. The model encapsulates interactions of transcription factors to explain 3 phenotypic attractor states, and identify missing and alternative pathways. Spatial distribution of cells and cell shape dynamics is not considered.

At the opposite extreme, internal cell signaling is omitted in favor of positions and spatiotemporal dynamics of cells. SPP “agent-based” models with forces of attraction, adhesion, and repulsion include [29], as reviewed in a previous section. There, nonlocal forces are considered in groups of hundreds of cells.

Modified agent-based models in [151] can be compared to the cell-pair study in [91]. Additional assumptions are included in [151], notably Vicsek-type alignment, cell membrane curvature, elasticity, actomyosin cable force, and a tendency to move in the outwards normal direction. The authors show the dependence of tissue wound-healing velocity on the width of a confining stripe. On a similar topic, CPM computations are used in [152] to investigate the effect of adhesive micro patterns on the motion of single cells and collective cell migration. Similarly, modified Viscek-type models are considered by [25] endowing each cell with a polarization direction. Various types of cell–cell polarity coordination mechanisms are considered.

Yet another modification of particle-based models is [153], where a cell is represented by a spring linking 2 discs, to depict the “cell body” and the “pseudopod.” The model demonstrates various outcomes of binary collisions in a 2D domain, as well as order–disorder transitions and velocity waves in large 2D cell groups with and without CIL. This paper explicitly bridges from 1 to many cells.

At the next level of complexity, we find representations of cells as spheroids or ellipsoids (Fig 5A1, [27, 154–156]). Drasdo computed a 2D growing monolayer [154], and later extended it to 3D [27]. Palsson built 3D simulations of deformable off-lattice ellipsoids, with cell division, chemotaxis, cell–cell adhesion, and volume exclusion. The platform has been used to model the dynamics of Dictyostelium aggregation and cyclic AMP (cAMP) signaling [155], as well as the posterior lateral line primordium, a cluster of about 100 migrating cells in zebrafish embryos that generate sensory structures on the surface of the fish [14].

Other recent papers add variations on these themes. For example, [157] provides a computational model for 3D off-lattice spheroid cells, but with greater level of detail for cell–cell and cell–matrix adhesion. The model is coarse grained from the biological levels of 10⁶ down to 100 sites per cell. After quantifying the stable distance between cell pairs (for a given set of forces, intrinsic cell properties, and adhesion sites), the authors go on to show the behavior of a larger group of cells.

Recently, these approaches have been extended to more finely resolve cell shapes in 3D [9]. Here, each cell is a triangulated 3D object, composed of discrete viscoelastic elements. Simulations track up to a thousand liver cells. While computationally very expensive, these approaches can provide valuable insights. For example, it is shown that highly deformable cells move more easily to heal a lesion than stiffer or more rigid cells [9]. They are also useful for specific systems where it is crucial to capture details accurately, as, for example, in drug design in silico for liver disease.

In [158, 159], a confluent epithelium in confinement is modeled by vertex-based polygonal cells, and the type of collective migration (directed flow, vortex chain, or turbulent) is related to a dimensionless quantity (“cell motility number”) that combines cell motility, cell density, and size of the confinement pattern.

In the phase-field category, the paper by [22] scales up from small numbers of cells in [91] to larger populations. Cell collisions are modeled by an energetic penalty for overlapping phase fields, and cell adhesion is described by a reaction–diffusion equation in [22]. Their model cells exhibit bistable shapes, (either symmetric or keratocyte-like), and both elastic and adherent collisions are predicted. Strong adhesion results in bands of dense closely packed cells that move as a collective. The authors compare their work with that of [160] (see Links with experiments) where a CPM model is used to describe large confluent cell rotations in a confined region. As noted by [22], the level of detail in their phase-field description is better suited for non-confluent cell models, but likely too detailed for confluent layers of cells where individual cell details average out. In [161], similar methods are extended to a 3D phase field computation for a single cell interacting with curved or grooved surface. This sequence illustrates the trade-off between geometric complexity (1D versus 2D versus 3D) and collective complexity (1 cell versus few versus many).

Overall, it emerges from many papers that the vertex-based models are good at describing epithelia, where fragmentation or EMT is absent or not important. Vertex-based models are not well suited to track cell death or fragmenting clusters, since edges and nodes are shared by neighboring cells. Breakage of a cell requires that shared edges and nodes be duplicated or reassigned, an inconvenience. To track tissue fragmentation or loss of cells, center-based models or CPM platforms have an advantage, since these represent each cell by an individual geometric object or by a set of pixels.

Continuum models are also commonly used at the level of tissue dynamics. Here, the cell identity is omitted altogether, in favor of cell densities, local flows, and material properties. Analogies are made with fluids, viscoelastic material, elastic sheets, or foams. Advantages include the ability to harness traditional methods of physics, mathematics, and fluid or engineering computational software. Examples of this approach are numerous, and we mention only a few.

One example of this class, [162], is a 3D physical description of an epithelial sheet. The authors show bistability in cell aspect ratio, bending and buckling instabilities of epithelia, and transitions that lead to formation of epithelial tubes and spheres. Experimentally testable scaling laws are given for such morphogenetic transitions. While there is no subcellular molecular detail, the connection between local cell shape and epithelial behavior is an important contribution.

In [163], the tissue is modeled as an active continuous medium, with a Maxwell viscoelastic constitutive law. A reaction–diffusion-transport equation accounts for actomyosin, whose contractility is coupled to cell motion and polarity. The authors explore predicted motion of tissue in confined regions, closure of wound gaps, and the relationship between traction forces and sizes of cohesive cell clusters. Mechanical waves, as well as tissue rigidity cycles, are observed, where fast fluidization is followed by a gradual period of stiffening.

The continuum model in [164] treats a tissue as a compressible fluid and includes both cell division and death in a type of Stefan free boundary problem. The authors show that behavior of the tissue in wound closure and in colony expansion depends on 3 parameters: 2 physical constants and the proliferation rate. They argue that all these can be estimated from limited experimental data, with examples calculated for IEC-6, a rat cancer cell line and MDCK cells, a canine kidney cell line. As the authors point out, continuum models are appropriate, provided the size of the tissue or the characteristic length of the wound is sufficiently large relative to the size of individual cells.

In some cases, an approach that combines features of both continuum and discrete cell identity is implemented. One example is [165]. Here, a 1D cell monolayer consists of a long contractile element with myosin creating a strain rate resulting in length changes. This element is flanked by cells in front and rear. Binding and unbinding of adhesion sites are also included. The model tissue exhibits durotaxis, movement directed up a stiffness gradient, under appropriate conditions on the myosin and adhesion parameters. The authors use this model to conclude that a monolayer is more effective at durotaxis than single cells.

Other examples of hybrid treatments of particle-based and continuum approaches include [166, 167]. These papers tackle the important question of how to derive appropriate continuum models from underlying SPP models. This kind of work and, in particular, expository papers that summarize the conclusions in ways that biological modelers can understand, could help link work in the literature that is currently underappreciated or not understood.

A recent work, [168], describes fingering at the front of an epithelial sheet using several of the above approaches. Cells are represented by pairs of points, as in [153], and also by Voronoi polygons for computations corresponding to experiments in [131]. The authors develop an active fluid model for the epithelium. Using these combined approaches, they demonstrate that stable fingering of the tissue edge requires leader cells. They characterize the wavenumber of the fingering instabilities (distance between stable fingers) using stability analysis of the continuum PDE model.

In summary, while the literature on tissue-scale collective migration is rapidly growing, it is still in stages of infancy as far as coherence, coordination, and clear direction are concerned. We see many exploratory steps in many individualistic directions. The list of core questions and the discovery of unifying principles is beyond the horizon, providing ample challenges for the community. We can draw an analogy between the current state of the art and the behavior of a growing but uncoordinated population of cells. There is no emerging unified front. Each is exploring independently and occasionally aggregating with a few others, but the global organizing principles are yet to emerge.

Links with experiments.

The links between molecular players such as Rho family GTPases and epithelial morphogenesis have been known for some time. (Experimental literature reviewed in [169]). The role of Rho GTPases in the leader-follower identities and in front–rear tissue polarity is reviewed in [170].

Over the years, investigators have been asking how external constraints, geometry, strain fields, gradients, and other factors affect collective behavior. How does collective cell migration emerge from transfer of mechanical information between cells? For a cell to be a “leader,” should it be more sensitive to stimulation than other cells? Should it have elevated or more responsive GTPase activity, for example?

These and similar questions are posed in the experimental work of [134]. Here, the authors investigate the roles of the GTPases RhoA and RhoC and their GEFs in collective cell migration. It is found that cells tend to crawl in the direction of maximal principal stress, a process called plithotaxis. In their scratch-wound assay of MDCK and human bronchial epithelial cells, wounding leads to coordinated waves of motion, with cells at the front edge reacting first, followed by those successively further back. The authors speculate that mechanical cues are induced by leader cells on followers behind then by normal strain, and alongside them by shear strain to coordinate motion.

In [134], the relationship of tissue speed to distance from the leading edge is quantified using particle image velocimetry. The authors employ shRNA to knock down GTPase Cdc42 or Rac1 and many of their GEFs (screening some 81 GEFs in total). They find that the speed and directionality of the cells drop everywhere. When RhoA is depleted, there is a reduction in the spatial gradient of cell speed. They suggest that the molecular mechanism may include signaling downstream of cadherin, as well as Merlin-Rac1 signaling.

The experiments in [129] elucidate the effect of pulling stress on Rac and Rho GTPases. The authors investigated which properties of underlying molecular machinery allow for coupling between mechanical forces and correlated cell motion. The correlation length scale of collective force transmission was determined experimentally in [171]. Observed behavior was then modeled using a thin elastic sheet of height h with some elasticity and isotropic contraction stress. The authors investigate the emergence of leader cells and found a typical length scale of about 170 μm.

A number of experimental studies have explored how epithelia behave in grooves [130, 132], various topographic surfaces, or confined settings [105], including arrays of posts [172], arrays of grooves [173], or adhesive [160] surfaces. In [132], the width of grooves and adhesive strips is varied to investigate how these affect collective migration of an epithelial sheet of MDCK cells. The authors map out the force and velocity fields in each case. Their narrowest grooves are 20 μm wide, so cells are in a single file, and the authors observe “push-pull” or “caterpillar” type motion, resembling the relaxation–contraction cycles in the model by [16]. Larger tissue widths exhibit vortices of cell motion. The authors deduce that the constraints of the geometry influence cell rearrangement, as well as junctional forces between cells. (Compare with [133], who observed sustained oscillations in epithelial sheets).

Micropillar arrays form the playing field in [172] to observe collective migration and EMT in mammary epithelia and breast cancer cell lines. Dispersal of single, highly motile mesenchymal cells from a spreading epithelial front was shown to agree quantitatively with a minimal physical model of binary mixture solidification. While such a model lacks cell detail, it has several advantages, including simplicity, basic summary statistics for comparison with experiments, and availability of an analytic solution. Methods from physics can be used to inform the link between the material properties and the overall macroscopic behavior.

Both single fibroblasts and epithelial monolayers were studied in [130], showing that cells tend to be more organized on grooved versus flat substrata. Mechanical exclusion interactions, rather than strength of junctions between cells, affect the distance over which the topographic guidance signal propagates between cells. The same paper proposes a CPM computational model to describe the observations. In the CPM model, following the the style of [15], each cell is assigned a phenomenological “polarity” vector that both guides and is affected by cell displacement. Aside from the customary Hamiltonian area constraint and cell–cell adhesion energy, there is also a term for a phenomenological motile force, implemented as a “migration energy” (dot product of the polarity vector and the cell centroid position). The authors investigate how this motile force affects spatial correlation of the velocity. They observe emergence of streaming patterns. The authors also describe the influence of “leader cells,” whose polarity vector is set to an external cue. Leaders are embedded in the tissue interior, and they coordinating neighbors over some “interaction distance.”

In a similar style of experiment-model study, MDCK cells are seeded on a circular adhesive domain in [160]. Once a critical density is attained, the confluent culture collectively rotates. The same behavior is then captured in a CPM computational model that includes a motile force and a polarity term that is reinforced by cell displacement, with some persistence time. The authors show that rotation takes place when the size of the circular domain is on the order of the correlation length of cells. That correlation length, in turn, depends on the cell—cell adhesion energy, the motile force magnitude, and the polarity persistence time. Note the comparison with the purely computational paper by [22] using phase-field methods, where cell collisions, rather than adhesive confluent cells, are the subject of focus. It is also interesting to compare the CPM model in [160] with the vertex-based treatment of a similar situation in [158].

A recent paper, [105], demonstrates the fact that single cell properties affect much of the collective behavior of a cell population. In this ground breaking work, the authors link single cell Rac1 polarity to the emergent rotations of confluent cells (from few to many) confined to a ring or closed curve. A simple mechanical model captures the balance of directed motility and contact forces to demonstrate the principles at work. This paper spans subcellular to multicellular scales and provides a great example of elegant experiments to be emulated by others.