Modeling Contact Inhibition of Locomotion of Colliding Cells Migrating on Micropatterned Substrates

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.

In Figure S2 we show the dependence of the sticking percentage on the simulation length for two different cases. We want to demonstrate that changes of the simulation length can a) change the collision behavior quantitatively near the reversal-sticking transition, but b) do in general not qualitatively change the outcome.
In the left panel we show the 50% contour line of the sticking/reversal transition for two different simulation lengths. Here, half of the of the collisions result in a pair of sticking cells and the other half ends in a reversal. With shorter simulation times the 50% contour line is shifted to a very slightly lower adhesion strength. The critical qualitative characteristics also stay the same: Increasing the propulsion strength α moves the sticking/reversal transition to larger adhesion strengths.
In the right panel of Figure S2 we show how changes in k CR affect the sticking percentage, and how this is affected by simulation length. This effect only occurs near the sticking-reversal transition. In agreement with our analysis of parameter variations (see Table 1 and Figure S3), increasing k CR can decrease the percentage of pairs that stick. Though the effect of longer simulation times is stronger here, the qualitative behavior remains the same for both simulation times. In longer simulations more cells detach and the sticking percentage decreases and for both simulation lengths, the sticking percentage increases for lower values of k CR .   These parameters are used throughout the paper; any deviation from them is explicitly noted.

Supplemental Methods
We solve the equations of motion on a large rectangular domain of size L×L (L = 100µm) with periodic boundary conditions. To save computation time, each cell has its own computation box with the size L small × L small (L small = 50µm) and periodic boundary conditions, in which we solve the phase field and reaction-diffusion equations. This box is recentered if the cell migrates close to the boundary. Then, the small box is shifted in the larger domain. The values of the newly created points are assigned by a periodic shift. Since the cell does not come close to the boundaries we assume that all fields (ρ,I,φ) are zero outside of the box. The computational boxes for the individual cells are always centered around the center of mass (COM) of a cell, where the COM is approximated to the nearest grid point of the larger grid. Both types of grid have the same properties, i.e. they are both regular, rectangular grids, with the same ∆x = ∆y. This has the consequence that each point of a computational boxes can be mapped clearly and without any interpolation/approximation to a grid point in the larger grid. When the equations of motion involve terms of two different cells (for example i =j φ (i) (r, t)φ (j) (r, t)), we transfer the required values from each smaller box to the larger grid from where they can be mapped back to every small box.
We solve the phase field equation with a semi-implicit spectral method to compute the first derivative for the advection and adhesion terms, the second derivative for the line tension and the fourth derivative for the bending term. All other terms and the reaction-diffusion equations are handled explicitly. Below, we introduce the substitutes I ≡ I(r, t) and φ ≡ φ(r, t). If those quantities are used at different locations or times, we will explicitly state the dependence.
In the multi-cell case, the phase field equation is supplemented by a term from the functional derivative of the cell-cell interaction energy, which is handled explicitly. In calculating this term, we use a 17-point finite difference stencil. This is to avoid lattice artifacts, like cell-cell interfaces favoring fixed directions with respect to the axis. Details of the derivation can be found in Section 4.
The Gaussian random variable Ξ(r, t) = t+∆t t dt ξ(r, t ) has zero mean and variance of Furthermore, we assume the noise to be uncorrelated between lattice sites and in time. To make sure that the strength of the noise is independent of the grid and time resolution we normalize it by √ ∆t/∆x. In our simulations we compute Ξ(r, ∆t) = Y η 2 ∆t/∆x −2 , with Y being a Gaussian distributed random number with zero mean and variance one.
Dividing by φ in Eq. (S7) is obviously a numerical issue when φ 1. Following the usual method for phase field equations Ref. [46] in the main text, we only apply this equation above a threshold value, φ ≥ 10 −4 . For φ < 10 −4 , we solve this equation, but multiplied by φ; this serves as a numerical convenience to reduce I and ρ where φ 1 and thus outside of the region where the reaction-diffusion equations would be solved in the sharp interface limit.
In particular, the diffusion term in Eq. (S7) ∇[φD I ∇I] is computed explicitly with a half-point stepping. In one dimension at time t and grid point i it has the form The equation for ρ(r, t) is solved identical, it only lacks the noise ξ(r, ∆t). We use OpenCL to solve the equations on GPUs. To compute the FFT needed for the spectral methods, we use clFFT in batch mode.
Initial conditions In our simulations, we confine the cells to stripes with width d = 26µm unless otherwise noted. The cells are initially separated by 40µm, but polarized toward one another. For each parameter combination we conducted 100 simulations, if not stated otherwise. To reduce the variability between runs with different parameters, we use the same seed to create the random numbers for the same run. Of course we never reuse a seed for runs within the same parameters. All default parameters are shown in S1 Table. Automatic detection of the collision outcome To automatically detect the outcome in our simulations we are defining the different events by the rules given in Table S2. It should be noted that the different possible results are checked in the order they are listed in Table S2. As soon as the conditions of one case apply the others are not checked anymore. Additionally, the simulations were aborted when the cells had a distance of 50µm between them (for reversal and walk-past cases) or they migrated 100µm in the same direction (chain). For the sticking case the simulation times are discussed in the relevant section. Table S2: Definitions of the different cases not included a cell turns around due to internal noise (distance between the cells is greater than 3 × r cell when a cell turns) rotation cells turn around and pass each other more than two times sticking y-velocity of both cells is below 0.005µm/s or the distance between the cells is smaller than 20m for the second half of the simulation and both cells turn around or the distance between the cells is smaller than 35m and the y-velocity of both cells is below 0.02µm/s walk-past no cell turns around or both cells turn around and pass each other two times over the course of the simulation ("double walk-past") and the distance between the cells is greater than 3 × r cell when the simulation ends. reversal both cells turn around and the distance between the cells is greater than 3 × r cell when the simulation ends and the cells have never passed each other. chain only one cell turns around and the distance between the cells is less than 3 × r cell when the simulation ends and the cells have never passed each other. A turn around is detected if the y-velocity of a cell changes its sign. For the latter three cases we only compare the sign of the start and the end velocity to detect a turn around. To compute the rates for each case we divide the count for a specific outcome by the total number of simulation, excluding the first entry of the table ("not included").

Derivation of higher order stencil for adhesion term
In this section we derive a higher order finite-differences stencil to compute the adhesion term for multiple cells Here, f = i =j |∇φ (j) | 2 , which is calculated using spectral methods. Eq. (S10) is solved with the finite-difference where h is the lattice spacing and is the centered finite-difference. For stability reasons Eq. (S10) is solved with a high order finite-differences scheme. One way to do this is to use the series expansion (F.B. Hildebrand's Numerical Analysis, Eq. 5.3.12): (S13) If we apply this to ∂ x (f ∂ x φ), we find when we omit all terms O(δ 5 ) or higher. The terms in Eq. (S14) in one dimension at time t and grid point i evaluate to: The points which are of the grid are computed with linear interpolation, yielding and Equation (S14) can be easily extended to two dimensions. We call the derived stencil G 1 . However, we can repeat the previous derivation on a different grid, along the two axes a 1 = (1, 1) and a 2 = (−1, 1). In this case, the lattice spacing is h √ 2. On this skewed grid, where δ 1 u i,j = u i+1/2,j+1/2 − u i−1/2,j−1/2 .
This resulting stencil is called G 2 . Then we can average the two derived representations for the derivative, to create a more isotropic stencil G. The weights taken from [1] are