Multiscale modelling of motility wave propagation in cell migration

The collective motion of cell monolayers within a tissue is a fundamental biological process that occurs during tissue formation, wound healing, cancerous invasion, and viral infection. Experiments have shown that at the onset of migration, the motility is self-generated as a polarisation wave starting from the leading edge of the monolayer and progressively propagates into the bulk. However, it is unclear how the propagation of this motility wave is influenced by cellular properties. Here, we investigate this question using a computational model based on the Potts model coupled to the dynamics of intracellular polarisation. The model captures the propagation of the polarisation wave and suggests that the cells cortex can regulate the migration modes: strongly contractile cells may depolarise the monolayer, whereas less contractile cells can form swirling movement. Cortical contractility is further found to limit the cells motility, which (i) decelerates the wave speed and the leading edge progression, and (ii) destabilises the leading edge. Together, our model describes how different mechanical properties of cells can contribute to the regulation of collective cell migration.

Collective cell migration is a fundamental process both during embryonic development and pathophysiology such as wound healing or cancer metastasis [1][2][3] . Conceptually, the dynamics of the migrating cells have been best understood in epithelial monolayers 4,5 which are the simplest tissues that line organs throughout the body 6 . As cells cooperatively move together 7,8 , each cell needs to establish its own polarity: confining certain biochemical processes to the front of the cell while others to the rear [9][10][11] . The most fundamental hallmark of cell polarisation manifests in cytoskeletal dynamics: polymerisation of F-actin at the leading edge, a process coordinated by a variety of intracellular signaling molecules 10 . Actin polymerisation and the activity of myosin II molecular motors produce cytoskeletal flows which translate to cell motility 1,7,[12][13][14][15][16] .
The synchronisation of cell polarity during collective migration as well as the guidance from mechanical environmental cues are poorly understood 17 . The process likely involves cell-cell junctions mediated by adhesion proteins (e.g., transmembrane receptor E-cadherin) coupled to the contractile actomyosin cytoskeleton through cytoplasmic adaptor proteins (e.g., α-catenin) 4 . Actomyosin contractility at the cell cortex can modulate both adhesion strength 18 , and cell shape, both exerting an important influence on tissue remodelling 1,6 . A common experimental system to study collective cell migration in vitro is the wound closure assay, in which a barrier, the "wound", divides a monolayer culture. After removal of the barrier, cells migrate into the empty region and eventually create or restore a continuous monolayer of cells. Live imaging of such assays [19][20][21] indicated that the number of cells in the cell-free region increases, mainly due to active cell migration and not proliferation [19][20][21][22] . Furthermore, the onset of migration is delayed for the cells deep in the bulk compared to those in the vicinity of the front boundary. In several cultures the onset of motility can be observed as a polarity wave propagating backward from the leading edge of epithelial monolayer [19][20][21] .
Several theoretical models have been proposed to explain the coordination between cells during collective migration 20,[22][23][24][25][26][27][28][29][30][31][32][33][34] . However, it remains elusive how the intercellular and intracellular mechanobiology regulates the initiation and propagation of the polarisation wave through a monolayer of cells. Recently, we developed a one-dimensional model of this mechanism, which involved mechanical forces and biomechanical feedback between cells. The model predicted a traveling wave that transmits polarisation information and initiates motility in the bulk of the monolayer 35 .
The one-dimensional model of the epithelial layer was based on a very simplified representation of the cell monolayer as a chain of active particles connected by elastic springs characterised by a single parameter. However, a realistic representation of the cell sheet should include a more detailed description of cellular mechanics, such as the contractility of the cell cortex, cell-cell adhesion, and cell-area extensibility.
Here we extend our model for the expansion of an epithelial monolayer to two dimensions in order to incorporate the dynamics of cell shape as well as the polarity driven active cell motility. The new two-dimensional (2D) model is based on the theoretical and computational framework of the Cellular Potts Model (CPM). We demonstrate that this model can capture the propagation of the motility wave through the monolayer, and also allows us to analyse how the properties and patterns of cell motility are affected by different components of cellular mechanics.
This paper is organised as follows. Section 2 presents the development of the model at intercellular and intracellular levels. At intercellular scale, the model represents the dynamics of cellular shapes and interactions between cells. At intracellular scale, the model describes how the self-generated cell forces are coupled to the cell polarisation. In the next section, results are presented and discussed. First, we characterise regions of the model parameters where the polarisation propagates through the monolayer sheet. This results in the characterisation of two behaviours (different than the sheet migration): swirling motion and transient migration. After establishing a phase diagram for the model parameters, we then focus on the parameter regime where the sheet migration occurs and analyse the propagation of the motility wave through the monolayer. Finally, a summary of the findings are presented and discussed with earlier observations and also future directions.

the model
A barrier-removal assay is often utilised to study collective motility of epithelial cells. Our model represents the migration of cells toward the cell-free region, after the barrier removal -and focuses on the role of intercellular interactions and intracellular mechanics in the process. The cell-cell interactions are represented using the CPM 23,36 . The CPM is a lattice model which is computationally and conceptionally simpler than most off-lattice models (e.g., vertex model), while it provides a realistic description of cell shapes 37,38 . The intracellular polarity dynamics was formulated as a set of ODEs, coupled to each model cell of the CPM. This representation of cell motility is a generalisation of the model we used in the context of a simple one-dimensional chain of cells for the onset of collective cell movement 35 . The model is implemented using the open-source software package CompuCell3D (CC3D) 39 . intercellular dynamics. The CPM represents the cells on a two-dimensional lattice, where each cell covers a set of connected lattice sites or pixels; each pixel can only be occupied by one cell at a time. In this paper, the lattice is a rectangular surface (1500 and 320 pixels in the xand y-dimensions, respectively). The expansion and retraction of the cell boundaries is determined by minimising a phenomenological energy or goal-function E, defined in terms of the area σ A and perimeter σ L of each cell σ of N cells (indices σ = … N 1, , ) 33,[40][41][42][43] , and the motile force → σ F modeled as: The first term of Eq. (1) models the compressibility of cells by penalising the deviation of cell areas from a preset value A 0 , 100 pixels (see Table 1). The second term represents the contractility of the cell cortex as a spring with zero equilibrium length. The penalty parameter λ cont represents cortical actomyosin contractility, around the lateral cell membrane 44 . The third term describes the cell-cell adhesion mediated by adhesion molecules, such as E-cadherin. J is the boundary energy cost at neighbouring lattice sites → i and → j . The Kronecker δ function prevents counting pixels that belong to the same cell. When both lattice sites → i and → j correspond to cells, adh , otherwise when one or both lattice sites represent empty space or boundary wall the boundary energy cost J is set to zero. Note that λ < 0 adh to represent that cells preferentially expand their boundaries shared with neighbouring cells. This is however balanced by the contractile tension along the cell cortex. The last term represents a motile force driving cells into the direction of  center of mass. The prefactors λ area , λ cont , λ adh , and F max reflect the relative importance of the corresponding cellular properties to set cellular morphology. The dynamics of the CPM is defined by a stochastic series of elementary steps, where a cell expands or shrinks accommodated by a corresponding area change in the adjacent cell (or empty area) 36,39 . The algorithm randomly selects two adjacent lattice sites → i and → j , occupied by different cells σ σ ≠ → → i j . The elementary step is an attempt to copy σ → i into the adjacent lattice site → j , which takes place with probability where ΔE is the change in functional (1) due to the elementary step considered, and the temperature parameter T is an arbitrary scaling factor. A Monte Carlo step (MCS) of the simulation, the natural unit of time in the model, is set to L elementary steps -where L is the total number of lattice sites in the simulated area 39 . Together, Eqs. (1 and 2) imply that cell configurations which increase the penalties in functional (1) are less likely to occur. Thus, the cell population evolves through stochastic rearrangements in accordance with the biological dynamics incorporated into the effective energy function E.
In the initial condition, the area of each cell was set to the equilibrium value A 0 , i.e., the size of each cell is 10 × 10 pixels. The total number of = N 3000 cells are placed in the domain: 100 × 30 cells in the xand y -dimensions, respectively. Then, the cell-free region is a 500-pixel-wide (in the x-dimension) empty region. This is long enough that the leading cells will not reach the end of the empty region before the end of the simulations. The surrounding wall cells are used to prevent the cells from sticking to the lattice boundaries. The barrier and wall cells (each is 10 × 10 pixels) have the CC3D "Freeze" attribute and thus, they are excluded from participating in the pixel copies of the Potts model 45 . intracellular dynamics. Active cell motility is powered by cytoskeletal dynamics which is regulated by cell polarity, a spatial imbalance of signaling molecules 1,17,46 . Following 35 , we represent cell polarity as a vector quantity → σ p , and we assume that the motile force is a nonlinear (Hill) function of polarity: n n n max with half-saturation constant α > 0, and maximal motile force F max . To describe the dynamics of cell polarisation, we adopt the earlier models 24,25,35 similar to the one recently used in 21 as: where 1/β is the characteristic persistence time of polarisation and the second term represents the reinforcement of polarisation through actual movement 15,47 , where → σ v and γ are the velocity of the center of mass of the cell σ and a characteristic scale of the velocity, respectively. This is qualitatively equivalent to earlier models in which cell polarity aligns with cell velocity due to the inherent asymmetry created in a moving cell 25,29,48 ; reviewed in 26 . Thus, the displacement of the cells is determined by the intracellular motile force → σ F in combination with other intercellular interactions represented in Eqs. (1)(2)(3)(4).
To mimic the initial migratory stimulus through the presence of the cell-free region to which the leading cells (i.e. the first row of cells) are exposed when the barrier is lifted, we assign an initial polarity to the leading cells to trigger their motility. This is consistent with earlier studies showing that the motile force production is initiated at the first few rows of the leading edge and travels backward into the monolayer 19,20,49,50 .
The magnitude of the initial polarity of the leading cells is set as the steady-state polarity of a migrating single cell; see Fig. 1(a). Due to the strong nonlinearity of Eq. (3), the model single cell motility exhibits a bistable behaviour as shown in 35 ; see also in Fig. 1(a). Such bistability has been also found in more detailed PDE models 51,52 . Recall that the cell velocity (and thus polarity) is determined by the motile force → σ F and intercellular interactions represented in Eqs. (1)(2)(3)(4). Thus, when a cell is initialised with polarity greater than the threshold parameter α, it will tend to move with a stable velocity and polarity. Conversely, if the initial polarity is lower than α the cell gradually looses its polarity and motility; see Fig. 1(a). Such bistable behavior has been experimentally observed in 53,54 . Thus, following the experiments 19,20,49,50 , at the onset of migration, when the barrier is removed, we set the polarity of the leading cells equal to the steady-state value; see Table 1. We set the maximum motile force = F 1000 max , at which according to the single cell simulations the cell velocity starts to saturate; see Fig. 1(b-d) and Table 1.
The single cell velocity can be estimated by considering the motility term in the energy function. A pixel change which moves the cell center in the direction of the cell polarization, i.e. reducing the energy, leads to a displacement of the cell center δ ∼ ∼ − r R N N / c 1/2 where = N A 0 is the number of pixels in the cell and ∼ R N 1/2 is the cell radius. The number of such energy reducing pixel changes in a single MCS timestep of the CPM model is proportional to the cell perimeter (~N 1/2 ). Thus, the estimated single cell velocity, i.e. displacement per unit time is: , which is comparable to the numerical values obtained in the simulations (~0.4-0.5). In order to estimate the required magnitude of the motile force coefficient F max , we have to take into account that individual pixel changes also lead to temporary changes (increase or decrease) of the cell area and perimeter. For the parameter values we used, the area constraint appears to be the stronger one. Therefore, for single pixel www.nature.com/scientificreports www.nature.com/scientificreports/ change attempts that move the center of mass in the direction of the cell polarization, to be successful, the decrease of energy due to the motile force need to be larger than the energy increase due to compression or expansion of the cell by one pixel relative to the reference value A 0 . Thus, we have . This is qualitatively in agreement with the range of parameter values identified in the numerical simulations.
In order to estimate the spatial and temporal units of our model, we can use a typical cell diameter of ~20 μm which means that the length unit, i.e. pixel size, of our model is ~2 μm. The typical value of single cell velocity is 20 μm/h corresponding to ~0.5 pixel/MCS. Thus, we can estimate the time unit MCS of the model as approximately 3 minutes. With this value the duration of our cell sheet simulations, in the range of 300-500 MCS, corresponds to 15-25 hours, which agrees with the typical timescale of scratch assay experiments.
We run each simulation in three stages. The first stage yields equilibrium cell shapes, starting from a uniform square grid initial configuration. The duration of this step is 2900 MCS, and the motile force is switched off = F 0 max . During this stage cell shapes are determined by the competition between cell-cell adhesion and the contractility of the cell cortex (Fig. 2) as previously reported 33,40,55 . In the "hard" regime, when contractility is strong, the cells tend to minimise their perimeter and exhibit shapes that are close to hexagons. Conversely, in the "soft" regime cell-cell adhesion dominates and the cells have irregular shapes with elongated boundaries.
The second stage (for 100 MCS) simulates the full dynamics of the cells, i.e., the cell polarisation coupled to the CPM through the motile force, while the cells are still confined behind the barrier. We denote the end of this stage as = t 0; see Fig. 3(a). Finally, at the onset of stage 3, the barrier is lifted, the first row of cells are polarised; see Fig. 3(b,c). The simulation ends when the polarisation wave traverses the domain.  Table 1.

Results and Discussion
We use the model described in the previous section to systematically investigate how the various model parameters representing cellular mechanics affect collective cell movement. We start with establishing parameter regimes where different collective behaviours are exhibited. This allows us not only to analyse the propagation of polarisation wave, but also to describe other collective cellular behaviours.
One interesting observation is that in a certain range of parameters, when the motile force in the CPM is switched on, cells spontaneously polarise and generate a swirling movement even in a confluent layer without a free edge (i.e., before the barrier is lifted). Typical cell trajectories corresponding to this type of characteristic swirling motion are shown in Fig. 4(a), top. The spontaneous polarisation and swirling motility happen when the contractility of the cell cortex is relatively weak and/or the polarisation threshold α is low. Otherwise, the cells remain unpolarised and almost stationary (Fig. 4(b), top). The corresponding phase diagram in the model parameter space is shown in Fig. 5(a).
After the barrier is lifted, we observe that the cells migrate into the cell-free region; see Fig. 4(a,b), bottom. In the regime where spontaneous polarisation occurs, the invasion and the swirling movements coexist -the latter is prominent in the bulk, further away from the moving edge; see Fig. 4(a), bottom. Such coexistence of the swirling and directed migration was indeed observed in experiments 56 .
To assess the effect of the swirling on the overall alignment of the migrating cells, we use an order parameter defined as: Figure 2. Cell shapes at a stochastic equilibrium state, for various regimes of cortex contractility λ cont and adhesion λ adh coefficients. Other simulation parameters are listed in Table 1. www.nature.com/scientificreports www.nature.com/scientificreports/ where Φ varies between 0 (uncorrelated random movement) and 1 (fully aligned migration). We find that the swirls reduce the overall alignment of the migration, captured by a decreased order parameter Φ; see Fig. 4(c). In the absence of the swirls, cell alignment gradually increases as the polarisation wave propagates into the monolayer; see Fig. 4(d). Such increase in the alignment of migrating cells with the propagation of the motility wave has been observed experimentally in 20 . Interestingly, the directed collective migration into the free space requires that the motile force generated by a polarised cell (F max ) is sufficiently strong. Otherwise either the polarisation cannot propagate backward into the bulk or the cell layer cannot expand. In either case the migration stops after a short transient period. The minimal value of the motile force required for the propagation of the polarisation wave into the bulk and for sustained directed migration increases with the cortical contractility of the cells; see Fig. 5(b).
We now focus on the parameter regime which keeps the cells stationary and unpolarised until the barrier is lifted. If the motile force is strong enough, such cells become motile and migrate into the free space. The progression of the average position of the cells is shown in Fig. 6 for various parameter values. To calculate the average position of cells, we bin the domain along the x-dimension, where each bin is five-cell wide, and calculate the x -component of the average position of cells over time. Figure 6 shows a progression of the leading edge with constant speed, in agreement with experiments 56 . The stationary cells are recruited into the collective migration at a constant rate by the propagating polarisation wave -as predicted by our previous one-dimensional (1D) model 35 . The slope of the dashed line shows the propagation speed of the polarisation wave, which decreases when the cell contractility is increased; compare the slopes of the dashed lines in Fig. 6(a,b).
In order to investigate how the speed of the polarisation wave and the expansion velocity of the edge cells change as a function of model parameters, we ran a series of simulations and the results are summarised in Fig. 7. The polarisation wave speed and the velocity of the leading edge are most strongly affected by the cortex contractility. Increasing the cortex contractility decelerates the propagation of the polarisation wave and also the www.nature.com/scientificreports www.nature.com/scientificreports/ velocity of the leading edge. Parameters representing cell-cell adhesion and compressibility have much weaker effects (Fig. 7). Note, that the overall value of the edge velocity is roughly similar to the single cell velocity as in the case of the 1D model 35 .
At low cortex contractility the edge velocity (~0.3) is not sensitive to the presence or absence of swirling motion (Fig. 7(b), filled symbols). This finding is consistent with the experimental observation that the average   Table 1 www.nature.com/scientificreports www.nature.com/scientificreports/ velocity of the migrating cell front did not show much variation in the presence of swirls 56 . Consistently with earlier reports 57 , Fig. 7(c-f) also indicates that cell-cell adhesion and the area compressibility have little effect on the polarisation wave speed and the edge velocity.
For certain parameter values the simulations indicate patterning instabilities at the free edge where some of the leading cells leave the monolayer. We find that such behavior develops when the cortex contractility is strong. When the cortex contractility is weak (e.g., λ = 2 cont ), the energy cost for the advancement of the leading cells is low, resulting in a smooth migration of the monolayer; see Fig. 8(a). However, when with the cortex contractility is strong (e.g., λ ≈ 6 cont or larger), the displacement of the leading cells is somewhat restricted, the stochastic fluctuations are amplified resulting in the development of instabilities at the free edge; see Fig. 8(b,c). This is analogous to the experimental observations 44 that the tension in the cortical actomyosin ring prevents the initiation of new leader cells.

Summary and outlook
In this paper, we studied the collective migration of epithelial cells in a confluent monolayer exposed to free space. Our simulations identified that the cortex contractility, not adhesion of the cells is a key model parameter that controls the transitions between swirling movements and well aligned cell sheet migration. Strong contractility can destabilise the advancement of the leading edge, block the propagation of the polarisation wave, and inhibit the collective coherent migration.
Our model also indicates that the cells cortex contractility plays an essential role in regulating the formation of the swirling movement, i.e. the swirls are more likely to form when the cortex contractility is weak; see Fig. 4(a,b). This is explained using Eqs. (1)(2)(3)(4). Cells with a weak cortex contractility (e.g., λ = . 0 3 cont ) can more easily change their shape which can also lead to the displacement of the center of mass of the cells. This displacement can spontaneously polarise the cell, due to the bistable behaviour of the cell polarity → σ p ; see Fig. 1(a). The spontaneous cell polarisation becomes more likely when the polarisation threshold α is lower. This resembles the experimental observations showing that the intensity of the polarisation promotes the appearance of the swirls 58 . Then, the motile force → σ F generated by a spontaneous polarisation is transmitted to the neighbouring cells leading to the formation of swirling movements. This suggests that sufficiently strengthening the cells cortical contractility would increase the overall alignment of the migration by transforming the swirling movements into a directed cell sheet migration; see Fig. 4(a,b).
On the other hand, excessive strengthening of the cells cortex contractility can block the collective migration in a monolayer; see Fig. 5(b) and Movie 5. A tight cortical contractility results in symmetric hexagon-like cell shapes corresponding to the "hard" regime ( Fig. 2) that blocks the cells polarisation by limiting their deformation and motility. The "hard" cell shapes can only have relatively little displacement of the center of mass, corresponding to weak cell polarity and motile force. This blockage of the active motility can also be caused by increasing the decay rate of the cell polarisation β. At high β, the cells do not possess non-zero steady-state for the polarity and , the cell-cell adhesion coefficient λ adh (c,d), and the area strain λ area (e,f). To calculate the edge velocity, we bin the domain along the x-dimension, where each bin is three-cell wide, and calculate the average velocity of cells in the front bin which is connected to the monolayer and the cell-free region. Each symbol is derived from a single simulation run and corresponds to mean ± SD. In (b), the filled symbols are at λ cont = 0.3 and 0.6 where the swirling movements are formed. Other simulation parameters are in Table 1 www.nature.com/scientificreports www.nature.com/scientificreports/ their polarity magnitude can only stabilise at zero; see Fig. 1(a) and our discussion in 35 . Then, starting with an initial polarity, the cells gradually loose their polarity and stop moving. Therefore, the polarisation wave becomes unstable and disappears gradually after the barrier is lifted; see Movie 6.
In the present study, we focused on the migration of cells in the absence of cell proliferation, and where all the cells had the same constant mechanical characteristics. Future works may consider the effects of the cell proliferation combined with active cell migration on the collective movement and the propagation of the polarisation wave. It is also interesting to consider the adaptive responses of cells to the environment, where the cellular properties can vary. For instance, the alignment of the cells' migration may be enhanced by coupling the cortical contractility to the polarisation threshold in order to inhibit the swirling movements; see Fig. 5(a,b). Likewise, recent studies on the efficiency of wound healing have shown that a monolayer can effectively cover the wound region and achieve a permanent gap closure when the instabilities of the leading edge are suppressed 59 . Accordingly, a coupling of the cortical contractility to the morphology of the leading edge in order to maintain a moderate cortical contractility could be studied in the context of a wound closure assay.

Data availability
The computer code is available to download from the website at https://github.com/hr-khataee/2DModel. The code has been developed with CC3D version 3.7.8.