Mathematical modelling of cell migration: stiffness dependent jump rates result in durotaxis

Durotaxis, the phenomena where cells migrate up a gradient in substrate stiffness, remains poorly understood. It has been proposed that durotaxis results from the reinforcement of focal adhesions on stiff substrates. In this paper we formulate a mathematical model of single cell migration on elastic substrates with spatially varying stiffness. We develop a stochastic model where the cell moves by updating the position of its adhesion sites at random times, and the rate of updates is determined by the local stiffness of the substrate. We investigate two physiologically motivated mechanisms of stiffness sensing. From the stochastic model of single cell migration we derive a population level description in the form of a partial differential equation for the time evolution of the density of cells. The equation is an advection–diffusion equation, where the advective velocity is proportional to the stiffness gradient. The model shows quantitative agreement with experimental results in which cells tend to cluster when seeded on a matrix with periodically varying stiffness.

to infected sites, or as part of the wound healing process (Alberts et al. 2002). Another example where cell migration is of great importance is in the growth and spread of cancer (Wang et al. 2005;Yamaguchi et al. 2005), both in local invasion of the surrounding tissue and in formation of metastases.
Two common mechanisms for cell locomotion are "swimming" and "crawling". Flaggelated cells such as certain bacteria move by rotating a flagella which act as a propeller to drive them forward. A typical example is that of E. coli, and its motion has been described as a "run and tumble" motion, where it alternates between two phases (Berg et al. 1972). During the running phase it moves in a more or less straight path, followed by a tumbling phase, where it re-orients itself with negligible change in position. The second type of locomotion, crawling, is the process we will focus on in this study. It is a complex phenomena which is often described as a cyclic process, consisting of four distinct phases (Kurosaka and Kashina 2008). The first phase is the polarization phase during which the cell defines its front end. The second is the protrusion phase, in which the cytoskeleton changes shape by extending a protrusion at the leading edge. The third phase is the attachment phase during which it adheres to the substrate on which it is crawling. The last phase is the retraction phase, where the cell pulls itself forward (Alberts et al. 2002;Kurosaka and Kashina 2008), and the trailing edge retracts.
Cell migration that occurs in the human body depends heavily on the properties of the microenvironment, in particular its mechanical properties, and has been the subject of extensive research. The extracellular matrix (ECM) making up the microenvironment is a complex fiber network which is made up of proteoglycans and fibrous proteins, mainly collagens, elastins, fibronectins and liminins (Frantz et al. 2010). It functions as a scaffolding for cells, and plays an important role in proliferation, differentiation and survival of cells (Keogh et al. 2010;Murphy et al. 2012). Abnormalities in the ECM can cause a range of syndromes such as osteogenesis imperfecta and Marfan's syndrome (Järveläinen et al. 2009). Important properties of the ECM related to cell migration include fiber density (Kaufman et al. 2005;Sander 2014), fiber orientation which can give rise to contact guidance (Schwarz and Bischofs 2005) and cell-substrate adhesiveness which can induce haptotaxis (Carter 1967). In this article we are interested in the impact of ECM elasticity on cell migration, in particular the phenomena known as durotaxis, where cells tend to move towards regions of higher ECM stiffness.

Durotaxis
Durotaxis was first observed in 2000 by Lo et al. (2000), and has since been observed repeatedly, for different cell types such as fibroblasts (Kuboki et al. 2014), vascular smooth muscle cells (Isenberg et al. 2009), endothelial cells, malignant mammary adenocarcinoma cells (Joaquin et al. 2016), on substrates with different mechanical properties. In the experiment performed by Lo et al. fibroblasts were cultured on a flexible polyacrylamide sheet coated with type I collagen. The sheet consisted of two regions of different rigidity, and cells were approaching the boundary from the soft and from the rigid side. The cell density was low enough for cell-cell interactions to be negligible. The experiment showed that when cells approached the boundary from the soft side, they moved across it into the stiffer region, whereas when they came from the stiffer side the protrusion crossing into the softer side stopped and cells did not cross the boundary. When there is a difference in compliance and the cell exerts forces on the substrate, mass is being pulled toward the stiffer side. However, the authors argue that the skewed mass distribution is not sufficient to explain the phenomena, because the displacements are too small. Instead they proposed that the mechanism is due to cells sensing small changes in stress and strain in the substrate, which is translated into increased traction forces causing a bias in movement direction.

Previous models of cell migration
Mathematical modelling has been used to study a range of different interactions between cells and their microenvironment. Some models focus on the subcellular processes such as dynamics of protrusions and stress fibers, and formation of focal adhesions (Harland et al. 2011;Kim et al. 2012Kim et al. , 2015. At a larger spatial scale, the entities of interest are often individual cells and individual fibers, such as in the model of Schlüter et al. (2012). They developed an individual-based model where cells exert forces on the ECM, and assumed that cells had the capacity to realign matrix fibers, modelled as thin cylinders. Their results showed that the cells in their model had a slight preference for stiffer regions. They also found that the cell speed was lower on very stiff matrices. van Oers et al. (2014) developed a hybrid cellular Potts and finite element computational model. The ECM was modelled as a linearly elastic and isotropic medium, and cells were assumed to exert forces on the ECM so it deformed. They also assumed that a strained ECM is stiffer along the orientation of strain than perpendicular to it, a phenomenon known as strain stiffening. Their model also captured durotaxis through an additional term in the energy function, which accounted for the increased probability of cell extension along local strain orientation, and reduced the probability of retraction in that direction.
A mathematical model to investigate durotaxis was developed by Stefanoni et al. (2011). They developed a 2D numerical model based on a modified version of the Langevin equation, where the stochastic force depends on the local stiffness. The cell is assumed to be able to probe the surrounding tissue and sense the local displacement. The stiffest direction is the one where the displacement is smallest. A probability distribution was then constructed for the angular distribution so that the stochastic force is more likely to point in the direction of higher stiffness. However, the stiffness only impacts the angular distribution and not the radial distribution. The authors investigate first the case of a homogeneous and isotropic substrate, where random motility is recovered, and go on to the case of a biphasic domain similar to the one used in Lo et al. and found good agreement with their results.

The role of adhesion dynamics in durotaxis
During the migration cycle, cells assemble and disassemble focal adhesions. A recent study by Fusco et al. (2017) suggests that the lifetime of adhesion sites depends on the stiffness. The authors investigated the lifetime of focal adhesions on three different substrates with elastic moduli of 30, 200 and 1000 kPa. The average lifetime were shown to be about 6.5, 8 and 12 min respectively. Joaquin et al. (2016) managed to create a 3D extracellular matrix with varying stiffness, while keeping the other mechanical properties, such as fiber density and protein concetration, constant. They investigate the relationship between substrate stiffness and cell migration velocity, as well as the preference for cells to migrate towards stiffer regions. Their study suggests that cell velocity is uncorrelated with the absolute stiffness, but depends on the stiffness gradient. Their results support a previously suggested theory by Plotnikov and Waterman (2013), which states that there may be a stiffness-dependent mechanism by which focal adhesion sites become reinforced, hence increasing their lifetime, which in turn could cause durotaxis.
A number of computational studies have investigated various aspects of adhesion dynamics on cell migration. For example the model by Ziebert and Aranson (2013), where they accounted for substrate stiffness on the dynamics of adhesion sites. Their model consists of two continuum fields, one for the moving cell boundary and one for the dynamics of the actin cytoskeleton. The dynamics of the number of adhesion contacts is governed by a reaction-diffusion equation and depends on the orientation of the actin cytoskeleton, substrate deformation, and accounts for excluded volume interactions. One of their findings was that their model produced predictions that cells tend to stay on a rigid region, when planted on a region with a discontinuity in stiffness. Harland et al. (2011) developed a model for sliding adhesion sites and formation and contraction of stress fibers, to investigate the impact of substrate stiffness on durotaxis. Their model showed that in the absence of a stiffness gradient, the expected drift of a cell was zero. When using a linearly increasing elastic modulus, their model predicted the cell drift to be proportional to the ratio of stiffness gradient to the square of the absolute stiffness. They also found that for most stiffness gradients, there exist an optimal stiffness for which the drift is maximized. They use two different assumptions in their model. One is that the number of stress-fibers (and hence adhesion sites) is constant, and that once a fiber contracts enough so its length is reduced below a certain threshold, it is removed and replaced by a new fiber with random orientation. The second assumption is that the formation of stress fibers is a stochastic process, where the formation intensity is stiffness dependent. Dallon et al. (2013b) introduced a 2D mathematical model for single cell migration, where the cell consists of a nucleus and a number of adhesion sites. The adhesion sites connect to the nucleus through elastic springs, and can be attached or detached to the substrate. Times for attachment and detachment are governed by a Poisson process. When the adhesion site attaches it chooses its position at random. The new position is chosen so that the distance from the cell nucleus is uniformly distributed between 10 µm and 15 µm, and the direction is chosen uniformly distributed between − π/2 and π/2 with probability 0.6, and between π/2 and 3π/2 with probability 0.4. They find that the cell speed is mainly influenced by the mean attachment time, and less by the mean detach time or the strength of the forces the cell exerts on the substrate.
In this work we aim to investigate the impact of stiffness-dependent adhesion lifetimes on the motility of cells, in particular on durotaxis. We do not model subcellular mechanisms, instead we use a model similar to the one introduced by Dallon et al. (2013b) to focus on the individual cell level, where a cell is assumed to consist of a nucleus and fixed number of sites which adhere to the substrate. The adhesion sites connect to the cell nucleus through elastic springs, and the cell moves by updating the position of adhesion sites at random times. To investigate the importance of adhesion lifetimes on durotaxis we assume that the rate at which an adhesion site updates its position depends on the local stiffness of the substrate.
The paper is organized as follows. In Sect. 2 we introduce the mathematical framework for modelling cell migration, in the form of a single-cell stochastic model. In Sect. 3 we analyze the model and derive a macroscopic continuous description for the time evolution of the cell density, and derive an advection-diffusion equation using a diffusion scaling. We then compare the stochastic and deterministic models. To investigate the impact of adhesion lifetimes on directed cell migration, we use biologically relevant parameters and compare our model predictions with experiments in Sect. 4. We summarize our results and discuss their relevance in Sect. 5, as well as possible extensions to our model.

Model formulation
Our mathematical model of a cell is a 1D model inspired by the model introduced by Dallon et al. (2013b). In the original model a cell was assumed to consist of a nucleus with position X and a number of adhesion sites at positions x i which could be either attached or detached to the substrate. A function Φ i (t) was used to indicate if adhesion site i is attached or not, taking value 1 if the site is attached, and 0 if it is detached. The position of a cell X , which is the position of the nucleus, satisfied the following equation derived from Newton's second law of motion under the assumption that the acceleration term can be ignored: where μ is a drag coefficient, α i the spring coefficient of site i, located at position x i . The model also allows for varying rest lengths of the springs, l i . In a later work by Dallon et al. (2013a), they analyzed a simplified version of this model, by informally considering the limit as the spring coefficients become very large. In that case the nucleus reaches its equilibrium position instantaneously, which can be obtained by setting dX / dt = 0. In our model we make a number of additional assumptions. We consider a cell in 1D with n adhesion sites, located at x 1 , x 2 , . . . , x n ∈ R. We assume that adhesion sites have detach time 0, meaning that each time an adhesion site detaches, it instantaneously attaches at a new randomly chosen position. At each site we assume that the rest length l i is zero, and that all spring coefficients are equal. These assumptions together give the position of the cell nucleus as the center of mass of the adhesion sites: Each adhesion site updates its position independently, where the waiting times between updates are exponentially distributed. The new position of an adhesion site is generated from a normal distribution centered around the current position of the nucleus, with variance σ 2 . When an adhesion site updates its position, the nucleus is assumed to change its position instantaneously. A cartoon figure of our mathematical representation of a cell is shown in Fig. 1.
To capture the notion that adhesion lifetimes are dependent on the stiffness of the substrate we assume that cells move on an elastic substrate with elastic modulus given by the function E = E(x) > 0, which is assumed to be at least twice continuously differentiable. The elastic modulus determines the mean waiting time until the next update of an adhesion site. For example, if site i is located at x i , the time until the next update is exponentially distributed with rate parameter λ i = G(E(x i )), where the functional form of G describes how the cells respond to the local stiffness E(x i ), which will be discussed in Sect. 2.1. This means that the average lifetime of site i, located at x i , is 1/λ i . We do not consider the case where the substrate becomes deformed under the force exerted by cells, but plan to investigate this in future work.
A stochastic simulation based on the Gillespie algorithm for the migration of a single cell is given in Algorithm 1. 2: while t < t end do 3: Compute time until next event τ ∼ Exp(Λ). Generate a uniform random number r ∼ U ([0, 1]).

Stiffness sensing
We investigate two different phenomenologically derived principles for how the stiffness impacts the lifetime of adhesion sites. Both are based on the hypothesis that cells probe the stiffness of the substrate at each adhesion site. If the cell exerts a force at each adhesion site, the site will experience a large displacement in regions where the substrate is soft, and smaller displacements in more rigid regions. We assume that the parameter governing the exponential waiting times is given by some basic rate β ≥ 0, and a contribution which is proportional to the relative difference of the magnitudes of displacement at site i, u i , and the average magnitude of displacements at all other sitesũ = 1 n−1 j =i |u j |, therefore we choose The nonnegative parameter γ ≥ 0 influences how much weight is given to the displacement difference. If the displacement |u i | is large compared to the average at the other sites,ũ, then the site updates more frequently. To compute u i we assume that the mechanical behaviour of the substrate is that of a Hookean spring of rest length 0 and spring constant given by E, so that if the cell exerts a force of magnitude F, the magnitude of the displacement is given by In the case of constant stiffness E, we have λ i = β for all i. Moreover, we assume that the force exerted is the same at every adhesion site so the force term F can be factored out and incorporated in the parameter γ . Therefore we simply write The second principle of stiffness sensing we consider is similar but uses absolute displacement instead of relative displacement, which takes the form We refer to Eq. (3) as using the relative relationship, and (4) as the absolute relationship. We will always choose parameters so that it is ensured that λ i ≥ 0.
In the next section we derive a population level description in the form of a partial differential equation (PDE), describing the time evolution of the density of non-interacting cells.

Model analysis
We consider the density of cells q(X , t) at spatial position X ∈ R at time t ≥ 0. By location of a cell we mean the location of its nucleus. Let the function f c (r | X ) be the conditional probability density function (pdf) for a cell making a jump of size r , when it is currently residing at position X , and assume that the waiting times between adhesion site jumps are exponentially distributed. The density of cells at time t + Δt is then given by where Λ is the total rate parameter at which the cell nucleus changes position, and is the sum of the individual rate parameters Λ = λ 1 + · · · + λ n . For both choices of λ i given by (3) and (4), the total rate is Λ = nβ independent of spatial position. The first term on the right hand side describe cells located at X at time t that did not jump in the interval [t, t + Δt], and the integral-term correspond to cells being located at X − r at time t, and performing a jump of size r . The term O(Δt 2 ) describe the gain and loss of density whenever a cell exhibit 2 or more jumps in the time interval of size Δt. By expanding e −ΛΔt in its Taylor series we obtain and rearranging gives We now take the limit Δt → 0 and obtain the partial integro-differential equation, often called a transport equation: Notice that the rate parameter Λ is the sum of all individual rate parameters λ i , as it describes the rate at which any event occurs. In what follows the rate parameters λ i will be spatially varying.

The distribution of nucleus jumps
We now simplify our analysis by considering only two adhesion sites, located at x 1 and x 2 . Since the nucleus is located at the center of mass of x 1 and x 2 , that is X = (x 1 + x 2 )/2, we can express the position of the two adhesion sites as where 2Y is the distance between the two adhesion sites, and Y ≥ 0. It can be shown (see "Appendix C" for details) that in the stochastic model of motion described in Sect. 2, the distance Y converges to a random variable distributed halfnormal with mean σ √ 2/3π and variance σ 2 (1 − 2/π ) /3, whenever the new position of an adhesion site is normally distributed around the current nucleus position with variance σ 2 . We denote the pdf of Y with f Y , given by and proceed to compute the pdf governing nucleus jumps through where f C (r | X , y) is the conditional probability density function of a nucleus jump, conditioned on the current nucleus position X , and the current distance y.
Assume that the current nucleus position X is known, as well as the current distance y. The nucleus then jumps as a result of either site i = 1 performing a jump, or the site i = 2 performing a jump. Whenever the left site makes a jump, the nucleus is on average displaced to the right, and whenever the right site jumps, the nucleus is on average displaced to the left. In both cases the average displacement of the nucleus is exactly half as much as the jump of an individual adhesion site. Therefore, the probability density f C (r | X , y) is a mixture-distribution of two components, given by In other words, it is a sum of the pdfs of two normal distributions with means y/2 and −y/2 respectively, and variance σ 2 /4. The first component in the mixture comes from the case that site i = 1 updates, and the second component from the case that site i = 2 updates. The weights of the mixture w 1 and w 2 are given by It is important to observe that the weights in the mixture are functions of the positions x 1 and x 2 , or equivalently as formulated above, functions of the current nucleus position X and the distance y between the sites. To summarize, the probability density function governing the size of a nucleus jump f C (r | X ), when the nucleus is located at X , is given as a mixture of two components, where each component comes from whether the first or second adhesion site updates. The weights depend on where the sites are located.
To proceed we use λ i defined in (3) and (4) described in Sect. 2.1, and we derive an advection-diffusion equation for the cell density using a large time diffusion scaling.

Large time diffusion approximation
We now consider the case where the rate of updates is given by Eq. (3), namely For two sites x 1 = X − y and x 2 = X + y this gives and the total rate is then given by We now substitute u i = 1/E(x i ) and simplify to get We now expand E in its Taylor series around X to obtain where the term T i is the ith term in the Taylor series of E, depending only on X , and the sum containing only odd indices because all even terms cancel in the series. Using (8) and (9) in Eq. (7), we thus obtain the conditional jump distribution as 2σ 2 /3 dy.
One can see that the last two terms of sums of integrals can be written as: That is, we obtain the sum of higher order non-central moments of the normal distribution with mean r /2 and variance σ 2 /4. We know that the ith moment of a normal distribution with mean r /2 and variance σ 2 /4 is a sum of the form where the coefficient c j is given by the product of a binomial coefficient and the double factorial: We therefore obtain the final result of the form where we have incorporated the constant √ 3 into the constant c j . To summarize, the jump distribution weighted by the total jump rate Λ(X ) = Λ = 2β, as formulated in the transport equation (6), is given by We now insert this expression into Eq. (6) and obtain We now proceed with the following rescaling for jump-size and time: for a small parameter ε. This corresponds to a cell which performs small jumps with high frequency. Writing q(X , t) = q X , 1 ε 2 τ =q (X , τ ) the equation is and we perform the change of variables, r = sε to obtain ε 2 ∂q ∂τ (X , τ ) = −Λq(X , τ ) We now proceed with Taylor expanding the functionq(X − sε, τ ) in space around X . Also notice that the last term containing the higher order moments of the normal distribution is of at least order ε 3 , and will therefore be neglected in the limit ε → 0. For the sake of clarity we keep the second integral, but write it compactly as O(ε 3 ).
Our equation now is We now assume that the stiffness E is varying slowly in comparison to the typical size of a cell, and expand the fraction E /E in its Taylor series around X : Computing the integrals of our equation, we obtain the final result Dividing through by ε 2 and taking the limit ε → 0, we obtain the equation after writing the first two terms as a single derivative of a product. We can perform the exact same analysis when using the absolute relationship (4), with the same assumption on a slowly varying function E, and the same diffusionscaling to obtain the equation which is different only in the advective velocity being 2E (X )/E 2 (X ) instead of E (X )/E(X ). For details regarding the derivation of Eq. (14), see "Appendix A".

Comparison between stochastic simulations and PDE model
We now compare the stochastic simulation to the solutions of the PDE models. We compute the trajectories of 1500 cells using the stochastic simulation, compute the occupancy in discretized intervals of the spatial domain, and normalize so that the total cell density is 1. The PDE model is solved using a Crank-Nicolson finite difference scheme (see "Appendix B" for details) using the same spatial discretization, and total density chosen to be 1. We compare the models for a linearly increasing stiffness function E(x) = 3+0.9x, on the domain Ω = [− 3, 3], with parameters σ = 0.05, β = 0.5, γ = 1 for t ∈ [0, 400]. In the stochastic simulations all cells are starting out at the origin, with both adhesion sites located at the origin, and in the PDE model we use a Dirac-type initial condition, located at the origin. In the stochastic simulation, cells are not able to migrate past the domain boundaries at ± 3. This is implemented by checking if a cell is attempting to move outside the domain. If it is, it simply does not move. In the PDE model we use zero-flux boundary conditions. Figure 2a shows the cell density at four different instances in time, the red dashed line corresponding to the stochastic simulation and the blue line to the solution of PDE (13), when using the relative relationship (3). A comparison between the mean position of a cell is shown in Fig. 2b. Figure 3a, b shows the same comparison when using the absolute relationship (4).
In both cases the stochastic simulation and the solution of the PDE show good agreement.

Periodically varying stiffness result in cell clusters
In this section we use our model to investigate if cells cluster on stiff regions, when seeded uniformly at random on a substrate with a sinusoidal variation in stiffness. This is motivated by the study performed by Joaquin et al. (2016), where the matrix was modified so that there was a spatial variation in the stiffness, in the form of a sine wave. In the experiments four different cell types were used, and all four cell types formed clusters on the stripes of high stiffness but had different final morphologies. From a correlation analysis they concluded that the stiffness gradient was the driving factor in this durotactic phenomena, and not the absolute stiffness.
We now want to investigate if our model will produce the same behaviour, namely that cells cluster on regions of high stiffness. To do so we need to find appropriate parameter values. Recall that the mean distance between two sites is a random variable with mean σ √ 2/3π . Now we wish to choose the cell sizes to be on average 20μm, which holds if we choose which requires us to set σ ≈ 40 µm = 0.04 mm. We choose the time to be measured in hours and the rates to have the unit updates per hour. The stiffness function is chosen to match the stiffness profile from the experiments, and is given by and the domain to be Ω = [− 3, 3] mm, so there are a total of 8 peaks of high stiffness in the domain. We use uniformly distributed initial positions of cells in the stochastic simulation, and constant cell density throughout the domain in the PDE model, normalized so the total density is 1. For our choice of σ , the nucleus jump is on average of size 9.21 µm, and the rate at which it jumps is 2β, in the absence of a stiffness gradient, so choosing β = 1 results in cells which move about 20 µm per hour on average, which is in a realistic range for many cell types. Figure 4 show the results of a simulation with 50 cells using the relative relationship (3) along with the solution to the PDE model (13) after 24 h, for four different values of γ ∈ {1, 2, 3, 4}, with unit h −1 in the case of using the relative relationship, and mm −1 h −1 in the case of using the absolute relationship. It can be seen that for larger values of γ , the cells tend to cluster at the peaks. This is expected since the parameter γ governs how big impact the difference in displacement (relative or absolute) has on the update rates. The stiffness function is shown in the plots to illustrate where the peaks are located, and the position of cells are plotted on this superimposed curve to better illustrate where they are located. Note however that the cells move on a 1-dimensional line. Figure 5 shows the same comparison when using the absolute relationship (4) and the corresponding PDE model (14).
To better understand how the lifetime of adhesion sites impact clustering we also investigate the difference between lifetimes of the two adhesion sites, for γ ∈ {0, 1, 2, 3, 4}. We do this by computing the minimum and maximum of λ 1 and λ 2 for each cell as it is moving around, and average over 2000 simulations. The results can be seen in the table below. This means that for example, in the case of using γ = 4 the site that updates most frequently does so on average 15% more often than the site that updates less frequently (Table 1).
We now introduce a way to measure to which degree cells cluster by defining the following quantities. We define H to be the part of our spatial domain where the  stiffness is higher than the average stiffness over the entire domain, and L to be the part where the stiffness is lower than the average. Then we define the total cell density on high stiffness and low stiffness regions, at time t, through respectively. We consider the ratio i.e. ratio of cell density located on a stiff region, to the total cell density. Notice that this measure ranges from 0 to 1, and takes the value 1/2 in the case that the cells are uniformly spread between the soft and stiff regions. Figure 6a, b show the evolution of the ratio R(t) for different values of γ , corresponding to the experiments shown in Figures 4 and 5. The ratio is computed using the PDE model. In order to compare our predictions with the experiments carried out by Joaquin et al. (2016), we analyzed a supplementary video published with that study. Images were extracted at 4, 8, 12, 16, 20 and 24 h and segmented using Otsu's method implemented in MATLAB (Otsu 1979). The total cell area was calculated in the high and low stiffness regions. From these quantities we calculated the fraction of the cell population being located on the stiff region for each image, corresponding to (18). The resulting data is shown in Fig. 6c where we have fit the parameter γ using least squares in our two PDE models and calculated the ratio defined in (18). It can be seen that the two PDE models provide near identical predictions, but use different values of γ . When fitting the PDE model based on the relative sensing we obtained γ = 13.28 and when using the model based on the absolute sensing we obtained γ = 6.5.

Conclusion
The mechanisms causing durotaxis remain poorly understood. A hypothesis has been proposed, that through repeated tugging, cells probe the stiffness of their surrounding and adhesion sites on stiff regions become reinforced (Joaquin et al. 2016;Plotnikov and Waterman 2013). In this work we developed a mathematical model for investigating directed cell migration, under the assumption that adhesion sites become reinforced on stiff regions. The model assumes that a cell migrates by updating the position of its adhesion sites at random times, and a new position is chosen randomly around its current nucleus position. The cells are assumed to be able to sense the local stiffness by exerting a force at every adhesion site. We propose two sensing mechanisms. The first is based on the relative displacement at adhesion sites, and the second based on the absolute displacement. To model reinforcement we assume that adhesion sites where the displacement is small update less frequently than sites where the displacement is large.
From the individual based model we derived a continuous description in the form of a partial differential equation for the density of cells. The resulting PDEs are advection diffusion equations, where the cells are predicted to move in the direction of increas-  Joaquin et al. (2016) and the two PDE models when γ is fitted using the least-squares method ing stiffness. In the case of assuming a relative sensing mechanism cells move with a velocity proportional to ratio of stiffness gradient to absolute stiffness, and in the case of assuming an absolute sensing mechanism, the velocity is proportional to the ratio of stiffness gradient to the square of absolute stiffness. The stochastic and deterministic models showed good agreement for the length and time scales typically observed in experiments.
Another type of taxis observed in biology is chemotaxis, where an organism responds to the concentration of some chemical to which it is chemosensitive. A continuum description is given by the classical chemotaxis equation (Erban and Othmer 2004) for the density p of cells where S(x, t) is the chemical concentration, and χ(S) describes the chemotactic sensitivity, and in general S may be prescribed to evolve according to some evolution equation. It is easily seen that if we assume that the chemical is independent of time, and choose the sensitivity to be of the form χ(S) = 1/S(X ) or χ(S) = 2/S 2 (X ), the chemotactic equation become Eqs. (13) and (14) respectively. To investigate how well the model is able to capture durotaxis, we chose biologically relevant parameters, and a sinusoidal stiffness function. This type of experiment has shown to produce clusters of cells on the "peaks" of high stiffness (Joaquin et al. 2016). Our models predicted cells to cluster, to different degrees depending on the parameter γ , which governs how much the adhesion lifetime is impacted by the stiffness. For example we showed that, to expect about 70-80% of the cells to be located in the stiffer regions, the lifetime of adhesion sites on the stiff side of the cell has to remain for 15-30% longer than the average lifetime of sites on the soft side of the cell. Although no experiments have yet investigated the lifetime of individual adhesions on a material with a stiffness gradient, we know from the study by Fusco et al. (2017) that the average lifetime of adhesions can be twice as long on stiff regions compared to soft regions. This suggests that a 15-30% difference in update frequency lies within a realistic range.
No visual difference between the two PDE models can be seen from Fig. 6c when γ is chosen so that the models best fit the data. There is a small discrepancy at early and late times between model and data. The early time difference could possibly be due to effect such as cell-cell adhesion or a higher overall cell speed. The faster decay for large times could possibly be the result of crowding, which our models does not account for. In solutions of the PDE models shown in Fig. 6 we used a basic motility rate β = 1. If it is increased to 2, corresponding to an overall cell speed of about 40 µm per hour, the fit becomes better. However, based on the result of the experimental study (Joaquin et al. 2016), a speed of 20 µm per hour is more realistic.
To be able to disentangle how much of the observed drift of cells is the result of reinforcement of adhesion sites, one would need to measure the lifetime of individual adhesion sites within a cell, as well as the cell velocity. This would in turn be sufficient to estimate our model parameter γ , which describes how much of the displacement difference that influences adhesion site lifetimes. To our knowledge such a study has not yet been conducted. However, the lifetime of individual adhesion sites has been measured (Fusco et al. 2017), and it is well known that matrices with a range of stiffness profiles can be designed. It remains to study the velocity of cells and the lifetime of adhesion sites on a matrix with a stiffness gradient.
The 1D model of migrating cells might appear unrealistic for modelling 3D cell migration in a fibrous ECM. However, there are environments in which cell migration takes places in a highly organized ECM, along fibers. Examples include glioma cell migration in white matter tracts (Bellail et al. 2004), and stem cell migration in rats, in both transplantation into the auditory nerve and after spinal cord injuries (Palmgren et al. 2012;Tysseling et al. 2010). It is also believed that 1D cell migration along fibers is more similar to migration in 3D matrices than cell migration on typical 2D planar surfaces (Doyle et al. 2009(Doyle et al. , 2013. This makes 1D models suitable for modelling; they are more analytically tractable than 3D models, and they can provide results which more accurately translate into in vivo cell behaviour.

Possible extensions
The purpose of our model is to investigate the impact of varying adhesion lifetimes on directed single cell migration. However, to more realistically model a large population of cells, one should incorporate interactions between cells, such as cell-cell adhesion or crowding effects, and interactions between cells and the ECM. For example, some cells have the ability to remodel the ECM, either through exerting physical forces, or through chemicals which break down the ECM. The former could be captured by explicitly modelling the ECM as an elastic material that deforms as cells exert forces on it. This could lead to interesting feedback between cells and the substrate. Another interesting property is the so-called strain-stiffening that can occur, which can alter the response of cells to the increased stiffness of the ECM close to large collections of cells e.g. close to tumours. Finally, should future experiments imply that reinforcements of adhesion sites play only a small role in the observed drift in durotaxis, our model can easily be extended to include for example asymmetric forces exerted by the cell or an asymmetric distribution governing new positions of adhesion sites, both of which would influence the drift velocity. One can also consider using a correlated random walk as the underlying mechanism of motility. In our model we assumed that the new position of an adhesion site was distributed normally around the current nucleus position, but one can make it more general by considering any other suitable distribution. However, in contrast to typical random walk models where only the first two moments of the jump kernel are of interest, using any other distribution in our model result in non-standard distributions of the nucleus jump distribution (data not shown). To the best of our knowledge, the normal distribution constitutes a special case where computations of all probability density functions are possible.
Understanding the interactions between cells and their environment is of fundamental importance and this paper will hopefully represent a step toward a deeper understanding of cell migration and its relation to substrate stiffness.

A Derivation of the advection-diffusion equation for the absolute relationship
We now assume that which can be rewritten in terms of the stiffness, using |u i | = 1/E(x i ): .
We expand E in its Taylor series around X and obtain where T i (X ) is the ith coefficient in the Taylor series. As we showed in Sect. 3.2, under the rescaling of space and time, all terms of order y 3 and higher will go to zero, as ε → 0, and is therefore not included in the proceeding calculation. We therefore have and Λ = λ 1 + λ 2 = 2β, along with mixture weights We then obtain Inside our transport equation we have the term Λ as well, so we therefore cancel the factor 1/2β. The transport equation is now ∂q ∂t (X , t) = −Λ(X )q(X , t) We now perform the same rescaling of jump-size and time: for a small parameter ε. Writing q(X , t) = q X , 1 ε 2 τ =q (X , τ ) our equation becomes ε 2 ∂q ∂τ (X , τ ) = −Λ(X )q(X , τ ) and we perform the change of variables sε = r and obtain ε 2 ∂q ∂τ (X , τ ) = −Λ(X )q(X , τ ) + Rq (X − sε, τ )e − s 2 2σ 2 /3 1 2πσ 2 /3 2β + γ sε 2E (X − sε) E 2 (X − sε) ds.
We now expandq in its Taylor series around X , as well as E and E in series around X and obtain ε 2 ∂q ∂τ (X , τ ) = −Λ(X )q(X , τ ) Dividing through by ε 2 , and taking the limit ε → 0 we finally obtain the equation the Crank-Nicolson method is the average of the forward Euler method at time n and the backward Euler method at time n + 1 on a discretization of the spatial domain x ∈ [−L, L] into M points of equal spacing h, x i = −L + hi, and time domain discretized into N points of equal spacing Δt, t n = nΔt.

C Distribution of adhesion site distances
In this section any lowercase f is a probability density function and any uppercase F a cumulative distribution function, with the corresponding random variable as a subscript. Assume there are two adhesion sites, x 0 1 and x 0 2 at time t = 0. After k updates, the positions are x k 1 and x k 2 . Denote the distance between the sites by Δx k = x k 1 − x k 2 , and the position of the nucleus by μ k = x k 1 + x k 2 /2. When site i updates, its new position is given by where W is a normal random variable with mean 0 and variance σ 2 , and i ∈ {1, 2}. The cumulative distribution function for Δx k is F Δx k (z) = P(Δx k ≤ z).
We now consider the cumulative distribution after k +1 steps. Here we assume that the position of the adhesion sites at time-step k is x 1 and x 2 , dropping the superscript k, and we denote by x + 1 and x + 2 the update position of site 1 and 2 respectively, depending on which one that updates. Assuming that site i = 1 updates with probability λ 1 and that site i = 2 updates with probability λ 2 , with λ 1 + λ 2 = 1 we obtain and re-writing x + i = μ k + W , and using μ k = (x 1 + x 2 )/2 we obtain We now investigate the first term, and notice that Notice that the second equality is obtained from the fact that Δx k and W are independent. If we differentiate with respect to z we obtain the pdf of distance between adhesion sites, and using the fact that W is independent of Δx k we obtain and similarly for the second term we get Together these two terms gives us the expression