Mathematical modelling of mechanotransduction via RhoA signalling pathways

TL;DR

This paper develops a mathematical model to study how cells convert mechanical signals into chemical ones through RhoA signaling, revealing how cell shape and stiffness affect this process.

Contribution

The novel contribution is a two-way coupled model of RhoA signaling and cell mechanics using bulk-surface finite elements, showing robustness in cell deformation.

Findings

Cell deformation is robust to changes in substrate stiffness due to signaling-mechanics coupling.

Cell shape strongly influences mechanotransduction dynamics.

A threshold-like response occurs with changes in substrate stiffness.

Abstract

We derive and simulate a mathematical model for mechanotransduction related to the Rho GTPase signalling pathway. The model addresses the bidirectional coupling between signalling processes and cell mechanics. A numerical method based on bulk-surface finite elements is proposed for the approximation of the coupled system of nonlinear reaction-diffusion equations, defined inside the cell and on the cell membrane, and the equations of elasticity. Our simulation results illustrate novel emergent features such as the strong dependence of the dynamics on cell shape, a threshold-like response to changes in substrate stiffness, and the fact that coupling mechanics and signalling can lead to the robustness of cell deformation to larger changes in substrate stiffness, ensuring mechanical homeostasis in agreement with experiments.

Mechanotransduction, a process by which cells convert mechanical stimuli into chemical signals, plays a crucial role in cell functions. To better understand this phenomenon, we need to analyse how signalling processes and cell mechanics work together. For this purpose we derive and simulate a mathematical model of mechanotransduction related to the Rho GTPase signalling pathways, central to almost all fundamental cellular processes including cell polarity, movement, division, and cytoskeleton reorganization. The model introduces a two-way coupling between the signalling processes and cell mechanics. We use a numerical method based on bulk-surface finite elements to solve model equations numerically. Our simulation results illustrate novel emergent features such as a strong dependence of the dynamics on cell shape, a threshold-like response to changes in substrate stiffness, and the fact that the two-way coupling between mechanics and signalling can lead to the robustness of cell deformation to larger changes in substrate stiffness, ensuring mechanical homeostasis in agreement with experiments. These interesting insights help us unravel the underlying mechanisms in mechanotransduction.

Figures

A sketch of the model interactions.The diagram shows a simplified overview of the interactions in the models presented in this paper, where the blue ellipses represent the signalling molecules and the green rectangles the mechanical signals. The black arrows represent the reactions of the biochemical model in Eq (1) based on [18], i.e. ECM stiffness activates FAK on the cell membrane which activates RhoA. The grey arrows represent the coupling between chemistry and mechanics introduced in this work and described by Eqs (3)-(6), i.e. activated RhoA affects the cytoplasmic elastic stress which activates FAK and the cytoplasmic stiffness is a function of activated FAK.

Numerical simulation results showing ρa, ϕd, ϕa and |u| for the model in Eqs (3)–(6) for the axisymmetric shape and in the case of 2xD stimulus at a steady state at T = 100 s.Four different scenarios are considered: (A) C1 = 0 (kPa s)−1 (σ↛ϕa) and Ec = 0.6 kPa (ϕa↛Ec); (B) C1 = 0 (kPa s)−1 (σ↛ϕa) and Ec=f(ϕa) (ϕa→Ec); (C) C1 = 0.1 (kPa s)−1 (σ→ϕa) and Ec = 0.6 kPa (ϕa↛Ec); (D) C1 = 0.1 (kPa s)−1 (σ→ϕa) and Ec=f(ϕa) (ϕa→Ec). Within each subfigure, the rows represent ρa, ϕd, ϕa and |u| on a cross-section of the plane x1 = 0 of the axisymmetric cell, and the columns represent E=0.1,5.7,7·106kPa. Parameter values as in Table 1.

Numerical simulation results showing ρa, ϕd, ϕa and |u| for the model in Eqs (3)–(6) for the lamellipodium shape and in the case of 2xD stimulus at a steady state at T = 100 s.Four different scenarios are considered: (A) C1 = 0 (kPa s)−1 (σ↛ϕa) and Ec = 0.6 kPa (ϕa↛Ec); (B) C1 = 0 (kPa s)−1 (σ↛ϕa) and Ec=f(ϕa) (ϕa→Ec); (C) C1 = 0.1 (kPa s)−1 (σ→ϕa) and Ec = 0.6 kPa (ϕa↛Ec); (D) C1 = 0.1 (kPa s)−1 (σ→ϕa) and Ec=f(ϕa) (ϕa→Ec). Within each subfigure, the rows represent ρa, ϕd, ϕa and |u| on the surface of the cell, and the columns represent E=0.1,5.7,7·106 kPa. Parameter values as in Table 1.

Simulation results showing the mean, 1|Ω|∫Ω·dx, min and max values of f(ϕa),  div(u), ϕa and ρa as functions of substrate stiffness E.We consider different couplings with four different values for C1 and two different shapes at T = 100 s by which time the results are at a steady state. All other parameter values as in Table 1.

Numerical simulation results showing ρa, ϕd, ϕa and |u| for the model in Eqs (3)–(6) for the axisymmetric shape and in the case of the 3D stimulus at a steady state at T = 100 s.Four different scenarios are considered: (A) C1 = 0 (kPa s)−1 (σ↛ϕa) and Ec = 0.6 kPa (ϕa↛Ec); (B) C1 = 0 (kPa s)−1 (σ↛ϕa) and Ec=f(ϕa) (ϕa→Ec); (C) C1 = 0.1 (kPa s)−1 (σ→ϕa) and Ec = 0.6 kPa (ϕa↛Ec); (D) C1 = 0.1 (kPa s)−1 (σ→ϕa) and Ec=f(ϕa) (ϕa→Ec). Within each subfigure, the rows represent ρa, ϕd, ϕa and |u| on a cross-section of the plane x1 = 0 of the axisymmetric cell, and the columns represent E=0.1,5.7,7·106kPa. Parameter values as in Table 1.

Simulation results showing the mean, 1|Ω|∫Ω·dx, min and max values of f(ϕa),  div(u), ϕa and ρa as functions of substrate stiffness E.We consider different couplings with four different values for C1 and two different shapes at T = 100 s by which time the results are at a steady state. All other parameter values as in Table 1.

Numerical simulation results showing ρa, ϕd, ϕa and |u| for the model in Eqs (3), (4), and (6) for the axisymmetric shape and in the case of the 3D stimulus at a steady state at T = 100 s.Four different scenarios are considered: (A) C1 = 0 (kPa s)−1 (σ↛ϕa) and Ec = 0.6 kPa (ϕa↛Ec); (B) C1 = 0 (kPa s)−1 (σ↛ϕa) and Ec=f(ϕa) (ϕa→Ec); (C) C1 = 0.1 (kPa s)−1 (σ→ϕa) and Ec = 0.6 kPa (ϕa↛Ec); (D) C1 = 0.1 (kPa s)−1 (σ→ϕa) and Ec=f(ϕa) (ϕa→Ec). Within each subfigure, the rows represent ρa, ϕd, ϕa and |u| on a cross-section of the plane x1 = 0 of the axisymmetric cell, and the columns represent E=0.1,5.7,7·106 kPa. Parameter values as in Table 1.

Simulation results showing the mean, 1|Ω|∫Ω·dx, min and max values of f(ϕa),  div(u), ϕa and ρa as functions of substrate stiffness E, in the case of the model in Eqs (3), (4) and (6) and 3D stimulus.We consider different couplings, four different values for C1, and two different shapes at T = 100 s by which time the results are at a steady state. All other parameter values as in Table 1.

Numerical simulation results showing ρa, ϕd, ϕa and |u| simulation results showing ϕa and ρa for the model in Eqs (3), (4), and (6) for the axisymmetric shape and in the case of 2xD stimulus at a steady state at T = 100 s.Four different scenarios are considered: (A) C1 = 0 (kPa s)−1 (σ↛ϕa) and Ec = 0.6 kPa (ϕa↛Ec); (B) C1 = 0 (kPa s)−1 (σ↛ϕa) and Ec=f(ϕa) (ϕa→Ec); (C) C1 = 0.1 (kPa s)−1 (σ→ϕa) and Ec = 0.6 kPa (ϕa↛Ec); (D) C1 = 0.1 (kPa s)−1 (σ→ϕa) and Ec=f(ϕa) (ϕa→Ec). Within each subfigure, the rows represent ρa, ϕd, ϕa and |u| on a cross-section of the plane x1 = 0 of the axisymmetric cell, and the columns represent E=0.1,5.7,7·106 kPa. Parameter values as in Table 1.

Equations

Full text

Introduction

Intercellular signalling processes constitute the mechanisms through which cells communicate with and respond to their environment. Hence, signalling pathways are important in all physiological activities of the cell, such as cell division, cell movement, the immune response, and tissue development [1]. Aberrant cell signalling can often result in the development of diseases [2]. It is therefore important to understand signalling phenomena. Recent studies have found that alongside biochemical reactions, mechanics plays an important role in many signalling pathways [3,4]. This phenomena is referred to as mechanotransduction which, broadly speaking, is any process by which cells convert mechanical stimuli into chemical signals [5,6].

A large number of recent works study the role of Rho GTPases, primarily RhoA, in mechanotransduction in relation to different mechanical cues: extracellular matrix (ECM) stiffness and viscoelasticity, tensile stress (stretching), compressive stress (compression), and shear stress (fluid flow shear), see e.g. [7,8] for a review. Moreover, the coupling between biochemistry and mechanics is bidirectional, i.e., chemical signals can also affect the mechanical properties of the cell, such as molecules like focal adhesion kinases (FAKs) that influence F-actin dynamics and therefore the stiffness of the cell [6,9–11].

The formidable complexity of the phenomena involved in mechanotransduction means that much about how the mechanics and the chemical processes of the cell communicate is not yet understood and mathematical modelling is crucial in this regard. Whilst the mathematical modelling of biochemical cell signalling processes is fairly well developed, e.g., [12–14], the study of mechanotransduction is comparatively more recent, see [15] for a review. Typically the modelling involves solving coupled systems of partial differential equations (PDEs) with reaction-diffusion equations modelling the biochemistry coupled to equations based on (visco)elastic constitutive laws for the mechanics. The progress of such efforts has been rapid, ranging from early models employing simplifications such as one-dimensional geometries [16,17] to full three-dimensional simulations [18] using advanced computational techniques. Alongside continuum models, a number of recent works have employed discrete approaches such as spring-based models [19], or models that employ a Potts formalism [20–22]. Despite this rapid progress, the existing models typically make major simplifying assumptions such as assuming a constant stiffness of the ECM [18,23,24], as well as neglecting the two-way coupling in which signalling pathways affect the mechanics alongside mechanical cues inducing signalling processes.

In the present work, we seek to develop, analyse and simulate a model for mechanotransduction through the Rho GTPase signalling pathway which allows for a two-way coupling between the mechanics and the biochemistry. The dynamics of the signalling molecules FAK and RhoA are modelled using reaction-diffusion equations, where the ECM stiffness and elastic stresses of the cell activate FAK. Under simplifying assumptions, i.e., assuming no dependence on the cell elastic stresses, the biochemical component of the model is derived as a reduction of the model proposed in [18]. For the cell’s mechanical properties, we assume an elastic constitutive relationship [25] and allow the material properties to depend on the concentrations of the signalling molecules. We propose a numerical method based on bulk and surface finite elements [26] for the approximation of the model equations.

The results presented here show that our model can reproduce the qualitative results of [18], i.e., the mass of activated FAK and RhoA depend on ECM stiffness, with the dependence captured by a Hill function. On the inclusion of the two-way coupling between signalling processes and cell mechanics, we observe novel dynamics, such as the conservation of cell deformation under different values of the ECM stiffness, which underlines the importance of including these more complex models of the mechanics. The role of mechanotransduction in homeostasis in biological processes has been discussed in a number of biological works, e.g., [4,5,9,10,27] and our work presents a concrete example of how modelling can help elucidate potential mechanisms that underlay the mechanical homeostasis. Homeostasis of cell deformation, as observed in simulations of our model, has been observed experimentally [28]. Our focus is on elastic constitutive assumptions for the mechanics of the cell to enhance clarity of exposition and to avoid unnecessary technical complexities. This can be extended to allow for other constitutive laws such as viscoelasticity of the cell and/or of the ECM as has been done elsewhere in the literature in simpler settings in 1D [16,29]. This work thus serves as a starting point in modelling and analysis of the two-way coupling between mechanics and chemistry.

The paper is organised as follows. We first derive the reduced model for the Rho GTPase signalling pathway, based on the model proposed in [18]. Next, the mathematical model for the mechanotransduction related to the Rho GTPase signalling pathway is derived. Simulations of the model are presented in the section after. We conclude the paper with a discussion of the results. Details on the numerical method applied to simulate the model equations are given in S1 Appendix, Sect A.3.

Methods

A mathematical model for the Rho GTPase signalling pathway

One of the main signaling pathways involved in mechanostranduction is the Rho GTPase pathway, responsible for many important cellular processes, e.g. motility, cell adhesion, polarisation, differentiation, remodelling of the exoskeleton, and the ECM [8]. The RhoA signalling pathway is activated through the activation of FAK in response to tension on integrins, which depends on ECM stiffness [1,11], see Fig 1 for an overview.

A sketch of the model interactions.The diagram shows a simplified overview of the interactions in the models presented in this paper, where the blue ellipses represent the signalling molecules and the green rectangles the mechanical signals. The black arrows represent the reactions of the biochemical model in Eq (1) based on [18], i.e. ECM stiffness activates FAK on the cell membrane which activates RhoA. The grey arrows represent the coupling between chemistry and mechanics introduced in this work and described by Eqs (3)-(6), i.e. activated RhoA affects the cytoplasmic elastic stress which activates FAK and the cytoplasmic stiffness is a function of activated FAK.

Our model for mechanotransduction related to the RhoA-mediated intercellular signalling pathway is based on models developed in [18] and [23,24]. To incorporate the interactions between mechanics and signalling processes, we extend the model proposed in [18] by considering elastic deformations of the cell. Activated FAK is downstream in the RhoA GTPase signalling pathway and hence the activation of RhoA is a function of activated FAK. The activation of RhoA results in ECM remodelling and deposition of new fibres, increasing ECM stiffness and hence activation of FAK [30]. FAK is expressed in the cytoplasm of the cell and is activated on the cell membrane. To simplify the model and focus only on the most significant aspects from the perspective of qualitative behaviour, we reduce the model for the RhoA signalling pathway of [18] that includes the dynamics of FAK, RhoA, ROCK, Myo, LIMK, mDia, Cofilin, F-actin and YAP/TAZ by considering only the dynamics of FAK and activated RhoA. Such a reduction is possible since other molecules considered in the full model of [18] do not influence the dynamics of FAK and RhoA. Our rationale behind considering a simplified model is to more clearly elucidate the emergent features that arise when mechanics is coupled with signalling. It is not challenging to incorporate other biochemical species or different reaction kinetics within the framework we propose.

We let $[eqn]$ denote the cytoplasm and $[eqn]$ the cell membrane. We denote by $[eqn]$ and $[eqn]$ the concentrations of inactive and active FAK, and by $[eqn]$ the concentration of active RhoA. We recall that FAK, both active and inactive, is assumed to be cytoplasm resident and activated RhoA membrane resident. Our reduced model for the biochemistry consists of the following system of reaction-diffusion equations

$$[eqn]$$

where $[eqn]$ is the Laplace Beltrami operator modelling diffusion on the surface $[eqn]$ , see e.g. [26], $[eqn]$ is the ratio between the volume of the cytoplasm and the area of the cell membrane, $[eqn]$ are deactivation and $[eqn]$ are activation constants, E is the substrate stiffness, $[eqn]$ are the diffusion constants, C, n and $[eqn]$ are positive constants, $[eqn]$ , $[eqn]$ and $[eqn]$ are bounded nonnegative functions, and $[eqn]$ is an approximation for deactivated RhoA ( $[eqn]$ ) on the surface with $[eqn]$ the total mass of RhoA, a quantity conserved in the full model of [18] and assumed to be conserved here. Simulations illustrating the agreement between results obtained using the reduced model (1) with those of [18] for the full model are presented in S1 Appendix, Sect A.1.

Mathematical model for mechanotransduction

As a starting point for the mechanics, we consider small deformations and hence, assume a linear elastic constitutive law for the mechanics of the cell. Although viscoelastic or poroelastic behaviour of cells is proposed in many works [31,32], linear elasticity is often chosen for modelling simplicity and it can yield results consistent with experimental observations [33–35]. Our focus is modelling a bidirectional coupling between cell stiffness and signalling processes, a minimal model assuming a linear elastic law for the cell mechanics is therefore sufficient for this work, as it avoids unnecessary complications that arise from the consideration of viscous stresses. To demonstrate that our results remain relevant under more complex assumptions, a linear viscoelastic model is presented and simulated in Sect A.7 in S1 Appendix. We note that for the corresponding simulation results, the qualitative behaviour in both cases, linear elastic and linear viscoelastic, is the same.

An important simplification that arises under the small deformations assumption inherent in this work is that the dynamics of the signalling molecules may be effectively considered on the reference configuration with no additional terms arising due to the deformation. Models where the assumption of small deformations is relaxed will be addressed in future studies.

The cell nucleus plays an important role in governing the mechanical properties of the cell [36,37], whilst we predominantly neglect this in the present work, in S1 Appendix, Sect A.5, we have included simulations of the model with a ‘passive’ nucleus that is considered to be more rigid than the cytoplasm.

It has been shown that the stiffness of the cell increases as F-actin increases [18]. F-actin is downstream from activated FAK and, as apparent from the results in [18], we can use activated FAK as a proxy for F-actin. Therefore, we assume that the Young’s modulus Ec of the cell is a function of the activated FAK concentration. Based on experimental observations [38] and numerical simulations [18], we propose

$$[eqn]$$

where k7, k8 and p are non-negative constants. Then for elastic deformations of the cell, we have

$$[eqn]$$

with

$$[eqn]$$

and the Lame constants $[eqn]$ and $[eqn]$ are given by

$$[eqn]$$

where $[eqn]$ is the Poisson ratio of the cell.

Activated RhoA regulates remodellling of stress fibres inside the cell and stabilisation of actin filaments [39–41,43]. This mechanism is modelled by the stress on the boundary being dependent on activated RhoA concentration

$$[eqn]$$

where k6 is a positive constant, $[eqn]$ is a projection on the space orthogonal to the space of rigid deformations, i.e. rotations and translations. Alongside models where the cell is allowed to deform freely, to model a typical experimental set-up where cells are placed on a rigid substrate, we consider

$$[eqn]$$

together with condition Eq (4) on $[eqn]$ , where $[eqn]$ denotes the tangential projection of vector w, and hence $[eqn]$ is the shear stress and Eq (5) specifies a shear stress free condition. In this work, we choose $[eqn]$ to simulate a rigid substrate.

It has been shown that an increased contractility is associated with increased activated FAK, see e.g. [39]. Thus we assume that FAK is activated by the stress of the cell and as a proxy for the cytosolic stress we use the positive part of trace of the Cauchy stress tensor $[eqn]$ , where $[eqn]$ is the first stress invariant and the positive part reflects the fact that extension rather than compression causes the activation of FAK. This modifies the system in Eq (1) to

$$[eqn]$$

where $[eqn]$ . We can prove existence, uniqueness and boundedness of solutions to the system in Eqs (3)–(6) which we intend to report on elsewhere.

Initial and boundary conditions and parametrization

Using the model in Eqs (3)–(6) we investigate different scenarios demonstrating the interactions between mechanics and signalling processes. First, we consider the impact of the cell Young’s modulus Ec and compare the dynamics when considering a constant Ec versus the case where Ec depends on activated FAK as defined in Eq (2). We also model the effect of the stress on the signalling molecules FAK and simulate equations Eq (6) for C1 = 0 $[eqn]$ and C1 = 0.1 $[eqn]$ , respectively. Additionally we consider two experimental scenarios: (i) the cell is placed on a rigid substrate, modelled by the boundary conditions Eq (4) on $[eqn]$ and Eq (5) on $[eqn]$ or (ii) the cell is embedded in an agar substrate and we apply the force boundary condition Eq (4) on the entire cell membrane. We also distinguish between two different stimuli, similar to [18], (i) the so called ‘2xD stimulus’, where the substrate stiffness is only applied to the bottom of the cell, i.e. E is nonzero only on $[eqn]$ , and (ii) the ‘3D stimulus’ where the cell is embedded in an agar (substrate) and the impact of the substrate stiffness on the signalling processes is considered on the whole cell membrane. To analyse the impact of the cell shape on the dynamics of signalling molecules and mechanical deformations we consider both axisymmetric cells and polarised cells with a lamellipodium like shape. The diameter of the cell is larger for the lamellipodium cells such that the volume is similar to the axisymmetric cells.

The parameters are chosen as in Table 1. For numerical simulations, we use a Finite Element Method to discretize in space and a semi-implicit Euler method to discretize in time, with the mesh size $[eqn]$ m and time step $[eqn]$ s. Details on the numerical scheme and benchmark computations demonstrating the accuracy of the approach for a problem with a known solution are given in S1 Appendix, Sect A.3.

Table 1: Parameter values for simulations of the model in Eqs (3)–(6).

Results

Numerical simulations with 2xD stimulus

First we look at the results that would most reflect a cell on a substrate in vitro. Here, the substrate stiffness appears as a stimulus only on the bottom boundary of the cell, i.e. E is nonzero only on $[eqn]$ , and deformation is restricted in the vertical direction at the bottom boundary of the cell. The results for the axisymmetric shape of the cell are found in Fig 2, whereas results for the lamellipodium shape are presented in Fig 3. Note that results for $[eqn]$ and $[eqn]$ when Ec = 0.6 kPa and C1 = 0 $[eqn]$ are identical to the one without mechanics in S1 Appendix, Sect A.1. In this case, we see that the magnitude of the deformation $[eqn]$ is largest at the edge of the cell. The cell expands axisymmetrically at the base. As expected, the expansion is larger for higher concentrations of $[eqn]$ . For a lower substrate stiffness, E = 0.1 kPa, the cell barely expands. When C1 = 0.1 $[eqn]$ , the concentrations of $[eqn]$ and $[eqn]$ and the magnitude of the deformation $[eqn]$ increase, with a bigger increase for lower substrate stiffness and a smaller increase for larger substrate stiffness. When comparing Ec = 0.6 kPa and $[eqn]$ , the deformations show similar patterns, expanding at the base of the cell, however, the magnitude of the deformation is much lower in the case $[eqn]$ . This is probably because $[eqn]$ kPa for a small substrate stiffness E, but is doubled in magnitude for larger substrate stiffness, see Fig 4. The larger cell Young’s modulus Ec means it is harder for the cell to deform, resulting in a lower magnitude of deformation. This difference illustrates that, unlike the constant Young’s modulus case, a concentration-dependent Young’s modulus allows for potential homeostasis and adaptation of cell mechanics to different values of the substrate stiffness [28].

Numerical simulation results showing ρa, ϕd, ϕa and |u| for the model in Eqs (3)–(6) for the axisymmetric shape and in the case of 2xD stimulus at a steady state at T = 100 s.Four different scenarios are considered: (A) C1 = 0 (kPa s)−1 (σ↛ϕa) and Ec = 0.6 kPa (ϕa↛Ec); (B) C1 = 0 (kPa s)−1 (σ↛ϕa) and Ec=f(ϕa) (ϕa→Ec); (C) C1 = 0.1 (kPa s)−1 (σ→ϕa) and Ec = 0.6 kPa (ϕa↛Ec); (D) C1 = 0.1 (kPa s)−1 (σ→ϕa) and Ec=f(ϕa) (ϕa→Ec). Within each subfigure, the rows represent ρa, ϕd, ϕa and |u| on a cross-section of the plane x1 = 0 of the axisymmetric cell, and the columns represent E=0.1,5.7,7·106kPa. Parameter values as in Table 1.

Numerical simulation results showing ρa, ϕd, ϕa and |u| for the model in Eqs (3)–(6) for the lamellipodium shape and in the case of 2xD stimulus at a steady state at T = 100 s.Four different scenarios are considered: (A) C1 = 0 (kPa s)−1 (σ↛ϕa) and Ec = 0.6 kPa (ϕa↛Ec); (B) C1 = 0 (kPa s)−1 (σ↛ϕa) and Ec=f(ϕa) (ϕa→Ec); (C) C1 = 0.1 (kPa s)−1 (σ→ϕa) and Ec = 0.6 kPa (ϕa↛Ec); (D) C1 = 0.1 (kPa s)−1 (σ→ϕa) and Ec=f(ϕa) (ϕa→Ec). Within each subfigure, the rows represent ρa, ϕd, ϕa and |u| on the surface of the cell, and the columns represent E=0.1,5.7,7·106 kPa. Parameter values as in Table 1.

Simulation results showing the mean, 1|Ω|∫Ω·dx, min and max values of f(ϕa),  div(u), ϕa and ρa as functions of substrate stiffness E.We consider different couplings with four different values for C1 and two different shapes at T = 100 s by which time the results are at a steady state. All other parameter values as in Table 1.

For the two-way couplings between the mechanics and chemistry, i.e. $[eqn]$ and C1 = 0.1 $[eqn]$ , we see similar results for the deformation as when $[eqn]$ and C1 = 0 $[eqn]$ . The main difference is that the deformation for E = 0.1 kPa is now at a similar magnitude as for the larger substrate stiffnesses, demonstrating the importance of the signalling processes in the adaptation of cell mechanics to changing environmental conditions.

Comparing the simulation results for the two different shapes in Figs 2 and 3, the concentration of activated RhoA, $[eqn]$ , is slightly lower for the lamellipodium shape. For the lamellipodium shape, we observe the largest deformations at the corners furthest from the nucleus.

Fig 4 summarises the results at time T = 100 s by plotting the mean, $[eqn]$ , of $[eqn]$ , the volume change $[eqn]$ , $[eqn]$ , and $[eqn]$ as functions of the substrate stiffness E, with the bars being the range of these variables, for different values of the constant C1 in the activation of FAK by the cell stress. As expected, an increase in C1 results in an increase in the concentration of activated FAK, $[eqn]$ . The dependence of $[eqn]$ on the substrate stiffness E, especially for C1 = 0 $[eqn]$ and C1 = 0.5 $[eqn]$ , resembles a Hill function representing a threshold response. This agrees with simulations in [18] which themselves fit experimental observations presented in [44]. For most of the cases, the results for the lamellipodium shape are very similar to the results for the axisymmetric shape. However, for C1 = 0.1 $[eqn]$ , the magnitude of the threshold-like response in all variables is bigger in the lamellipodium case. In terms of the Young’s modulus, when Ec = 0.6 kPa we observe much larger volume changes than when $[eqn]$ in all the numerical experiments.

Numerical simulations for the 3D stimulus case on a rigid substrate

In numerical simulations for a 3D stimulus on a rigid substrate, the substrate stiffness affects the whole cell membrane and we consider the boundary conditions Eq (4) on $[eqn]$ and Eq (5) on $[eqn]$ . The results for numerical experiments can be found in Figs 5 and 6. Overall, the concentrations $[eqn]$ and $[eqn]$ are larger than in the case of the 2xD stimulus, which is in line with the results in [18]. The higher concentrations of $[eqn]$ results in larger deformations, where the maximum magnitude of the deformation in the case of the 2xD stimulus was $[eqn]$ m, see Fig 2, while the maximum magnitude of the deformation in the case of the 3D stimulus is $[eqn]$ m, see Fig 2. Similar behaviour is observed for the lamellipodium shape, see Figs 3 and S1. Another difference between two cases are larger variations in concentration and a larger difference between maximal and minimal values in the case of the 2xD stimulus than in the case of 3D stimulus, see Figs 4 and 6. Similar behaviour is observed also in the model for the signalling processes without mechanics, see S1 Appendix, Sect A.1.

Numerical simulation results showing ρa, ϕd, ϕa and |u| for the model in Eqs (3)–(6) for the axisymmetric shape and in the case of the 3D stimulus at a steady state at T = 100 s.Four different scenarios are considered: (A) C1 = 0 (kPa s)−1 (σ↛ϕa) and Ec = 0.6 kPa (ϕa↛Ec); (B) C1 = 0 (kPa s)−1 (σ↛ϕa) and Ec=f(ϕa) (ϕa→Ec); (C) C1 = 0.1 (kPa s)−1 (σ→ϕa) and Ec = 0.6 kPa (ϕa↛Ec); (D) C1 = 0.1 (kPa s)−1 (σ→ϕa) and Ec=f(ϕa) (ϕa→Ec). Within each subfigure, the rows represent ρa, ϕd, ϕa and |u| on a cross-section of the plane x1 = 0 of the axisymmetric cell, and the columns represent E=0.1,5.7,7·106kPa. Parameter values as in Table 1.

Simulation results showing the mean, 1|Ω|∫Ω·dx, min and max values of f(ϕa),  div(u), ϕa and ρa as functions of substrate stiffness E.We consider different couplings with four different values for C1 and two different shapes at T = 100 s by which time the results are at a steady state. All other parameter values as in Table 1.

Numerical simulations for the model in Eqs (3), (4), and (6).

To investigate a setting more close to a cell in vivo, we consider the coupled model in Eqs (3), (4), (6) with force boundary conditions on the whole cell membrane, without restricting the deformation on the bottom of the cell.

Numerical simulations in the case of 3D stimulus.

Simulation results for a 3D stimulus that models a cell surrounded by the extracellular matrix are presented in Figs 7 and 8. Corresponding results of the evolution of the mean of $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ can be found in S1 Appendix, Sect A.2. The results show the same differences between the different couplings as in Figs 5 and 6. Comparing Figs 5 and 7, the results for the concentrations $[eqn]$ and $[eqn]$ are indistinguishable, however there is a clear difference in deformation of the bottom of the cell and in the case of the fixed vertical deformations the magnitude of the deformation at the base of the cell is slightly lower than in the case of force boundary conditions. The same differences are observed for the lamellipodium shape case, see S1 and S2 Figs. The concentrations of the signalling molecules are also less sensitive to parameter changes than the volume change, as can be seen from the parameter analysis in S1 Appendix, Sect A.6. Comparing Figs 6 and 8, the main difference is in the behaviour of $[eqn]$ as a function of E. Even though the average volume change is the same, we see differences in the maximum and minimum values of the local volume change across the domain. In particular, the maximum local volume change when considering the model with a partially fixed boundary is larger and is located on the base of the cell, while the maximum local volume change when considering the model with the force boundary conditions is smaller, but the cell deforms more evenly in all directions.

Numerical simulation results showing ρa, ϕd, ϕa and |u| for the model in Eqs (3), (4), and (6) for the axisymmetric shape and in the case of the 3D stimulus at a steady state at T = 100 s.Four different scenarios are considered: (A) C1 = 0 (kPa s)−1 (σ↛ϕa) and Ec = 0.6 kPa (ϕa↛Ec); (B) C1 = 0 (kPa s)−1 (σ↛ϕa) and Ec=f(ϕa) (ϕa→Ec); (C) C1 = 0.1 (kPa s)−1 (σ→ϕa) and Ec = 0.6 kPa (ϕa↛Ec); (D) C1 = 0.1 (kPa s)−1 (σ→ϕa) and Ec=f(ϕa) (ϕa→Ec). Within each subfigure, the rows represent ρa, ϕd, ϕa and |u| on a cross-section of the plane x1 = 0 of the axisymmetric cell, and the columns represent E=0.1,5.7,7·106 kPa. Parameter values as in Table 1.

Simulation results showing the mean, 1|Ω|∫Ω·dx, min and max values of f(ϕa),  div(u), ϕa and ρa as functions of substrate stiffness E, in the case of the model in Eqs (3), (4) and (6) and 3D stimulus.We consider different couplings, four different values for C1, and two different shapes at T = 100 s by which time the results are at a steady state. All other parameter values as in Table 1.

Numerical simulations in the case of 2xD stimulus.

In Figs 9, S3 and S4 we report on simulation results in the case of 2xD stimulus and force boundary conditions applied to the entire boundary. For the concentrations, the results are similar to the results in the case of 2xD stimulus and no vertical deformation on the bottom of the cell, see Figs 2 and 9. However, the results for the deformation are different compared to the previous results. In Fig 9, the cell does not just expand but changes shape as the edges of the cell deform upwards, which is not possible in the case of the partially fixed boundary as we assume no vertical deformation at the base. The deformation of the cell upwards can also be observed in the case of the 3D stimulus, but it is smaller due to the impact of the ECM surrounding the cell, see Fig 7. We observe that for C1 = 0 $[eqn]$ the cell deforms upwards a little more than for C1 = 0.1 $[eqn]$ . This is due to the larger variation in the concentration $[eqn]$ for C1 = 0 $[eqn]$ compared to C1 = 0.1 $[eqn]$ . The same features are observed for the lamellipodium shape, see Figs 3, S2, S3 and S4.

Numerical simulation results showing ρa, ϕd, ϕa and |u| simulation results showing ϕa and ρa for the model in Eqs (3), (4), and (6) for the axisymmetric shape and in the case of 2xD stimulus at a steady state at T = 100 s.Four different scenarios are considered: (A) C1 = 0 (kPa s)−1 (σ↛ϕa) and Ec = 0.6 kPa (ϕa↛Ec); (B) C1 = 0 (kPa s)−1 (σ↛ϕa) and Ec=f(ϕa) (ϕa→Ec); (C) C1 = 0.1 (kPa s)−1 (σ→ϕa) and Ec = 0.6 kPa (ϕa↛Ec); (D) C1 = 0.1 (kPa s)−1 (σ→ϕa) and Ec=f(ϕa) (ϕa→Ec). Within each subfigure, the rows represent ρa, ϕd, ϕa and |u| on a cross-section of the plane x1 = 0 of the axisymmetric cell, and the columns represent E=0.1,5.7,7·106 kPa. Parameter values as in Table 1.

Discussion and conclusion

We have derived a model for mechanotransduction via the RhoA signalling pathway with ECM stiffness and intracellular mechanical properties serving as the mechanical cues. The modelling extends the work of [18] incorporating the explicit modelling of cell deformation based on an elastic constitutive assumption. We have extended on [18,23,24] and introduced a two-way coupling between the mechanics of the cell and biochemical signalling processes. This two-way coupling appears to be central to mechanical homeostasis observed in biological experiments [28]. We propose a robust numerical method, based on the bulk-surface finite element method (FEM), see e.g. [26], for the approximation of the model and report on simulation results for different scenarios, validating the results by comparison with simulations presented in [18] and experimental observations in [44]. Namely, we considered different levels of substrate stiffness for cells of different shapes that either sit on a rigid flat substrate or are embedded in a three-dimensional substrate.

Our broad conclusions are that cell shape strongly influences the dynamics of the signalling molecules and the deformation of the cell, as seen in all figures comparing the axisymmetric and lamellipodium shape, where the emergent patterns differ, which is in line with experimental observations [45,46]. Cell shape also affects experimentally observed features such as the threshold-like response to changes in substrate stiffness [44] which is reproduced by the model. In Figs 4, 6, 8 and S4, we see that for certain parameters (C1 = 0.1 $[eqn]$ and low substrate stiffness), the cell shape affects the mean concentrations of the signalling molecules and the mean volume change of the cell, and thus changes the threshold-like response.

Our simulations exhibit novel emergent features, that are inaccessible without the framework we propose, such as the bidirectional coupling between mechanics and signalling processes through allowing the Young’s modulus of the cell to depend on protein concentration that can allow for robustness in terms of the magnitude of deformation in response to differences in substrate stiffness. This is an example of a mechanical homeostasis mechanism that emerges only at this level of modelling complexity which is of relevance to biology [28]. Other instances of mechanical homeostasis are the stress being maintained in the cardiovascular system under mechanical perturbations [47] and the tensional homeostasis by the RhoA signalling pathway at the level of multiple cells [48,49], which is known to be governed by cellular stiffness sensing [50]. Another mechanism that experiences homeostatic response to substrate stiffness is that of the mechanical memory of the cell, describing the phenomenon of a cell responding less to substrates with lower stiffness if they have been cultured on stiff substrates [49,51]. Due to the bidirectional coupling between the mechanics and the chemistry in our modelling framework, an extension of this work by changing the chosen couplings could be used to model these other mechanical homeostasis phenomena.

Based on previous biological studies [18,38], we considered cases in which the mechanical properties of the cell (cell stiffness) depend on the concentration of signalling molecules. This coupling yields less sensitivity of total deformation to substrate stiffness whilst leaving the dynamics of the signalling molecules themselves broadly unchanged, see Figs 4, 6, 8 and S4. The insensitivity of the dynamics of the signalling molecules to deformation levels arises since in the model proposed here they are influenced by the local stress rather than deformation. We note that the above constitutes another emergent homeostasis mechanism that the modelling framework allows us to explore. We stress that our work serves as an example of how mechanotransduction may be modelled and more complicated models for the mechanics, biochemistry and couplings therefore are warranted based on the remarkable emergent features we observe even in our relatively simple setting. We expect such models to be particularly fruitful avenues for future work.

One such example of more complicated models for the mechanics could include a viscoelastic or poroelastic constitutive law. As presented in Sect A.7 in S1 Appendix, the assumption of a (linear) viscoelastic constitutive law leads to qualitatively similar results to those presented in this work for a purely (linear) elastic constitutive law. Our current assumption of a linear elastic constitutive law for the mechanics of the cell is limiting as it assumes small deformations. This framework needs to be extended to study the effect of large deformations and shape changes, which would include the effect these deformations have on the signalling molecules. This would be especially interesting, as this study shows that cell shape is one of the determinants affecting the mean concentrations of the signalling molecules.

The boundary conditions for the deformation we consider correspond to simple models of a cell in vitro (flat 2D substrate) or in vivo (homogeneous 3D matrix). We see that the cell on a 2D substrate appears to spread radially with minimal deformation orthogonal to the substrate while the latter exhibits a more uniform although smaller in total magnitude 3D deformation. Differences in deformation for different environments are in line with the literature as the effect of the substrate stiffness on cells is known to vary in 2D and 3D substrates [52]. An interesting extension that could be included in the above framework would be spatial variations in substrate stiffness or more complicated models for the substrate mechanics both of which are of much biological relevance [53–55].

This work shows how mechanistic modelling of mechanotransduction can reveal remarkable emergent properties. It lays the groundwork for future studies where further complexity can be added as required to model specific signalling pathways or to reflect other mechanical models derived from different constitutive assumptions. We anticipate that choosing a viscoelastic or poroelastic constitutive law for the mechanics of the cell is an interesting direction for future studies, as this is in line with recent experimental observations [31,32]. Given the fact that cell shape greatly influences the dynamics of the cell, as shown in this work, other reference geometries are also of interest as a subject for future work. Extending the signalling model of [18] further, we intend to couple the model of this work with a similar biomechanical model for the deformation of the nucleus coupled with the dynamics of signalling molecules within the nucleus, such as the YAP/TAZ pathway [56].

Supporting information

S1 AppendixNumerical scheme and simulations.The appendix contains a comparison of the reduced model and the full model of [18] in Sect A.1, temporal statistics in Sect A.2, the numerical scheme in Sect A.3, the conversion from $[eqn]$ to $[eqn]$ in Sect A.4, simulations for the model with nucleus in Sect A.5, a parameter sensitivity analysis in Sect A.6, and simulations for viscoelastic model in Sect A.7.(PDF)

S1 FigNumerical simulation results showing $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ for the model in Eqs (3)–(6) for the lamellipodium shape and in the case of the 3D stimulus at a steady state at T = 100 s.Four different scenarios are considered: (A) C1 = 0 $[eqn]$ and Ec = 0.6 kPa $[eqn]$ ; (B) C1 = 0 $[eqn]$ and $[eqn]$ ; (C) C1 = 0.1 $[eqn]$ and Ec = 0.6 kPa $[eqn]$ ; (D) C1 = 0.1 $[eqn]$ and $[eqn]$ . Within each subfigure, the rows represent $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ on the surface of the cell, and the columns represent $[eqn]$ kPa. Parameter values as in Table 1.(TIF)

S2 FigNumerical simulation results showing $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ for the model in Eqs (3), (4), and (6) for the lamellipodium shape and in the case of the 3D stimulus at a steady state at T = 100 s.Four different scenarios are considered: (A) C1 = 0 $[eqn]$ and Ec = 0.6 kPa $[eqn]$ ; (B) C1 = 0 $[eqn]$ and $[eqn]$ ; (C) C1 = 0.1 $[eqn]$ and Ec = 0.6 kPa $[eqn]$ ; (D) C1 = 0.1 $[eqn]$ and $[eqn]$ . Within each subfigure, the rows represent $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ on the surface of the cell, and the columns represent $[eqn]$ kPa. Parameter values as in Table 1.(TIF)

S3 FigNumerical simulation results showing $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ for the model in Eqs (3), (4), and (6) for the lamellipodium shape and in the case of 2xD stimulus at a steady state at T = 100 s.Four different scenarios are considered: (A) C1 = 0 $[eqn]$ and Ec = 0.6 kPa $[eqn]$ ; (B) C1 = 0 $[eqn]$ and $[eqn]$ ; (C) C1 = 0.1 $[eqn]$ and Ec = 0.6 kPa $[eqn]$ ; (D) C1 = 0.1 $[eqn]$ and $[eqn]$ . Within each subfigure, the rows represent $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ on the surface of the cell, and the columns represent $[eqn]$ kPa. Parameter values as in Table 1.(TIF)

S4 FigSimulation results showing the mean, $[eqn]$ , min and max values of $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ as functions of substrate stiffness E, in the case of the model in Eqs (3), (4) and (6) and 2xD stimulus.We consider different couplings, four different values for C1, and two different shapes at T = 100 s by which time the results are at a steady state. All other parameter values as in Table 1.(TIF)