Adhesion modulates cell morphology and migration within dense fibrous networks

Our simulation is based on the mathematical model for multicellular systems introduced in reference [81] and explored further in [82]. First we define a free energy functional F[ϕ, ψ] that will determine the system dynamics. Similarly to [56], we use two order parameters $\phi \left(\mathbf{r},t\right)$ and ψ(r, t). The cell is identified in space by the regions where ϕ ≈ 1 and the ECM is defined by the domain where ψ ≈ 1 (see figure 1). The interstitial space is defined as the region where both functions are close to zero. This free energy functional is the sum of a functional related to the cell's energy F_cell[ϕ], with another functional describing the cell's interaction with other biological entities F_int[ϕ, ψ], such as the ECM:

$$\begin{equation}F\left[\phi ,\psi \right]={F}_{\text{cell}}\left[\phi \left(\mathbf{r},t\right)\right]+{F}_{\mathrm{int}}\left[\phi \left(\mathbf{r},t\right),\psi \left(\mathbf{r},t\right)\right].\end{equation} \tag{ 1 }$$

Zoom In Zoom Out Reset image size

The first component of this functional, written in equation (2), is given by the Ginzburg–Landau free energy with a double well potential, with two stable solutions ϕ = {0, 1}, which correspond respectively to the exterior and interior of the cell (see figure 2),

$$\begin{align}\hfill {F}_{\text{cell}}\left[\phi \right]& ={\int }_{{\Omega}}\left[\frac{{\varepsilon }^{2}}{2}{\vert \nabla \phi \vert }^{2}+\frac{1}{4}{\phi }^{2}{\left(1-\phi \right)}^{2}\right]{\mathrm{d}}^{3}\mathbf{r}\hfill \\ \hfill & \quad +\frac{{\alpha }_{v}}{12}{\left({V}_{\text{target}}-V\left[\phi \right]\right)}^{2},\hfill \end{align} \tag{ 2 }$$

where Ω denotes the total volume of the system. The coefficient ɛ is a positive constant that is related to the width of the interface and to the surface tension σ (see SI section S2) (stacks.iop.org/JPhysCM/32/314001/mmedia). The squared gradient term acts as an energetic cost for maintaining the interface: more interface costs more energy. The coefficient α_v in the last term is the Lagrange multiplier for the volume $V\left[\phi \right]={\int }_{{\Omega}}h\left(\phi \right){\mathrm{d}}^{3}\mathbf{r}$, where h(ϕ) = ϕ²(3 − 2ϕ), α_v is the reinforcement coefficient for the volume constraint and V_target is the target volume of the cell: if $V\left[\phi \right]{ >}{V}_{\text{target}}$ this term removes material from the surface and when $V\left[\phi \right]{< }{V}_{\text{target}}$ it introduces material in the interface in order to maintain the volume constant. The function h(ϕ) reinforces the order parameter values to remain between the interval [0, 1], since h'(ϕ) = 0 for ϕ = 0 and ϕ = 1 (see reference [81] for a detailed description).

Zoom In Zoom Out Reset image size

The interaction component of the free energy is given by

$$\begin{equation}{F}_{\mathrm{int}}\left[\phi ,\psi \right]={\int }_{{\Omega}}\left[\eta \nabla h\left(\phi \right)\cdot \nabla \psi +\gamma h\left(\phi \right)h\left(\psi \right)\right]{\mathrm{d}}^{3}\mathbf{r}.\end{equation} \tag{ 3 }$$

The first term characterises the adhesion between the cell and the matrix (controlled by the parameter η), and the second penalises energetically the overlap between cells and matrix fibres (controlled by the parameter γ). Given the free energy functional, the equation for the cell dynamics is derived as a gradient descent equation with an advective term. Therefore, in this case we consider the ϕ change rate to be proportional to the gradient of the free energy functional [83, 84]

$$\begin{equation}{\partial }_{t}\phi +\mathbf{v}\cdot \nabla \phi =-\frac{\delta F}{\delta \phi }.\end{equation} \tag{ 4 }$$

Taking the free-energy functional derivative, we finally have the equation for the cell dynamics

$$\begin{align}\hfill {\partial }_{t}\phi +\mathbf{v}\cdot \nabla \phi & ={\varepsilon }^{2}\nabla \phi +\phi \left(1-\phi \right)\left(\phi -\frac{1}{2}-\gamma h\left(\psi \right)\right.\hfill \\ \hfill & \quad \left.+\eta {\nabla }^{2}\psi +{\alpha }_{V}\left({V}_{\text{target}}-V\left[\phi \right]\right)\right),\hfill \end{align} \tag{ 5 }$$

where the velocity v is written as the gradient of the chemical field c, which for metastatic cells can be assumed as the oxygen concentration in the tissue guiding the cell to a blood vessel during invasion:

$$\begin{equation}\mathbf{v}=\chi \nabla c,\end{equation} \tag{ 6 }$$

where χ is chemotactic parameter. In this study we assumed that this velocity is constant and directed along the x-direction, which allowed us to isolate the influence of the mechanical constraints in the movement. The details regarding the numerical implementation are in the supporting information (SI) and the software is available under the DOI 10.5281/zenodo.35342 [85, 86].

We propose a model for the system described in the previous sections using a DPD approach. Moreover, the DPD allows us to explore different perspectives that would be extremely complicated using only the PFM, such as the introduction of elastic fibres. Briefly stated, in DPD we solve Newton's equation of motion based on the total force that acts on each particle i [77, 87]:

$$\begin{equation}\sum _{i\ne j}{\mathbf{F}}_{ij}=\sum _{i\ne j}\left({\mathbf{F}}_{ij}^{\mathrm{C}}+{\mathbf{F}}_{ij}^{\mathrm{D}}+{\mathbf{F}}_{ij}^{\mathrm{R}}+{\mathbf{F}}_{ij}^{\mathrm{H}}\right).\end{equation} \tag{ 7 }$$

As usual in DPD models, each bead is interpreted through a coarse-grained perspective, where distinct numbers of particles are combined into a single bead. Each bead species has distinct properties, defined by their interaction potentials. The first term in this equation is the conservative force, which models the interaction between the different bead types, as function of the bead's types. Here, our conservative term is

$$\begin{equation}{\mathbf{F}}_{ij}^{\mathrm{C}}\left({r}_{ij}\right)={a}_{ij}{\omega }^{\mathrm{C}\mathrm{R}}\left({r}_{ij}\right){\hat{\mathbf{r}}}_{ij}-{b}_{ij}{\omega }^{\mathrm{C}\mathrm{A}}\left(r\right){\hat{\mathbf{r}}}_{ij},\end{equation} \tag{ 8 }$$

where the first term is the repulsive force, with a_ij being the maximum repulsion between the particles i and j, and the second term the attractive force with maximum intensity b_ij. In equation (8) r_ij = r_i − r_j is the distance between the two particles, and ω^CR(r), known as the weight function, usually assumes the form

$$\begin{equation}{\omega }^{\mathrm{C}\mathrm{R}}\left({r}_{ij}\right)=\begin{cases}\left(1-\frac{{r}_{ij}}{{r}_{1}}\right),\quad \hfill & \mathrm{i}\mathrm{f}\quad {r}_{ij}{< }{r}_{1}\hfill \\ 0,\quad \hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\;\hfill \end{cases},\end{equation} \tag{ 9 }$$

with r₁ = 1.25  μm. This term is responsible for the repulsion due the excluded volume of the beads. To model the adhesion of cells to the ECM fibre we added the second term of equation (8). This attractive term in the conservative force was proposed recently for the study of cell deformation in a surface [88]. The attractive weight function is given by:

$$\begin{equation}{\omega }^{\mathrm{C}\mathrm{A}}\left({r}_{ij}\right)=\begin{cases}\left(1-\frac{{r}_{ij}-{r}_{1}}{{r}_{2}-{r}_{1}}\right),\quad \hfill & \mathrm{i}\mathrm{f}\,{r}_{1}{< }{r}_{ij}{< }{r}_{2}\hfill \\ \quad 0,\quad \hfill & \text{otherwise}\hfill \end{cases},\end{equation} \tag{ 10 }$$

where r₂ = 2.50  μm. This term acts only near the fibre's surface and describes the biochemical adhesion forces between the cell and the fibres.

The dissipative force, ${\mathbf{F}}_{ij}^{D}$ in equation (7) represents the effect of viscosity. It uses the relative velocity between two particles to slow down the motion with respect to each other and is given by

$$\begin{equation}{\mathbf{F}}_{ij}^{\mathrm{D}}=-\gamma {\omega }^{\mathrm{D}}\left({r}_{ij}\right)\left({\hat{\mathbf{r}}}_{\mathbf{ij}}\cdot {\mathbf{v}}_{ij}\right){\hat{\mathbf{r}}}_{\mathbf{ij}},\end{equation} \tag{ 11 }$$

where γ is a coefficient, v_ij = v_i − v_j is the relative velocity between the particles and ω^D(r_ij) is a weight function that depends on the distance between the particles:

$$\begin{equation}{\omega }^{\mathrm{D}}\left({r}_{ij}\right)=\begin{cases}\left(1-\frac{{r}_{ij}}{{r}_{2}}\right),\quad \hfill & \mathrm{i}\mathrm{f}\quad {r}_{ij}{< }{r}_{2}\hfill \\ 0,\quad \hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\;\hfill \end{cases}.\end{equation} \tag{ 12 }$$

Finally, the random force ${\mathbf{F}}_{ij}^{R}$ represents the thermal or vibrational energy of the system and includes fluctuation effects,

$$\begin{equation}{\mathbf{F}}_{ij}^{\mathrm{R}}=\sigma {\omega }^{\mathrm{R}}\left({r}_{ij}\right){\xi }_{ij}{\hat{\mathbf{r}}}_{ij},\end{equation} \tag{ 13 }$$

where σ is a coefficient, ω^R(r_ij) the weight function and ξ_ij a random constant obtained from a Gaussian distribution.

The dissipative and random forces are related by their coefficients and weight functions,

$$\begin{equation}{\omega }^{\mathrm{D}}\left(r\right)={\left[{\omega }^{R}\left(r\right)\right]}^{2}\;,\quad {\sigma }^{2}=\frac{2\gamma {k}_{\text{B}}T}{m},\end{equation} \tag{ 14 }$$

where k_B is the Boltzmann constant and T the equilibrium temperature. These relations turn the DPD into a thermostat. Since the algorithm depends on the relative velocities between the beads and since all interactions are spherically symmetric, the DPD is an isotropic Galilean invariant thermostat which preserves the hydrodynamics [89].

Here, we propose a model to study the relation between cell shape and dynamics. Our cell consists of 452 beads distributed around a sphere of radius a_c = 6 μm, the cell radius. Cell surface beads c repel themselves by the conservative force with an intensity of a_cc = 12.5 × 10⁻¹¹ N. These surface cell beads are attached by a standard harmonic force F^H to a central ghost bead g, as depicted in figure 3(a):

$$\begin{equation}{\mathbf{F}}_{\mathrm{c}\mathrm{g}}^{\mathrm{H}}=-{k}_{\mathrm{c}\mathrm{g}}\left({r}_{\mathrm{c}\mathrm{g}}-{a}_{\mathrm{c}}\right){\hat{\mathbf{r}}}_{\mathbf{cg}},\end{equation} \tag{ 15 }$$

where r_cg is the distance between a c bead and the central bead g, and k_cg = 30 pN μm⁻¹. Figure 3(b) shows a snapshot of the cell in a bulk simulation with small surface beads so the ghost bead can be seen. In figure 3(c) we show the same snapshot but with the surface beads depicted in their regular size. As we can see, the cell remains approximately spherical—obviously less spherical than the raspberry model [90]—and, as we will demonstrate, it can deform to pass through obstacles.

Zoom In Zoom Out Reset image size

The other two species in the DPD system are the solvent and the fibre matrix. Solvent beads s have a number density ρ_N = 4.0, consistent with previous works [91, 92], and are repelled by c beads with a strength a_cs = 65 × 10⁻¹¹ N with a cut off radius of 2.5 μm. The fibres are modelled as polymers with beads of type f, and the interaction parameters between c and f beads are the same as between c and s, a_cs = a_cf. The only attraction interaction in the system is between c and f beads. In this case, the attractive parameter plays the role of the adhesion defined by the relation b_cf = η_DPDa_cf. Here η_DPD provides the strength of the adhesion between the cell and the polymer network.

In the case of a rigid fibre network, the polymers are placed randomly inside the simulation box and remain fixed, i.e., the corresponding beads are not integrated in time. In the case of flexible fibre network only the beads placed at the ends of the polymeric chains are fixed, and the rest of the polymer can oscillate similarly to a rope with fixed ends. In this flexible case, the beads in the same polymer are bound by the harmonic bond, equation (15), with k_ff = k_cg = 30 pN μm⁻¹. The case of rigid polymers can be imagined as the limit case of k_ff → ∞. This description is motivated by the fibrous structures of collagen present in the ECM, which may have different stiffness depending on the collagen concentration, its origin (skin, tendon, vasculature, organs, bone) and/or polymerization process of collagen-based scaffolds for in vitro experiments [93–95].

The cell is placed in the left buffer zone of the simulation box and a constant force in the longitudinal direction is applied to the cell mimicking the chemical gradient. Then we analyse the system properties during the cell migration in the fibre network region. The simulation ends when the cell reaches the right buffer. The SI presents a detailed description of the setup. We show in figure 4 a schematic depiction of the simulation box.

Zoom In Zoom Out Reset image size

For the cell to overcome the obstacles in dense fibre networks it needs to deform. In figure 5 we show some snapshots from the cell migration. Row (A) refers to cell modelled with DPD when the fibres are rigid, while row (B) shows examples for the flexible matrix scenario. For better visualisation, the fibres are smaller than the cell beads, and the solvent beads are not shown. Comparing (A.I) and (B.I) snapshots, we can see that even for a simple obstacle, which is to bypass a single fibre, the cell is more deformed by the rigid fibre. The deformation becomes larger when the cell is blocked by more than one rigid fibre, as the (A.II) and (A.III) snapshots show. On the other hand, the cell is only slightly deformed by the flexible matrix, as the (B.II) and (B.III) snapshots show. We expect that the time needed by the cell to escape from a fibre trap is directly related with the time necessary to deform its spherical shape and cross through the available space. The PFM presents similar results to the DPD with rigid fibres, which we can see when comparing the rows (A) and (C) in figure 5.

Zoom In Zoom Out Reset image size

To understand the velocity dependence on the fibre density ρ for rigid and flexible networks we present in figure 6 typical examples of the x component of the cell center of mass (CM) of the DPD cell as function of time for three different densities ρ = {0.20, 0.55, 0.70}. We present similar results for the PFM in figure 2 of section 3 in SI.

Zoom In Zoom Out Reset image size

In the examples for the lowest fibre density, the cell may be blocked a few times for both rigid and flexible fibre matrices, but it diffuses almost without obstacles for most of the time. In these low density matrices, the time required for the cell to leave the rigid or flexible network is similar. However, as ρ increases the scenario changes. For ρ = 0.55, figure 6(b), we can see that both curves have two trap regions, but the time spent by the cell to leave those regions is much longer in the rigid network. The difference in the stop time for the rigid and flexible matrices becomes even more pronounced for ρ = 0.70, figure 6(c).

Another interesting feature in these systems is the diffusion regime. We quantify the diffusion by scaling the MSD and the time, namely 〈|r(t)−r(0)|²〉 ≈ t^β. If the MSD increases linearly with time, β = 1.0, the system is in the normal or Fickian regime. If β < 1.0 the diffusion is anomalous and called subdiffusive, while in the case of β > 1.0 the regime is superdiffusive. The limit of β = 2.0 is called ballistic. To understand how the rigid or flexible matrix can influence the diffusion regime we evaluated the exponent β as function of the fibre density ρ. The value of β used here is the obtained as t → ∞, and not the short time exponent.

We can observe in figure 7 that a superdiffusive regime was observed at lower values of ρ for the three models, in agreement with experimental studies [96]. However, at higher network densities, ρ ⩾ 0.50, we observe a change in the regime for the DPD model with rigid fibres and for the PFM. There is a transition from a superdiffusive to a subdiffusive regime. This transition can be directly related to the longer times the cell gets blocked the fibre matrix. At low densities, since there is force pulling the cell, the movement is almost ballistic, with occasional collisions with the matrix. However, as the density increases, the cell MSD is strongly affected, with long periods of no mobility. On the other hand, the migration inside flexible matrix is only slightly affected. The movement changes from a near-ballistic to a superdiffusive regime. Here, the change is small and even at high densities the time that the cell remains blocked is small, as we have showed in figure 6(c). Therefore, the flexibility of the fibre network plays an important role not only in the velocity of migration, but also in cell deformation and in diffusion regime.

Zoom In Zoom Out Reset image size

We then explored how adhesion affects cell migration in both PFM and DPD models. The mean velocity is measured as the quotient of the distance traveled by the center of mass of the cell and its respective time interval. In figure 8 we present how the fibre density ρ affects migration for three different values of adhesion η. The overall behaviour is similar and independent of the adhesion coefficient: for lower fibre densities we obtain the highest velocities, while as the density increases, the velocity becomes smaller the cell gets trapped in the matrix. Nonetheless, by comparing the migration velocity for different values of adhesion at the higher densities ρ = {0.6, 0.7}, we observe that the velocity is higher for η = 1 than for the lower values η = {0.0, 0.5}. This dynamic is anticipated by experimental results of cell migration, but it is striking to observe this behaviour in this simple model.

Zoom In Zoom Out Reset image size

For further investigation concerning the dependence of the migration velocity on the adhesion we fixed the density of fibres at ρ = 0.6 and ran simulations for several values of adhesion. In figure 9 the mean velocity is plotted as a function of the adhesion coefficient for our computational models and it includes the in vitro results adapted from [97] for fibroblast migration as a function of the fibronectin concentration. We can see that the cell migration velocity increases with the adhesion coefficient until it reaches a maximum, after which it decreases. This suggests that there is an optimal value of adhesion for migratory events [98, 99]. Indeed, we have superposed the experimental data for the cell velocity as a function of the fibronectin concentration obtained by Maheshwari and colleagues [97] and the same qualitative behaviour is verified (see figure 9). Fibronectin appears as protagonist in cell adhesion and its concentration is a proxy for altering adhesiveness—more fibronectin present in the ECM leads to higher adhesiveness.

Zoom In Zoom Out Reset image size

Experimentally the adhesiveness is linked to the cell capability to exert traction forces in the fibres in order to pull its body forward. Thus this is the main reason by which the in vitro assay exhibits an increase in velocity with the adhesiveness and then decreases after reaching a critical point where adhesion is so strong that it suppresses migration. However, in our mathematical model we have a slightly different circumstance that leads to the same qualitative results. The first main difference we should stress here is that adhesion is a passive attribute while in nature adhesion is an active dynamical event performed by the cell. The migration of our droplet is affected by adhesion in two different manners: (i) changing cell plasticity and (ii) introducing an artificial pulling force. A higher adhesion increases the cell deformity (see figure 10) thus facilitating the migration through small pores. Moreover, the adhesion acts as a pulling force due to the leading edge of the cell that feels the fibres upfront and is attracted to them by the adhesive force. This effect mimics some of the features observed in real cells during migration.

Zoom In Zoom Out Reset image size

The study realised by Wolf et al used an in vitro model of tumour cells and immune cells (T-cells and neutrophils) migrating in 3D collagen scaffolds of different porosities to study the physical limits of migration [9]. In order to compare our results to their experimental data, calculations were performed to determinate the mean pore cross-section of the fibre networks. Each matrix was sliced in several planes orthogonal to the migration direction and for each plane the connected component algorithm (CCA) was applied to obtain the area of the pores (see section 5 of SI). Then the procedure was repeated for several random generated matrices of the same density and the mean values and their respective unbiased standard deviation were obtained. On the other hand, for the DPD matrices we have used a method described in [100] that estimates the area of the pores by the effective space available [101]. The pore sizes is inversely proportional to the density of fibres. The computational data is depicted in figure 11 together with experimental results adapted from [9] (figure 2(b)). The experimental data was measured for human fibrosarcoma cells of the line HT1080 migrating in rat tail collagen matrices treated with an MMP inhibitor (GM6001) and for untreated matrices (with medium).

Zoom In Zoom Out Reset image size

As shown in figure 11 the cell velocity increases linearly with the pore size for both experimental and computational models. For small pore sizes the cells start to have MMP-dependent migration, degrading the fibres and opening space for migration. In our model the cell becomes trapped for higher pore sizes than in the experimental setup due to the lack of this mechanism (figure 11). For the DPD model with flexible fibres (figure 11(b)) we observe a higher velocity than in the rigid case, which is expected since matrix elasticity favours migration.

In the experimental assay the authors had also repeated the migration study isolating the role of MMPs. For this purpose, the cells were treated with the potent broad spectrum MMP inhibitor GM6001, also known as ilomastat, which mitigates the influence of MMP activity during migration. In figure 12 we present the migration velocity as a function of the pore cross section for (a) PFM and (b) DPD results and compare with the GM6001 experimental data. As before, we can observe for the three cases a linear dependency with the pore cross section.

Zoom In Zoom Out Reset image size

Comparing the two lines in figure 12(a) we can notice that one is shifted in relation to another, as well when comparing the DPD rigid in figure 12(b). This difference is the result of a larger plasticity of the cancer cells used in the experiments and their capability of adhering to the fibres, which facilitates migration, when compared with the capacity of the simulated cells to deform. The fibrosarcoma cells are able to migrate under extreme confinement conditions as seen in the transwell filter assay where cells can pass through apertures of only 3 μm wide (figure 7(b) of [9]). The fact that the line presented in figure 12(b) for the DPD flexible is shifted to the left corroborate this hypothesis and reveals some level of equivalence between matrix elasticity and cell plasticity.

We then increased the adhesiveness in the PFM and ran simulations for η = 0.5 and η = 1.0 in order to evaluate the contribution of adhesion to migration. We can see in figure 13 that increasing adhesion from η = 0 to η = 0.5 the velocity dependency with the pore size rotates, becomes more similar to the result obtained experimentally with MMP inhibition (RT GM6001). When increasing even further the adhesiveness (η = 1.0) the PFM results became analogous to the experimental results without the inhibitor (RT Medium). These results indicate that MMP inhibitors are able to affect significantly adhesiveness by suppressing the cleavage of molecules that act as adhesion promoters. Moreover this brings adhesiveness as a key modulator of migration in crowded media. In addition, we suggest that more experimental assays should be carried to unveil how MMPs activity acts in different mechanisms of cell migration, using specific inhibitors to target isolated members of the MMP family. Some MMP isoforms are known to affect mainly collagen degradation and remodelling while others act regulation cell adhesion and migration.

Zoom In Zoom Out Reset image size