Abstract

We consider a generic model for cell motility. Even if a comprehensive understanding of cell motility remains elusive, progress has been achieved in its modelling using a whole-cell physical model. The model takes into account the main mechanisms of cell motility, actin polymerization, actin–myosin dynamics and substrate mediated adhesion (if applicable), and combines them with steric cell–cell and hydrodynamic interactions. The model predicts the onset of collective cell migration, which emerges spontaneously as a result of inelastic collisions of neighbouring cells. Each cell here modelled as an active polar gel is accomplished with two vortices if it moves. Upon collision of two cells, the two vortices which come close to each other annihilate. This leads to a rotation of the cells and together with the deformation and the reorientation of the actin filaments in each cell induces alignment of these cells and leads to persistent translational collective migration. The effect for low Reynolds numbers is as strong as in the non-hydrodynamic model, but it decreases with increasing Reynolds number.

1. Introduction

Substrate-based cell motility is a well-studied process for eukaryotic cells, such as keratocytes, fibroblasts and neutrophils. It plays a fundamental role in tissue growth, wound healing and immune response. Less explored are motility mechanisms where local adhesion is less evident, such as for cells moving in martigels or freely swimming microorganisms. We here consider a generic model, which accounts for a wide range for motility mechanisms by considering: (i) the generation of a propulsive force by actin polymerization, which acts against the cell's membrane, (ii) the formation of adhesive contact to the substrate, transferring this force to the substrate (if present), and (iii) a contractile action of actin–myosin complexes determining the cell polarity and being responsible for retraction of the cell's rear (e.g. [1,2] for a review on the forces involved in cell movement). Several experimental studies for fish keratocyte (e.g. [3–5]) indicate a self-organization process behind the motility mechanism, which has been adapted in various theoretical approaches [6,7]. Similar models (without the adhesive contact) have been proposed as generic models for cell motility in three-dimensional environments [8–10]. They all apply an active polar gel theory [11–13]. If considered in a confinement, a splayed polarization of the actin filaments can occur, which models the contractile stress due to the interaction of myosin and actin. If combined with the treadmilling process of polymerization and depolymerization of actin filaments (e.g. the one considered in [14–16]) and if applicable an effective treatment of the adhesive contact, a whole-cell physical model for moving cells can be constructed [10,17]. Such models have been established for single cells and are used to analyse motility of various cell types [10,18]. The results strongly support the physical view on cellular motility, which exploits autonomous physical mechanisms whose operation does not need continuous regulatory effort. Recently, such models have also been considered for collective migration [19]. Here, each cell is considered as an active polar gel and interactions between the cells are specified. The model predicts that collective migration emerges spontaneously as a result of inelastic collisions between neighbouring cells. These collisions lead to mutual alignment of the cell velocities and to the formation of coherently moving multi-cellular clusters. These results essentially confirm simpler agent-based modelling approaches of Vicsek type [20] with inelastic behaviour in the interaction rules [21], recent mesoscopic simulations based on active phase field crystal models [21] and continuum approaches which only consider the emerging macroscopic behaviour [22,23] using Cahn–Hilliard type models. All these approaches for collective migration neglect hydrodynamic interactions. The effect of these interactions on collective migration is controversially discussed. In the related problem of motility-induced phase separation [24], where clustering results from a reduction of the propulsion speed due to cell collisions in environments with high local density [25,26], a suppression of cluster formation is observed if hydrodynamic interaction is taken into account, while the hydrodynamic active Cahn–Hilliard model in [27] leads to arrested phase separation.

We here consider the hydrodynamic active polar gel model, which was used in [10] for a single cell, for multiple cells. Each cell is thereby described by a phase-field variable, which defines the confinement of the field variables of the active polar gel model for each cell. The interaction between the cells only considers steric interactions. Short-range repulsion is realized by a Gaussian potential using the phase-field variables [28]. Using a multi-mesh approach [29], which allows for an efficient numerical treatment by considering differently refined meshes for each variable, allows one to significantly reduce the computational cost and to consider numbers of cells, which are sufficient to observe collective migration.

The paper is organized as follows. In §2, we introduce the mathematical model and compare it with the non-hydrodynamic model in [19]. We further discuss numerical aspects. In §3, we first perform several computations for binary collisions before the onset of collective migration is studied for larger systems. The simulations do not indicate a suppression of collective motion if hydrodynamic interactions are considered. However, the number of cells which can be considered in the numerical study is limited and it remains open, if long-range polar order can also be established for much larger systems.

2. Mathematical model for cell motility

The mathematical model is based on physical phenomena and results from energy minimization, conservation laws and active components, taking into account the filament network, the cell membrane, cell–cell and cell–substrate interactions, as well as fluid properties.

2.1. Energy

Following [8,10], we consider the free energy of a single cell i2.1which consists of the energy of the filament network , described by an orientation field , which is the mesoscopic average of the actin filaments and the surface energy of the cell membrane . Each cell is described by a phase-field variable , defined as , where characterizes the thickness of the diffuse interface and denotes the signed-distance function between , in the considered case a bounded domain in , and its nearest point on . Depending on , we label cell i with and the outside with . The cell membrane is then implicitly defined by the zero-level set of . In [8], the cell has been considered as a droplet for which the surface energy reads
2.2where denotes the double-well potential and is the membrane tension. In [10], also a bending energy of the cell membrane was taken into account using the Helfrich energy in a phase-field approximation [30,31]
2.3with denoting the bending rigidity and the derivative of the double-well potential with the spontaneous curvature . The surface energy thus results as a combination of both energies
2.4

In the following, we will consider , and for simplicity. The energy of the filament network of cell i is given by
2.5

The gradient term with the positive Frank constant is a simplification of a general distortion energy formulation from the theory of liquid crystals, with the assumption of the same value of the stiffness associated with splay and bend deformations (e.g. [32]). Linking to the second term allows restricting to the cytoplasm: if the minimum is obtained for and thus the term does not contribute to the energy, and for , the term forms a double-well with two minima with and the form specified by the parameter . The last term in equation (2.5) guarantees for that points outwards in normal direction to the cell boundary. This is required to account for the effect of polymerization of actin filaments [33]. We will again only consider the case , and .

The overall energy for N cells and their interaction in a fluid environment is given by
with the kinetic energy and the velocity . For the sake of simplicity, we consider in the derivation equal density and viscosity for and the fluid outside , which is considered as an isotropic Newtonian fluid, so that
2.6with and . We further introduce the phase field containing all cells. Figure 1 provides a schematic description for two cells.

Schematic description for two moving cells. Shown are the splayed orientation field , as well as the streamlines of the velocity profile and the phase fields with the cell membranes corresponding to the zero-level sets of . (Online version in colour.)

The cell–cell interaction energy requires a coupling of all surrounding phase fields with . We here consider only steric interactions and model a short-range repulsion by a Gaussian potential. Following [28], we use the definition of to compute the signed-distance function , which is used to link cell i and cell j. Within the diffuse interface region, we obtain
2.7and thus we can write the Gaussian interaction potential within the phase-field description as
2.8with being non-zero only within the diffuse interface around , the interaction function
2.9and the strength of the repulsive interaction between cell i and cell j with respect to the evolution of cell i. Here, we consider a constant repulsive interaction strength, hence . The approach circumvents any non-local terms which are typically required for cell–cell interactions and has been analysed in detail in [28].

2.2. Non-dimensional form

Before we introduce the governing equations, we consider the energies in a non-dimensional form. We consider the characteristic values for space , velocity and energy , with characteristic length L, characteristic velocity V and fluid viscosity . This yields a time scale and a pressure . We further define the constants and , and the dimensionless quantities:
which are Reynolds, capillary, bending capillary, polarity and interaction number, respectively. Dropping the notation, we obtain the energies in a non-dimensional form
and again , and .

2.3. Governing equations

The hydrodynamic model is an extension of the model in [8,10]. The governing equations look similar, but now have to be considered for each cell with the additional contribution from the interaction terms. We denote the variational derivative or chemical potential of the orientation fields and phase fields by and .

The evolution equations for the phase-field variables are regularized advection equations with the advected velocity given by the fluid velocity . The introduced diffusion term is scaled with a small mobility coefficient . The equations read
2.10and are coupled with each other through the fluid velocity and the interaction terms, which are contained in the chemical potentials , which read
for .

The orientation field equations for each are the same as for the single-cell case and read
2.11for i = 1, … , N, where the left-hand side is the co-moving and co-rotational derivative where the vorticity tensor defined as takes rotational effects from the flow field into account. The first term on the right-hand side describes the alignment of with the flow field, with the deformation tensor . and are non-dimensional material parameters. The evolution equations are defined in , but due to the coupling with we have outside of cell i. The non-dimensional chemical potentials read

The flow field and pressure p are defined through the incompressible Navier–Stokes equations, which read
2.12and
2.13with friction coefficient , modelling substrate adhesion, hydrodynamic stress tensor , consisting of passive and active components, and a forcing term . The viscous stress is
2.14with and if the outer fluid and the cells have the same viscosity and a quotient if they differ. The active stress due to actin–myosin complexes is
2.15with the active force number and . The stress coming from the distortions of the filaments reads
2.16and for the Ericksen stress we consider the divergence to be defined through
2.17

The forcing term accounts for actin polymerization and reads , with the non-dimensional self-propulsion velocity . We again only consider the case .

If we set , we obtain the system considered in [10] with two additional terms in the Navier–Stokes equations. The first is the friction term , which has not been considered as the focus in [10] is on motility in environments without local adhesion, and the second is the forcing term , as actin polymerization is not taken into account in [10]. However, both terms had already been considered in [8]. Considering friction as a modelling approach for substrate adhesion requires strong approximations. More detailed approaches can be found in [18].

2.4. Non-hydrodynamic model

For comparison, we consider also a non-hydrodynamic model. As all stress and forcing terms have been considered in the Navier–Stokes equations, we cannot simply neglect the hydrodynamic interactions. Instead, we consider
2.18and
2.19with the advections only due to the self-propelled velocity . The chemical potentials and are defined as before. This model can be related to the model used for collective migration in [19]. However, several differences should be pointed out. We here neglect the treatment of adhesion bonds and the viscoelastic properties of the substrate. Furthermore, the cell–cell interaction is considered differently. We do only consider steric interactions and no cell–cell adhesion. However, the strongest difference is the treatment of the orientation fields . In [19], only one variable is used for all cells. As the equation contains diffusion/elasticity of the orientation field this induces an unphysical coupling of the actin filaments over cell boundaries.

2.5. Numerical approach and implementation

The system of partial differential equations is discretized using the parallel adaptive finite-element toolbox AMDiS [34,35]. We use a semi-implicit time discretization and an operator splitting approach that allows us to decouple all subproblems, similar to [10,28]. We further conduct a shared memory OPENMP parallelization to solve the phase-field equations and the orientation field equations via a parallel splitting method. Each linear system of equations is solved using the direct solver UMFPACK. As the computational mesh has to be fine along the interface, adaptive mesh refinement is heavily used. However, using a single mesh for all variables is not appropriate in this case as, for example, the phase-field variable only requires a fine resolution close to the zero-level set of but not at the zero-level sets of with . The efficiency would go down if the number of cells increases if a single mesh would be used. The multi-mesh strategy, considered in [36] for two meshes, overcomes these numerical problems and assigns a mesh to each phase-field variable, which can be independently refined. In [28,29], the approach is extended to arbitrary meshes and validated for related problems.

3. Simulations and results

3.1. Binary collisions of cells

We first study binary collisions of cells within a symmetric set-up with a fixed incidence angle of 45°. Figure 2 shows snapshots of the cell shapes and orientation fields together with the flow field if appropriate. The cells deform at collision, the deformation influences the orientation fields which set the new directions for cell motion. For the hydrodynamic model, each cell is accomplished with two vortices. Upon collision the two vortices, which come close to each other, annihilate. This leads to a rotation of the cells and together with the deformation and reorientation of the orientation fields set the new directions for cell motion. In both cases, the non-hydrodynamic and the hydrodynamic cases, the coupling between the involved fields leads to partly inelastic collisions and alignment. However, the strength of the alignment strongly depends on various parameters. Figure 3 shows the centre of mass trajectories for the non-hydrodynamic model and for the hydrodynamic model for different Reynolds numbers Re. The results show a tendency from more inelastic towards more elastic collisions for increasing Re.

(a) Non-hydrodynamic model. Shown are the cell shapes and the orientation fields. The parameters used are Ca = 0.0281, Be = 0, Pa = 0.1, In = 0.1125, , , , , and . (b) Hydrodynamic model. Shown are the cell shapes and the orientation fields, together with the flow field. The parameters used are Ca = 0.025, Be = 0, Pa = 0.1, In = 0.1, Fa = 1, Re = 0.001, , , , , , , and . The results for smaller values of the friction coefficient θ are similar (results not shown). The time instances for both cases are t = 3, 17, 30 and 45.

Centre of mass trajectories for binary collision for the cases considered in figure 2 and Re = 1. (Online version in colour.)

All simulations are performed within a two-dimensional computational domain of size . Each cell has a size, corresponding to a circle with radius . We apply periodic boundary conditions in each direction. A systematic study of the influence of various parameters on alignment (not shown) reveals mainly the same qualitative dependencies for the hydrodynamic and the non-hydrodynamic model, even if the mechanism behind alignment significantly differs. The alignment is more efficient at small incidence angles and it is stronger for higher capillary numbers (Ca) and smaller polarity numbers (Pa). Only the strength of the self-propulsion seems to have the opposite effect. While a larger value for leads to more elastic collisions in the non-hydrodynamic model, it leads to more inelastic behaviour in the hydrodynamic model. However, the effect is small if compared with the influence of the other parameters. The influence of the bending capillary number (Be) is negligible. All other parameters are kept fixed. Clearly, the binary interaction behaviour is beyond simple particle-based models, even if elastic deformations and/or hydrodynamic interactions are considered. The strength of alignment in the considered models is a result of the complex interplay between the cell shapes, viscosity, passive and active stresses, as well as actin polarizations and adhesion. The results further indicate the effect of the hydrodynamic interactions, with a tendency towards more elastic collisions for increasing Reynolds number Re.

3.2. Collective motion

We now investigate collective motion. For low cell densities, collective motion is dominated by binary collisions. So from the previous results we might guess the onset of collective motion also within the hydrodynamic model, at least for low Reynolds numbers Re. To quantify the effect, we introduce an order parameter
with the velocity vector of the ith cell. The parameter is 1 if all cells move in the same direction and 0 if no correlation of the directions exists. Figure 4 shows snapshots of the evolution for 23 identical cells, which initially move in random directions. The cell size now corresponds to a circle with radius . The domain sizes as well as all other parameters are as in the previous section with Reynolds number Re = 0.001.

Snapshots of the cell shapes, orientation fields and fluid velocity, if appropriate. (a) Non-hydrodynamic model and (b) hydrodynamic model. The snapshots correspond to the same times, shown in non-dimensional units. The parameters are the same as in figure 2. See also electronic supplementary material movies S1 and S2.

The result is quantified in figure 5, which shows the evolution of for the non-hydrodynamic model and the hydrodynamic model for two different Reynolds numbers Re. These results for the non-hydrodynamic model confirm the findings in [19]: without hydrodynamic interactions, collision of deformable cells can lead to collective migration if the collisions are inelastic. This is even true if for each cell a separate orientation field is used and thus any diffusion/elastic interaction between these fields is impossible. The situation with hydrodynamics has not been analysed before. The results indicate that also for low Reynolds numbers Re = 0.001, which essentially corresponds to the Stokes regime and is the most relevant situation for cell motility, collective migration can be observed. The time to reach collective motion is longer, but all simulations within this regime lead to persistent translational collective migration. Even if the mechanism is different, the analogy between inelastic binary collisions and collective migration seems to hold also for the hydrodynamic model with low Re. For Re = 1, the situation changes. The binary collision was more elastic and thus does not suggest collective migration. However, the more elastic collisions cannot suppress collective migration only the time to reach this state is significantly increased.

The diagram shows the temporal evolution of ω for the non-hydrodynamic model and the hydrodynamic model for two different Reynolds numbers Re.

Increasing the viscosity of the cells relative to the viscosity of the surrounding fluid (results not shown) has qualitatively no influence on these results. In both cases, Re = 0.001 and Re = 1 and collective migrations are reached faster as for and the fluctuations in before reaching collective motion are reduced.

These simulations indicate collective migration for deformable cells even under the influence of hydrodynamic interactions. In the low Reynolds number regime, all performed simulations result in collective migrations. The effect seems to be as stable as without hydrodynamic interactions. Only for Re = 1 is the time to reach collective migration significantly increased and even larger Re might be able to suppress the formation of collective motion.

In all these simulations, the system size and number of cells is relatively small. For larger systems the resulting collective behaviour might be more complex. We expect the formation of clusters, as for example, observed experimentally for myxobacteria [37] and found computationally for self-propelled rods [38,39] and within a microscopic field theoretical approach [21].

4. Conclusion

We have developed a computational model for the collective migration of cells. On a single cell level, the model is based on the well-established mechanisms of cell motility accounting for actin polymerization, motor-induced contractility and substrate adhesion (if applicable). The model uses the hydrodynamic active polar gel theory [11–13] and is comparable with the approaches in [6–8,10]. Each cell is treated individually using one phase-field variable per cell. Cell–cell interaction is considered through an additional potential with a short-range repulsive force as used and validated in [28,29]. The overall model only uses physical mechanisms, which do not need continuous regulatory effort. It describes details of the motility mechanism which allows one to study the influence of many parameters on the dynamic behaviour. The related non-hydrodynamic model [19] could already reproduce many experimentally observed phenomena. The overall question to answer is if these phenomena persist under the influence of hydrodynamic interactions, which is controversially discussed [25–27]. On the level of detail, which is considered in this paper, the effect of hydrodynamic interactions has not been studied before. Our results on the collision of two cells lead qualitatively to the same results as in the non-hydrodynamic model [19]. These binary cell interactions may be quantified in terms of inelastic or elastic collisions. In the hydrodynamic model, the variation of various parameters shows the same tendency to one or the other as in the non-hydrodynamic case. However, with a stronger deformation of the cells and a more elastic behaviour if the Reynolds number Re increases. As inelastic collisions have been reported as one indicator for collective migration [19], these results suggest the onset of collective migration also if hydrodynamic interactions are taken into account, at least for low Re. The simulations with 23 cells confirm this. All considered cases lead to persistent translational collective migration. Only the time to reach it differs and increases significantly with increasing Re. The considered parameters are Re = 0.001 and Re = 1. Even larger Re, which might be able to suppress collective migration, are irrelevant for typical situations of cell motility. These results provide valuable insights into the physics behind the biological processes in collective cell migration. It answers fundamental questions on collective motion for self-propelled particles and suggests some experimentally testable predictions. Can collective migration be found without cell–cell adhesion; is the effect stronger for cells with smaller membrane tension and larger elastic properties, as all predicted by our simulations; and can the effect of viscosity on collective migration be observed?

Competing interests

We declare we have no competing interests.

Funding

W.M. and A.V. acknowledge support from the German Science Foundation through Vo899/11. We further acknowledge computing resources at JSC through grant no. HDR06. We also thank the Isaac Newton Institute for Mathematical Sciences for its hospitality during the programme ‘Coupling geometric partial differential equations with physics for cell morphology, motility and pattern formation’ supported by EPSRC grant no. EP/K032208/1.