article: In cancer metastasis, embryonic development, and wound healing, cells can coordinate their motion, leading to collective motility. To characterize these cell-cell interactions, which include contact inhibition of locomotion (CIL), micropatterned substrates are often used to restrict cell migration to linear, quasi-one-dimensional paths. In these assays, collisions between polarized cells occur frequently with only a few possible outcomes, such as cells reversing direction, sticking to one another, or walking past one another. Using a computational phase field model of collective cell motility that includes the mechanics of cell shape and a minimal chemical model for CIL, we are able to reproduce all cases seen in two-cell collisions. A subtle balance between the internal cell polarization, CIL and cell-cell adhesion governs the collision outcome. We identify the parameters that control transitions between the different cases, including cell-cell adhesion, propulsion strength, and the rates of CIL. These parameters suggest hypotheses for why different cell types have different collision behavior and the effect of interventions that modulate collision outcomes. To reproduce the heterogeneity in cell-cell collision outcomes observed experimentally in neural crest cells, we must either carefully tune our parameters or assume that there is significant cell-to-cell variation in key parameters like cell-cell adhesion.
  Collective cell migration, in which cells crawl together in a coordinated way, is important for a wide range of biological functions, including the immune system [1], the early development of embryos [2] and cancer invasion [3]. The cell-cell interactions that drive collective behavior can also create new dynamics not seen in isolated crawling cells. This includes fingering instabilities [4–6], cell clusters that exhibit emergent chemotaxis in the absence of individual cell chemotaxis [7–10] or groups of cells which can be governed by a few leader cells [11]. How this collective behavior emerges is an ongoing and active area of research [12,13]. One well-studied way in which cells interact with one another to create collective motion is contact inhibition of locomotion (CIL), which was named and characterized by Abercombie more than 50 years ago [14], with observations dating back as far as 1921 [15]. CIL describes the tendency of a cell to change its direction of motion after contact with another cell. CIL can play a central role in the coordination of collective migration of cells, including the cluster-level polarization of neural crest cells in vivo [16] and the dispersal of hemocytes in the early stages of Drosophila development [17]. CIL may also play a role in cancer [18]. A powerful and commonly used technique to analyze CIL in vitro is the collision assay, where the velocity of pairs of cells before and after collision is measured [16]. However, these assays can have low efficiency, since cell-cell collisions are rare. Recently, several groups have developed 1D collision assays, where cell motion is restricted to an adhesive micropatterned stripe, increasing efficiency and ensuring a reproducible collision geometry [19–22]. These assays, originally used to study cell motility in the presence of confinement [23,24], can be used to study outcomes of cell-cell collision and to identify critical molecular mediators of CIL [20,22,25,26]. The experiments show that head-on collision of two cells can result in four possible outcomes: [19,20]: In the case of Xenopus cranial neural crest cells, Scarpa et al. were able to analyze a large number of cell-cell collisions and to generate quantitative statistics for the possible outcomes [20]. These experiments reveal that the majority of cell collisions resulted in reversals, a smaller fraction of collisions resulted in sticking, walk-past was uncommon and chaining was not observed (see Table 1). However, chaining-like behavior (cells following one another on contact) was observed in vivo in chick cranial neural crest cells [27]. In this work, we will use a computational model to evaluate the chemical and mechanical factors that control the interactions of motile, eukaryotic cells exhibiting CIL on narrow micropatterned stripes. Our approach, which extends the phase field technique presented in [28], models the mechanics of a changing cell shape as well as the biochemical polarity of the cell. The primary biochemistry we are interested in is cell polarity—i. e. what determines the cell “front” and thus the direction in which the cell is migrating. To do this, we use a minimal model of Rho GTPase kinetics [29] to describe the dynamics of a polarity protein, which we will assume is Rac, which is activated at the front of the cell. CIL is incorporated through the assumption that cell-cell contact leads to local generation or activation of a Rac inhibitor. We will simulate a large number of cell-cell collisions on narrow micropatterned stripes. Our aim is to probe which properties in the cell-cell interactions are responsible for the different outcomes observed in the experiments.
  The interface of a cell i is tracked by an auxiliary phase field ϕ (i) (r, t), which varies smoothly between ϕ = 0 (outside of the cell) and ϕ = 1 (inside) over a length scale ϵ; the cell interface is given by the contour ϕ = 1/2. Assuming any fluid flow can be neglected and that the interface is only driven by local forces, the motion of the cell interface is given by ∂ t ϕ (i) (r, t) = 1 τ α ρ (r, t) χ (r) - β | ∇ ϕ (i) | - 1 τ ϵ δ H δ ϕ (i) (1) from a force balance argument [28,30]. Here, τ is a friction coefficient. A full set of parameters and their values is given in S1 Table. We note that many groups have recently modeled both single [30–36] and collective [28,37–39] cell motility with phase fields. The first term on the right hand side of Eq 1 describes the active motion of the cell, arising from forces caused by actin polymerization at the leading edge and myosin-driven contraction of the cytoskeleton at the cell rear [40]. This arises because the first term of Eq 1 pushes the cell front outward where ρ is large (αρ > β) and contracts at the rear where ρ is low (β > αρ). To model the effect of the adhesive micropattern, we assume that cells only create protrusion if they are able to adhere to the underlaying substrate. This is implemented by including χ (r), which takes on values between χ (r) = 0 (cell cannot adhere) and χ (r) = 1 (cell can fully adhere). A more in-depth motivation can be found in [41]. In absence of the active motion term in Eq 1, the phase field ϕ will minimize a Hamiltonian H = Hsingle + Hcell−cell. The single cell Hamiltonian is H s i n g l e = γ (P) ∫ d 2 r ϵ 2 | ∇ ϕ | 2 + G (ϕ) ϵ + ∫ d 2 r κ 2 ϵ ϵ ∇ 2 ϕ - G ′ (ϕ) ϵ 2, (2) where γ (P) is the interface tension and κ the bending modulus. The double-well potential G (ϕ) = 18ϕ2 (1 − ϕ) 2 stabilizes the two phases ϕ = 0 (outside of the cell) and ϕ = 1 (inside). In the sharp interface limit ϵ → 0 and with a perimeter-independent interface tension, it is known that Hsingle is equivalent to the Canham-Helfrich Hamiltonian [42,43] (see discussion in [28,31]). γ (P) depends on the cell perimeter P and has the form γ (P) = γ 0 · 1 + γ p e r × (P - P c r), if P ≥ P c r 1, otherwise. The perimeter is calculated as P = ∫d2 r|∇ϕ|. Here, Pcr is a critical perimeter, and for perimeter values above this parameter cells have a component to their perimeter energy that behaves as an elastic membrane with an associated elastic energy Hel ∼ (P − Pcr) 2. For cell perimeters below Pcr the line tension is constant as is appropriate for a fluid membrane [42,44]. One reason we have added this aspect to our model is that when the cell-cell adhesion is very strong it can overcome interface tension, leading to a situation where it is energetically favorable for a pair of cells to increase their perimeter without limitation. Throughout this work we use γper = 0. 5μm−1 and Pcr = 58μm. Here, Pcr is slightly larger than the unperturbed perimeter of a moving single cell, which is P ≈ 56. 5μm for our default parameters. Note that if P increases two microns above Pcr, γ (P) doubles. We have not found any qualitative changes in the collision outcomes for different values of these parameters, as long as the growth of P without limitation is prevented. However, we did not conduct systematic variations of these parameters. The cell-cell interaction part of the Hamiltonian includes two physical interactions, volume exclusion and cell-cell adhesion: H c e l l - c e l l = ∑ i ≠ j ∫ d 2 r g 2 ϕ (i) (r, t) ϕ (j) (r, t) - σ ϵ 3 4 | ∇ ϕ (i) | 2 | ∇ ϕ (j) | 2. (3) The first term excludes volume by penalizing overlap between different cells with strength g. The second term, which models adhesion, favors contact between the membranes of different cells. The strength of this interaction is set by σ. We note that as a consequence of our phase field description of the cells, altering the strength of repulsion g or adhesion σ can also change the structure of the interface where cells overlap, i. e. how sharply the interface transitions from ϕ = 0 to ϕ = 1. This effect would not appear in a sharp-interface model. The chemical concentrations within the cell are modeled with reaction-diffusion equations of the type: ∂ t u (r, t) = D u ∇ 2 u (r, t) + f u (u (r, t) ). (4) As shown in [45,46], these equations can be solved in a complex geometry characterized by a phase field ϕ. They then have the form ∂ t (ϕ (r, t) u (r, t) ) = ∇ [ ϕ D u ∇ u ] + ϕ f u (u). (5) In the sharp interface limit, Eq 5 is equivalent to Eq 4 with Neumann (no-flux) boundary conditions. Accordingly, our equations for the membrane-bound state ρ (r, t) and the inhibitor I (r, t) in each cell i are ∂ t (ϕ (i) ρ (i) ) = ∇ [ ϕ (i) D ρ ∇ ρ (i) ] + ϕ (i) f ρ (ρ (i), ρ c y t (i), I (i) ), ∂ t (ϕ (i) I (i) ) = ∇ [ ϕ (i) D I ∇ I (i) ] + ϕ (i) f I (I (i), { ϕ }, { ρ }), (6) with Dρ, I being diffusion coefficients and fρ, I reaction terms. fρ describes the exchange between the active membrane bound state ρ (r, t) and an inactive uniform cytosol pool ρcyt (t) of the Rho GTPase Rac, in a modification of the wave-pinning scheme developed by Mori et al. [29]: f ρ (ρ (i), ρ c y t (i), I (i) ) = k b (ρ (i) ) 2 K a 2 + (ρ (i) ) 2 + k a ρ c y t (i) - k c [ 1 + I (i) (r, t) I 0 ] ρ (i) (7) ρ c y t (i) = N t o t - ∫ d 2 r ϕ (i) ρ (i) ∫ d 2 r ϕ (i) (8) with I0 = 1μm−2. The rate in Eq 7 models three basic processes: 1) cytosolic ρ binding to the membrane at base rate ka kb, 2) cooperative recruitment of ρ from the cytosol to the membrane, and 3) ρ detachment from the membrane to the cytosol with the rate kc (1 + I (r, t) /I0); the inhibitor I (r) specifies the difference in the base detachment rate kc. These processes conserve the total number of Rac molecules Ntot = ∫d2 r (ρ (r) + ρcyt) ϕ (r). Using this conservation, we find Eq 8 by assuming that the cytosolic state diffuses quickly on the time scales we study, and treating it as spatially constant. The reaction term for the inhibitor, fI, is written as: f I (I (i), { ϕ }, { ρ }) = - k - I I (i) (r, t) + f I ({ ρ }, { ϕ }) c e l l - c e l l + η ξ (i) (r, t) (9) 〈 ξ (r, t) ξ (r ′, t ′) 〉 = η 2 δ 2 (r - r ′) δ (t - t ′). (10) Here, the first term describes inhibitor decay with rate k−I and the second term models CIL through inhibitor generation by cell-cell contact (discussed extensively in the next section). The third term describes inhibitor generation by a fluctuating Gaussian Langevin noise ξ (r, t) arising from the intrinsic stochasticity of complex biochemical processes in the cell when there are small numbers of molecules involved [47,48]. We note that I (r, t) can be negative. Since |I|/I0 is small, positive (negative) I corresponds to a small increase (decrease) in the base decay rate of kc in Eq 7; the effective decay rate never becomes negative. While having a negative concentration appears to be unphysical, this is equivalent to writing an equation for Y ≡ I + I0, which is always positive. This can be seen by looking at the sharp-interface equation for the reaction-diffusion model, ∂tI = DI∇2 I − k−I I + ⋯, which can be transformed into ∂tY = DI∇2 Y − k−I (Y − I0) + ⋯, i. e. a reaction-diffusion equation with decay and a basal rate. We have written our model in terms of I to show more easily the change from kc, which can be very significant, even if the percentage change is small. The dynamics of these reaction-diffusion equations lead to a stable profile with a high value of ρ (r, t) defining the front of the cell and a low value indicating the back. The inhibitor’s effect on ρ will generally reorient the front of the cell away from high I. The resulting distribution of ρ (r) and I (r) can be seen in Fig 1 where we show a cell moving down a 1D stripe. The amplitude of the stochastic noise η controls the persistence of a single cell’s crawling motion [28]. For η = 0 the cell will crawl persistently in one direction and with increased levels of η the movement is more erratic. As mentioned above, and suggested by experimental observations [25], we implement CIL by assuming that additional inhibitor is produced where the phase fields of the two cells overlap. The production of additional inhibitor, coupled with diffusion in the cell interior, can result in cell repolarization away from contact areas [28]. As in earlier computational work [28] and based on experiments using NRK-52E cells [19] we model two distinct mechanisms for CIL, contact repolarization (CR) and front repolarization (FR), and vary their relative importance by the constants kCR and kFR, respectively. We sketch these mechanisms in Fig 2. Mathematically, the two mechanisms are described by including the following reaction terms for the i-th cell: f I i ({ ρ }, { ϕ }) c e l l - c e l l = k C R S ∑ j ≠ i N c e l l ϕ (j) ︸ CR + k F R S ∑ j ≠ i N c e l l ϕ (j) ρ (j) - O c r i t ρ c h a r ︸ FR, (11) where S (x) is a sigmoidal function, S (x) = max (0, tanh[x/x0]) with x0 = 0. 1 and Ncell the total number of cells. The first term of Eq 11 describes CR. The production of I will start as soon as there is cell-cell overlap and will cease when this overlap becomes negligible. The second term incorporates FR and, contrary to CR, also depends on the polarization of the neighboring cells: it only leads to additional inhibitor production if the i-th cell is in contact with another cell and if that other cell has elevated levels of Rac at the point of cell-cell overlap. We consider the possibility that the contact repolarization and front repolarization interactions occur by different signaling pathways, which set in at different times, e. g. if detection of cell-cell contact and cell front contact occur through different mechanisms. For instance, we could envision that cell-cell contact leads to mechanical force, leading to generation of an inhibitor at contact (i. e. CR), but that FR requires ligand-receptor binding to specifically recognize the front of another cell. These could have different rates, force requirements, or other kinetic details. We treat this possibility generically by allowing the rates kFR and kCR to differ, as well as introducing a critical overlap Ocrit in Eq 11, which sets a threshold that has to be exceeded for the FR to set in. Unless the overlap between the cells exceeds this threshold, the FR will not generate any inhibitor. When cells approach one another the values of ϕρ in the contact region are small and below this threshold. Thus, one of the effects of the critical overlap is that the FR sets in later, when the cells overlap sufficiently and not upon first contact. Therefore, the mechanical interactions between the cells are stronger when the FR sets in. If Ocrit is too large the threshold will not be exceeded and the FR does not produce any inhibitor at all. At lower values of Ocrit, FR only produces inhibitor in regions where the cells overlap reasonably and ϕρ of the other cell is sufficiently high. Hence, the distribution of the inhibitor production is different and will be concentrated to these regions. Additionally, as only the overlap which exceeds the threshold yields inhibitor production, the critical overlap will reduce the total amount of produced inhibitor for a given FR strength kFR. Typical values of ρ (r) at the cell front are around 1. 4μm−2. This is much larger than the value at the back, which is typically around 0. 01μm−2. Usually we choose values of Ocrit between 0 and 0. 2μm−2. This ensures that ρϕ crosses the threshold and inhibitor is generated when cells are in full front-front contact (i. e. when the contours ϕ (1) = 1/2 and ϕ (2) = 1/2 are close), but delays the inhibitor production when they are approaching. We emphasize that these two biochemical interactions have essential differences. For CR, cells only cease the production of inhibitor when they are not in contact anymore. In contrast, with FR a cell stops producing inhibitor when the other cells starts to repolarize. If Ocrit > 0, the two interactions will also start their onset at different times in the cell-cell contact. We note that, for simplicity, we have only included a critical overlap Ocrit in the FR interaction. Additional behaviors could potentially be found by including a similar parameter for the CR interaction. We have made this assumption both to reduce the parameter space slightly, and because we view the contact repulsion as a more fundamental effect: detecting cell-cell contact may be much simpler than detecting whether contact is specific to the front of the cell.
  Our simulations can reproduce all experimental observed cell-cell collision cases (reversal, sticking, walk past, and chaining) by varying only four parameters: adhesion strength σ, critical overlap Ocrit, and the CR and FR strengths kCR and kFR. In Fig 3, we present representative snapshots of all four cases along with the values of the four key parameters. To gain a better understanding of the relative importance of the critical parameters we carry out several parameter sweeps. These sweeps are performed by varying one or two parameters, while keeping the others fixed to their default values (which are given in S1 Table), or otherwise as noted. Most of the parameters we have varied in our parameter sweeps are not readily accessible in experiments. One experimental parameter that can easily be varied, however, is the width of the adhesive stripe, d. Consequently, we tested how narrower and wider stripes change the outcomes of our simulations. Our modeling results predict that varying d strongly affects the walk-past transition, as is shown in Fig 8 where we plot the percentage of walk-past events as a function of kCR for various values of d. Decreasing the stripe width d can strongly decrease the walk-past rate. Intuitively, this is clear: if the stripe is too narrow the cells cannot pass each other. Our simulations reveal that this transition is sharp; for example, for kCR = 0. 004s−1 the maximum walk-past rate drops from ∼60% (d = 24μm) to ∼4% (d = 23μm) to 0% (d = 22μm). Fig 8 also reveals that the dependence of walk-past probability on the stripe width is not necessarily monotonic. For instance, for kCR = 0. 004s−1 we see that the walk-past probability increases for increasing values of stripe width, reaches a maximum for d = 26μm and then decreases again. The reason for this decrease in walk-past behavior is the lack of confinement. For larger stripe widths cells sometimes begin the walk-past movement, but repolarize too quickly to complete it, resulting in a reversal instead of a walk-past. Comparable to Fig 6 cells can walk past each other several times for very small kCR, kCR = 0. 002s−1. We note that if a cell pair performs two walk-pasts and then separates, we also count this as a walk-past. For more details on the exact classification, please see S1 Text. We next ask the question whether our model can reproduce the quantitative statistics of cell-cell collision found in experiments of neural crest cells which revealed that ∼25% of cell pairs stick, ∼73% of cell pairs reverse, and ∼2% walk past [20]. As is evident in Fig 4, the coexistence region between sticking and reversal requires careful tuning. We manually varied all four relevant parameters and find that we can generate reversal statistics that quantitatively agree with the experiments. Specifically, we find that for kFR = 0. 04s−1, kCR = 0. 01s−1, Ocrit = 0. 15μm−2, and σ = 3. 85σ0 our simulations result in ∼69% reversal, ∼27% sticking, and ∼4% walk-past (see Table 1). Small deviations in these parameter values can cause large changes in outcomes; we will use this to screen for parameters which are critical to the experimental outcomes below. Our model can be used to predict and test the outcome of various biochemical interventions in the experiments. For example, the CIL response in neural crest cells can be blocked with the application of two treatments (DshDEP+ and Y-27632) [16], yielding a reduction in reversal and an increased number of sticking events—an increase to 44% sticking with DshDEP+ and 67% with Y-27632 [20]. In our model, a natural starting assumption for modeling the effect of these treatments in reducing CIL would be decreasing kCR. Can reducing kCR increase the number of sticking events, or do we need to invoke other effects of these treatments on, e. g. cell-cell adhesion? To test this hypothesis, we start with the parameter set that produce results consistent with the control experiments of [20] (Table 1) and determine the required change in model parameter values to increase the sticking probability to above 50%, as observed with their intervention experiments. The parameters that are more critical to controlling this outcome are likely candidates for the most relevant effects of treatments that increase sticking. We summarize these results in Table 2 by showing the required change in each parameter values (size and sign). The full results, including the percentage of each outcome as parameters are varied, are shown in S3 Fig. Consistent with our hypothesis, we can reproduce the effects of CIL-inhibiting treatments by altering the behavior of the inhibitor generation, either by reducing kCR or by increasing kFR. For lower values of kCR we see more sticking and walk-past outcomes and fewer reversals (S3 Fig). This can be understood by realizing that a weaker CR causes the repolarization to be slower, resulting in more cells that stick together. Surprisingly, an increased FR strength yields more sticking outcomes. Why does a stronger biochemical interaction lead to more sticking cases? We can observe two different mechanisms that lead to a pair of sticking cells. A stronger FR results in more cells that start the walk-past movement. The relatively strong CR prevents most of these from completing and the cells repolarize before the walk-past is carried out. These aborted walk-pasts often result in sticking cases. Additionally, more chains form at larger FR. These are unstable and their breakup often yields the sticking outcome. If kFR is increased further more chains form which become increasingly stable. Hence, the sticking rate decreases again (S3 Fig). Our simulations also show that other parameters can result in significantly more sticking events. Most of these parameter changes can be explained at an intuitive level. For example, decreasing the propulsion or increasing the contraction strength will make it harder for cells to separate (see Fig 4) and will thus increase the sticking probability. In addition, a stronger line-tension γ makes it harder for the cells to push the membrane outwards and will result in more sticking. More surprising, however, is the effect of increasing the activation rate for Rac, kb. Even though a higher activation rate yields a stronger recruitment of ρ from the cytosol to the membrane, the maximum level of ρ at the cell front does not rise. It is observable that the area of a single, freely migrating cell grows and the velocity drops slightly. The lower velocity even indicates a lower value ρ at the cell front and thus a weaker protrusion and more sticking. Taken together, our results suggests that our hypothesis, that CIL-inhibiting treatments strongly regulate the CIL response in [20] through the modification of kCR and kFR, may need to be extended; we address this further in the Discussion. Our model is easily extendable to more than two cells to address experimental studies of collective migration on stripes [19]. In Fig 9 we show outcomes of simulations of multiple cells on a single, long stripe for four parameter combinations that each yield a high percentage of one of the possible two-cell collision outcomes. Parameters that promote cell-cell reversal lead to cells that robustly repolarize upon contact with another cell and remain well distributed over the stripe (Fig 9i, S8 Movie). Parameters that create sticking cell pairs lead to cells that repolarize upon collision and mostly stick with the colliding cell. However, if a single cell collides with a pair of already sticking cells, we can observe that one cell of the pair may detach, resulting in non-motile cells that are mainly in pairs (Fig 9ii, S9 Movie). For larger adhesion strengths, we can also observe larger clusters. For walk-past parameters, cells begin to form trains but instead of having persistent motility as in [19], the trains break up (Fig 9iii, S10 Movie). Finally, for chain parameters, we observe the formation of persistent cell trains (Fig 9iv, S11 Movie). The trains form in a similar manner to the experiment of Ref. [19]: isolated cells collide and form moving chains. Subsequent collision with other isolated cells result in even larger chains.
  We have presented a model to describe the crawling and collision dynamics of cells on a micropatterned stripe. We applied the modeling approach from [28], including modeling contact-related signaling events known to occur in neural crest cells. This model can reproduce the four distinct phenotypes of cell-cell collisions observed in either neural crest cells [20] or healthy and metastatic epithelial cells [19,22]. Tuning only the four parameters of cell-cell adhesion σ, critical overlap Ocrit and the strengths of contact/front repolarization kCR and kFR is enough to have a very high rate of each outcome. We also have identified a set of parameters that yields all of the cases seen in neural crest cells, with rates in good agreement with those seen by [20]. We are also able to reproduce additional qualitative aspects of experimental collisions. This includes the peak in the velocity for the reversal case (S1 Fig), and the presence of “cell trains” (Fig 9), which Desai et al. reported in [19] for NRK-52E cells. The transition between these states in our model is summarized in Fig 10, which emphasizes that multiple parameters must be changed to control the different outcomes. If contact repolarization is strong (e. g. kCR ≥ 0. 03s−1, Fig 7) together with a weaker adhesion strength (σ ≤ 3. 84σ0, Fig 4), we can create a 100% reversal rate. Increasing the adhesion above this value yields a sharp transition to a high sticking rate. If instead of increasing adhesion from the 100% reversal rate parameters, we reduce the contact repolarization and tune the front repolarization to a medium value (kFR ∼ 0. 05s−1, Fig 7), we can robustly create chain events. Increasing the adhesion and introducing a nonzero critical overlap then yields a high walk-past rate—though this requires some degree of tuning (Figs 5 and 6). We also make predictions for how the rates of walk-past depend on the stripe width. We find that the walk-past rate is at its maximum of 99% for d = 26μm (Fig 8). For narrower channels it drops rapidly to zero, as the cells can no longer squeeze past each other. Surprisingly, wider stripes (d = 30μm) can also reduce the walk-past rate for some parameter combinations. This occurs when the absence of confinement allows the cells an additional route to separate prior to undergoing walk-past, leading to more reversals and fewer walk-past cases. Interestingly, [22] have also found that walk-past rates can be non-monotonic in stripe width for metastatic MDA-MB-231 mammary epithelial cells, though not for nontransformed MCF-10A epithelial cells. However, they see an opposite behavior to our simulations—with stripes of minimum width leading to a minimum in cell walk-past. We report a set of parameters that creates a mix of outcomes similar to those seen experimentally by [20], kFR = 0. 04s−1, kCR = 0. 01s−1, Ocrit = 0. 15μm−2 and σ = 3. 85σ0. With these parameters we observe 69% reversals, 27% sticking, 4% walk-past, and no chaining outcomes. These results, though, can be highly sensitive to small variations in parameters (Table 2, S3 Fig), similar to the sensitivity to adhesion that occurs near the sticking-reversal transition (Fig 4). However, the most sensitive parameters are not necessarily the expected ones. Variations in parameters that control the biochemical interactions (kCR and kFR) must be large (≥ 50%) to have an effect. Contrarily, the cells are very sensitive to changes in mechanical parameters that affect the sticking transition (α, σ, contraction strength β and line tension γ). Here, even small variations of 1% or less can switch collisions from majority reversal to majority sticking. This high sensitivity of the cells tells us that either a) cells are very carefully tuned towards the sticking transition, or b) cell-to-cell variations or variability of effects on contact play a major role in the coexistence of reversal and sticking—those cells that reverse may have lower cadherin expression levels or other biochemical differences. We argue that b) is more likely: typical scales of variation in, e. g. protein expression levels are 10% and higher [50,51], and eukaryotic cell responses to signals may also vary significantly from cell to cell [52,53]. Cell-to-cell variations would be expected to significantly increase the robustness of certain types of responses. For instance, even if the sticking-reversal transition is sharp, as in Fig 4, if there is a population of high-adhesion cells and a population of low-adhesion cells, sticking will occur when two high-adhesion cells stick; this will not depend on the precise value of σ for the high-adhesion cells. Our simulations use, as a central hypothesis, the idea that cell-cell contact leads to the generation of a chemical that inhibits a polarity protein, taken to be Rac, which we view as a reasonable first hypothesis to model CIL in neural crest [25]. What do we gain from modeling cell interactions to this level of biochemical detail, rather than using more generic mechanisms like flocking [5,54], “velocity alignment” [55–57] or others [9,37,58–61]? Some of our results are qualitatively surprising in comparison to those that would be expected from minimal models. With a minimal model (e. g. our simplified model of CIL in [9]), we would intuitively expect cell-cell adhesion to monotonically reduce walk-past. However, we find that increasing adhesion can promote walk-past events (Fig 5). This arises from the interaction of adhesion, which influences the cell-cell contact area, with the inhibitor generated at contact. We also note the importance of the different thresholds at which contact and front repolarization arise (the critical overlap). Neither of these effects would arise in a simplified model, and we argue that, for a full understanding of cell-cell interactions, modeling biochemical interactions like CIL and cell contact dynamics are both necessary. Our results show that the coupling of detailed cell shape dynamics with biochemistry can lead to unintuitive behaviors, like the nonmonotonic effects of adhesion. This illustrates the importance of complementing minimal models of cells as self-propelled particles [5,9] with increasingly detailed models for cell shape and mechanics. The influence of cell shape and adhesion on collective cell migration has also recently been studied in the context of the jamming transition in epithelial monolayers [62]; our results suggest that combining signaling with shape may lead to new effects. In partial contrast with our results, recent experiments studying the interaction of pairs of mammary epithelial cells colliding on micropatterns have found that increasing E-cadherin expression decreases the odds of walk-past [22]. These data suggest that either modeling the interaction of epithelial cells requires more than our straightforward CIL mechanism, or that only one portion of a non-monotonic dependence on E-cadherin expression has been explored for these systems. Many extensions of our modeling may be interesting. One natural choice would be to develop a more detailed description of the statistics of cell protrusions in neural crest, e. g. including stochastic protrusion and retraction [63] as has been modeled for Dictyostelium [64,65]. Inclusion of hydrodynamic flow within the cells, as we have previously studied for single cells [31,41], could potentially be relevant as well, especially if we were to extend our study from collisions in micropatterns to microchannels [66]. Another feature that we have not yet included in our modeling is the secretion of chemoattractant and subsequent chemotaxis to this signal, or “co-attraction, ” which is known to play a role in promoting the cohesion of neural crest cells [67]; this could be included in our model by techniques along the lines of recent approaches [10,60,61,68]. However, the graded response to a chemoattractant across a cell is itself an area of a great deal of interesting research [64,65,69–72], and including this feature could lead to significant additional complexity. What can we understand about experiments studying collisions of neural crest cells [20] in the context of our simulations? Scarpa et al. did not report any chaining events. We can eliminate chains by either making the front repolarization rate high, or including a contact repolarization effect. Even a weak contact repolarization can dominate a collision and prevent chains, especially when front repolarization is also present (Fig 7). The walk-past case shows that a nonzero CR is necessary for the cells to separate after passing each other. However, the CR cannot be too strong without suppressing walk-past (Fig 6). Taken together, this suggest that in modeling neural crest cell collisions, we should assume both FR and CR are present, with FR stronger than the CR. From these results, we reach three broad conclusions. First, based on our model, we would propose that treatments that increase sticking are likely to have effects either on cell-cell adhesion (e. g. cadherin expression levels) or on single-cell properties (e. g. changing cell speed or contractility). Secondly, because in order to create the experimental mix of outcomes from a simulation, we needed to carefully tune parameters, we suggest that the range of outcomes in [20] may arise more from cell-to-cell variation than it does from the stochastic motion of individual cells. Third, we predict that altering either cell-cell adhesion or micropattern size can lead to non-monotonic changes in outcome frequency.
summary: Many cells cooperate with their neighbors to move as a group. However, the mechanisms of these cell-cell interactions are not well understood. One experimental tool to analyze interactions is to allow cells to collide with one another, and see what happens. In order to better understand what features these experiments measure, we develop a computational model of cell-cell collisions, and identify the biochemical and mechanical parameters that lead to different outcomes of collisions. We can recreate all known types of collisions seen in experiments, including cells reversing on contact, sticking, or walking past each other. Our model suggests that what happens in a collision may depend strongly on the mechanical forces between the two cells.
keywords: cell physiology
  cell motility
  medicine and health sciences
  neuroscience
  biomechanics
  biological locomotion
  epithelial cells
  cell polarity
  developmental biology
  stem cells
  neuronal tuning
  developmental neuroscience
  animal cells
  biological tissue
  neural crest
  biochemistry
  biochemical simulations
  cellular neuroscience
  neural stem cells
  cell biology
  anatomy
  physiology
  cell migration
  epithelium
  biology and life sciences
  cellular types
  computational biology
title: Modeling Contact Inhibition of Locomotion of Colliding Cells Migrating on Micropatterned Substrates