Modeling collective cell migration in geometric confinement

Monolayer expansion has generated great interest as a model system to study collective cell migration. During such an expansion the culture front often develops ‘fingers’, which we have recently modeled using a proposed feedback between the curvature of the monolayer’s leading edge and the outward motility of the edge cells. We show that this model is able to explain the puzzling observed increase of collective cellular migration speed of a monolayer expanding into thin stripes, as well as describe the behavior within different confining geometries that were recently observed in experiments. These comparisons give support to the model and emphasize the role played by the edge cells and the edge shape during collective cell motion.


V Tarle et al
To eliminate spurious forces arising near the entrance to the stripe, we do not calculate the contour forces for the cells that are along the vertical part of the confinement (between stripes), as described in the SI (stacks.iop.org/ PhysBio/14/035001/mmedia). The parameters used in the simulations are given in table S1, and were used in [8,14] to fit experimental data of the same MDCK cells used in [16]. The simulations start with randomized initial conditions (random locations of the cells within the initial area, figure S1), as explained in [8,20] and SI. The curvature of the edge of the monolayer is calculated as in [8]. In the SI we give more details regarding the model, while a complete technical description appears in [8].

Stripe width effect on expansion velocity
In figure 1(D) we show a snap-shot from the simulation, where cells expanded from an initial straight configuration into several parallel stripes of different widths. In these simulations, we used values for the positive feedback of F border with curvature, such that we are in the regime where the edge is unstable to the formation of fingers [8]. In the thinner stripes it is apparent that the geometry is defining the formation of a single leader cell, while in the wider stripes several leader-cells can initiate. In figure 1(C), we plot the corresponding motile forces of the edge cells (F border ), which shows the main effect of the confinement: in wide stripes there are usually two leader-cells where the leading edge contour of the monolayer curves due to the confining borders (see where F border is maximal). In the thinner stripes the confinement naturally selects a highly curved cell at the front of the column, and due to the curvature-motility feedback this cell becomes a highly motile leader-cell. We identify leader-cells as particles at the leading edge that have a convex curvature that is larger than some threshold value [8].
Recently the formation of leader-cells at the edges of the confining stripe were identified [21], as our model predicts.  [8], as seen here after 22 h (scales are in microns, pink particles denote leading edge contour). ((C)-(E)) Snap-shot at time 22 h from a simulation of the expanding monolayer confined into stripes of different widths. Note that the simulations were done for a single channel, or pairs of channels, as indicated by the black horizontal lines in (D), and we placed all the results next to each other for presentation purposes, as was done in [16]. In (C), we plot the vectors of F border . (D) shows the particle positions and indicates the leading-edge cells in pink, and leaders by bold purple circles. (E) shows the velocity field, where large vortices are observed in the wide channels. Note that in our particle-based simulations each point-particle represents a cell of unknown shape. The leading edge cells are especially highly spread and have large lamellipodia, which is why the stripes are more fully covered by the cells in the experimental images [16] compared to the images of the particles produced by the simulations (figures 1(C)-(E)). In the simulations we suppressed cell divisions in order to extract the effects of geometric confinement alone. With divisions restored (figure S6) the channels get slightly filled more fully, but qualitatively the behavior is unchanged, highlighting that the width-dependence of the migration is not driven by cell proliferation.
We next quantified the effect of the formation of leader-cells on the dependence of the collective migration on the stripe width. We denote the leading edge of the monolayer in these images (figure 1(D)), which allow us to define the average location of the front and calculate its rate of propagation (see SI and figure S2). This is shown in figure 2, for different values of F min and F max . We find that when F max is reduced (figure 2(A), (D) and (G)), or F min is increased (figures 2(C) ,(F) and (I)), such that the force-curvature slope is too low to trigger the instability transition, the increase in velocity for thin stripes almost disappears. When both F min and F max are increased (figures 2(B), (E) and (H)), the force-curvature slope remains constant and therefore leader-cells are produced, but their effect on the overall migration decreases when F min is large. The reason is that the maximal forces applied at the leader cell are saturated, and therefore the relative difference between the force on the leader and the cells near it decreases.
Note that even without the appearance of leadercells, i.e. no curvature-motility feedback (figures 2(A) and (D) dashed black line), there is an increase in the average migration speed in the narrower stripes simply due to the confinement that aligns the cellular motions and traction forces [22]. We calculated the orientation order-parameter, which is the mean cosine angle of the cellular motion inside the stripes θ = S cos ⟨ ( )⟩ [12], and find: ∼ S 0.9 for stripes of width 20-µ 50 m, while it decreases to ∼ S 0.7 for stripes of width 125-µ 400 m. These numbers are very similar to those observed in the experiments [16,23], but the alignment order param- eter becomes very noisy in the wider channels, due to the large scale vortices in the cellular flows (figure 1(E)).
When cells are treated with blebbistatin [16], it was observed that the cellular migration becomes faster for wide stripes and slower for the thin stripes, and therefore overall the migration speed becomes largely independent of the width. This is similar to the behavior we find in the simulations for increasing F min in figures 2(E) and (F): compare the nominal behavior Indeed it was observed previously that blebbistatin treatment gives rise to faster expansion of cell monolayers, with decrease in the number of leader-cells and finger formation [24,25].
Comparing to our simulations we conclude that the blebbistatin mainly affects the cells that are either flat or concave along the edge, causing these cells to have larger outwards motility [8]. When F min is increased the gradient in the outward forces between the leader (highly curved tip) cell and the surrounding cells is smaller, and the column of cells is thicker compared to the case of smaller F min (inset of figure 2(F)). The thinner column of cells obtained for the larger curvature-force slope (stronger tendency to form fingers) maintains a higher curvature at the tip and therefore a faster forwards migration. This is the essence of the mechanism behind the width-dependent migration velocity, driven by the curvature-traction force feedback. Note that the concave and flat edge cells are often observed to form acto-myosin contractile cables, and when blebbistatin interferes with the ability to produce contractile forces, the cytoskeleton in these cells converts into migratory lamellipodia. The loss of the actomyosin cables slightly decreases the migration speed in the thinnest stripes (figure S5), further contributing to the slowing down observed in the µ 20 m stripes when blebbistatin is applied [16].
The above observations from our simulations, namely that the velocity depends on the stripe width due to leader cells that are initiated by the geometric confinement, can be approximately captured by the following simple analysis. The mean velocity of the monolayer front v can be approximated as arising from a simple average of the total motile forces acting along the edge, by the leader and non-leader cells. This gives an average front velocity, in the form (1) where w is the width of the front (stripe), a is the size of the leader along the edge, α = > v v 1 leader 0 / is the ratio of the contribution to the front progression velocity of the leader and non-leader cells. If the the average length-scale of the most unstable mode driven by the curvature motility feedback (average distance between leader cells) is λ [9], then we expect the total number of leader cells along the front n leader to have the form: λ + n w 1 leader / , where the 1 arises from the confinement-induced leader cell that exists even in the limit of w 0 → . The width dependence of the mean number and density of leaders, as well as the average pulling force density exerted by the leading edge cells (F border ), are shown in figure S7. In the limit of the thinnest stripes we expect a single leader cell, simply due to the confinement-induced high curvature effect. Our model's prediction of a constant density of leader cells per unit length of monolayer front (for wide channels) was recently confirmed experimentally [21]. Comparing the simulated density of leaders ( figure  S7(b)), and the migration speed ( figure 2(D)), we find that indeed the migration speed seems to be dominated by the density of leader cells, validating the assumption leading to equation (1) (figure 3), and as proposed in [21].
In figure 3 we plot the dependence of the mean front velocity on the stripe width from equation (1). We find that when there is a strong difference between the propulsive force of the leader and the non-leader cells, we obtain a steep increase of the velocity for think stripes. When either the leaders produce weak forces, or the non-leaders produce large forces, the overall width dependence of the front velocity disappears. The case of Phys. Biol. 14 (2017) 035001 fast non-leader cells fits qualitatively the observations when cells are treated with blebbistatin [16].
Note that the concave regions of the monolayer edge contribute a forward pushing force due to the contractile forces of acto-myosin cables (F cables , figure S4). Stronger contractility of these cables induce faster propagation of the leading edge in our model [8]. However, we found that even when the cables are removed we obtain very similar width dependence of the migration velocity (see figure S5). The main mechanism for the width-dependent velocity is therefore attributed to the formation of leader cells.

Saltatory motion in thin stripes
In [16] it was noted that the motion along the thinnest stripes does not advance smoothly, but rather through periods of fast and slow (almost stopping) motion forward. The experiments were described in detail in [16], and a brief description is given in the SI. In figure 4, we compare simulation and experimental kymographs of the cells' velocity as the column advances along the stripe. In both cases we find this saltatory motion. In the simulations this occurs due to the leading cell experiencing a strong forward traction force F border , due to the curvature-motility feedback ( figure 1(A)). The cells behind this leader do not experience this force directly, and are following the leader due to their mutual attractive forces and the tendency to align their motion (Vicsek interaction). However these processes occur with some time delay during which the leader has already moved forward faster than its following cells. The leader is therefore subsequently stopped by the attractive interaction with the cell behind it, until these follower cells catch up and the cycle renews. The result is the observed saltatory motion.
We explored systematically the effects of the cellcell interactions on the frequency of the events of slowing down and re-accelerations that give rise to the saltatory motion inside the narrowest stripes. The results are given in figure S8. As expected, increasing the cell-cell adhesion strength (U 3 in our model, see SI, equation (S1)) leads to lower frequency of these events ( figure S8(a)). Increasing the strength of the alignment interaction (β in our model, see SI, equation (S4)), similarly increases the correlation in the cell movements and decreases the rate of saltatory events ( figure S8(b)). Finally, the rate of saltatory events is found to be maximal when the adhesion range between cells is between nearest-neighbors but long enough to prevent the leading cell from detaching ( figure S8(c)). The tendency to detach from the column of cells will be explored elsewhere, and here we restricted the simulations to regimes where this behavior was minimal.

Expansion into converging and diverging stripes
Recently the collective migration of cells in converging and diverging stripe geometries was investigated [21], and particular attention was given to the possible role of leader cells. It was found that in the diverging stripe the migration rate is much slower compared to the converging case, while the rectangular geometry leads to an intermediate value.
We simulated, using our model, these geometries ( figure 5(A)) and found a similar trend (figures 5(C) and S1(B)). We next utilized our simulation to explore the role of the leader cells in this phenomena. When the model includes the curvature-motility feedback, we indeed find leader cells produced at the two edges of the stripe, and appearing spontaneously along the leading front (figures 5(B) and S2). However, when we used a constant F border which is independent of curvature, we found that the differences in the migration rate between the three geometries were preserved ( figure  5(D)). This suggests that the appearance of the leader cells is not the main driving force for these differences in migration speed. We also checked the contribution of the contractile forces of acto-myosin cables at the leading edge to the migration speed (F cables ). We found that the cables give a minor effect, mainly to the converging geometry ( figure 5(E)) where there are stronger concave curvatures ( figure 5(A2)).
The main mechanism that is driving the differences between the migration speed in these three geometries, within our model, turns out to be simply the Bernoulli effect of flow within a slab with varying cross-section. This is demonstrated in figures 5(F)-(H), where we plot the cell density ρ ((F1)-(F3)) and velocity v x along the channel axis ((G1)-(G3)), averaged over the channel cross-section at different positions along the channel (w(x)). We finally plot the normalized cell flux (figures 5(H1)-(H3)). The  (a) Model prediction that the two leader cells moving up the sides of the stripe (black circles and arrows, connected by black curve that represents the leading edge), induced by the confinement and the curvature-motility feedback, will continue into the expanded region and cause highly anisotropic invasion that is elongated along the boundary. The high curvature of these edge cells is indicated by the red arcs. (b) Scheme of the geometry of the stripes used in the experiments. (c) An example of a simulation using our model, where we kept the edge cells in contact with the boundaries by using a rounded corner (radius µ 30 m), and with the addition of an attractive force at the edge of the adhesive region (force strength µ 800 m h −2 , acting up to a distance of µ 9 m from the edge of the boundary). The arrows indicate the force acting on the leading edge cells, which provides the curvature-motility feedback. Parameters as in table S1, except for: T 0 = 6 h, curvature-motility feedback slope of µ 15 000 m 2 h −2 . (d) Snap-shots of the experimental results (at the times indicated on the images) of the cells emerging from the µ 84 m stripe into the wide square region (see in (b)), confirming the schematic prediction made in (a), and the model simulations of (c). The experimental conditions are as in [16], and are described in detail in the SI. Scale bars in (c) are µ 100 m. last quantity measures the overall cell flow through each cross-section of the channel. The flux of cells is given by where all the quantities depend on the position x along the channel. As this overall flux is found to be constant along the channel, it clearly demonstrates that the increase (decrease) in the flow for the converging (diverging) geometry arises mainly due to the conservation of mass condition.

Conclusion
We have shown in this work that our model for expanding monolayers [8], which incorporates the curvature-motility feedback for cells at the culture edge [9], can explain the observed dependence of collective migration speed on the geometry of the confining boundaries. The higher migration speed for thinner stripes, which remained a puzzle since its observation [16], is shown to be naturally related to the curvaturemotility feedback. This provides a strong support to our model of collective cell migration.
Using the insight gained by our model, that the confining stripe induces two leader-cells where the monolayer is curved by the boundary, we reasoned that these two leaders should give rise to a very anisotropic expansion of the cells if the channel is suddenly expanded in width. This is schematically illustrated in figure 6(a), where the high curvature of the cells at the edges of the adhesive strip are indicated. Unfortunately we found that the particle-based simulations we used above, do not allow the leader cells to follow a sharp L-turn, as the model does not account for the deformable, amoeboid-like shape of the cells which do engulf the corner. The simulations with the sharp corner showed the cells continuing to move forward, which is the persistent motion expected due to the alignment interactions. In order to make the particles in the simulations maintain contact with the edge of the adhesive region, we made the corner round (figure 6(c)). We also added a force acting to pull cells that are close to the boundary (closer than µ 9 m) towards the boundary. This force captures the observed tendency of cells to move along the edge of an adhesive region, by forming enhanced adhesions along this edge (Micahel Sixt, private communication). Note that this force is normal to the edge, so does not in itself contribute to faster cell motility, which is induced for the leader cells by our curvature-motility feedback model. With these changes implemented, we find that the simulations indeed follow qualitatively our schematic prediction, with leading edge cells at the boundary remaining highly curved and leader-like (figures 6(c) and S8). Remarkably, we found that cells in an experiment (overall geometry shown in figure 6(b)) follow closely our theoretical prediction, as shown in figure 6(d).
We therefore demonstrate that our model captures the collective migration patterns in a variety of confining geometries, thereby demonstrating that curvaturemotility feedback plays an important role in collective cell migration phenomena.