Physical Confinement and Cell Proximity Increase Cell Migration Rates and Invasiveness: A Mathematical Model of Cancer Cell Invasion through Flexible Channels

Cancer cell migration between different body parts is the driving force behind cancer metastasis, which is the main cause of mortality of patients. Migration of cancer cells often proceeds by penetration through narrow cavities in locally stiff, yet flexible tissues. In our previous work, we developed a model for cell geometry evolution during invasion, which we extend here to investigate whether leader and follower (cancer) cells that only interact mechanically can benefit from sequential transmigration through narrow micro-channels and cavities. We consider two cases of cells sequentially migrating through a flexible channel: leader and follower cells being closely adjacent or distant. Using Wilcoxon's signed-rank test on the data collected from Monte Carlo simulations, we conclude that the modelled transmigration speed for the follower cell is significantly larger than for the leader cell when cells are distant, i.e. follower cells transmigrate after the leader has completed the crossing. Furthermore, it appears that there exists an optimum with respect to the width of the channel such that cell moves fastest. On the other hand, in the case of closely adjacent cells, effectively performing collective migration, the leader cell moves $12\%$ faster since the follower cell pushes it. This work shows that mechanical interactions between cells can increase the net transmigration speed of cancer cells, resulting in increased invasiveness. In other words, interaction between cancer cells can accelerate metastatic invasion.


Introduction
One of the main reasons of death in cancer patients is metastasis [28].When the cell invades through dense tissues, it often applies forces on its immediate surroundings to migrate through narrow channels and cavities [39].Cancer cells that are more invasive, tend to be very pliable and dynamic, both internally [11] and externally [7,14,33,45], and hence they are able to adjust their cytoskeleton and geometry.In addition to that, cancer cells have been shown to apply a significantly larger traction force on substrates, compared to benign cells [28]; Importantly, the forceful, invasive interactions of tumor cells with their environments have been used to provide rapid, clinically relevant cancer diagnosis and prognosis [20,30].
Cells are complex viscoelastic objects, and their shapes are determined by a local balance between retraction and protrusion of its boundaries [10].The driving force for cell deformation is typically triggered by external stimuli, for example, cell-substrate adhesion, biochemical signaling and forces exerted on the cell [32,33].Furthermore, cell adhesion and motility are affected by the mechanics of their environment [53] and, for example, cells appear to be more rounded on a softer substrate and more elongated on a stiffer substrate [24,28,54].In other words, cell shape is also an outcome of mechanical equilibrium as a result of external forces [33].However, the mechanisms that determine global cell morphology in relation to its function remain poorly understood [15,18,32].
Directed cell migration is an active process involving cell polarization that is often driven by chemotaxis or mechanotaxis [15].During single cell migration, the cell must first polarize for all modes of migration [24,25,39].Cell polarization into well-defined front (leading edge) and rear parts [8,15], is one of the most significant responses of animal cells to their environment [24].Indeed, there are many possible causes of cell polarization: external forces [49], Rho family small guanosine triphosphate (GTP)-binding proteins (Rho GTPases) [8,24], chemokines and cytokines [17,25] etc. Subsequently, the cell migrates in the direction of the polarity axis, that is, the longitudinal axis of the cell.
Mathematical modelling has been acknowledged as an important tool to help turn blurred concepts and ideas into testable, quantifiable and rigorous hypotheses and to reveal the correlations between various factors, which are otherwise difficult to determine in complex biological phenomena and microscopic experiments.In general, mathematical models in biology can be categorized as agent-based models (for small-scale materials) and continuum models (for large scale materials).Agent-based modelling is extremely beneficial to model cellular activities like cell division, differentiation and migration etc., since the model treats every cell as an individual and hence, it is capable of tracking cell positions and cell-substrate interactions.A further advantage of agent-based models is that most input parameters are expressed in terms of directly measurable quantities.Since in this manuscript, we will focus on cell geometry and single cell metastasis, we utilize agent-based modelling.Agent-based models have been widely used to investigate the evolution of cell geometry.Rens and Edelstein-Keshet [40] developed an algorithm approach to model the impact of cellular forces with the cellular Potts model, which is based on the assumption of minimizing the energy configuration; the work of Zhao et al. [56] mainly focuses on the intracellular environment where the finite-element method is applied to describe the total energy of the cell; the model in Cusseddu et al. [8] emphasizes the impact of Rho GTPases.
The current work is an extension of our previous work [36], where the Poisson effect, cell polarization and focal adhesion were not considered, yet still provides a basic model to depict the cell shape evolution under external stimuli and dynamic processes that occur during cell invasion.We use the finite-element method to approximate the solutions to all partial differential equations (PDEs) in the bulk domain.As energy is consumed during cell migration and deformation, energetic effects may also have an impact on the cell migration and invasion; those aspects will be studied in future studies.Furthermore, to the best of our knowledge there are currently no in-vitro results for cells that are migrating through deformable channels, which makes impedes experimental validation of the model.Our simulation study may therefore serve as a precursor for the development of in-vitro experiments of cells transmigrating through deformable channels.
The manuscript is structured as follows: The biological assumptions and the resulting mathematical model are presented in Section 2. In Section 3, numerical results from the simulations with one set of input parameter values are displayed.Since the model contains randomness and we aim at generalizing the conclusions in Section 3, Monte Carlo simulations are presented in Section 4. Finally, in Section 5, we summarize our findings and indicate the possible directions to further this research.

Mathematical Modelling
In this manuscript, we extend the phenomenological model in Peng et al. [36], where the Poisson effect was not taken into account in the mechanics of the cell.In that study, we were only modeling the interplay between the cell and its immediate extracellular environment.We are aware that to simulate the Poisson effect, the intracellular environment should be modelled as in Zhao et al. [57,56].However, in the work of Zhao et al. [56,57], the interaction between the cell and the substrate, in particular, the traction forces exerted by the cell to the substrate were not incorporated.On the other hand, from a computational point of view, the model will be relatively complicated if we model the intracellular and extracellular environment separately, as well as the impact of the traction forces on the substrate.Hence, we are interested in continuing the phenomenological model [36] to depict the Poisson's effect, in particular, when the cell needs to transmigrate through a small pore or channel.In this problem, we define the intracellular mechanics by means of a spring model with nodal masses on the cell boundary that are connected to their immediate neighbours and to the midpoint of the cell.

Biological Assumptions
We are aware that it is impossible and infeasible to include all the cellular activities of viable cells in the model.Hence we simplify reality regarding cell polarization and migration.Cells need to be polarized before migration, which results into leading edges and non-leading edges of cells.According to Zhao et al. [56], leading edges are characterized by having a positive inner product of the outward unit vector, n c , and velocity vector, v c , i.e. (n c , v c ) > 0. Note that (•, •) represents the inner product of two vectors.We consider the combination of cell migration and deformation by dividing the cell boundary into a set of nodal points, of which the coordinates are represented by x j .We consider a two-dimensional model, hence we consider the projection of cells onto a substrate layer.The boundary segments are line segments with two neighboring nodal points as vertices.Since the nodal points define the cell membrane, we use a weighted average to approximate the outward normal unit vector at a nodal point by the two line segments that contain the point as a vertex, hence, where e i,j stands for the line vector, that is, e i,j = x j − x i , which connects x i and x j , and n i,j c is the outward normal unit vector of the line segment e i,j , hence (n i,j c , e i,j ) = 0. We assume that the leading edges are more sensitive to the signalling molecules.
Since the cytoplasm of cells contains water and polymers, one often models cells as viscoelastic objects.Therewith they are deformable and compressible only up to a certain extent.Therefore, we assume that the difference between the area of the equilibrium status and the area at time t cannot exceed 10%, that is, where A(Ω C (t)) is the area of the cell at time t and A(Ω eq C ) represents the area of the cell when it is in its equilibrium shape.
We introduce two indicators, namely circularity and aspect ratio to evaluate the evolution of the cell shape quantitatively [15,28].For the two-dimensional case, the circularity is a measure of how circular a two dimensional object is.The circularity is defined by where A(Ω C (t)) represents the area of cell Ω C at time t and l(∂Ω C (t)) is the circumference of the cell boundary.The circularity value ranges 0 and 1, and the circularity value of a circle is 1, objects degenerating to lines and curves provide circularities tending to zero.Objects with circularity value between 0 and 1 may have an elliptic shape, or have an irregular shape that can contain many protrusions and lamellipodia.The aspect ratio is defined by: where d 1 (Ω C ) and d 2 (Ω C (t)) represent the length and width of Ω C at time t, respectively.The aspect ratio indicates the internal polarity axes of the cell, and illustrates how elongated a cell is; the aspect ratio of a circle is 1.Due to inhomogeneities in the extracellular domain, the migration of cells is subject to randomness.To model the random (uncertain) part of cellular displacement, a memorized random walk is embedded in the model [46], which is developed from the generalized Langevin model.Let dW (t) represent a vector Wiener process in which the contributions to the different coordinate directions are treated as independent events, and a normal distribution with zero mean and variance δt for all components, with δt denoting a time interval.From statistical principles, we have the following link between random walk described by the Wiener process for the particle position and the evolution of the probability density function of the position of the particle: , where D c represents the diffusivity of the cell phenotype in this extracellular matrix.This term models the random displacements caused by unpredictable localized microstructural inhomogeneities.Then, the spontaneous velocity of the cell v s (t) and the memory of the velocity V s (t) are described by [43] where β s (v s (t)) is the velocity decay rate with the velocity dependency, α s is the memory rate and γ s is the memory decay rate.Since Takagi et al. [46]  2 , where S(Ω C (t)) is the aspect ratio of the cell defined in Equation (2.3).The assumption is based on an experimental observation that elongated cells migrate faster, and therefore, we use inverse proportionality of the velocity decay rate on the aspect ratio of the cell [18].

Cellular Traction Forces
It has been documented that cancer cells exert forces on their direct environment [28], which results in the deformation of the extracellular environment, which , in turn, has an impact on the cell as well.Furthermore, cells, in particular cancer cells, are known to be able to permanently change their immediate environment.These permanent changes could give rise to shrinkage or growth processes, and to mechanical displacements as a result of forces that are exerted by the cells themselves.Similar to Peng et al. [36], we therefore use a morphoelastic approach to model passive convection of the cell, and (permanent) displacements of the extracellular matrix.Morphoelasticity is commonly used to combine the modeling of mechanical stresses and processes like growth and shrinkage [3,13,19,42].The morphoelastic model reads as [19]: where ρ is the density of the extracellular matrix, is the local (Eulerian) effective strain tensor to be solved, L = ∇v is the deformation rate gradient, and α is a non-negative constant that accounts for the amount of plastic deformation proportional to local (Eulerian) effective strain.This parameter α could have an electro-chemical origin, however, in the current study, we treat it as a constant for the sake of simplicity.Note that if α = 0, then as soon as the force f = 0, the tissue will gradually recover to its original shape and volume.Furthermore, let y be a scalar function, then the material derivative is given by Dy Dt = ∂y ∂t + v∇ • y where v is the migration velocity of any point within the domain of computation.In order to have a fixed boundary of computational domain Ω, we use a homogeneous Dirichlet boundary condition for the velocity.From a mechanical point of view, we treat the computational domain Ω as a linear and isotropic medium.Further, as a result of the presence of liquid phases in the tissue, the mechanical balance is also subject to viscous, that is friction, effects.Therefore, we use Kelvin-Voigt's viscoelastic dashpot model, of which the stress tensor reads as where E s is the stiffness of the substrate, ν s is the Possion ratio of the substrate, is the strain tensor, µ visco is the weight of the viscosity in viscoelasticity, µ 1 and µ 2 are the dimensionless shear and bulk viscosity, respectively.The morphoelasticity model provides a set of nonlinear partial differential equations for the displacement velocity and the strain tensor.Since the initial displacement is zero everywhere, the displacement of the domain can be approximated by integrating the velocity over time: u(x(t), t) We assume that the traction forces applied by the cell are modelled by Dirac Delta distributions.Let Ω ⊂ R d be an open region, then the Dirac Delta distribution is given by We currently model (cancer) cells that are migrating through narrow channels and we incorporate the pushing forces that are exerted by the cancer cells on their immediate environment.These forces are directed in the outward normal direction of the cell boundary and they are only applied on the cell boundary segments that are in mechanical contact with the extracellular obstacles.These obstacles can be other cells or the wall of the narrow channel.Following the same form of traction forces as in our previous work [36], we denote N i m = {j 1 , . . ., j m } as the set of line segments of the cell membrane that are in mechanical contact with the obstacles, then the traction forces for cell i are expressed as where n(x(t)) is the unit outward pointing normal vector at x(t), s i j (t) is the midpoint of line segment j of cell i, and ∆Γ i,j is the length of line segment j.Note that the cell index is marked as a superscript in the notations.Most of the time, for the simplification, we neglect it unless we explicitly describe the case for multiple cells.Here, Q(d(x), t) is defined analogously as in [36]: Here, E * is the total equivalent Young's modulus defined by [38] 1 where E c and E s are the Young's modulus of the cell and the ECM, respectively, and ν c and ν s are the Poisson ratio of the cell and the ECM, respectively.Furthermore, d(x) is the invagination depth of the cell due to the compression by the obstacle and l m is the total curve length of the portion of the membrane of the cell that is in mechanical contact with the obstacles.Note that the division by l m is necessary for normalisation since Q represents the force per unit of peripheral length of the cell.The total force follows from integration of Q over the boundary portion that is in physical contact with the obstacle.

Cell Cytoskeleton
In Peng et al. [36], the cell boundary is modeled by a set of points on the cell boundary where all the points are connected to the cell center by springs.In the current study, we extend this principle by connecting adjacent cell boundary points by springs.Using this formalism, it becomes possible to accommodate for surface tension.Here, we assume that the elasticity of the springs that connect the node to the center and to its neighboring nodes can have different values.The cell cytoskeleton is schematically represented in Figure 2.1.Similar to Peng et al. [36], Chen et al. [6], Vermolen and Gefen [50], without any extracellular stimulus (i.e.here we only present the elasticity part of the cell deformation; the spontaneous displacement of the cell that is described in Equation (2.4) will be incorporated later into the equation of motion), the position of nodal point x j over the cell boundary is determined by ) where E c and E m represent the stiffness of the spring connecting x j to the centre point and neighbouring points respectively; xj = xj (t) − x c (t), xj,j±1 = xj (t) − xj±1 (t) are the vectors connecting the nodal point j to the centre point and the neighbouring points when the cell is in equilibrium shape, and xj (t) represents the equilibrium position of the nodal point j relative to the centre x c (t).
To maintain the right orientation of the cell, we consider the rotation matrix defined by such that the angle of orientation, φ, can be computed from φ = arg min (2.10) The orientation of the cell is represented by the angle of the vector connecting the 'front and tail' of the cell.The overall displacement of the nodes of the cell boundary are determined  by translation and rotation.This matrix B(φ) monitors the angle of rotation of the cell with respect to the cell position (and hence boundary nodes) at the previous timestep.This angle of rotation is important for the determination of the equilibrium points of the cell boundary nodes.The equilibrium points reflect the position and shape to which the cell boundary nodes will converge to if the cell is not subject to any external cue for migration, that is, the cell no longer migrates.Without any external cues, the cell will always return to its initial orientation.Subsequently, Equation (2.8) is rewritten as Note that here we have not yet incorporated the memorized random walk of the cell that is described in Equation (2.4).The memorized random walk will be added to the equation of motion in the next section.

Concentration of Generic Signal
For cancer cells, examples of such cues for migration are oxygen and nutrients (such as sugarlike chemicals) [21,41].We consider a generic chemical, which originates from an artificial point source.Since we only want to scrutinize the model's ability to predict collective cell migration that is impacted by external mechanical factors, we keep this cue the same for all cells (both leader cells and follower cells).Therefore, contrary to our previous work, we use a steady-state diffusion equation (i.e.∂c ∂t = 0): where c s (x) is the concentration of the signalling molecule, D is the diffusion coefficient which has been taken constant in the current study, k is the secretion rate of the signal source, x s is the position of the source.Furthermore, Ω is an open bounded simply connected subset in R 2 with boundary ∂Ω.The steady-state problem is closed by the following Robin condition which deals with a balance between the diffusive flux across the boundary and the flux between the boundary and the region far away from the domain of computation.The symbol κ s , which is non-negative, represents the mass transfer coefficient.Note that as κ s → 0 then the Robin condition tends to a homogeneous Neumann condition, which represents no flux (hence isolation), which in the current steady-state case, yields an ill-posed problem in terms of existence (and uniqueness).Whereas κ s → ∞ represents the case that c s → 0 on the boundary, which, physically, is reminiscent to having an infinite mass flow rate at the boundary into the surroundings.The Robin condition, also referred to as a mixing boundary condition, is able to deal with both these two limits and all cases between these limits.From a mathematical perspective, the solution c s to Equation (2.12) does not live in the right finite element space (H 1 (Ω)) [35].Therefore one observes issues regarding accuracy and convergence of the finite element approximation in regions close to the point where the Dirac distribution is acting, that is, x s .However, away from x s , the accuracy of the finite element solution is as good as the classical accuracy in the L 2 -norm of the solution [23].We note that in our application, we are only using the solution away from the point of action x s .

Focal Adhesions
Focal adhesions play an important role in cell migration.The movement of a cell is driven by continuous remodeling of the cytoskeleton and is mediated by the lamellipodia.In general, a cell moves in three steps: (1) Protrusion: the rear part of the cell is attached on the substrate while the membrane of the front (leading) part is extended; (2) Adhesion: new adhesion is generated at the newly extended membrane of the front part; (3) De-adhesion and retraction: the rear part of the cell is detached from the substrate and due to the cytoskeleton, the rear part retracts to maintain the cell geometry.Once cell contraction has completed, the cycle of movement turns back to Step (1).A more detailed description of cell migration, or crawling, over the substrate can be found in Ananthakrishnan and Ehrlicher [1], Bershadsky and Kozlov [4], Mitchison and Cramer [31].
For the implementation, we simplified the aforementioned three steps into two steps, that is, the leading part firstly detaches from the substrate and moves forward while the non-leading part is attached, then the leading part is attached to the substrate and the non-leading part detaches and migrate due to cell contraction.For each nodal point on the cell membrane, we consider a stochastic Markov Chain model, such that the impact of the previous time step (i.e. the previous status) is also accounted for.We assume that the probability that the nodal point is detached from the substrate is related to the local shear force and the local gradient of the concentration of the signalling molecules.In other words, if there is larger local shear force and/or higher gradient of concentration of the signalling molecules, then it is more likely that the nodal point detaches from the substrate.Denoting 1 for the attached status and 0 for the detached status of the nodal point, respectively, then the probability of the occurrence of status s when the previous status is s is given by where X n is the next status and X n−1 is the present status; λ 1 = λ 1 (σ 12 (x, t), c s (x, t)) = |σ 12 (x)| + 100 ∇c s (x) + 5, and λ 2 = λ 2 (σ 12 (x)) = |σ 12 (x)| + 5 respectively, and σ 12 is given by Equation (2.6), where the matrix is symmetric and the non-diagonal element is related to the shear stress.Equation (2.13) can be rephrased as a Markov Chain matrix: The initial status is randomly assigned to the nodal points on the cell membrane.To determine the status of such a point, a random value from a uniform distribution between 0 and 1 is drawn and tested with the transition probabilities.

Cell Geometry Evolution
The cell shape is determined by the positions of the nodal points on the cell boundary.The movement of the nodal points is determined by various processes, which are the retraction towards the equilibrium shape, fluid flow, spontaneous displacement (random walk), chemotaxis and possible external mechanical forces.In Section 2.5, the focal adhesion part of cell migration was formulated in terms of a stochastic model where the front and rear parts of the cell detach and attach onto the substrate by chance.Nodal points on the cell boundary that adhere onto the solid substrate are only subject to movement that is determined by the local displacement of the substrate.On the contrary, if a nodal point on the cell boundary is detached, then its movement is determined by all the other aforementioned processes, and hence where γ is a small positive constant to prevent the denominator being zero and v is the velocity of the substrate determined by Equation (2.5).Here, we take E c and E m as constant.Furthermore, µ s (S(Ω C ), x j ) is the weight of chemotaxis/mechanotaxis that is related to the aspect ratio of the cell and whether the nodal point in on the leading part.It has been found that the receptors are distributed differently of the polarized cell, that is, the leading edges and non-leading edges have a different response to the chemotaxis [9,51].Hence, we assume µ s (S(Ω C ), x) is defined by where v c is the cell velocity at time t that is defined by the centre position of the cell: and here, d 1 (Ω C ) is the length of the cell and the initial velocity is zero.The first term in Equation (2.16) indicates that the points on the leading edge react more sensitively to the signalling molecules.Since it is assumed that the cell migrates in the direction of its longitudeaxis, ( in the second term computes the projection of the vector (x j −x c ) onto the cell velocity vector.We divide this by the diameter of the cell for the sake of normalization.
The second term in Equation (2.16) is meant to emphasize that the more to the front the point is located, the more sensitive it responses to the signalling molecules.Furthermore, when a cell is elongated, it is more sensitive to the signalling molecules.In summary, the cell geometry is obtained by connecting these nodal points on the cell membrane in the predefined (anti-)clockwise order.Cells that penetrate through narrow channels and narrow blood vessels are subject to friction with the channel or vessel walls.This friction causes a reduction of the migration speed.Denote ∂Ω ob as the boundary of the obstacle (for a microtube, this represents the walls).Then if a nodal point collides with an obstacle, then both the repelling force, exerted by the cell, and the force that results by the compression of the springs are taken into account.Note that the repelling force does not contribute to the friction if the springs are stretched or in their equilibrium length.For every nodal point on the cell membrane, there are three springs connected with the centre and the two adjacent nodal points.Hence, the force that is caused by these springs are given by where B( φ) xj is the equilibrium position of x j if there is only orientation but no other mechanics, represent the unit vector of the springs, respectively when they are compressed.In this modeling framework, we consider the force that is exerted by the cell in order to retain the equilibrium shape configuration.One may compare this spring force to the relaxation of skin after a force has been applied.The model also considers a second force, which is the force that the cell actively exerts in order to open the channel further so that the cell can continue its migration.The current formalism superposes these two forces.Then, the friction force orthogonal to the obstacle that has an impact on the migration of the cell is determined by where f s (x; t) is the force from the springs defined in Equation (2.17) and f (x; t) is the repelling forces actively exerted by the cell that is defined in Equations (2.5) and (2.7).It is assumed that the magnitude of the friction is proportional to the repelling force.Furthermore, n ob is the unit orthogonal vector pointing outward the cell centre, hence the vector points into the substrate.In our model, we simply subtract a friction part of the velocity in the tangential direction of the obstacle.Hence, the displacement of the nodal point which collides the wall of the microtube is given by where µ f is the cell friction coefficient, f (x j (t)) is the repelling force exerted by the cell and τ ob (x) is the tangential direction (also unit vector) of the obstacle boundary ∂Ω ob .Furthermore, the cell is not allowed to penetrate into the wall itself and hence velocity components normal to the wall are subtracted.In other words, the displacement of the nodal point that collides the obstacle is rephrased as (2.20)

Numerical Results
In this section, we investigate how cells influence each other's migration when they migrate through a narrow pore or channel.In the in-vitro experiments which contain microtubes, the width of the microtube is often set a bit smaller than the cell diameter, which experimentally models the transmigration of cells through narrow channels in real tissues.However, as we mentioned earlier, in the in-vitro experiments, the microtubes are undeformable [26,27,55], which is contrary to real tissues, where the immediate environment of cells, consists of viscoelastic material.Therefore, inspired by the microtube experiments, we adjust the experimental setting such that the channel can deform due to the forces exerted by cells in our simulations.We consider two different settings, which are distinguished by how many cells are migrating through the channel at the same time.In Case (1) (see Figure 3.1(a)), there is only one cell in the computational domain penetrating through the channel.In other words, the next cell is only allowed to enter the channel when the previous one has completely exited.We are aware that this setting does not probably reflect reality since cancer cells mostly migrate collectively.However, this setting is helpful to determine and to quantify whether the follower cell benefits from the leader cell, which possesses the ability to permanently widen the channel.The permanent nature of the deformation has been incorporated by the use of the morphoelastic component in the solid mechanical part of the model.In Case (2) (see Figure 3.3(a)), we consider two cells that simultaneously migrate through the narrow channel.This case reflects the mechanism of collective migration somewhat better.In this case, we consider the forces that the cell exert in order to repel and push each other.Note that here we assume that there is only mechanical contact between the two cells, as likely occurs in cadherin deficient cells, such as the commonly used metastatic, breast cancer cell line MDA-MB-231 [55].
In both cases, we have similar experimental settings and in this manuscript, all the simulations are carried out in two dimensions.In the computational domain Ω = (−100, 100) × (50, 50), a sinusoidal deformable channel that starts at X l = −64 and ends at X r = 60, is given by the following top and bottom part of the channel walls y = ±(b + A sin(wx)), b, A, w > 0. (3.1) Note that 2(b + A) is set to be smaller than the cell diameter to ensure that the cell is confined in the channel.The advantages of a sinusoidal channel are: (1) it is a well-defined geometry, where the channel width ranged between 2(b − A) and 2(b + A); (2) it resembles the walls of a channel in a tissue, where the wall is formed by cells [16].We assume that the cells are located initially at the left side of the channel and at the other side of the channel, there is a point source of signalling molecule.The parameter values that are used in the simulations in this section are shown in Table 3.1.Note that we also define the initial positions of the most left and right point of the channel, since the channel will deform due to the forces exerted by the cells.As a result, the position of the channel is altering in the same time: we define {X channel } as the set of initial positions of the channel, then the new positions of the channel {x channel } can be tracked by where u(X, t) is the displacement of the computational domain, post processed by solving the morphoelasticity model in Equation (2.5), and X is the initial position.We say that the cell enters the channel when there is more than one nodal point of the cell membrane in the channel (the moment is denoted by t i in for cell i) and the cell leaves the channel completely when there is no nodal point of cell membrane in the channel (the moment is denoted by t i out for cell i).
To rephrase it, the entering time and exiting time are defined by t i in = arg min t 0 {x i j ∈ ∂Ω C i : ∃x i j enters the channel}; for the cell i and ∂Ω C i , and we remind the reader that x i j denotes the cell membrane and j-th nodal point on the cell membrane of the cell i.Subsequently, the transmigration time of cell i is given by T i = t i out − t i in , and the (average) speed of the cell migration reads as vi = , where x i c (t) represents the central position of cell i at time t.

Case (1): Single Cell Migration
The initial setting of the simulation is shown in Figure 3.1(a).We consider a deformable sinusoidal micro-channel that separates the left and right parts of the domain.The cells can migrate through the channel.To ensure most settings are the same for every cell, the initial positions of cells are the same, that is, all the cells start moving from location where the centre position is (−75, 0).Several snapshots at consecutive times of the simulation are shown in Figure 3.1.Furthermore, to have a better visualization of the simulation, a short animation is attached (see Supplementary Material Video 1). Figure 3.2(a) shows the cells' speed when they are migrating through the channel.Note that here t = 0 min in the figure is the moment when the cell enters the channel.In general, it can be seen that the curve of the first cell has the largest amplitude of the oscillation, while the migration speed profiles of the other cells exhibit smaller amplitudes.Furthermore, it is clear that the first cell moves the slowest and takes the longest time in the channel compared to the follower cells.However, regarding the cell speed, it is hard to see the difference between the second and the third cell, except that the third cell spends less time in the channel (since the black curve disappears the earliest in the figure).Furthermore, in Figure 3.2(c) and (e), cells are in an elongated shape when they migrate, as the circularity and the aspect ratio never recover to 1, which is the case for the equilibrium (and initial) shape.However, it seems that the difference in the circularity is not that significant, while regarding the aspect ratio, the follower cells are more elongated, which accelerated their velocity according to our model assumptions.
For the sake of quantification, we show the transmigration time and average speed of three cells in Table 3.2, and the change of the width of the channel is shown in Table 3.3.Since cells are migrating to the right-hand side, the channel is stretched to the right as well, as it is seen in Figure 3.1.Furthermore, we also noticed that the sinusoidal channel becomes more and more flat after the transmigration of several cells, which is mainly caused by the fact that the repulsion force is proportional to the penetration depth of the cell.From Table 3.2, it can be concluded that, in general, the follower cells move faster and move over a shorter distance, compared to the leader cell.Furthermore, every cell contributes to making the channel wider and after all three cells penetrating through the channel, even the minimal width of the channel is comparable with the cell size; see Table 3.3.As a result, the follower cells benefit from the leader cell that widens the channel such that the follower cells deform less and there is also less friction from the wall of the channel.However, we observed that the speed of Cell 3 is 1.53%After the second cell exits and before the third cell enters the channel 6.81 14.08 After the third cell exits the channel 8. 47 14.26 smaller than the speed of Cell 2, which may imply an optimum channel width to maximize the cell speed.Such an optimum could easily exist since the wall is also expanded by the cell in order to migrate.Furthermore, there is a significant difference of the cell aspect ratio between Cell 2 and Cell 3, which results from the fact that before Cell 3 enters the channel, the channel has been expanded significantly such that Cell 3 is not required to change shape as much; hence, cell polarization and the cell aspect ratio are less compared with Cell 2. Consequently, Cell 3 transmigrates more slowly than Cell 2.

Case (2): Cell Doublet Migration
We consider the second case, where the follower cell follows the leader cell immediately entering the channel.That is, there are more than one cell in the computational domain at the same time, which is so-called cell doublets in in-vitro experiments [55].The setting of cell doublets is a simplified type of collective migration.Further, this experimental setting (see Figure 3.3(a)) results into the possibility that cells collide against each other.We assume here that cells can also exert force on the other cell as long as they are compressed.For a better visualization, we refer to the Supplementary Material Video 2, attached to this manuscript.Figure 3.3 shows several snapshots at consecutive times of the simulation, and Figure 3.2(b) illustrates the speed of the two cells inside the channel.The results of the simulations are summarized in Tables 3.4 and 3.5.The difference between the average speeds of the follower cell and the leader cell is very small, and due to the randomness involved in the model, these differences are not significant from a statistical point of view.Another, more important, effect that we see is that the average speeds of both the follower and leader cells are about the same as the average of the follower cell in Case (1), which has increased compared to the average speed of the original leader cell.This fact clearly shows that collective cell migration is beneficial for the cells.This also demonstrates that cancer metastasis benefits from collective cell migration.Similarly to the conclusion we drew in the previous section, Figure 3.2(d) and (f) display the circularity and the aspect ratio of both cells.Different from what we saw in Case (1), there is a significant difference regarding the circularity and the aspect ratio of the follower cell and the leader cell: the follower cell alters its equilibrium shape much less than the leader cell, as both indices of the follower cell are closer to 1 than for the leader cell.Although, we do not  explicitly take energy into account in the current model, we expect from these simulations that the follower cell consumes less energy than the leader cell.

Comparison between the Two Cases
From Table 3.2 and 3.4, we observe that there is a significant increase in the average speed of Cell 1 in Case (2).To have a clearer view, we plot the average of the time series of the velocity of Cell 1 and Cell 2 that are collected from 50 simulations with the same input parameters in Table 3.1; see Figure 3.4 for the results.It can be concluded that for the selected parameter set, the leader cell (Cell 1) is moving faster in Case (2) than Case (1), as the follower cell in Case (2) exerts the repelling force, while for the follower cell (Cell 2), there is hardly any difference with respect to the speed, however, in Case (2), Cell 2 needs more time to transmigrate through the channel completely since the leader cell is blocking the channel.However, the most important conclusion so far is that in Case (2), both the leading and follower cells transmigrate faster than the leader cell in Case (1).This shows that collective cell movement can increase metastatic migration rates, where we illustrated this behavior for two cell-interaction modes: (1) the leader cell 'paves' the way through for the follower cell (Case (1)), and (2) the cells interact by pushing each other thereby producing a net larger joint force onto the channel walls (Case (2)).
Regarding the invasiveness, we focus on the channel width after the leader cell (Cell 1 in Case ( 1)) or the cell doublet (both Cell 1 and Cell 2 in Case (2)) have exited the channel completely.The results can be found in Table 3.3 and Table 3.5.Both minimal and maximal width of the channel in Case (2) (5.82 and 13.72 µm, respectively) are wider than in Case (1) (5.75 and 13.09 µm, respectively), since in Case (2), both cells exert forces on the channel walls at the same time and within a certain distance.This implies that cell doublets enhance the invasiveness and increase its speed as compared with single cell, which is in line with experimental observations by Merkher and Weihs [29], Tulchinsky and Weihs [48].
The aforementioned observations hold for the current set of input values.Next, we will vary the input parameters according to statistical distributions and see whether the above conclusions are general.

Monte Carlo Simulations
In Section 3, we considered series of 50 simulations with the same input parameters.To investigate whether the conclusions (the potential existence of an optimal width of the channel that increases migration speed under for interacting cells) from Section 3 hold in general, Monte  Carlo simulations are performed where the input parameters are subject to the statistical distributions that are given in Table 4.1 and the other constant parameters are the same as in Table 3.1.In this section, we will show box-plots of the average cell speeds, and conduct the Wilcoxon test to statistically investigate whether the follower cell and leader cells benefit from each other for both Case (1) and Case (2) for the parametric ranges that we use.Since the output parameters do not necessarily follow normal distributions, we compute the Spearman correlation coefficients [52] to determine the correlation between input-and output parameters.
In this section, we compute the Spearman correlation [52] as it does not request the data to obey the normal distribution.Furthermore, we did the Shapiro test [44] , which evidenced that most output data does not follow a normal distribution.

Case (1): Single Cell Migration
In total, 1019 samples are collected with all the input and output variables mentioned in Table 4.1.As the channel is deforming continuously, the average speed of the cell during the penetration is more appropriate to investigate.Figure 4.1 illustrates the average speed of all the cells and the scatter plot between the first and second cells.It is hard to conclude from the box-plot (Figure 4.1) whether there is a significant difference between cells.Therefore, Wilcoxon's test is used and the results are shown in Table ??.As the leader cell, Cell 1 does help the follower cells (Cell 2 and Cell 3) move faster since the p-value from the Wilcoxon's test is almost 0.However, Cell 3 moves more slowly than Cell 2, which can be explained by the energy consumption and the confinement: from Figure 4.2, it can be seen that after Cell 2 has transmigrated through the channel, the channel is expanded more, which results in less confinement for the cell and then the aspect ratio of the cell is smaller; we speculate from the perspective of the energy consumption, that a cell may prefer to stay mostly in its original equilibrium rather than migrate faster, to some extent.Furthermore, we see a relatively strong correlation between the speeds of each cell; see Figure 4.1.As a result, to investigate how the properties of the substrate and the flexible channel influence the cell displacement, we will mainly focus on Cell 1.In Figure 4.2, we investigate how the width of the flexible channel before the cell enters, affects the average velocity of the cell respectively.As it is expected, once every cell leaves the Outputs with respect to cell i, i = 1, 2, 3 †

T i
The transmigration time of cell i through the channel vi The average speed of cell i penetrating through the channel cell dis(i) The length of the channel that cell i has migrated, as the channel is deformed due to the cellular traction forces cell max as(i) The maximal aspect ratio of cell i during the penetration cell min c(i) The minimal circularity of cell i during the penetration Outputs with respect to width of the channel, where i = 1, 2, 3 † min width 0 The minimal width of the initial channel min width(i) The minimal width of the channel after cell i completely exits the channel, that is, before cell i + 1 enters the channel max width 0 The maximal width of the initial channel max width(i) The maximal width of the channel after cell i completely exits the channel, that is, before cell i + 1 enters the channel channel, it contributes to widening the channel in general.In Figure ???? -??, the correlation for Cell 1 is opposite for the other two cells.Before Cell 1 enters the channel, the minimal width of the channel can be much smaller than the cell size.According to our model assumptions that the cell area can only alter at most 10%, the cell has to become more elongated, which results into a larger aspect ratio.Therefore, the cell should move faster.However, in the meantime, a narrow channel also increases the friction from the channel on the cell, and it requires more energy from the cell to deform more, which actually slows down the cell -that is why Figure ???? shows a parabolic behavior and it shows an opposite tendency compared to Figure ???? and ??.Once the leader cell (Cell 1) expands the channel to a width that is comparable to the cell size (see Figure 4.2), the follower cells (Cell 2 and Cell 3) do not need to deform as much as before, thus, they stay more akin the equilibrium shape and less elongated since cells are not confined.Subsequently, a decreasing tendency appears in Figure ???? and ??.In other words, there should exist an optimum of the width of the channel, such that the cell moves fastest with a combination of a relatively large aspect ratio and a relatively small deformation-this is verified by the fact that Cell 2 moves significantly faster than Cell 3 in the Wilcoxon's test (p-value = 8.28 × 10 −47 ) in Table ??.

Case (2): Cell Doublet Migration
For Case (2), since the leader cell (Cell 1) and the follower cell (Cell 2) are entering the channel at about the same time, the data regarding the channel width after the leader cell exits is less interesting.There are 1044 samples collected from the Monte Carlo simulations of Case (2) and similarly, we start with analyzing the average speed of the cells.From Figure 4.3, there is hardly any difference between the average speed of the cells, which is similar to the preliminary results in the previous section.The Wilcoxon's test verifies that two cells are moving at more or less the same velocity with p-value = 0.9999 ; see Table ??.In the meantime, there exists a strong positive correlation between the average velocity of two cells (corr = 0.86).Two cells are moving together through the channel and they will definitely collide, since from Case (1), we have concluded that the follower cell benefits from the expanded channel that is caused by the leader cell.Because of the collisions, even though the follower cell is able to move faster, the leader cell blocks the way, hence, the average velocity is more or less the same between two cells.

Cell Doublets enhance Cancer Cell Invasiveness and Migration
As mentioned in Section 3, we are interested in knowing whether the model predicts that the leader cell can also benefit from the the follower cell, and vice versa, or whether the follower cells is slowed down by the leader cell.To this extent, we analyze the model setting in Case (2).To achieve this, we collect the data from the Monte Carlo simulation such that every dataset is paired, and each paired dataset has exactly the same input parameters for the Monte Carlo simulations for both cases.Subsequently, the results collected from the both cases are comparable.The distributions of the input parameters are shown in Table 4.1, whereas the main interesting outputs are the average velocities of Cell 1 and Cell 2 in the both cases.In this section, 1000 paired samples are collected.
In Figure 4.4(a), we show the box-plot of the average speeds of the two cells in both cases.It is not straightforward to draw any conclusion for Cell 1, while there is a dominant difference in Cell 2 for different cases, that is, the average velocity of Cell 2 in Case (2) is less than the average speed of Cell 2 in Case (1).We conduct the Wilcoxon's test between the paired samples and the results are shown in Table ??.As we have expected and similar to the conclusion in Section 3, Cell 1 transmigrates through the flexible channel significantly faster in Case (2) than in Case (1) (with p-value 1.52 × 10 −17 ) due to the cellular forces exerted by the follower cell in Case (2).As a result, Cell 2 cannot pass through Cell 1 in Case (2) even though Cell 2 is migrating faster thanks to the expanded channel (otherwise, Cell 2 will no collide with Cell 1 when there is distance between them initially), and Cell 1 also exerts the repelling force on Cell 2, which also slows down the transmigration of Cell 2 significantly (p-value = 4.81 × 10 −154 ).In summary, under the setting of collective migration, the leader cell also benefits from the follower cell such that the leader cell transmigrates faster and easier through the narrow channel.Hence, collective cell migration accelerates cell migration and may be responsible for increased metastatic rate.   1) and the cell doublet exits the microchannel completely in Case (2).As each dataset is paired, we conduct the Wilcoxon's test in Table ?? and we conclude that both the minimal and maximal width in Case (2) are significantly wider than in Case (1) when all the input parameter values are the same.

Conclusions and Discussions
In this manuscript, we attempt to answer the following research questions: 1. Can a follower cell benefit in transmigration speed through a narrow flexible from a leading cell that may have widened the channel?
2. Can cell doublets (i.e. a pair of cells that are migrating close to each other) benefit in transmigration speed by the mechanical repelling forces that they exert on each other?
We designed these questions as two categories of simulations, so-called Case (1) and Case (2), respectively.
Case (1) answered the question that indeed, the follower cells move significantly faster than the leader cell.However, there exists an optimum regarding the confinement on the cell, such that the cell moves fastest through the channel or the pore.If the confinement is too small, there is barely any polarization of the cell, hence, the cell would prioritize maintaining the equilibrium shape instead of other cellular activities; on the other hand, if the confinement is too large, besides the friction that burdens the migration, most energy is spent on the deformation of the cell, then the cell cannot migrate optimally; furthermore, the cell is under a higher risk of death due to the large mechanical stress [12].This behavior is verified in the experiment conducted by Mak et al. [27], where they show that cells take much longer time to transmigrate through the constrictions and only part of the cell population survive.Another possible explanation is from the perspective of minimizing the energy configuration: to some extent, a smaller deformation costs less energy to deform and to recover to its equilibrium shape, hence, most energy can  be spent in the migration; however, if there is hardly any confinement on the cell, the cell prefers to stay in its (nearly) equilibrium shape rather than consume any energy on any cellular activities such as migration or deformation.We remark that this energetic view is speculative since energy is not considered in the current formalism.
Case ( 2) is developed to mimic collective migration of cells in a simplified setting, where we simulate migration of cell doublets.As we assumed that once a cell is compressed due to encountering any obstacle, a cell exerts repelling forces on the obstacle to obtain more spaces.Therefore, in the simulations, we observed that the leader cell moves significantly faster than in Case (1), whereas the average velocity of the follower cell is much smaller.In fact, the follower cell is moving faster before it touches the leader cell, since the two cells are initially distanced at the entrance of the channel.However, the leader cell blocks the way of the follower cell, which results in a slower average velocity of the follower cell.The experimental results of Merkher and Weihs [29], show that invasiveness of cancer cells increases when in close proximity, which likely results from additive and synergistic contributions, see [48].This is qualitatively reproduced by our model (see Figure 4.4(b)).In summary, adjacent cells do benefit from each other in various aspects.
There are many possibilities to extend our current model, in terms of parameters considered in the model as well as in mathematical analysis perspectives.The model can be extended, for instance, to include a limit on deformation as large deformations may lead to irreversible cell damage [12] and death.Specifically large deformations may damage the cytoskeleton and the plasma membrane, whereas the springs that maintain the cell intact break.With a smaller deformation, cells may remain in a "softened" structure, which could result in the experimentally observed increased cell speed.We would like to attempt to answer the question in Mak and Erickson [26], that after the cell penetrates through the first constriction, why the following-up constrictions cost less time for the cell when invading serial constrictions.On the other hand, instead of splitting the cell membrane into finite line segments, to keep the smoothness of the cell boundary, we can convert it into a free boundary problem that incorporates a PDE for the cell mechanics, which ensures no self-intersecting of the cell membrane, and probably increases the accuracy of the modelling framework for the fluid-structure interaction between cell and fluid.To further consider limits in deformation, the mechanics of the stiff nucleus can be implemented in the model, as it plays an important role in the invasiveness [48] and the velocity of cells [22].In parallel, the material properties of the channel walls, such as their stiffness, are likely to affect cell migration and can be considered.In addition, to evaluate combined effects of biomechanical and biochemical cell-matrix interactions, biochemical degradation of the cell environment, e.g. by the secretion of matrix metalloproteinase (MMP), can be added.
In summary, the current manuscript demonstrates the bi-directional interactions of cells with their environment and how that is affected by cell proximity.We have extended our earlier cell invasion model to describe transmigration of cells through narrow, flexible channels, such as those existing in various tissues.Furthermore, we have described the bi-directional impact of cells interactions with adjacent neighbors and with their direct environment through the forces that are exerted by the cells.The cell environment may be a pore (through the cross-linked fibrous structure in the extracellular matrix) or any other narrow channel (a small blood vessel).In order to transmigrate, cancer cells deform and exert forces to widen the channel, and two closely adjacent cells are able to interact to accelerate this process.The current model reveals mechanisms underlying the observations of increased invasiveness and the cancer cell metastatic rate upon collective cell migration, when the deformation of the immediate environment is considered.
Cell collides with the other cell

Figure 2 . 1 :
Figure 2.1: A schematic of the cell cytoskeleton and the repelling (pushing) forces when the cell has mechanical contact with an obstacle.(a) Cell cytoskeleton is built up by a series of elastic springs.For every nodal point on the cell boundary, it is connected to the cell centre and two neighbouring points by a elastic spring separately.Here, the arrows indicate the vectors.(b) Repelling forces are exerted by the two cells when they collide with each other and they are compressed.Dashed curves represent the equilibrium shape of the cells, and the black curve between two cells is the contacting surface.Note that the forces, indicated by the arrows, are perpendicular to the black curve.

Figure 3 . 1 :Figure 3 . 2 :
Figure 3.1: The screenshots of the simulation for Case (1), where the next cell appears only when the previous cell leaves the channel completely.

Figure 3 . 3 :
Figure 3.3: The screenshots of the simulation for Case (2), where cells follow each other to enter the channel together.

Figure 3 . 4 :
Figure 3.4: The time series of the average cell speed of Cell 1 and Cell 2 in Case (1) and Case (2).Black solid curve and green dashed curve represent Case (1) and Case (2), respectively.The translucent envelops around the curves are the corresponding 95% confidence interval.

Figure 4 . 1 :
Figure 4.1: The figure shows the Monte Carlo results of Case (1), in particular, the results of the average velocity of the three cells.The box-plot of the average velocity of the two cells, where the median, 25% (Q1) and 75% (Q3) percentile (lower and upper bound of the box), the minimum and the maximum of the data excluding the outliers and the outliers are shown.

Figure 4 . 2 :
Figure 4.2: The figures show the Monte Carlo results of Case (1), in particular, the results regarding the maximal and the minimal width of the channel after each cell exits the channel completely.The maximal and minimal width of the flexible channel are shown as orange and green boxes, respectively, at certain moments.

Figure 4 . 3 :
Figure 4.3: The figure shows the Monte Carlo results of Case (2), in particular, the results of the average speed of the three cells.The box-plot of the average velocity of all three cells, where the median, 25% (Q1) and 75% (Q3) percentile (lower and upper bound of the box), the minimum and the maximum of the data excluding the outliers and the outliers are shown.

Figure 4 .
Figure 4.4(b) shows the box-plot of the microchannel width after the leader cell (Cell 1) exits the microchannel completely in Case (1) and the cell doublet exits the microchannel completely in Case(2).As each dataset is paired, we conduct the Wilcoxon's test in Table?? and we conclude that both the minimal and maximal width in Case (2) are significantly wider than in Case (1) when all the input parameter values are the same.

Figure 4 . 4 :
Figure 4.4: Using the same input values in the Monte Carlo simulations in both cases, the box-plot show cell speed and the channel width after the leader cell and the cell doublets exit the channel in both cases.The orange and green box represent Case (1) and Case (2), respectively.(a) The average speed of Cell 1 and Cell 2 in different cases.(b) The maximal and minimal width after the leader cell (Case (1)) and the cell doublet (Case (2)) exit the channel completely.

Table 3 .
1: Parameter values of the cell invasion model used in Section 3

Table 3 .
(1)Numerical results (i.e.transmigration time and average speed of three cells) of the simulation shown in Figure3.1 for Case(1).

Table 3 .
3: The maximal and minimal width of the channel after different cells exit the channel of the simulation shown in Figure 3.1 for Case (1).

Table 3 .
4: Numerical results (i.e.transmigration time and average speed of two cells) of the simulation shown in Figure 3.3 for Case (2).

Table 3 .
5: The maximal and minimal width of the channel after different cells exit the channel of the simulation shown in Figure 3.3 for Case (2).

Table 4 .
1: Input parameters and outputs of the Monte Carlo simulations for both cases