## Abstract

We propose a chemo-mechanical model based on stress-dependent recruitment of myosin motors to describe how the contractility, polarization and strain in cells vary with the stiffness of their surroundings and their shape. A contractility tensor, which depends on the distribution of myosin motors, is introduced to describe the chemical free energy of the cell due to myosin recruitment. We explicitly include the contributions to the free energy that arise from mechanosensitive signalling pathways (such as the SFX, Rho-Rock and MLCK pathways) through chemo-mechanical coupling parameters. Taking the variations of the total free energy, which consists of the chemical and mechanical components, in accordance with the second law of thermodynamics provides equations for the temporal evolution of the active stress and the contractility tensor. Following this approach, we are able to recover the well-known Hill relation for active stresses, based on the fundamental principles of irreversible thermodynamics rather than phenomenology. We have numerically implemented our free energy-based approach to model spatial distribution of strain and contractility in (i) cells supported by flexible microposts, (ii) cells on two-dimensional substrates, and (iii) cells in three-dimensional matrices. We demonstrate how the polarization of the cells and the orientation of stress fibres can be deduced from the eigenvalues and eigenvectors of the contractility tensor. Our calculations suggest that the chemical free energy of the cell decreases with the stiffness of the extracellular environment as the cytoskeleton polarizes in response to stress-dependent recruitment of molecular motors. The mechanical energy, which includes the strain energy and motor potential energy, however, increases with stiffness, but the overall energy is lower for cells in stiffer environments. This provides a thermodynamic basis for durotaxis, whereby cells preferentially migrate towards stiffer regions of the extracellular environment. Our models also explain, from an energetic perspective, why the shape of the cells can change in response to stiffness of the surroundings. The effect of the stiffness of the nucleus on its shape and the orientation of the stress fibres is also studied for all the above geometries. Along with making testable predictions, we have estimated the magnitudes of the chemo-mechanical coupling parameters for myofibroblasts based on data reported in the literature.

## 1. Introduction

Cytoskeletal tension is fundamental to many cellular processes such as cell motility [1], cytokinesis [2], tissue morphogenesis [3–6] and remodelling, as well as pathological processes such as tumour growth, metastasis [7] and fibrosis [8]. A number of cells such as fibroblasts, cardiomyocytes, epithelial cells, endothelial cells and smooth muscle cells develop significant contractile forces as part of their physiological function [9–12]. One important function of cellular forces is to act on the surrounding extracellular matrix (ECM), align the matrix and reorganize tissues as they form and develop. Cells modulate their contractility and adapt to different tissue-specific and extra-cellular environments in both normal and pathological settings. Cells must reorganize their shape in a highly dynamic and orchestrated manner during both homeostasis and migration. Cell contractility changes rely upon spatial and temporal coordination of biochemical and physical processes at the molecular, cellular and tissue scale. Stiffness of the extracellular environment plays a crucial role in determining the contractility of cells. Molecules such as Rho-family GTPases act as biochemical switches that couple cytoskeletal organization to distinct environmental signals and mechanics to regulate the changes in contractility. Because of the coupled nature of the mechanical and biochemical signalling processes, progress in understanding how they interact to control cell contractility and, more importantly, cell shape has proved challenging.

Predictive quantitative models that link biochemical signalling pathways that control recruitment of molecular motors leading to changes in contractility with the mechanics of the extracellular environment can provide a means to interpret and predict cell functions in normal and pathological conditions. However, most of the current models do not treat these chemo-mechanical processes in a consistent and unified manner. For a given shape, the polarization of cells in three-dimensional matrices and cells on two-dimensional substrates has been studied by accounting for the contribution of mechanical stresses to contractility [13]. Deshpande *et al*. [14] presented a general model for the contractility of cells to predict experimental observations including the dependence of cellular contractility on the substrate stiffness and the stress fibre concentration at the focal adhesions. A similar mathematical model for a single cardiac cell is used to investigate the cardiac mechanics by accounting for the dynamic reorganization of the cytoskeleton [15]. These microscopic approaches based on tension-dependent assembly of stress-fibres and cross-bridge cycling have also been implemented, but this work is also based on phenomenological approaches rather than on principles of non-equilibrium thermodynamics alone [14,15]. Macroscale models for chemo-mechanical coupling in tissues are also available [16–18], but these models do not address how mechanics of the ECM influence cell contractility and shape. A key drawback of these and other models reviewed in [19,20] is that they do not consider the chemical, mechanical and interfacial contributions to the free energy of a given shape and polarization of a cell, and thus are not able to predict which configuration (e.g. spindle shaped or spherical) is favoured from an energetic perspective. It is also well known that cells migrate towards regions of larger mechanical stiffness, a process referred to as durotaxis. However, the energetic driving forces for this process have not been quantitatively studied to date. Given the key role that cytoskeletal tension plays in physiological processes, it is important to develop a mathematical description that treats the mechanics of the cell and the ECM and the biochemistry of stress fibre formation on an equal footing.

In this paper, we develop a unified framework to study the contractility and polarization of cells as a function of the stiffness of their surroundings by considering the chemical free energy available from the assembly of stress fibres and the mechanical energy due to the deformation of cells and the ECM. The key kinematic variables that describe the cell are the strain and the contractility (the cell is polarized when this tensor is anisotropic) tensors and the free energy of the cell-matrix system is expressed in terms of these variables. We explicitly include a term that couples the stress and contractility tensors, and thus are able to capture how the contractility of the cell can change in response to mechanosensitive signalling pathways (e.g. the Rho-ROCK and the MLCK pathways). By using the principles of non-equilibrium thermodynamics, we relate the rate of myosin recruitment and cross-bridge sliding to variations of the free energy with respect to contractility and the active strain. We show that the rate of change of contractile strain is linearly related to stress, providing derivation of the Hill law [21] based on considerations of free energy, rather than phenomenology. We then show that the overall free energy of cells is lowered in stiffer environments due to recruitment of myosin motors, providing a means to evaluate the thermodynamic driving force for durotaxis. By comparing the predictions of our simulations with experiments on cells deformed using the atomic force microscope we obtain estimates of the chemo-mechanical coupling parameters. We also study how the free energy and polarization of cells depend on their shape in three-dimensional matrices, on two-dimensional substrates and on micropost arrays. The chemo-mechanical parameters we have derived can be used to quantitatively model cell sorting, motility and morphogenesis in complex three-dimensional environments.

## 2. A one-dimensional model for the chemo-mechanical free energy of the cell-matrix system

To illustrate the key ideas that underpin our chemo-mechanical formulation, we start with a simple one-dimensional example that contains all the relevant features of the full three-dimensional formulation to be presented in §3. In this model, we start with a contractile element (e.g. a cell) in its quiescent state, represented by an active component denoting the actomyosin network and a passive component that accounts for the stiffness of the elastic components (e.g. the cytoskeleton). We now apply a force (or stress) to this contractile element, which leads to changes in both the overall strain, , and the contractility of molecular motors, treated as force dipoles whose magnitude per unit volume is *ρ*_{0} in the quiescent (stress-free) state (refer to figure 1). All the symbols in the chemo-mechanical model are given in appendix A. In the one-dimensional model, *ρ* also characterizes the motor density (refer to appendix B). Our goal is to determine the change in strain and contractility, and d*ρ*, as a function of the applied stress, *σ*. To this end, we invoke the second law of thermodynamics, which states that the overall free energy should decrease
2.1
where the first two terms denote the changes in the mechanical energy and chemical energies per unit volume of the cell, respectively, and the last term is the mechanical work done by the molecular motors. While the first term depends on the strain and the second term depends on the density of molecular motors, the third term provides the coupling between the strain and contractility. The mechanical energy consists of the elastic energy of the passive elements and the mechanical work done by the external stress
2.2
where *K* denotes the elastic modulus of the passive components. To write the chemical contribution, we note that, in the quiescent state, the density of the motors attached to the cytoskeleton, characterized by *ρ*_{0}, is determined by a balance of binding and unbinding processes and minimizes the net chemical free energy. Perturbations of the contractility from this value should lead to an increase in the chemical free energy. On the other hand, with applied stress, force activated signalling pathways lead to additional recruitment of molecular motors due to an increase in the cytosolic concentration of Ca^{2+} or phosphorylated Rho through the pathways sketched in figure 1; this leads to a decrease in free energy. In other words, stress alters the equilibrium density of attached motors to a higher value; the change in free energy due to this change can be written as
2.3
where the first term ensures that the motor density in the absence of stress is the quiescent value, *ρ*_{0}, and the second term represents chemo-mechanical coupling, wherein stress favours the recruitment of molecular motors leading to a decrease in free energy. The chemo-mechanical coupling parameters *β* and *α* are related to the molecular mechanisms that regulate the engagement of motors and stress-dependent signalling pathways, respectively.

As the molecular motors are modelled as force dipoles, the work done by the motors as they deform the cytoskeleton and the ECM can be written as (refer to appendix B for a derivation of this expression and for its generalization to three dimensions). Through this term, the chemical energy associated with myosin recruitment that involves ATP hydrolysis can be converted into mechanical work. A schematic of how the total free energy is partitioned and converted between three categories is given in figure 2. Using equations (2.1)–(2.3) and the second law of thermodynamics, equation (2.1) can be cast in the form where the total free energy is 2.4

To study the *kinetics* of motor recruitment and the *dynamics* of contraction through sliding of motors, it is useful to consider the time dependence of the variation of the total free energy, written as
2.5

By choosing the rate of change of the contractility and strain in the form
2.6
where *k _{ρ}* > 0 and

*k*> 0 are the kinetic constants that govern the rates of motor recruitment and cell contraction, we ensure that the rate of change of free energy is always negative as required by the second law of thermodynamics.

_{ɛ}Next, we show that the kinetic law we derived for cell contraction is indeed the linear version of the Hill relation [21], when the rate of motor recruitment is fast, i.e. In this limit, we can equate the terms in the brackets in the first equation of (2.6) to relate the contractility to applied stress:
2.7
Using this relation in the second equation in (2.6), we get
2.8
where
2.9
where *σ*_{m} can be identified as the ‘stall stress’ and is the maximum rate of contraction. *Thus, our analysis shows that the Hill law follows as a direct consequence of the second law* and that the standard form only holds when motor recruitment is fast; in the more general case, coupled equations (equation (2.6)) govern the rate of contraction and motor recruitment.

### 2.1. Limits on the magnitudes of the chemo-mechanical coupling parameters: the stability criterion

It is well known from feedback theory of coupled systems that large values of a feedback parameter, *α*, can lead to instabilities (in the present case, divergences in the strain or contractility). In appendix C, we determine the upper and lower bounds on the chemo-mechanical feedback parameter, *α*, as well as the relation between the chemical ‘stiffness’, *β*, and the mechanical stiffness, *K*, namely *β* > *α* > 1/*K*. If the feedback parameter exceeds the chemical stiffness (note that both have the same dimension), the contractile strains increase when an external stress is applied, signalling a negative mechanical stiffness, which is unphysical. Similarly, larger mechanical feedback leads to increasing contractile strain and, unless the mechanical stiffness is large enough, the system is unstable against mechanical collapse.

## 3. A cell adherent to microposts analysed using the one-dimensional model

To illustrate the ability of our chemo-mechanical model to capture the key features of stress-driven increases in contractility, we studied the stress developed and contractility in a cell (modelled as a one-dimensional contractile element) adhered to microposts (figure 3*a*). We find that the contractility increases with increasing stiffness and for large values of chemo-mechanical feedback parameters and with decreasing values of chemical stiffness. While a larger value of the feedback provides the driving force for recruitment of motors, a larger chemical stiffness makes recruitment more difficult. Furthermore, the overall free energy of the cell–micropost system decreases with increasing post stiffness, provided that the feedback parameter is in the range identified by stability arguments. The decrease in free energy with stiffness provides the basis for durotaxis (preference of cells for stiffer surroundings) in a simple one-dimensional setting.

We solve the force balance equation: along with the geometric boundary condition: (indicating that the total strain between posts is zero), where *K*_{p} is the elastic modulus of the post, *ɛ* is the equilibrium strain of the cell and *ɛ*_{p} is the equilibrium strain of the post. This results in the relation: Substituting this relation into the one-dimensional constitutive law (equation (C2) in appendix C), we can predict the equilibrium stress, contractility in the cell and strain:
3.1
Using this solution in equation (3.1), the overall free energy of the cell–micropost system, including the energy of the cell (equation (2.4)) and the elastic energy of these two posts, can be expressed as
3.2

We have plotted these quantities as a function of the stiffness of the posts in figure 3. The key predictions of our model are discussed below.

### 3.1. Cell contractility and active stress increase with post stiffness

For the parameters that satisfy the stability criterion: *β* > *α* > 1/*K*, the equilibrium stress, *σ*, and motor density, *ρ*, increase with post stiffness, *K*_{p}, and, finally, saturate at (equation (3.3) and figure 3*b*–*c*), which is in agreement with experimental results for 3T3 fibroblasts on micro-textured surfaces or flexible cantilevers [11,12]:
3.3
Note that cell contractility stress, *σ*, is approximately proportional to post stiffness, *K*_{p}, in the low stiffness range (equation (3.3) and figure 3*b*), consistent with the previous experimental observations [10–12,27].

### 3.2. Cells tend to migrate toward stiffer surroundings

For the parameters that satisfy the stability criterion: *β* > *α* > 1/*K*, the overall free energy, *U,* decreases with post stiffness, *K*_{p} (equation (3.4) and figure 3*d*). These predictions agree with the observations that cells prefer to migrate towards stiffer regions on a substrate with a stiffness gradient, which is known as durotaxis [28–30]:
3.4

### 3.3. Effect of chemical stiffness *β* on stress, contractility and free energy

For the parameters that satisfy the stability criterion: *β* > *α* > 1/*K*, the stress and contractility decrease with *β*, while the free energy increases with *β* (equation (3.5) and figure 3*b*–*d*):
3.5
These results are consistent with the fact that small values of chemical stiffness lead to more efficient recruitment of molecular motors, which in turn results in higher contractility and smaller free energy.

### 3.4. Effect of the chemo-mechanical feedback parameter *α* on stress, contractility and free energy

For the parameters that satisfy the stability criterion: *β* > *α* > 1/*K*, the stress and contractility increase with *α*, while the free energy decreases with *α* (equation (3.6) and figure 3*b*–*d*). Unlike the chemical stiffness, large values of the chemo-mechanical feedback parameter lead to more efficient recruitment of molecular motors, which in turn results in higher contractility and smaller free energy:
3.6

### 3.5. Determination of the chemo-mechanical coupling parameters from experimental data

Recently, Mitrossilis *et al.* [12] studied the force exerted by a myoblast held between a rigid plate and a deformable cantilever as sketched in figure 4. The steady-state force *F* exerted by the myoblasts as a function of the stiffness of the cantilever, *k* [12], is shown in figure 4. Using equation (3.1), our one-dimensional model predicts that
3.7
describes the ratio between increases in cell force and increases in the cell length *rε*, where *r* is the cell size. It is an effective stiffness that accounts for contributions from both active and passive elements. (equation (3.7)) is the stall force that the cell exerts on the rigid cantilever. Fitting this relation with the experimental results, we obtain and . From these two parameters, we can obtain the chemo-mechanical coupling parameters:
3.8
For myoblasts, we estimate the passive modulus of the cytoskeleton to be *K* = 1 kPa. The quiescent contractility has been estimated in our previous work to be *ρ*_{0} = 0.5 kPa [31]. The size of a myoblast is taken to be *r* = 10 µm [32]. Using these values in the above equation, we obtain *α* = 1.91 kPa^{−1} and *β* = 2.26 kPa^{−1}. These experiments can also be analysed using our three-dimensional model, but, as we show in appendix D, they give similar estimates for the chemo-mechanical coupling parameters.

## 4. A three-dimensional model for cell polarization

The one-dimensional model can be generalized to obtain the total free energy *U* for a cell in three-dimensional environments, which includes the mechanical energy *U*_{mech}, the chemical energy *U*_{chem} and the work done by the motors *U*_{motor-work}:
4.1

We assume that, in the stress-free state, the contractility tensor is isotropic, with magnitude, *ρ*_{0}. Here *σ _{ij}*,

*ɛ*and

_{ij}*ρ*denote, respectively, the stress, strain and contractility tensors (all of which are symmetric),

_{ij}*K*and

*μ*denote the bulk and the shear moduli of the cell,

*β*denotes the chemical stiffness,

*α*

_{v}and

*α*

_{d}denote the volumetric and deviatoric chemo-mechanical feedback parameters, which naturally arise in the three-dimensional model (just like the bulk and shear moduli). To understand their origin, note that the stress

*σ*, strain

_{ij}*ɛ*and contractility

_{ij}*ρ*can be written as the sum of the volumetric and deviatoric tensors, respectively: 4.2

_{ij}The volumetric and deviatoric parts of the chemo-mechanical feedback parameters relate the volumetric and deviatoric parts of the contractility to the corresponding components of the stress tensor. To study the *kinetics* of motor recruitment and the *dynamics* of contraction, the time dependence of the variation of the total free energy can be written as:
4.3
Similar to equation (2.6), the rates of change of the contractility and strain tensors are given by
4.4
where and are the kinetic constants that govern the rates of volumetric and deviatoric motor recruitment, and and are the kinetic constants that govern the rates of volumetric and deviatoric cell contraction, respectively.

When the rate of change of contractile strain and contractility vanish (), a steady condition is achieved. In this configuration, the stress and contractility are related to the strain through relations: 4.5 where the effective contractility, effective bulk modulus, effective shear modulus, effective modulus for motor density, and effective modulus for polarization, are defined as: 4.6

These relations illustrate how mechanical stresses can lead to the polarization of cells: while the contractility tensor is isotropic in the stress-free state, the chemo-mechanical coupling parameters lead to polarization of cells, as evidenced by the presence of deviatoric components in the presence of multi-axial stress states.

Similar to the one-dimensional model, the upper and lower bounds of chemical stiffness, *β*, and the volumetric and deviatoric chemo-mechanical feedback parameters, *α*_{v} and *α*_{d}, can be determined from feedback theory as shown in appendix C. The volumetric or deviatoric contractility should increase when a volumetric or deviatoric stress is applied. Similarly, a volumetric or deviatoric stress cannot lead to an increase in the volumetric or deviatoric contraction. Using the criteria that and in equation (4.6), the chemical stiffness and chemo-mechanical feedback parameters must satisfy the stability condition
4.7

## 5. Examples to validate the three-dimensional chemo-mechanical model

Next, we consider three examples to validate our three-dimensional chemo-mechanical model, namely cells on a bed of microposts, cells on two-dimensional substrates and cells in three-dimensional matrices. In all of these examples, our focus is on the elucidation of the variations of the chemical and mechanical free energies as a function of the stiffness of the posts, substrates and matrices.

### 5.1. Cell supported by deformable microposts

Motivated by previous experimental and computational work [14], we modelled the steady-state stress developed and the contractility, in a cell on a bed of microposts (modelled using plane stress boundary conditions; refer to appendix E for details), illustrated in figure 5. A cell, of length *l* = 50 µm and thickness *t* = 10 µm, is supported at the four corners by 64 linear elastic springs (eight springs per side). Elastic springs are assumed to be of length *l*_{s} = 10 µm and can rotate with the cell edges. The effect of the post stiffness is investigated by varying post stiffness *K*_{s} over the range accessible to experiments: 0.01 kPa < *K*_{s} < 0.5 kPa. The circular nucleus occupies 10% of the initial cell volume [33]. The nucleus is represented as a linear elastic material with elastic modulus *E*_{n} = 1.2 kPa [34] and Poisson's ratio *v*_{n} = 0.3. In the remaining sections, we use the same volume ratio and elastic constants for the nucleus unless otherwise specified.

The cytoskeletal deformation and contractility are interpolated using planar 4-node bilinear quadrilateral elements. The constitutive equations (equations (E 2)–(E 3) in appendix E), the equilibrium condition, and the boundary conditions constitute a well-posed boundary value problem. We implemented the constitutive equation in a user material model in the finite-element package ABAQUS. To visualize the distribution of motor density, we define the average motor density, , and the magnitude of motor polarization, *η*, as:
5.1
where *ρ*_{1}, *ρ*_{2} and *ρ*_{3} are the three principal contractility values with *ρ*_{1} > *ρ*_{2} > *ρ*_{3}. The motor polarization orientation, *ϕ*, is defined as the eigenvector direction of maximum principal contractility, *ρ*_{1}.

The steady-state distributions of the stress and motor density for selected *K*_{s} show that the volumetric stress (figure 5*a,e*) and von Mises stress (figure 5*c*,*g*) are proportional to (figure 5*b*,*f*) and *η* (figure 5*d*,*h*), respectively. For soft supports, figure 5*a* indicates that the volumetric stress is small except in the regions adjacent to the supports. As we increase support stiffness, the deflection of posts decreases and the stress in cells increases (figure 5*e*) in good agreement with the previous experiments [10–12,27] and the prediction of our one-dimensional model (§3). The average motor density, , always peaks near the constrained corners and decreases with increasing distance from these supports (figure 5*b*,*f*). For small *K*_{s} (figure 5*b*), the supports are not stiff enough to sustain high such that except near the corners. For large *K*_{s} (figure 5*f*), there is only a small deflection of posts with the significant cell shape change, which results in higher motor density, especially near the constrained corners. We can understand this distribution by noting that stiffer posts lead to larger tensile stress on the corners of the cell (figure 5*e*), which further results in an increase in motor density (figure 5*f*) through the linear relation equation (C1) in appendix C. The distribution of *η* (figure 5*d*,*h*) is very similar to that of (figure 5*b*,*f*), except that there is high *η* around the nucleus, which shows that the stress fibres are aligned around the nucleus. The stress concentration (figure 5*c*,*g*) and stress fibre alignment (figure 5*d*,*h*) in the cytoskeleton surrounding the nucleus are consistent with the recent work by Dowling *et al*. [35].

The distributions of stress and motor density in figure 5 are qualitatively the same as in the previous work by Deshpande *et al.* [14] (who did not consider the nucleus), which provides evidence for the validity of our chemo-mechanical model. More significantly, unlike previous approaches, our model can predict the free energy distributions in the cell (figure 6*a–f*); the total free energy, *U*_{total}, is the sum of chemical energy, *U*_{chem}, the strain energy and the motor potential, *U*_{strain+motor}. All the distributions show similar trends: a minimum adjacent to the supports and an increase with increasing distance from the posts. The influence of the post stiffness on the free energies of the cell and the microposts is plotted in figure 6*g*. Our simulations show that *U*_{chem} and *U*_{total} generally decrease with an increase in support stiffness, while *U*_{strain+motor} shows the opposite trend. The dependence of support stiffness can be understood by comparing the free energies of two limit cases, that of an unconstrained cell and a fully constrained cell. Using equations (4.1) and (4.5) for the soft supports, we find
5.2
while, for stiff supports, we find
5.3

For the chemo-mechanical parameters in the range required for stability (3*Kα*_{v} > 1), we conclude that, compared with the unconstrained cell, the fully constrained cell has lower *U*_{chem} and *U*_{total}, and higher *U*_{strain+motor}. Indeed, large stiffness leads to an increase in the density of recruited motors and contractility and hence a reduction in the chemical free energy in *proportion* to the chemo-mechanical feedback parameter, *α*_{v}. On the other hand, the cells cannot mechanically deform in response to the increased contractility, leading to vanishing mechanical energy in this limit. Our work also shows that, when presented with posts of varying stiffness, the cells will move towards stiffer regions. We will further explore this phenomenon of durotaxis using the next example.

### 5.2. A hemi-spherical cell on an elastic substrate

Motivated by previous experimental and computational work [35], we modelled the steady-state stress developed, *σ*, and the contractility, *ρ*, in a cell on a substrate, illustrated in figure 7. The cylindrical substrate and the spherical nucleus are modelled as linear elastic materials with Poisson's ratio 0.3. The radius and thickness of the substrate are taken to be *R*_{s} = 5*R*_{c} and *H*_{s} = 2*R*_{c}, where *R*_{c} = 10 µm is the cell radius. The bottom and lateral boundaries of the substrate are fixed. The effect of the substrate and nuclear stiffnesses are investigated by varying *E*_{s} and *E*_{n} over the range: 1.2 kPa < *E*_{s} < 30 KP, 1.2 kPa < *E*_{n} < 30 kPa. While much of the work on mechanotransduction has focused on the mechanics of the cytoskeleton and the cellular microenvironment, it is now emerging that the mechanical properties of the cell nucleus and its connection to the cytoskeleton may play a major role in cancer metastasis, as deformation of the large and stiff nucleus presents a substantial obstacle during the passage through the dense interstitial space and narrow capillaries and can also have an impact on the contractility of cells [36–38]. In what follows, we explore the consequences of varying the stiffness of the nucleus on the contractility of the cell and the shape of the nucleus.

The distributions of the stress and motor density for selected substrate and nuclear stiffnesses *E*_{s} and *E*_{n} indicate that the von Mises stress (figure 7*a*,*e*,*i*) and volumetric stress (figure 7*c*,*g*,*k*) always shows the same variations as *η* (figure 7*b*,*f*,*j*) and (figure 7*d*,*h*,*l*), respectively. As can been seen in figure 7*f* and *h*, and *η* are highest between the substrate and nucleus, as adhesion of the cell to the stiff substrate prevents significant contraction in this region of the cytoskeleton. Similarly, a stiff nucleus also prevents the contraction in the surrounding cytoskeleton, which leads to high values of and *η* (figure 7*j*,*l*). This suggests that tensile stress is concentrated in the cytoskeleton between the substrate and nucleus (figure 7*e*,*g*), or surrounding the nucleus (figure 7*i*,*k*), which further leads to the high *η* and in these regions (figure 7*f*,*h*,*j*,*l*). Nagayama *et al*. [39] found that two types of stress fibres exist on both the apical side and the basal side of adherent cells on a substrate in vascular smooth muscle cells cultured on a two-dimensional substrate. Apical stress fibres run across the top surface of the nucleus and basal stress fibres underneath the nucleus, which is in good agreement with the predicted high motor density region (figure 7*f*,*h*,*j*,*l*). Our simulations show that stiffer nuclei lead to more density of apical stress fibres (figure 7*f*,*h*). Moreover, basal stress fibres increase with the increasing stiffness of substrate (figure 7*j*,*l*). Our simulations thus provide a method to control the ratio of these two types of stress fibres by tuning stiffnesses of the nucleus and the substrate.

It is interesting to note that the free energy distribution (figure 8*a*–*c*) has trends opposite to that of volumetric stress (figure 7*c*,*g*,*k*) and average motor density (figure 7*d*,*h*,*l*). As in the case of the cell supported by microposts, a stiffer environment results in large tensile stress, high motor density and small deflection, which further leads to the lowering of the total free energy through recruitment of motors. Heat maps of the free energies of the cell and substrate as a function of substrate and nuclear stiffnesses (*E*_{s}, *E*_{n}) are given in figure 8*d*–*f*. We find that, whereas either increased substrate or nuclear stiffness leads to a decrease in *U*_{chem} or *U*_{total} (figure 8*e*–*f*), the effect is significantly amplified when both these factors are present simultaneously. An inverse trend is predicted for *U*_{strain+motor} in figure 8*d*: *U*_{strain+motor} generally increases with an increase in substrate or nuclear stiffness. This work clearly provides a thermodynamic basis for durotaxis, where cells move towards stiffer substrates by remodelling their cytoskeletal network as evidenced by an increase in the average motor density and polarization in figure 7*f*,*h* with increasing substrate stiffness.

### 5.3 Nuclear deformation and cell polarization in elongated cells

Next we consider the polarization of cells and the shape of their nuclei in three-dimensional matrices [40] as a function of the aspect ratio of the cells. For a given stiffness of the matrix, the state of mechanical stress in the cell will depend on its shape and therefore the contractility and polarization of the cell and the nuclear shape, which through the chemo-mechanical coupling parameter will also be shape dependent. Recent experiments show that the orientation and deformation of the nucleus can be regulated by the cell elongation and spreading [41]. Here we consider spherical or elongated cells (with spherical nuclei) in a three-dimensional matrix, and study how cell elongation leads to cell polarization and nuclear elongation. We model cells as prolate spheroids described by the shape (*x*/*a*)^{2} + (*y*/*a*)^{2} + (*z*/*b*)^{2} = 1, where *a* and *b* represent the lengths of the semi-minor and semi-major axes of the prolate spheroid, respectively. In order to compare cell shapes with different aspect ratios, *γ* = *b*/*a*, we assume that the reference volumes of the cells are the same in all the cases. The cell is embedded in a cylindrical matrix whose radius is 10 times *R*_{c}, where *R*_{c} = 10 µm is the radius of the spherical cell. The matrix is modelled as a linear elastic material with elastic modulus *E*_{s} = 30 kPa and Poisson's ratio *v*_{s} = 0.3. All the material parameters for the cytoskeleton are the same as in the previous examples, except that the quiescent contractility is taken to be *ρ*_{0} = 1 kPa.

Similar to the previous two examples, the distributions of the stress and motor density for selected cell aspect ratio *γ* indicate that the von Mises stress (figure 9*a*,*e*) and volumetric stress (figure 9*c*,*g*) always show the same variations as *η* (figure 9*b*,*f*) and (figure 9*d*,*h*), respectively. The high motor polarization *η*, localized on the side of the nucleus (figure 9*b*,*f*), indicates a high degree of stress fibre alignment in this region, which is consistent with the observed spatial distribution of actin filaments in an elongated cell [41]. We find that the average density motor is largest along the long axes of the cells and the minimum along the short axes of the cells (figure 9*d*,*h*). We can understand this from the fact that the shape anisotropies lead to a concentration of tensile strains along the long axes of the cells (figure 9*c*,*g*).

The influence of the cell aspect ratio, *γ*, on the nuclear aspect ratio and the total free energy of the cell-matrix system is plotted in figure 10. Our simulations show that the nucleus is more elongated within elongated cells, while the total free energy also decreases with increasing aspect ratios. We find that the nuclear aspect ratio can be as high as 1.63 for cell aspect ratio *γ* = 5.2, which is qualitatively consistent with experimental observations [41]. The reduction in total free energy of elongated cells is also consistent with shape transitions of contractile cells in three-dimensional matrices [40]; experiments show that cells in stiffer matrices tend to be more polarized and elongated, consistent with the fact the total chemical free energy is lower.

## 6. Summary and future work

In summary, we have developed a model that explicitly includes stress-dependent recruitment of molecular motors to develop a chemo-mechanical description of polarization and strain in contractile cells. The basic variables that describe the model are the strain tensor and the contractility tensor that depends on the spatial distribution of the molecular motors. In contrast with previous models, we do not keep track of individual fibres, but capture the effect in a coarse-grained sense; the principal direction of the contractility tensor with the largest principal value gives the orientation of the stress fibres and the difference between the largest and smallest principal values of contractility gives the polarization of the cells. We have derived the free energy of the cell in terms of these variables. The free energy includes the elastic strain energy of the cytoskeleton, the mechanical work done by the molecular motors and the chemical free energy of motor recruitment. The dynamic laws for the evolution of the active strain and the contractility tensor are derived based on the second law of thermodynamics. In this manner, we are able to recover the well-known Hill relation for active stresses in tissues, based on fundamental principles of irreversible thermodynamics rather than phenomenology.

To illustrate the power of our method, we have calculated the steady-state spatial distribution of strain and contractility tensors in (i) cells on microposts, (ii) cells on substrates, and (iii) three-dimensional matrices. Our work shows that the chemical free energy of the cell decreases with the stiffness of the extracellular environment as the cytoskeleton polarizes in response to stress-dependent recruitment of molecular motors. This observation provides the basis for durotaxis, where cells preferentially migrate towards stiffer regions of the extracellular environment. The effect of cell shape on the shape of the nucleus has also been studied. For simplicity, all the parameters in this model are assumed to be constant, and this assumption is sufficient to capture the steady-state mechanical behaviour of cells. The evolution of stress and motor density distributions can be modelled using the kinetics of motor recruitment and dynamics of contraction in equation (4.4). However, on a time scale of tens of minutes to hours, myosin motors, actin filaments and other associated proteins self-assemble to form well-organized acto-myosin stress fibres [42]. During this process, all the material parameters could depend on strain (*ɛ*) and motor density (*ρ*); incorporating these material nonlinearities will be the subject of future work. If the dynamics of cell-matrix adhesions can be included in our model, it would be possible to study cell shape changes and orientation in response to time varying loads and motility of cells. We hope to address these issues in forthcoming publications.

## Competing interests

We declare we have no competing interests.

## Funding

This work is supported by grant nos. NIH-R01EB017753, NIH-U54CA193417 and NSF-CMMI131239.

## Acknowledgement

We thank Dr Rebecca Wells and Dr Abhilash Nair for discussions and comments on the manuscript.

## Appendix A. List of symbols and definitions in the chemo-mechanical model

See table 1.

## Appendix B. Work done by molecular motors in terms of the strain and contractility tensors

Figure 2 in the main text illustrates a region of the stress-fibre network subject to internal forces that are exerted by myosin II motors bound to actin filaments. Each motor contains two head domains that walk along the actin filaments. This process generates a force couple, i.e. a pair of equal but oppositely directed forces and that act on the stress-fibre network. Here and are coordinates of the head domains for the motor labelled *k*, and is the distance between them (which is approximately 200 nm) [43], is the component of the unit vectors that denotes the orientation of the motor and is the magnitude of the force exerted by the motor. The order of magnitude of the force is approximately 1 pN [44] and the unit vector denotes the direction of the force.

Denoting the displacement field in the cytoskeleton by *u _{i}*(

*x*), the mechanical work done by the force dipole is Let there be

_{j}*N*bound motors in the region of interest that has volume

*V*. The total work done by these motors is B 1 To relate the work to the contractility tensor, note that which when used in the above equation gives the work done by motors per unit volume, i.e. the motor potential energy: B 2 where is the motor density or contractility tensor. In the above derivation, we have used the definition of strain and have assumed that the contractility tensor is symmetric. Indeed, for the total moments of the force dipoles to vanish, the relation which ensures that the contractility tensor is symmetric.

## Appendix C. Limits on the chemo-mechanical feedback parameter (one-dimensional model)

When the rate of change of contractile strain and contractility vanish (), an equilibrium or steady-state condition is achieved from equation (2.6) in the main text. In this configuration, the contractility and strain are related to the stress through relations:
C 1
These relations allow us to determine the upper and lower bounds on the chemo-mechanical feedback parameter, *α*, also the relation between the chemical ‘stiffness', *β*, and the mechanical stiffness, *K*. It is well known from feedback theory that large values of a feedback parameter can lead to instabilities (in this case divergences in the strain or contractility). First, we consider the quiescent case with *σ* = 0, where both the contractile strain and the contractility diverge when *Kβ* = 1. Thus for stability, the mechanical stiffness should satisfy the condition *K* > 1/*β*; we can understand this by noting that when the chemical stiffness *β* is large, it is difficult to change the contractility (from the value *ρ*_{0}), while this value can be readily increased when *β* is small. However, increasing contractility leads to increasing contractile strain and, unless the mechanical stiffness is large enough, the system is unstable against mechanical collapse. Thus, the mechanical stiffness should be larger than 1/*β*. Next, we consider the case when an external force is applied. On physical grounds, it is clear that both the contractility and strain should increase with increasing force, which from equation (C1) gives us the conditions *β* > *α* and *K* > 1/*α*. These conditions can also be understood in a manner similar to the arguments that provide the relationship between chemical and mechanical stiffnesses. If the feedback parameter exceeds the chemical stiffness (note that both have the same dimensions), the contractile strains increase when an external stress is applied, signalling a negative mechanical stiffness, which is unphysical. Similarly, the mechanical stiffness should satisfy *K* > 1/*α*; if this were not true, the contractility decreases when an external stress is applied, which results in an unphysical negative chemical stiffness. These conditions can also be derived by writing the stress and contractility in terms of strain:
C 2
where the effective contractility effective modulus effective modulus for motor density are given by
C 3

Using the criteria that we obtain the condition that the chemo-mechanical feedback parameter must satisfy.

## Appendix D. Determination of the chemo-mechanical coupling parameters

The estimation for the chemo-mechanical coupling parameters can be improved using a more realistic three-dimensional model. Reflecting the experimental set-up, we assume that the cell adheres between one compliant post and a rigid substrate and develops uniaxial contraction without lateral constraints. Using equation (4.5) in the main text, the relation between cell force *F* and post stiffness *k* is the same as equation (3.7) in the main text, but with a modified expression for the effective cell stiffness:
D 1
Using the fitting parameters, and and estimates *E* = 1 kPa and and Poisson's ratio of the passive element to be *v* = 0.3, we have In the remaining sections, we use these fitting parameters for the three-dimensional constitutive law of the cytoskeleton unless otherwise specified. We obtain different values of *α* and *β* using one- and three-dimensional models because the three-dimensional model considers lateral contraction of the cells that is absent in the one-dimensional model. Note, however, that the magnitudes of these parameters are comparable in the two cases.

## Appendix E. Constitutive equations for the plane stress problem

Since the lateral surfaces are traction-free, E 1 where subscript 1, 2 denotes an in-plane variable, and subscript 3 denotes an out-of-plane variable.

Substitute equation (E1) into governing equation (4.5) in the main text, E 2 E 3

## Footnotes

One contribution of 19 to a theme issue ‘Integrated multiscale biomaterials experiment and modelling: towards function and pathology’.

- © 2015 The Author(s)