Research articleA biochemo-mechano coupled, computational model combining membrane transport and pericellular proteolysis in tissue mechanics
Abstract
We present a computational model for the interaction of surface- and volume-bound scalar transport and reaction processes with a deformable porous medium. The application in mind is pericellular proteolysis, i.e. the dissolution of the solid phase of the extracellular matrix (ECM) as a response to the activation of certain chemical species at the cell membrane and in the vicinity of the cell. A poroelastic medium model represents the extra cellular scaffold and the interstitial fluid flow, while a surface-bound transport model accounts for the diffusion and reaction of membrane-bound chemical species. By further modelling the volume-bound transport, we consider the advection, diffusion and reaction of sequestered chemical species within the extracellular scaffold. The chemo-mechanical coupling is established by introducing a continuum formulation for the interplay of reaction rates and the mechanical state of the ECM. It is based on known experimental insights and theoretical work on the thermodynamics of porous media and degradation kinetics of collagen fibres on the one hand and a damage-like effect of the fibre dissolution on the mechanical integrity of the ECM on the other hand. The resulting system of partial differential equations is solved via the finite-element method. To the best of our knowledge, it is the first computational model including contemporaneously the coupling between (i) advection–diffusion–reaction processes, (ii) interstitial flow and deformation of a porous medium, and (iii) the chemo-mechanical interaction impelled by the dissolution of the ECM. Our numerical examples show good agreement with experimental data. Furthermore, we outline the capability of the methodology to extend existing numerical approaches towards a more comprehensive model for cellular biochemo-mechanics.
1. Introduction
A deeper understanding of the mechanics of living cells interacting with the extracellular matrix (ECM) is essential for various physiological and pathological processes, e.g. embryonic growth, wound healing, inflammatory response, angiogenesis and tumour metastasis. Despite this outstanding importance, the detailed processes involved are not yet fully understood. Especially, the influence of the ECM for regulating cell migration is still topic of current research [1–3]. Therefore, in our overall goal, we aim to establish a comprehensive computational model enabling the study of individual cell migration with a specific focus on the influence of physical factors, such as matrix stiffness and interstitial flow.
In this publication, we address advection–diffusion–reaction processes within the ECM and on the cell membrane and its coupling to the properties of the mechanical system. A vast number of chemical species and processes interplay on the cellular level [4,5]. During proteolysis, the ECM is dissolved due to biochemical reactions initiated by the cell. It is widely known from experimental data that the mechanical state of the ECM, i.e. the stress and strain state, influences the dissolution rate [6–8]. ECM fibrils often show significantly different reaction kinetics when exposed to loading or strain. However, most mathematical models do not yet consider this specific interaction of mechanics and biochemistry. We propose a fully coupled model for such transport and dissolution processes within a deforming ECM including interstitial flow.
The ECM is represented by a two-phase poroelastic medium model, which is similar to porous medium approaches used in other biomechanical applications (e.g. [9–11]). This homogenized approach is able to represent the mechanical resistance of the ECM as well as the influence of interstitial flow. The solid phase of the porous medium is modelled by a hyperelastic constitutive law. The flow is described by the Brinkman equation. Therein, considering a variable porosity and its impact on the geometrical structure and mechanical behaviour of the ECM is of major importance when aiming for a realistic model, as the porosity is altered significantly during cell migration either through pure mechanical deformation or through biochemical proteolysis. Proteolysis is modelled as a dissolution process of the porous medium. Models for dissolution of porous media are known and often applied in soil mechanics [12,13]. We will adapt those concepts for our degradation model. The kinetics of the chemical reactions are linked to the mechanical state of the collagen matrix via a change of chemical potential. The cell membrane is represented by a diffusion–reaction model on a curved surface. We also consider the coupling between species within the ECM volume and on the cell membrane due to chemical reactions. In recent years, transport models for cell-interior mechanisms have been developed [2,14,15] and also advection–diffusion models for tissue engineering [16,17] and surface-volume coupled models for ECM degradation [18–20] have been published. Therein, detailed models for the expression of enzymes at the cell membrane and the advection–diffusion–reaction processes leading to proteolysis are given. For a review of computational approaches for modelling cellular systems to study anti-angiogenic therapeutics in cancer the interested reader is referred to [21]. However, despite the growing interest in mathematical models for cellular processes, the influence of the matrix deformation in this context is rarely considered. Very advanced approaches can be found in [22,23]. Therein coupled models of poroelasticity and reactive transport are presented with application to tissue mechanics. In addition, interaction between mechanical and biochemical processes are included, for instance, for modelling growth. We use a similar model extended by a formulation for surface-volume reactions and matrix degradation specifically designed for proteolysis in cell mechanics. Recently, an extensive numerical model for cell invasion dynamics into an ECM was presented in [24]. In contrast to our approach, the fibre network was resolved explicitly there. The degradability is modelled by cross-linking and uncross-linking events, i.e. the varying load bearing behaviour of the dissolving matrix is also considered. It further includes similar transport equations like our model and the references given before, although interstitial flow, i.e. advective transport, is not yet taken into account in this published model. And as a most important difference to the model presented in this contribution, in [24] the reaction kinetics are assumed to be independent of the mechanical loading. We present a model which is capable of resolving this missing mechano-chemo coupling. To the best of our knowledge, this is the first time a reactive surface-volume transport model with interstitial flow, extended by including a deformable macroscopic solid phase as well as dissolution processes depending on the structure’s thermodynamical state and its influences on the mechanical characteristics of the ECM, is proposed.
The outline of this publication is as follows. In §2, the model formulation will be presented. This includes the advection–diffusion–reaction equations for scalar quantities within the porous medium and on a curved surface as well as the model for the mechano-chemo coupling, i.e. the dissolution process and the deformation-dependent reaction kinetics. Two numerical examples are analysed in §3. The first example is a reproduction of an experiment, where the degradation kinetics of strained collagen networks is analysed. The second example extends an existing numerical model for transport and reactions within the ECM by the mechanical coupling of the dissolution process. Finally, a conclusion is drawn in §4.
2. Model and governing equations
(a) Extracellular matrix as porous medium
Simply spoken, the ECM is composed of a scaffold of mainly collagen fibres and interstitial fluid. Thus, from a continuum mechanical, macroscopic point of view, it can be described as a porous medium. In this publication, we model incompressible flow through a saturated, deformable porous medium. In particular, we consider large deformations, varying porosities and nonlinear material behaviour. A detailed derivation and further theory can be found in [25,26]. The governing equations are reviewed in the electronic supplementary material, section S1. In the main part of this publication, we only present the constitutive equation, as they are important for the dissolution model given later. It can be derived [26] as
(b) Proteolysis as advection–diffusion–reaction process within extracellular matrix
(i) Scalar transport within the extracellular matrix volume
For describing scalar transport within a porous ECM, the governing equations on the microscale are homogenized to the macroscale [27]. For a scalar transported within the fluid, one obtains
(ii) Transport processes on curved cell boundaries
For the diffusion equation, the mathematical theory in this case is based on the Laplace–Beltrami operator, being the equivalent of the Laplace operator on curved surfaces. Here, the focus lies on the integration of the transport equations into a coupled problem. Numerical analysis for transport processes on curved surfaces have been performed extensively in recent years (e.g. [28,29]). Also in the context of biomechanics and cell modelling, computational transport models on curved surfaces representing the cell membrane have been used extensively [20,30,31]. In this publication, we will only restate the equations in a matter-of-fact manner without going into detail about the finite-element model.
The advection–diffusion–reaction equations for N scalars ci, i=1,…,N on a curved surface reads
(iii) Chemical reactions
For modelling chemical processes in cell mechanics, there is the need of being able to incorporate a large number of reaction equations into the formulation. The reactions, more precisely the reaction kinetics, determine the form of the reaction terms σi(c) in the transport equations (2.2)–(2.4).
The coupling of the surface and the volume transport problem is established by reactions on the cell membrane, which include both chemical species bound to the membrane and chemical species transported within the ECM. One has to keep in mind, that a reaction between surface- and volume-bound species is a heterogeneous chemical reaction. Therefore, in those cases a characteristic length λ is used to scale surface concentrations to volume concentrations [18]. The specific reaction formulations will be elucidated in the respective examples.
(iv) Reactive dissolving porosity model
The idea is to formulate a model in which the porosity will change due to chemical reactions. This allows to reproduce the dissolution process of the ECM scaffold. The basic theory of dissolving porosity models presented here can be found in more detail in [32]. For more extensive reviews about general thermodynamics in porous media, we refer to [33,34].
We introduce the internal variable Δϕsol representing the additional change of reference porosity (and therefore solid mass) due to the dissolution process. The Hemholtz free energy of the structure Ψs(E,Δϕsol) thus also depends on Δϕsol (we neglect the dependency of the free energy on the temperature for simplicity here). We propose the following form for the Helmholtz free energy by adding the additional factor (1−ϕsol/(1−ϕ0)):
Now, we will focus on the chemo-mechanical coupling specifically for proteolysis of the ECM. In order to establish a link to the chemical model, a relationship between the structural density and the collagen concentration is needed. Providing that collagen is the dominant constituent of the ECM, the reference porosity is related to the molar concentration [C1] of collagen type I by
While the detailed mechanisms of collagen–enzyme interaction are still unclear, the effect observed in experiments is commonly an increasing or decreasing reaction rate. In this publication, we will use the following law for a decreasing reaction rate coefficient k:
3. Numerical examples
We use the finite-element method for solving all involved fields. The numerical approach is based on schemes presented in [25,26,42]. Detailed explanations can be found in the electronic supplementary material, section S2, and a validation example of the numerical scheme for surface-volume reactions in the electronic supplementary material, section S3.
(a) Dissolution of strained and unstrained extracellular matrix
This example is motivated by the experimental studies in [6]. There, the strain dependency of the dissolution process was examined. A collagen micronetwork was pre-strained by displacing two micropipettes and exposed to bacterial collagenase (Clostridium histolyticum). The amount of collagen fibres was estimated from optical images. An edge-intensity image processing method was used in order to obtain a measure for the collagen concentration. Representative results of the experiments are depicted in the electronic supplementary material, figure S2. It is stated there that the elapsed time between 5% and 90% edge loss is ΔTS = 1113 s ± 260 s for strained fibrils and ΔTU = 654 s ± 149 s for unstrained fibrils.
We will compare those results with our computational model. Four chemical species are considered for the transport model: collagen type I C1 composing the ECM, the degrading enzyme E, the collagen-enzyme complex C1E and the denatured collagen C1*. The chemical reactions involved can be schematically represented by a Michaelis–Menten mechanism, see e.g. [7],
The results for the Green–Lagrange strain and the Cauchy stress are well in the range, which were estimated in [6], see the electronic supplementary material, section S4, for details. The resulting change of chemical potential is illustrated in figure 2a. As expected, it shows a similar spatial distribution as the strain and stress field. The evolution over time of the collagen I concentrations at three different points are shown in figure 2b. The point A is placed within the strained region, while the points B and C are located in a certain distance of the centre of the domain (figure 1). As expected from the model formulation, the degradation process is slowed down significantly within the strained region. The point of time, where the degradation process starts, is different for each spatial point observed due to the diffusion time of the enzyme. As the point A is the closest to the pipettes the collagen concentration decreases there first, only after a delay there is significant change of concentration in B and C. The degradation dynamics in point B and C do not differ significantly, as the strained region is very localized to the area between the pipettes. Thus, the main difference between time evolution at these two points is due to the diffusion delay. The kinetics fit well to the experimental results from [6] in the electronic supplementary material, figure S2, with respect to the qualitative behaviour as well as to the time scales of collagen dissolution. Note, that only the unstrained reaction rate, i.e. at point C, has been fitted, the deceleration of the kinetics in point A follows directly from the physical parameters and the proposed model without any further fitting.
Figure 1. Dissolution example: dimensions of computational domain. The whole domain is a square region with 200 μm edge length. The pipette domains are two square regions with 20 μm edge length. The points A (centre of domain), B and C denote the points where the collagen concentration will be evaluated. Figure 2. Dissolution example: solution field for chemical potential and evolution of C1-concentration over time. The point A is located in the strained region between the pipettes, the points B and C further away, see figure 1. (a) Chemical potential (dm2 s−2) and (b) C1-concentration over time.
However, one has to keep in mind that there are still some uncertainties in the model. The material of the solid phase is assumed to be of neo-Hookean type, which is a simple elastic nonlinear constitutive law. It is known, that collagen fibres do not only exhibit nonlinear elastic, but also viscoelastic behaviour [6,44]. As the material model primarily determines the change of reaction kinetics, the use of a well-fitted viscoelastic model could improve the results. This could induce additional stabilization of the collagen at early time instances, as the viscous stresses are expected to be high in the beginning before relaxing over time. Actually, one could hypothesize that this produces an effect which can be observed when looking at the early behaviour of the collagen concentration in the experiments in [6] and the electronic supplementary material, figure S2. The collagen concentration is higher in the strained region than in the unstrained region at all times, whereas in our computational results, the collagen in point A is smaller than in the unstrained points at early times, because it is closer to the enzyme source and therefore more quickly subjected to proteolysis. Not until some diffusion time, the stabilizing effect starts to be visible in the concentration being higher in point A than in point B and C. However, the experimental results could also have other reasons like local matrix compaction or that unpolymerized collagen is attacked first [6]. Still, we are confident of our model as it gives good results without the need of unphysical parameters or assumptions. Besides, it is formulated in a general way based on the internal energy, such that other, possibly more suitable, material models can be included, further improving the results without changing the methodology.
Computational models can grant insights into the sensitivity of systems with respect to many parameters, which are not easily accessible for experimental set-ups. In the following, we will study the influence of certain parameters on the degradation time. The displacements and therefore the resulting strain, which were applied in the experiment in [6] are not easily controlled, however, in our simulation it is a given input parameter. Also, the mechanical material behaviour of collagen matrices is hard to measure. The degradation enzyme concentration can be controlled in experiments, but the dynamics and the diffusion processes might be of interest. Furthermore, the density of the ECM fibres can only be estimated. In figure 3, the applied displacement, the stiffness, the degradation enzyme concentration and the structural density are varied and their influence on the degradation time in strained and unstrained regions is evaluated. If one quantity is varied, the others are fixed at the respective value given in the electronic supplementary material, table S1. In [6] the time between 95% and 10% of initial collagen concentration was used as a measure for the degradation time. As this value is not always reached in our simulations for every set-up with every parameter setting, subsequently we will define the degradation time as the time span until 80% of the collagen is degraded, starting from the point in time when 95% of the initial collagen concentration is present. In general, all four plots in figure 3 show expected behaviour. Higher displacements and hence higher strains stabilize the collagen within the ECM (figure 3a). If the applied strain is low, the degradation time in strained regions is even lower than in the unstrained region, which is due to the diffusion delay. Concerning collagen matrix stiffness (figure 3b), the degradation time increases with a stiffer matrix. The reason for this is that under the same strain, stiffer materials exhibit a higher mechanical energy state, which results in higher chemical potential and stability. Similarly, a lower structural mass density leads to a higher chemical potential under the same loading state, also inducing higher degradation times (figure 3c). These three parameters are material characteristics of the mechanical behaviour of the ECM. Those parameters do not influence the degradation kinetics of the unstrained region. Clearly, if the enzyme concentration is varied the degradation times in the whole domain change accordingly (figure 3d). Still, the stabilizing effect of strain is visible for every analysed setting. The effect of different ECM porosities, i.e. different matrix densities, is analysed in figure 4. To obtain a similar setting, for all simulations in which the porosity was altered, the initial collagen concentration, the stiffness and the maximum enzyme concentration was adapted accordingly. With rising solid volume, the macroscopic stiffness rises and the amount of degradation enzyme needs to be increased as well to assure the same ratio between collagen and degradation enzyme. For instance, Young’s modulus and the maximum enzyme concentration were increased by a factor of 10 compared with the former values, if the porosity was changed from 0.99 to 0.9 as this corresponds to a change of solid volume fraction from 0.01 to 0.1. Note that this is a rough approximation concerning the stiffness, as the constitutive behaviour is nonlinear. Those adaptions lead to the results depicted in figure 4a. There, the stabilizing effect of the applied strain reduces with higher solid volume fraction. However, if the reaction coefficient is decreased in the same way as the solid volume fraction is increased, the stabilization of collagen is evident for all considered porosities, see figure 4b. If the latter results could be verified in experiments, this may indicate that there is an additional dependency of the reaction rates ( or others) on the porosity, which is not yet included in the model. For instance, it might be necessary to incorporate the change of collagen fibre surface, where the degradation takes place, due to a changing collagen fibre volume. In a volume averaged sense, this would lead to an altered effective reaction rate.
Figure 3. Dissolution example: parameter study. The degradation time within an unstrained and a strained region (point A and C, see figure 1) is compared for varying model parameters. Degradation time is defined as time span between 95% and 20% of local initial collagen concentration. (a) Varied displacement , (b) Varied stiffness Es, (c) Varied density ρs0 and (d) Varied maximal enzyme concentration [E]max. Figure 4. Dissolution example: parameter study. The degradation time at an unstrained and a strained region (point A and C, see figure 1) is compared for varying porosity. In the left subfigure, Young’s modulus and the maximum enzyme concentration were changed by the same factor as the solid volume fraction 1−ϕ. In the right subfigure also the reaction rate coefficient was adapted. Degradation time is defined as time span between 95% and 20% of local initial collagen concentration. (a) Varied porosity ϕ and (b) Varied porosity ϕ with adapted reaction rate coefficient .
(b) Cell proteolysis within extracellular matrix
The second example is motivated by the numerical example presented in [18]. Therein, a complex transport model coupled with Brinkman flow for ECM proteolysis is solved. We will extend this example by explicitly resolving the cell membrane including diffusion on the membrane and production of MT1-MMP and other chemical species localized at the root of the pseudopod [24]. Also, we incorporate a deformable, dissolving structural phase into the model, i.e. the collagen matrix, and the chemo-mechanical coupling, this time induced by a pulling force applied on the ECM by the cell. The geometry of the domain is depicted in figure 5. Details on the simulation set-up can be found in the electronic supplementary material, section S5. The cell itself is not modelled explicitly, only the transport processes on the cell membrane are considered. We include seven species transported within the ECM and five chemical species which are bound to the cell surface, similar to [18]. Additionally, the diffusion processes of the respective scalar on the membrane are accounted for. The scalars within the ECM are collagen type I C1, the matrix metalloprotease MMP2, the proenzyme pro-MMP2, the inhibitor TIMP2, the complex MMP2⋅TIMP2, the denatured collagen C1* and the degradation product MMP2⋅C1. On the cell membrane, the modelled chemical species are the degradation enzyme MT1-MMP, its complex MT1-MMP⋅TIMP2* with inactivated TIMP2 and the complexes MT1⋅TIMP2⋅MMP2, MT1⋅TIMP2⋅MMP2⋅MT1 and MT1⋅C1. An overview over the modelled species is given in the electronic supplementary material, table S2. In this table, short notations are introduced that will be used in the following.
Figure 5. Cell example: dimensions of computational domain. The whole domain is a rectangular region (30 μm high, 40 μm width). The cell geometry is composed by the cell body and a protrusion. At the root to the protrusion, MT1-MMP, TIMP2 and pro-MMP2 are produced [24]. The cell is 10 μm high and 15.66 μm long. At the leading edge, a pulling force is applied. (Online version in colour.)
In the following, the modelled reaction system is described. Both enzymes M1 and M2 exert proteolytic activity, i.e. they can degrade the collagen matrix. M1 is acting on the cell surface, while M2 is diffusing within the ECM. The inhibitor T2 can bind both degradation enzymes M2 and M1. The inhibition of M1 is at the same time the first step towards the activation of pM2, which is secreted by the cell at the cell membrane. Here, M1 and M1T2 act as a catalyst for the production of M2. The corresponding reaction equations and reaction terms are given the electronic supplementary material, section S5, and in equations (S20)–(S39).
The matrix of the ECM is modelled in a similar way as in the previous example. The same neo-Hookean material model (S17) is used for the structural phase The degradation model for the ECM is included as in the first example, i.e. the collagen degradation depends on the chemical potential as follows
We will compare the results obtained by degradation of the ECM with and without the presence of a pulling force applied by the cell. Regarding flow and deformation within the ECM, the results are nearly indistinguishable for those two cases. They are described in detail in the electronic supplementary material, section S5. The solutions for the M2 concentration and the porosity (which corresponds to the collagen concentration) are depicted in figure 6. The activated M2 is transported with the flow showing the highest concentration at the cell rear. The porosity is the highest near the cell membrane, as the degradation process is mainly invoked by M1 and M2 near the cell. In the vicinity, a very homogeneous porosity field is observed, due to steady M2 degradation. For analysis of the kinetics of the collagenolysis, the evolutions over time of the C1 and C1* concentration averaged over the domain are depicted in figure 7. There, a difference between both simulations can be observed. The concentrations reach a nearly steady state; however, with the applied force produces a slightly higher C1 concentration and a lower C1* concentration level. This is expected due to the stabilizing effect of the applied force. Still, the difference between loaded and unloaded configuration is very small. This has several reasons. First, the induced strain within the ECM by the cell is very small, which in turn produces only small changes in the chemical potential and the reaction rate. Second, only the degradation kinetics of M2 are modelled as load dependent. However, the main degradation here is due to M1. As already mentioned before, the dependency of the degradation kinetics due to M1 is not included although in [38] an increased degradation speed of loaded collagen due to M1 is observed in experiments. We did not include this dependency for the sake of simplicity and clearance of the interpretation of the results, yet it can be included into the model formulation in a straightforward way. In order to amplify the effect of M2 degradation the initial concentration of M2 is increased to 0.1 μmol dm−3 and also the same constant M2 concentration enters at the inflow. The evolutions over time of the C1 and C1* concentration averaged over the domain for this set-up are depicted in figure 8. There, the collagen dissolution due to M2 is emphasized and the difference between the loaded and unloaded case is more evident compared to a zero initial M2 concentration, however still small.
Figure 6. Cell example: spatial distributions of M2 concentration and porosity at the end of the simulation (t=5400 s). Results for simulation without force applied by the cell (simulation with force gives nearly indistinguishable results). (a) M2 concentration (μmol dm−3) and (b) porosity (−). Figure 7. Cell example: temporal evolution of C1 and C1* concentration averaged over the domain. Comparison of results with and without applied force. (a) C1 concentration (μmol dm−3) and (b) C1* concentration (μmol dm−3). Figure 8. Cell example: temporal evolution of C1 and C1* concentration averaged over the domain with initial presence of M2. Comparison of results with and without applied force. (a) C1 concentration (μmol dm−3) and (b) C1* concentration (μmol dm−3).
This example demonstrates the capability of the method to handle complex coupled flow and reaction processes. The stabilizing effect of mechanical loading is negligible for this scenario. However, this conclusion might be premature, as the model might not yet be complete to capture all relevant processes. More complex reactions with different dependencies of the reaction rates on the loading might be necessary. Additionally, the geometry of the cell and the flow pattern can have large effects on the results and an inhomogeneous ECM with spatially varying porosities and stiffness might amplify the importance of the biochemical–mechanical coupling.
4. Conclusion
In this contribution, we presented a computational model for the interaction of surface- and volume-bound scalar transport and reaction processes with a deformable porous medium. Flow and pressure induce displacements and alter the porous structure leading to different load bearing and transport behaviour. We used a strongly coupled algorithm to solve the system with the finite-element method. For the porous medium, a mixed finite-element approach suitable for high porosity gradients was used. The model was applied to specific phenomena in cell and tissue mechanics. Therefore, we developed a continuum formulation which incorporates the bi-directional coupling of the chemical degradation process and the material behaviour of the ECM. The destabilizing effect of the dissolution process on the structural integrity was considered via a changing reference porosity, which is similar to a damage parameter often applied in solid continuum mechanics. Changes of the mechanical energy induce changes in chemical potential, which in turn affects the reaction kinetics. We postulated an exponential relationship between reaction rate coefficient and chemical potential, which is based on known experimental and theoretical results. The whole model does not include any unphysical parameters. To the best of our knowledge, this is the first contribution extending a coupled transport-porous-medium model for the ECM by surface-volume reactions on the cell membrane, a dissolving collagen matrix and deformation-dependent dissolution kinetics. Thus, it gives a computational tool to investigate the interplay between biochemical dynamics and mechanical properties in cellular systems. The first numerical example reproduces an experiment, where different dissolution times for strained and unstrained fibres were observed. Our computational results produce good agreement with experimental data, even though the dissolution model is still quite simple. A parameter study indicates the influence of the most important material parameters. The second example comprises flow and transport processes within the ECM and the cell membrane. It includes matrix-tethered and free chemical substances involved in the degradation process. The mechano-biochemical coupling is triggered by a pulling force applied by the cell. The capability of the proposed methodology to model complex interactions and to analyse specific settings and parameters is demonstrated.
The presented method is capable of including further chemical interactions and other constitutive behaviour, like visco-elasticity, if necessary. Future developments may include the extension of our model to incorporate the intracellular signalling processes and their feedback to the cell environment. Also, more complex material behaviour, e.g. anisotropic fibre models and fibre remodelling in the ECM might be of interest [47]. For instance, this could give insights into the dynamics and the persistence of directional proteolytic migration of cells due to durotaxis.
Data accessibility
The datasets supporting this article have been uploaded as part of the electronic supplementary material.
Authors' contributions
A.-T.V. designed the model, performed implementation and simulations and drafted the manuscript. A.D.R. helped design the model, performed simulations and co-drafted the manuscript. W.A.W. helped with stating the problem and with drafting the manuscript. All authors gave final approval for publication.
Competing interests
We declare we have no competing interests.
Funding
This work was partly funded via the International Graduate School of Science and Engineering (IGSSE) of the Technical University of Munich, Germany, within the focus area Biomaterials.
