Phenomenlogical Modeling of Durotaxis

TL;DR

This paper presents a phenomenological model of cell durotaxis based on substrate stiffness-dependent focal adhesion distribution, reproducing experimental persistence and predicting directional movement towards stiffer substrates.

Contribution

It introduces a novel 2D model linking focal adhesion distribution to cell motion and validates it with a Fokker-Planck equation, highlighting the role of adhesion in durotaxis.

Findings

Model reproduces increased persistence on stiffer substrates

Predicts cell movement towards stiffer regions (durotaxis)

Derives a validated Fokker-Planck equation for the model

Abstract

Cells exhibit qualitatively different behaviors on substrates with different rigidities. The fact that cells are more polarized on the stiffer substrate motivates us to construct a two-dimensional cell with the distribution of focal adhesions dependent on substrate rigidities. This distribution affects the forces exerted by the cell and thereby determines its motion. Our model reproduces the experimental observation that the persistence time is higher on the stiffer substrate. This stiffness dependent persistence will lead to durotaxis, the preference in moving towards stiffer substrates. This propensity is characterized by the durotaxis index first defined in experiments. We derive and validate a 2D corresponding Fokker-Planck equation associated with our model. Our approach highlights the possible role of the focal adhesion arrangement in durotaxis.

Figures

b). In Fig. 2c and 2d, we show trajectories of cells simulated on both uniform soft and hard substrates. Consistent with the experimental observation [[27](#bib.bib27)], cells crawl more efficiently on stiffer substrates.

a shows the initial condition in the all cases discussed in the following; in particular we apply a narrow gaussian distribution to approximate the $\delta$ in $p(x,\Phi;t=0)=\frac{1}{2\pi}\delta(x),\Phi\in(0,2\pi)$. The solution shows a peak in $x$ which varies from being at positive values (for $\Phi\simeq 0$) to negative ones (for $\Phi\simeq\pi$); the peak heights are independent of $\Phi$ as expected via rotational symmetry.

c). We also define $\hat{p}(x,t)=\intop_{0}^{2\pi}p(x,\Phi,t)d\Phi$ and compare our PDE result for this quantity with direct simulations (Fig 5d). The very good agreement between PDE and direct simulation results validates the Fokker-Planck equation approach. Note that the above trends continue with the population continuing to break up into a peak at positive $x$ and a straggler peak at negative $x$ corresponding to cells that have wandered out into the uniform less stiff side of the gradient profile.

Equations

$$\sum _{i=1}^{N_f}{r}_i\sin ({\theta }_i-\Delta \Phi )=0$$

$$\tan \Delta \Phi =\frac{{\sum }_{i=1}^{N_f}{r}_i\sin ({\theta }_i)}{{\sum }_{i=1}^{N_f}{r}_i\cos ({\theta }_i)}$$

$$\int \cdots {\int }_{-{\theta }_{max}}^{{\theta }_{max}}\left(\prod _{i=1}^{N_f}\frac{d{\theta }_i}{2{\theta }_{max}}\right){\arctan }^2\left(\frac{{\sum }_{i=1}^{N_f}\sin ({\theta }_i)}{{\sum }_{i=1}^{N_f}\cos ({\theta }_i)}\right)$$

$$\sigma (k)=\frac{1}{\alpha k+\beta }$$

$$<x^2>=v^2{\tau }_p^2(\frac{t}{{\tau }_p}+e^{-\frac{t}{{\tau }_p}}-1)$$

$$k(x)=\{\begin{array}{ll}{k}_{left}=1& -1000<x<-L\\ k_{left}+\frac{(k_{right}-k_{left})}{2L}(X+L)& -L\le X\le L\\ k_{right}=5& L<x<1000\end{array}$$

$$DurotaxisIndex(t_i)=\frac{N_{right}-N_{left}}{N_{right}+N_{left}}$$

$$\begin{array}{rl}p(x,\Phi ;t+dt)& =\underset{0}{\overset{2\pi }{\int }}p(x-v\cos ({\Phi }_0)dt,{\Phi }_0;t)\\ & f(x-v\cos ({\Phi }_0)dt,{\Phi }_0-\Phi )d{\Phi }_0\end{array}$$

$$f(x,{\Phi }_0-\Phi )=a(x)e^{-\frac{({\Phi }_0-\Phi {)}^2}{2\sigma (x{)}^{2}dt}}$$

$$\frac{\partial p}{\partial t}=-\frac{\partial p}{\partial x}vcos\Phi +\frac{\sigma (x{)}^2}{2}\frac{{\partial }^{2}p}{\partial {\Phi }^2}+vsin\Phi \frac{\partial }{\partial x}[\sigma (x{)}^{2}dt\frac{\partial p}{\partial \Phi }]$$

$$\frac{\partial p}{\partial t}=-\frac{\partial p}{\partial x}vcos(\Phi )+\frac{\sigma (x{)}^2}{2}\frac{{\partial }^{2}p}{\partial {\Phi }^2}$$

$$\begin{array}{rl}& \frac{dx}{dt}=vcos(\Phi )\\ & \frac{d\Phi }{dt}=\eta (t)\end{array}$$

Full text

Cells are capable of sensing and responding to the mechanical properties of their external environment. For example, cytoskeletal stiffness [1], cellular differentiation [2, 3, 4, 5] and cell morphology and motility [6, 7, 8] are all strongly influenced by ECM stiffness. In particular, it has been shown experimentally that cells prefer crawling towards the stiffer parts on substrates with spatially varying rigidity, a property which is referred to as durotaxis. Durotaxis is a universal property of motile cells, despite the diverse shapes and structures among different cell types. It has been proposed that durotaxis is critical for fine-tuning cell path-finding and wound healing [9, 10]. Also, there is increasing evidence showing that durotaxis is involved in cancer metastasis, since tumors are usually stiffer than the surrounding materials [11, 12].

A standard approach to modeling cell motility is to assume that cells execute a persistent random walk [13, 14, 15]; sometimes Lévy walks are used instead [16, 17]. Recently, Elizavata and colleagues applied persistent random walk ideas to understand durotaxis by relating persistence to substrate stiffness [18]. Their approach did show how this assumption could lead to durotaxis, but did not propose any direct mechanical reason for this correspondence; also they did not fully analyze their model in the relevant case of a two-dimensional spatial domain. In this study, we propose a simple intracellular mechanism that naturally leads to stiffness dependent persistence which, in agreement with the above findings, results in durotaxis. Our approach combines direct simulations with the derivation of a quantitatively accurate 2D Fokker-Planck equation, for which the numerical solution matches well with simulation data.

Our basic hypothesis is built on the fact that cells are observed to be more polarized when they move on stiffer materials. Cells have sophisticated mechanisms to sense stiffness, involving various cellular components and subsystems including the plasma membrane [19, 20], actin filaments [21, 22], actomyosin-based contractility, integrin-based focal adhesions [23, 24], etc. Once cells sense a stiffer substrate, they take on a more elongated shape [25, 26] as a response. Now, cells move by protrusion which occurs with the help of focal adhesions which allow force transmission to the substrate. We will assume that the change in shape to being more polarized implies that focal adhesions (FA) are formed within a narrower wedge on the cell front. It is also possible that the total number of FAs present at some fixed time increases on stiffer substrates, as FA are observed to be more stable on stiffer substrates [6]. In our model, both the distribution and the total number of FA directly control the variance of deflection angles in cell motion over a short time interval. This mechanism will create the necessary relationship between stiffness and persistence.

In experiments, the locations of cells moving on a 2D surface are typically recorded at fixed time intervals. Accordingly, we model the cell as a rigid object moving with velocity $v$ and rotating its motion direction $\Phi$ (its polarization) by an angle $\Delta\Phi$ at fixed time intervals $\Delta t=t_{i+1}-t_{i}$ which we take to be our unit of time. To determine $\Delta\Phi$, we assume there are a number $N_{f}$ of focal adhesions which are positioned at distances $r_{i}$ from the cell center and angles $\theta_{i}$ relative to the current direction of motion; these are chosen randomly from uniform distributions with ranges $(r_{min},r_{max})$, $(-\theta_{max},\theta_{max})$ respectively. We assume in line with the previous arguments that $\theta_{max}$ is determined by local substrate stiffness $k$ as $\theta_{max}=A/k$, where $A$ is a constant factor. The basic picture of our cell is given in Fig. 1, Our calculations will assume that $N_{f}$ remains constant. The driving force from each focal adhesion is assumed to have a constant magnitude and to point in the current moving direction. The net driving force is canceled by the friction acting on the cell, thereby determining the velocity. At each time step, the dynamical formation and disruption of FAs cause a possible imbalance in the driving torque. With fast relaxation, the cell will rotate by an angle $\Delta\Phi$ at each time step to satisfy zero net torque:

[TABLE]

whose solution is:

[TABLE]

For the purpose of illustration, we typically set $r_{i}=1$ for all $i$ and $N_{f}=12$ in our model.

Clearly the variance of the induced distribution for $\Delta\Phi$ determines the persistence of the motion. Here we use a Monte Carlo sampling method to evaluate this variance. We use $10^{6}$ sampling steps and have checked that this gives us an accurate evaluation for the range of parameters we have investigated. For the case of fixed radii, we obtain

[TABLE]

One can also compute the variance for the more general situation with a distribution for the radii as well. Typical results of this calculation are shown in Fig 2a. For use later on, we have fitted the data for the case $N_{f}=12$ with fixed radii to a simple function of $k$.

[TABLE]

As expected, increasing $N_{f}$ or decreasing $\theta_{max}$ reduces the variance. Thus, rigidity-dependent changes in the focal adhesion dynamics can indeed be used to model the mechanism underpinning the persistence-stiffness correlation.

Since focal adhesions are dynamically formed and destroyed, in our model at each time step the locations of all focal adhesions $\theta_{i}$ are reselected with no correlation to their previous value, hence $<\Delta\Phi(t_{i})\Delta\Phi(t_{i+1})>=0$. Thus on a uniform substrate, approximating the distribution of $\Delta\Phi$ to be Gaussian with the calculated width reduces our model to a version of the worm-like chain, where the mean squared displacement is:

[TABLE]

Here the persistence time is defined as $\tau_{p}=-\frac{1}{v\ln<\cos\Delta\Phi>}$. Since $<\cos\Delta\Phi>=e^{-\frac{\sigma^{2}}{2}}$, where $\sigma^{2}=var(\Phi)$, $\tau_{p}=\frac{2}{v\sigma^{2}}$. We simulate cell trajectories on uniform substrates with different stiffness and verify the previous results for $\sigma^{2}$ (see Fig 2b). In Fig. 2c and 2d, we show trajectories of cells simulated on both uniform soft and hard substrates. Consistent with the experimental observation [27], cells crawl more efficiently on stiffer substrates.

Next, we study the effect of stiffness gradients on cell motility. We set impose a constant stiffness gradient in the central region with constant low stiffness $k_{left}$ on the left side and high stiffness $k_{right}$ on the right side. We fix both $k_{left}$ and $k_{right}$ and vary the width of central region.

[TABLE]

Initially all our cells are placed at the origin and given a random initial direction. For small width, at a time when half of the cells go into the stiff region on the right, the other half are still hovering within the central gradient region (Fig. 3a). As the width increases, fewer and fewer cells enter the soft region on the left (Fig. 3b and c). This is caused by the fact that larger width allows more moving steps inside the gradient region and cells have more time to adapt to the direction of stiffness gradient. We further characterize these results by the durotactic index (DI) ([28]). We calculate DI defined below every ten time steps in our simulation:

[TABLE]

where $N_{right}$ and $N_{left}$ are the number of cells instantaneously moving to the right and to the left respectively. This index ranges from $[-1,1]$. Larger Durotaxis Indices indicate more cells are moving towards the ascending gradient direction.

We find that the curve can be divided into three sections (Fig. 3d). In the first section, nearly all cells are still in the gradient region. The index increases rapidly, which suggests that cells start being guided by the stiffness gradient. Consistent with experiment observation [26], the magnitude of DI is highly correlated with the magnitude of the gradient. In the second section, part of the cells are in the gradient region while the other have entered the uniform stiffness region. In the last section, the index starts decreasing because all cells move into the uniform rigid region and begin to execute random walks. The DI curve in the first section elucidates the role of gradient stiffness on cell motility.

To facilitate understanding of our simulation data, we now develop a Fokker-Plank equation for the probability distribution $P(x,y,\Phi;t)$ governing a population of particles in our model. We will be specifically interested in cases with a stiffness gradient, which we choose to lie along the $x$ direction. We focus on the variation with $x$ and $\Phi$ and introduce $p=\int dyP$ as a two-dimensional density. For any single cell, the next position $x(t+dt)$ depends on the current position and angle via $x(t)+v\cos(\Phi(t))dt$. We can therefore represent a single step in our stochastic process via

[TABLE]

where we will use the aforementioned Gaussian approximation

[TABLE]

Here $a(x)=\sqrt{\frac{1}{2\pi\sigma(x)^{2}\mathrm{d}t}}$ is the normalization coefficient as long as the width is significantly smaller than $2\pi$. Note that now the variance depends on $x$ through an $x$ dependence in the stiffness $k$.

In standard manner we can assume small $dt$ and expand $p$ around the current values of its arguments. After some simplification, we obtain

[TABLE]

The horizontal location $x$ and moving direction $\Phi$ are directly coupled in the last term on the right side, which is of the order of $dt$. We have checked that this third term can be neglected for the discrete update steps in out simulation, Consequently, Eq(Phenomenlogical Modeling of Durotaxis) can be simplified to

[TABLE]

This is of course the same as the Langevin equation for single cells:

[TABLE]

where $<\eta(t)>=0$,$<\eta(t)\eta(t^{\prime})>=\delta(t-t^{\prime})\sigma(x)^{2}$.

We then solve Eq (9) for both the uniform and stiffness gradient substrate cases. Fig 4a shows the initial condition in the all cases discussed in the following; in particular we apply a narrow gaussian distribution to approximate the $\delta$ in $p(x,\Phi;t=0)=\frac{1}{2\pi}\delta(x),\Phi\in(0,2\pi)$. The solution shows a peak in $x$ which varies from being at positive values (for $\Phi\simeq 0$) to negative ones (for $\Phi\simeq\pi$); the peak heights are independent of $\Phi$ as expected via rotational symmetry.

For the stiffness gradient case, the stiffness distribution is described by Eq (5) with $L=400$. The full distributions are shown at several times in Figs. 5a and 5b. Now there is a clear peak as a function of the direction; this can be seen most directly by plotting the integrated distribution $\tilde{p}(\Phi,t)=\intop_{-\infty}^{\infty}p(x,\Phi,t)dx$. Most cells adapt their moving directions from their initial directions $\Phi(t=0)$ to $\Phi$ near zero, exhibiting durotaxis (Fig 5c). We also define $\hat{p}(x,t)=\intop_{0}^{2\pi}p(x,\Phi,t)d\Phi$ and compare our PDE result for this quantity with direct simulations (Fig 5d). The very good agreement between PDE and direct simulation results validates the Fokker-Planck equation approach. Note that the above trends continue with the population continuing to break up into a peak at positive $x$ and a straggler peak at negative $x$ corresponding to cells that have wandered out into the uniform less stiff side of the gradient profile.

In our calculations so far, we fixed the number of focal adhesions and varied their distributions. In fact, FA’s are observed to be more stable on a stiffer substrate[6], which suggests a larger number of effective FA’s at any specific time. An alternative model would have a fixed FA distribution and vary the number of FA on substrates with different stiffness. In addition, it is possible that cells move with faster speed on stiffer substrates. Our model is capable of including such effects, all of which can contribute to durotaxis.

In this study, we discussed a possible underlying mechanism for durotaxis, namely a stiffness dependence of FA formation. It is known that FAs can dynamically sample rigidity to act as mechanosensors [24], but it remains elusive how FAs formation can directly control cell motility. We show that stiffness dependent FA formation causes a positive correlation between persistence time and substrate stiffness, which leads to durotaxis [18]. Also, we derive the corresponding 2D Fokker-Planck equation associated with our model and show that it gives consistent numerical agreement with our simulations. Our work can potentially help in predicting cell motility in more complex physiological environments such as those we arise during cancer metastasis.

Our model implicitly assumes that cells are incompetent in sensing rigidity gradients without moving around. For chemotaxis, a close analog of durotaxis, a eukaryotic cell is capable of comparing chemical concentration between its two ends, even though a typical bacterium bacteria is not [29]. It is technically hard to test such an ability in durotaxis, mainly because the cytoskeleton is essential for both cell motility and mechanosensing. Recently it has been shown that some cells can exhibit durotaxis as a cluster even if isolated constituent cells are ineffective [30]; in this case motion is not necessary.

Acknowledgements. This work was supported by the National Science Foundation Center for Theoretical Biological Physics (Grant NSF PHY-1427654). HL was also supported by the CPRIT Scholar program of the State of Texas.