Friction-limited cell motility in confluent monolayer tissue

Mechanical forces are important factors in the development, coordination and collective motion of cells. Based on a continuum-scale model, we consider the influence of substrate friction on cell motility in confluent living tissue. We test our model on the experimental data of endothelial and cancer cells. In contrast to the commonly used drag friction, we find that solid friction best captures the cell speed distribution. From our model, we quantify a number of measurable physical tissue parameters, such as the ratio between the viscosity and substrate friction.


Introduction
Various models of collective cell motion exist in the literature, including particle-based models [1][2][3][4], cellular Potts models [5][6][7], vertex models [8], phase field models [9] and continuum-scale models with [10][11][12][13][14][15][16] and without cell polarization [17,18]. Often, continuum-scale models have the advantage of being analytically more tractable than the corresponding particle-based models. However, a general problem with most models is how to choose the right constitutive relation between quantities, such as cell strain and stress or to decide on the right functional form of various dissipative processes. Observations at the level of individual cells suggest that cell deformation is viscoelastic [19,20], i.e. at short time scales, any deformation responds elastically to a mechanical loading, whereas loads applied over a longer period will be followed by a viscous relaxation and thus an irreversible deformation.
Several constitutive relations have been employed to describe the deformation of single cells, see, for example, [21]. Examples of such relations include a Maxwell fluid ( figure 1(A)), consisting of a dashpot in series with a spring, and a Kelvin-Voigt solid ( figure  1(B)) consisting of a dashpot in parallel with a spring. Depending on the coupling of the basic building blocks (dashpots and springs) a broad range of fluidlike and solid-like behaviors can be obtained. One prominent example of a basic rheological model is the standard linear solid, where the dashpot in the Kelvin-Voigt element is coupled in series with a spring. When subjected to a sudden strain, the stress in the standard linear solid model decays exponentially towards a constant non-zero value.
Aggregates of cells in tissues also behave viscoelastically [22]. Based on stress relaxation experiments on freely suspended cell monolayers, Harris et al [23] used the standard linear solid model to estimate an apparent viscosity for Madin-Darby canine kidney cell monolayers. Guevorkian et al [24] obtained an apparent viscosity for murine sarcoma cell aggregates using a micropipette aspiration technique and the standard linear solid model in series with a dashpot. Forgacs et al [25] performed parallel plate compression experiments on chicken cell aggregates and estimated the viscosity using a rheological model of two parallel-coupled Maxwell fluid elements. In addition, the viscosity of breast cancer tumors has been measured using ultrasonic shear-wave imaging experiments combined with a Maxwell model [26].

Methods
We shall here consider the Oldroyd-B rheological model (see, for example, [27]), which follows from the immersion of a Maxwell fluid of viscosity η 1 and shear modulus G in a Newtonian fluid of viscosity η 2 (figure 1(C)). Defining the total viscosity η 0 = η 1 + η 2 , the relaxation time λ 1 = η 1 /G and the retardation time λ 2 = (η 2 /η 0 )λ 1 , the Oldroyd-B constitutive relation is: where ∇ σ and ∇ γ are the upper convected derivatives of the deviatoric stress tensor σ, γ = 1 2 (∇v + (∇v) T ) is the strain rate tensor and v is the local mean velocity. Experiments measuring the stress relaxation of tissue after a sudden compression have found that the stress relaxes exponentially with one or two characteristic time scales depending on the external loading conditions [23,25]. The Oldroyd-B model only captures a single characteristic relaxation time scale λ 1 under the same experimental conditions. While more rheological elements could be added to account for more time scales, we shall here aim at simplicity and stick with a description based on a single time scale.
Moreover, the upper convected derivative can often be greatly simplified in slowly moving tissues. Defining a characteristic time t γ on which γ changes, a characteristic flow speed U v and length scale L v , the upper convected derivative can be cast in dimensionless form: (2) In the limit where the time scale, L v /U v , used for a fluid element to traverse a length scale L v , is large compared to the characteristic (relaxation) time t γ , the quantity t γ U v /L v (the Deborah number) is small and the upper convected derivative reduces to a partial derivative with respect to time, In the experimental data, we shall consider below that the used cell lines have a characteristic flow velocity U v ∼ 1 µm min −1 and a characteristic size of L v ∼ 20 μm (see appendix A). The typical relaxation time is of the order of a minute. The prefactor in the upper convected derivative is therefore relatively small t γ U v /L v ∼ 0.05.

Substrate friction
The friction between motile cells and the substrate is often considered to be linear in the cell velocity [1,10,17,18,28,29]. That is, the friction is similar to the viscous drag on a fluid flowing over a surface at low Reynolds numbers. In contrast, measurements on individual cells suggest that the traction force exerted by a cell is non-trivially related to the filmamentous actin speed [30]. The traction force seems to follow an increasing trend for low actin speeds, crossing over to a decreasing trend for larger actin speeds.
It remains an open question as to whether these measurements can be generalized to the larger scale dynamics of, for example, a confluent cell layer. From a modelling point of view, various forms of friction can be derived from the dynamics of discrete contact points between cells and the substrate, see, for example [31]. Here, we consider the dynamical consequences of various friction terms using, as a starting point, a simple stochastic differential equation. In conclusion, we find that a term which crosses over to solid friction (or socalled Coulomb friction) for larger velocities offers the best quantitative agreement with our data and is consistent with the experimental observation that the speed distribution of the cells have exponential tails (see figure 3(A)). The Coulomb friction does not depend on the speed, only on the velocity direction v = v/|v| scaled with a friction coefficient α. For simplicity, let us for a moment disregard the interaction between cells and consider only the balance between the motility force and friction force of an individual cell. We first consider the motility force to be a basic stochastic process similar to the Ornstein-Uhlenbeck process, which describes a viscously damped and stochastically-driven particle. This process has also been commonly used to describe cell motion [32,33]. The process is formulated in the cell velocity v, where the motility force ξ is a Gaussian white-noise field of strength ξ i (t)ξ j (t ) = 2δ(t − t )δ ij and the friction (or damping) force is given by ψ(v). Note that we have scaled the damping term and time with a general amplitude of the noise field. If the friction force is assumed to be isotropic and always acts in a direction opposite to the local velocity, it will assume a form Under sudden stress, the spring of elastic modulus G deforms instantaneously whereas the dashpot deforms at a constant rate like a fluid of viscosity η. When the Maxwell element is released, the spring regains its original length, but the dashpot irreversibly maintains its deformation. (B) Kelvin-Voigt solid. Under sudden stress, the Kelvin-Voigt solid deforms with a characteristic time scale η/G. The deformation is reversible, and, when released, the Kelvin-Voigt solid regains its original shape. (C) Oldroyd-B fluid. When subjected to a sudden strain 0 , the stress decays exponentially with a timescale η 1 /G towards zero.
From the stochastic differential equation (4), one obtains the corresponding Fokker-Planck equation for the speed distribution, Because of the rotational symmetry of the damping term, the equation only contains the radial term in velocity space, i.e.
We find the stationary distribution P(v) by demanding that the time derivative vanishes, i.e. we have the equation For a given ψ(v), the solution to the steady-state equation can be written in the form where K is a normalization constant. In the case of Coulomb friction, where for all v = 0 the term ψ(v) = −α is a constant, the steady-state distribution becomes an exponential A stochastic differential equation similar to the one considered here has been suggested to describe the dynamics in granular systems [34]. For completeness, we write the contrasting Gaussian tail achieved when a linear damping term is used, −αv , At a microscopic level, a velocity-independent friction arising from contact with the substrate would require the contact area of individual cells to be roughly constant, or similarly require the number of bonds formed between a cell and the substrate to be independent of the cell velocity. For this to be the case, the rate of bond formation and breakage would have to be equal. In dry friction [35], however, the friction force is often observed to decrease at larger slipvelocities, suggesting that the bond formation cannot catch up. A decreasing friction force would change the tail of the velocity distribution of equation (9). We emphasize that the data we shall consider here does not reveal a consistent departure from an exponential tail at larger velocities.
As suggested from observations of individual cells [30], the singular transition of Coulomb fric-tion when going from no motion to motion might not be fully representative of the substrate friction. In the following, we therefore consider an approximately linear friction for small velocities crossing exponentially over to a velocity-independent friction for larger velocites by using a friction term of the The characteristic speed w 0 determines the point of crossover to velocityindependent friction. Note that this friction would still result in an exponential speed distribution for larger speeds.

Equations of motion
The tissues, we consider, are confluent monolayers of cells residing on a substrate. The layers are taken to be incompressible such that the divergence of the local mean velocity vanishes and the projected area of each cell is conserved. We model the tissue dynamic by a momentum balance equation where cell motility is introduced via a persistent random forcing. From the considerations in the previous section, the interaction between tissue and substrate will be given by a friction term of the form Here, ρ is the mean density, p is pressure and σ is the deviatoric stress tensor. The term m(x, t) is the force driving the cell motility, λ m is a persistence time and φ is filtered white Gaussian noise. Note that in contrast to the considerations above, the noise term, in the case of non-interacting motile cells, would lead to persistent random motion. Similar persistent random motion has been considered, for example, in [36,37]. As each cell is a coherent body, there is a minimum length scale-similar to the cell size-below which the velocity field is approximately constant. We therefore impose a length scale m on the random forcing by filtering the noise field ξ with a Gaussian function of The strength of the white Gaussian noise field was taken where the indices i, j run over the spatial directions x, y. It follows that the filtered noise field will have an exponentially decaying spatial correlation: We note that the cell motion is strongly over-damped (the flow occurs at negligible Reynolds numbers) and it is therefore safe to neglect the inertial terms on the left-hand side of equation (12).

Model parameters
The model described above has eight parameters, which we have listed in table 1. We can then use the viscosity η 0 , along with the density ρ, the relaxation time λ m and a characteristic length m , to cast the model equations in a dimensionless form, where we have defined the dimensionless control parameters: That is, we have five parameters to fit by the numerical simulation of equations (17)- (21). From these five parameters, it is not possible to extract all the eight parameters in table 1. For example, we can only determine the ratio between the friction constant and the viscosity as well as the ratio between the viscosity and the noisy driving force. We note that we cannot explicitly determine the viscosity. In figure 2, we have computed the distribution of speed for different values of w 0 used in the friction term for typical tissue parameters. We see that for large w 0 , we recover the Gaussian tail characteristic for the Ornstein-Uhlenbeck process, whereas, when w 0 is smaller than typical speeds in the tissue, we have an exponential tail in the speed distribution. We apply our model to the experimental data on epithelial and endothelial tissue from [14,38]. The model is simulated numerically using a pseudospectral method on a two dimensional (2D) periodic domain. The Fourier and inverse Fourier transforms in the spectral formulation are performed using the fast Fourier transform algorithm, and the non-linear terms are evaluated in real space. Time integration of the stress tensor and the noise term are performed using an exponential time differencing scheme [39]. In each time step, the velocity field and the pressure is found by a relaxation procedure. The periodic 2D computational domain consists of 256 × 256 grid points corresponding to a box of length ∼ 200 μm in physical units. The time step was ∼0.01 min in physical units. From the model parameters, we recover parameters with physical units for each cell type by matching the  (17)- (21). For the large cross-over value, w 0 , in the friction term, we go from a distribution with an exponential tail to a distribution with a Gaussian tail.
simulation mean speed v 0 with the experimental mean speeds (listed in appendix A) and the simulation velocity correlation length 0 with the experimental velocity correlation length. Finally, the main constituent of biological tissue is water; we therefore do not consider the density ρ as a free parameter of the model.

Results
We fit the model by performing a parameter grid search and by choosing the parameters simultaneously resulting in the smallest relative squared error between the experimental P e (v) and simulated speed distribution P s (v), the spatial velocity correlation function (simulated C s vv (r) and experimental C e vv (r)) and the temporal velocity correlation function (simulated C s vv (t) and experimental C e vv (t)), i.e. the parameters that minimize the sum where the averages are taken over the limited ranges shown in figure 3. We note that the choice of error function is somewhat arbitrary and that individual parts could have been weighted differently. Nonetheless, we believe that similar results would be obtained from related error functions. The set of best fit parameters for each cell type is listed in table 2 and follows from the limit of pure solid friction also considered in [38]. The model reproduces, fairly accurately, the spatial and temporal velocity correlations-see figures 3(A)-(C) and where the results for both the epithelial and endothelial tissue is consistent with the findings in [38]. The choice of friction force in the model leads to an exponentially decaying tail of the speed distribution. Note that the only slight difference between the model and experiments is in the spatial correlation function, where the model has a slightly more pronounced dip than in the experimental data (see figure 3(B)). The spatial correlation in the model predominantly enters through the correlation length scale introduced by the structured noise (see equation (16)). The length scale of the noise dominates over the one introduced by the viscous forces. Below the length scale of the structured noise, our model slightly overestimates the spatial correlation, which is also apparent in figure 4 where indeed, locally, the correlation in the models seems stronger than in the experiments. Figure 4 shows an example of a velocity field obtained from a best fit together with a random zoom on an experimentally measured velocity field. In the experimental data, however, the spatial correlation is long ranged, probably due to stronger interaction forces between the neighboring cells (seen by the slightly slower decay in the velocity field in the experimental panel of figure 4). To model the cross-over between the small and large scale correlations more accurately, we would have to generalize further our visco-elastic terms.
Finally, we should mention that numerical simulations where the friction term is replaced by a drag term −αv do not lead to a significant change in the correlation functions. The choice of a drag or a Coulomb friction mainly affects the speed distribution. The velocity correlation functions can be calculated analytically when advection of the noise field is neglected and when a drag term is used instead of a friction term (see appendix B). The resulting analytical correlation functions are displayed as black lines in figure 3 and closely resemble the simulations.

Conclusions and discussion
In a continuum model of tissue dynamics, we find that a solid friction term between cells and a substrate in a mono-layered tissue describes the exponential tail in the cell speed distribution fairly well. This is in contrast to many existing continuum and particle models, which typically consider a drag term linear in the velocity field and therefore a Gaussian tailed speed distribution. In addition to the experimental data considered here, a similar exponential tail has been observed, for example, for dilute suspensions of the MDA-MB-231 cell line [40]. Non-Gaussian tails have also been observed experimentally for other types of tissue [36,37,41]. To our knowledge, little is known about the cell-substrate dissipation process on the tissue scale, and it would be interesting to obtain firmer experimental knowledge on the appropriateness of either a linear drag or a Coulomb type friction. By modifying the frictional force to make a cross over between a drag-like friction at small velocities to Coulumb friction at larger velocities, our model can  further capture speed distributions with a Gaussianlike intermediate range, consistent with experimental observations [1]. The self-propulsion of cells in the proposed model incorporates a length scale and a persistence time of the local motility force field, which accounts for the finite extent of a single cell and for the tendency of a single cell to change its velocity in a certain time scale. Several papers have used related but more elaborate approaches, where the motility force field evolves in time due to some specified dynamics. Basan et al [42] proposed that the local motility force tends to align with the tissue velocity, and, combined with an assumption of cell locomotion being suppressed by neighbors, Zimmermann et al [43] were able to model the traction patterns observed experimentally in spreading epithelial tissue [44,45]. Other authors have envisioned the polarization field as a 2D nematic liquid crystal [10,13,29], which allows for complex flow patterns with vortices and jets. In contrast to these existing models, the proposed motility evolution contains no assumptions of velocitymotility alignment or nematic behavior.
Including a motility force explicitly, as we have done in this paper, has the advantage that the local velocity need not be aligned with the local motility force. This behavior has been observed in an expanding monolayer, where the velocity and the local motility (traction) force under certain circumstances were anti parallel [45]. A natural next step is therefore to incorporate tissue boundaries in the proposed model, such that monolayer expansion and the generated tractions can be studied. This would also allow the model to be compared with, and used to study, the classical scratchwound assay experiment [46][47][48], the observed fingering of tissue edges [48,49] and the propagation of strain rate waves in spreading tissue [50].
where Γ indicates the integration path in the complex plane for the inverse Laplace transform, and we have defined: The velocity component correlation is proportional to: We note that the term e −k 2 is a result of the noise being filtered. If φ had been white Gaussian noise, the term would not have been present.
The velocity correlation function C vv (r, τ ) of r = |x − x | and τ = t − t > 0 is then proportional to: where the Fourier transform has been turned into a Hankel transform and J 0 (kr) is the zeroth order Bessel function of the first kind. The function f (k, t) can be found analytically by applying the inverse Laplace transform to equation (B.2) and also the integral over t 1 can be performed analytically. Neglecting the transient behavior (t, t 0), we find the spatial correlation function for equal times τ = t − t = 0: