## 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
** S**. Therein

**=**

*C*

*F*^{T}

**depends of the deformation gradient**

*F***and is used for the pull-back operation of the fluid pressure. The Green–Lagrange strains**

*F***are defined as**

*E***=0.5(**

*E***−**

*C***1**) and

*Ψ*

^{s}denotes the macroscopic strain energy function of the skeleton. This formulation can be very convenient, as

*Ψ*

^{s}can be chosen as a classical strain energy function known from solid mechanics.

### (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 *v*^{c}=*v*^{f}−*v*^{s} denotes the convective velocity. All scalar quantities are summarized within the vector ** c** indicating that the reaction term

*σ*

_{i}(

**) potentially depends on all involved scalars. Equation (2.2) represents the convective formulation of the macroscopic instationary advection–diffusion–reaction problem within a deforming porous continuum. As the porous medium is deforming, the equation is formulated with respect to a moving observer, similar to an arbitrary Lagrange–Eulerian (ALE) formulation. Equation (2.2) holds for a scalar transported within the fluid. For a scalar**

*c*#### (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 *c*_{i}, *i*=1,…,*N* on a curved surface reads
** v** is the transport velocity and the operator

**∇**

_{Γ}denotes the tangential surface gradient. The latter can be calculated by the normal projection of the spatial gradient along the outward pointing unit normal

**, i.e.**

*n***∇**

_{Γ}=(

**1**−

**⊗**

*n***)**

*n***∇**.

#### (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})):

*μ*

^{s}

_{0}and the Helmholtz free energy

*ϕ*

_{ref}, which corresponds to the stress free reference state

*ϕ*

_{ref}=

*ϕ*

_{0}+Δ

*ϕ*

_{sol}. For Δ

*ϕ*

_{sol}=0 and

*μ*

^{s}

_{0}=0, the original energy function

*ϕ*

_{sol}/(1−

*ϕ*

_{0}). Using the above presented approach for the stress (2.1), one then obtains from the strain energy function (2.5)

*ϕ*

_{sol}. If Δ

*ϕ*

_{sol}=1−

*ϕ*

_{0}the solid phase is dissolved completely and thus does not any longer contribute to the mechanical stress. Equation (2.6) represents the key equation for the mechanical damage model. It should also be noted, that, as the reference porosity is allowed to change, the reference mass of the solid phase

*m*

^{s}

_{ref}changes equivalently:

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 [*C*1] of collagen type I by
*ρ*_{0,C} in [mol m^{−3}]. The factor 1/*J* represents a pull-back operation, as the molar concentration [*C*1] is written with respect to a spatial volume, whereas the reference porosity refers to a material volume. It is widely known that proteolysis is often influenced by mechanical stress and strain. Considering the dissolution process due to matrix metalloproteinase (MMP), experimental evidence suggest that mechanical strain stabilizes the ECM against proteolysis [6,7,36]. Also, there is theoretical work proposing that mechanical strain provides thermodynamical stabilization by loss of configurational entropy of the collagen fibres [37]. However, there are also cases where the opposite is observed. In [38], the preferential degradation due to MMP-1 of loaded collagen is reported. For a review of the current knowledge about collagen degradation mechanisms and the influence of mechanical load on the degradation rate, see [39].

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*:
*k*_{0} in the unloaded case, the molecular weight of collagen *m*_{mol,C1}, the gas constant *R* and the temperature *T*. The reaction rate changes due to the change Δ*μ*^{s} of chemical potential of the structural phase. Equation (2.9) is comparable to the reaction rate equation fitted from experimental data in [40]. There, the reaction rate is directly related to the applied force via an exponential function. We will use the chemical potential instead as it is more suitable for thermodynamical analysis and closer to the theoretical results from [37]. Also, note that equation (2.9) is very similar to the Arrhenius equation for reaction rates depending on temperature and activation energy, e.g. in [41], and therefore Δ*μ*^{s} can be interpreted as a change of the Gibbs energy of activation. Furthermore, using an analysis of the thermodynamics on the microscale and assuming only small changes of the chemical potential of the solid *μ*^{s} along the microscopic dissolution front, it can be shown [32], ch. 3.6.3 that the following identity holds:

## 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 Δ*T*_{S} = 1113 s ± 260 s for strained fibrils and Δ*T*_{U} = 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 C1*E* and the denatured collagen C1*. The chemical reactions involved can be schematically represented by a Michaelis–Menten mechanism, see e.g. [7],
*k*^{+}, *k*^{−} and *k*^{cat} denote the reaction rates. The kinetics are described by the set of reaction terms given in the electronic supplementary material, equations (S9)–(S12). The enzyme is transported within the porous fluid, while all other scalar quantities are bound to the structural phase of the ECM. As convective effects are very small, the flow of the fluid within the porous medium is not considered in this example. Instead, the focus lies on the mechano-chemical reaction model. As the matrix is expected to exhibit nonlinear constitutive behaviour [6], we assume a compressible Neo-Hookean material for the mechanical model of the ECM. As the dissolution process is inhibited due to mechanical loading, we model the reaction rate coefficient *k*^{+} as presented in equation (2.9) to be variable. The other reaction rate coefficients are assumed to be constant. The rate constant *T*_{U}≈700*s*. The rate constant *k*^{−} is then chosen as such, that the Michaelis constant *ρ*_{0,C} of collagen was estimated to be in same order of magnitude as water. All necessary parameters for this example are summarized in the electronic supplementary material, table S1. Further details on the set-up of the example and the numerical scheme can be found in the electronic supplementary material, section S4.

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 2*a*. 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 2*b*. 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.

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 3*a*). 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 3*b*), 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 3*c*). 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 3*d*). 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 4*a*. There, the stabilizing effect of the applied strain reduces with higher solid volume fraction. However, if the reaction coefficient *b*. If the latter results could be verified in experiments, this may indicate that there is an additional dependency of the reaction rates (

### (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 *MT*1-*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 *MMP*2, the proenzyme *pro*-*MMP*2, the inhibitor *TIMP*2, the complex *MMP*2⋅*TIMP*2, the denatured collagen C1* and the degradation product *MMP*2⋅C1. On the cell membrane, the modelled chemical species are the degradation enzyme *MT*1-*MMP*, its complex *MT*1-*MMP*⋅*TIMP*2* with inactivated *TIMP*2 and the complexes *MT*1⋅*TIMP*2⋅*MMP*2, *MT*1⋅*TIMP*2⋅*MMP*2⋅*MT*1 and *MT*1⋅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.

In the following, the modelled reaction system is described. Both enzymes *M*1 and *M*2 exert proteolytic activity, i.e. they can degrade the collagen matrix. *M*1 is acting on the cell surface, while *M*2 is diffusing within the ECM. The inhibitor *T*2 can bind both degradation enzymes *M*2 and *M*1. The inhibition of *M*1 is at the same time the first step towards the activation of *pM*2, which is secreted by the cell at the cell membrane. Here, *M*1 and *M*1*T*2 act as a catalyst for the production of *M*2. 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
*ρ*_{0,C} have been approximated in the same way as in the previous example. The reaction rates of the chemical reactions are listed in the electronic supplementary material, table S4. They are all taken from [18,46] and references therein. The diffusion coefficients for the respective transport model in the electronic supplementary material, table S5, are also taken from [18] for the species transported in the ECM. All diffusion coefficients on the cell membrane are chosen to be in a realistic order of magnitude [5].

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 *M*2 concentration and the porosity (which corresponds to the collagen concentration) are depicted in figure 6. The activated *M*2 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 *M*1 and *M*2 near the cell. In the vicinity, a very homogeneous porosity field is observed, due to steady *M*2 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 *M*2 are modelled as load dependent. However, the main degradation here is due to *M*1. As already mentioned before, the dependency of the degradation kinetics due to *M*1 is not included although in [38] an increased degradation speed of loaded collagen due to *M*1 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 *M*2 degradation the initial concentration of *M*2 is increased to 0.1 μmol dm^{−3} and also the same constant *M*2 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 *M*2 is emphasized and the difference between the loaded and unloaded case is more evident compared to a zero initial *M*2 concentration, however still small.

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.

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3699064.

- Received November 2, 2016.
- Accepted February 3, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.