## Abstract

The constitutive model presented in this article aims to describe the main bio-chemo-mechanical features involved in the contractile response of smooth muscle cells, in which the biochemical response is modelled by extending the four-state Hai–Murphy model to isotonic contraction of the cells and the mechanical response is mainly modelled based on the phosphorylation-dependent hyperbolic relation between isotonic shortening strain rate and tension. The one-dimensional version of the model is used to simulate shortening-induced deactivation with good agreement with selected experimental measurements. The results suggest that the Hai–Murphy biochemical model neglects the strain rate effect on the kinetics of cross-bridge interactions with actin filaments in the isotonic contractions. The two-dimensional version and three-dimensional versions of the model are developed using the homogenization method under finite strain continuum mechanics framework. The two-dimensional constitutive model is used to simulate swine carotid media strips under electrical field stimulation, experimentally investigated by Singer and Murphy, and contraction of a hollow airway and a hollow arteriole buried in a soft matrix subjected to multiple calcium ion stimulations. It is found that the transverse deformation may have significant influence on the response of the swine carotid medium. In both cases, the orientation of the maximal value of attached myosin is aligned with the orientation of maximum principal stress.

## 1. Introduction

As the contractile component of hollow organs such as the intestines, the airways and blood vessels, smooth muscle cells react to stimulations, such as electrical field and calcium ion (Ca^{2+}) transients, by contracting the hollow organs. Abnormal contractility of smooth muscle cells is an important cause of many diseases, such as asthma, incontinence and hypertension. For example, airway hyper-responsiveness is a characteristic of asthma and generally ascribed to the excessive contraction of airway smooth muscle cells. In past decades, significant progress has been made in developing experimental techniques to measure mechanical responses, mainly contraction, of smooth muscle cells or smooth muscle cell-based tissues under various conditions. In a recent development, Tan *et al.* [1] have measured the distribution of force exerted by a smooth muscle cell by seeding the smooth muscle cell on a bed of poly(dimethylsiloxane) (PDMS) microposts and determining the deflections of the posts. Alford *et al.* [2] have used vascular muscular thin film method to measure the dynamic stress generation during contraction of microfabricated tissues of vascular smooth muscle with differing tissue and cell-level architecture. Owing to the complexity of the problem, most experimental investigations focus on details of individual elements regulating the contractility of smooth muscle cells, but an integrative understanding of how the different regulatory elements function together remains elusive. Hence, interpretation of experimental measurements at cell or tissue level remains a challenge.

A number of constitutive/mathematical models have been developed to simulate contractile response of a cell, aiming to interpret experimental observations. Appropriate constitutive models should be capable of integrating key biochemical features into the mechanical response of a cell. Earlier attempts to model cell contractility have simply prescribed a thermal strain to cell by regarding a cell as an isotropic elastic continuum [3] or a discrete set of elastic filaments [4]. Biochemical features of the cell under stimulation are neglected. Deshpande *et al.* [5] have developed a generalized bio-chemo-mechanical modelling framework that accounts for the dynamic formulation of stress fibres and the cross-bridge cycling between the actin and the myosin filament that generates the tension. The capability of the modelling framework has been demonstrated through reasonable correlation between experimental measurements and predications for the contractile responses of smooth muscle cells, fibroblasts and mesenchymal stem cells on a bed of PDMS microposts [6], the formation of stress fibres perpendicular to the direction of cyclic stretching [7] as well as the response of cells to the stiffness of an elastic substrate [8]. The responses and remodelling of a three-dimensional cytoskeletal network in response to compression [9], shear [10] and *in situ* dynamic loading [11] have been studied numerically, based on the modelling framework.

To develop bio-chemo-mechanical-based constitutive models for smooth muscle cells, a number of attempts have been made using the four-state Hai–Murphy model [12] to describe the kinetics of myosin phosphorylation and cross-bridge interactions with thin filaments. Notably among these, the model of Bursztyn *et al.* [13] describes the intracellular Ca^{2+} concentration and stress produced by a cell in response to the depolarization of the cell membrane by considering three Ca^{2+} control mechanisms: voltage-operated Ca^{2+} channels, Ca^{2+} pumps and Na^{+}/Ca^{2+} exchangers. More recently, Murtada *et al.* [14–16] have developed a model considering the relative sliding between myosin and actin filaments in force production based on Hill's three-element model as well as the dispersion effects of contractile fibres. Stalhand *et al.* [17] have developed a three-dimensional modelling framework in which a free-energy function is proposed as the sum of mechanical energy, chemical kinetics as well as the Ca^{2+} concentration and the Hai–Murphy model is included in the developed evolution law as a special case. As the Hai–Murphy model is a biochemical model for isometric contraction of smooth muscles, these models may mainly be applicable for simulating isometric responses. The applicability of these models to isotonic contractions remains questionable. In this paper, we show that the Hai–Murphy model neglects the effect of strain rate on the kinetics of cross-bridge interactions with actin filaments for isotonic contraction of smooth muscles.

In addition to the above-mentioned attempts, Wang *et al.* [18] and Maggio *et al.* [19] have developed a modified Hai–Murphy model to simulate non-isometric contractions by integrating Huxley's sliding filament theory of muscle contraction with the Hai–Murphy model. In the model, a variable representing the distance along the thin filament to a binding site from the cross-bridge as in a Huxley cross-bridge model has been incorporated into the four-state Hai–Murphy model to form a group of partial differential equations. The model has been used to simulate airway smooth muscles and uterine smooth muscles, respectively, with reasonable correlation to experiments.

The aim of this article was to develop a constitutive model that is capable of describing the bio-chemo-mechanical features in contractility of smooth muscle cells. The article is organized as follows. Section 2 gives brief description on the four-state Hai–Murphy model. Section 3 presents a constitutive model for one-dimensional idealized smooth muscle cells, which is used to simulate shortening-induced deactivation in §5. The model is extended to two-dimensional smooth muscle cells in §4. The two-dimensional constitutive model is used to simulate swine carotid media strips under electrical field stimulation and contraction of a hollow airway and a hollow arteriole buried in a soft matrix subjected to multiple Ca^{2+} stimulations in §6.

## 2. Synopsis of the Hai–Murphy model

The key bio-chemo-mechanical process involved in contractile response of smooth muscles is described in the electronic supplementary material, appendix A. Smooth muscles take 30 times longer to contract and relax than does skeletal muscle and can maintain the same contractile tension for prolonged periods at less than 1% of the energy cost. Part of the striking energy economy of smooth muscles is that the steady-state force is maintained at relatively low levels of myosin light chain phosphorylation and cross-bridge cycling rates, which has been termed the ‘latch states’ [12]. Hai & Murphy [12] developed a four-state model for isometric contraction of a smooth muscle to describe the kinetics of the fractions of myosin in the different states, as shown in figure 1. The four-state model postulates the coexistence of latch bridge and cross-bridge, in which the latch bridge is generated by the dephosphorylation of attached phosphorylated cross-bridges. The model could be described via a set of ordinary differential equations as
2.1
2.2
2.3
and
2.4where, as shown in figure 1, [*M*], [*MP*], [*AMP*] and [*AM*] are non-negative and non-dimensional quantities with [*M*]+[*MP*]+[*AMP*]+[*AM*]=1, representing fractions of free unphosphorylated, phosphorylated, attached dephosphorylated and attached phosphorylated myosin, respectively, and *k*_{1},…,*k*_{7} are the rate constants which can be obtained by fitting the model behaviour against experimental data [12]. The value of attached cross-bridges, *η*=[*AM*]+[*AMP*], can be related to isometric contraction stress, *σ*_{o}, of a smooth muscle via
2.5where *σ*_{max} is the maximal tension corresponding to the maximal number of attached cross-bridges that is permitted by biochemistry [20]. Among all rate constants, *k*_{1} and *k*_{6} are the only Ca^{2+}-regulated rate constants. Based on the assumption that the affinities of myosin light chain kinase and myosin light chain phosphatase for detached and attached cross-bridges are similar, *k*_{1}=*k*_{6} and *k*_{2}=*k*_{5} were assumed for data fitting [12].

## 3. Constitutive model for one-dimensional smooth muscle cells

We first analyse an idealized one-dimensional smooth muscle cell, i.e. all acto-myosin stress fibres are aligned along a same direction, contracted owing to the rise in the intracellular calcium ion level.

### (a) Relation of active contraction stress versus strain rate

Numerous studies, including Murphy [21] and Bates *et al.* [22], have demonstrated that the relationship of isotonic shortening velocity and tension of smooth muscle can be described by a hyperbolic curve as the ‘characteristic equation’ of Hill [23]. However, unlike skeletal muscle, the stress–velocity relation of smooth muscle is not unique and exhibits a family of phosphorylation-dependent curves [24]. Following Hill's equation and considering the effect of phosphorylation of myosin (see the electronic supplementary material, appendix B), the relation between active contraction stress *σ* and strain rate of the smooth muscle cell under isotonic contraction can be given as
3.1where *k*_{10} is a positive rate constant, a reference strain rate, and *ξ* is the phosphorylation level of myosin, *ξ*=[*MP*]+[*AMP*]. The sign convention of solid mechanics is adopted here, i.e. negative for shortening and compression and positive for lengthening and tension. The maximum shortening strain rate, , i.e. the strain rate corresponding to unloaded contraction velocity, can be calculated by letting *σ*=0 in equation (3.1), i.e.
3.2which dictates linear relationship between the maximum shortening strain rate and myosin phosphorylation. The maximum shortening strain rate is an intrinsic property of a smooth muscle [25] and the measurement of cross-bridge cycling rate [26]. Experimental studies have demonstrated that the maximum shortening strain rate is proportional to cross-bridge phosphorylation, *ξ*, in many types of smooth muscles with different types of activation and at different temperatures [26]. Equation (3.2) is consistent with the general trend observed in experiments.

As examples, figure 2*a*,*b* shows a comparison of experimental results and predictions via equations (3.1) and (3.2) on the relationship of active contraction stress and shortening strain rate, at selected values of *ξ*, for swine carotid media [21] and bovine tracheal smooth muscle [27], respectively, with constants *k*_{10} and *ν*_{o} calibrated against experimental results. The corresponding predictions on the relation of maximum contraction strain rate and phosphorylation of myosin are compared with experimental results in figure 2*c*. The agreement between predictions via equations (3.1) and (3.2) and experimental results is very good in all these examples.

When the shortening strain rate is equal to or exceeds the maximum shortening strain rate, the cell is under unloaded state, i.e.
3.3When the smooth muscle cell is under eccentric contraction, i.e. following the experimental findings of Hanks & Stephens [28], the relationship of active contraction stress and contraction strain rate is described as a straight line with the same slope as the hyperbolic relationship at a strain rate of zero, given as (see electronic supplementary material, appendix B),
3.4The load capacity of a smooth muscle under tension can be obtained via a quick stretch test. According to Dillon *et al.* [29], the yield point for smooth muscle strips obtained from media of swine carotid arteries is 1.6 times the isometric stress *σ*_{o} under steady-state conditions. However, as the focus of this paper is related to isotonic contraction of smooth muscle cells, yielding of smooth muscle cell under tension is not included in the model.

### (b) Biochemical model

As mentioned in §2, the four-state Hai–Murphy model describes the biochemical process for isometric contraction. As shortening velocity of the smooth muscle cell increases from isometric state, it is possible that more attached cross-bridges become detached in order to reduce force production, i.e. detachment rates increase for both latch bridge and cross-bridge, as schematically shown in figure 1. Therefore, higher detachment rates need to be used in modelling in order to reflect the influence of isotonic shortening. Barany [25] has investigated the relation between ATPase activity of myosin and maximum velocity of muscle shortening via experiments on a large variety of muscles, including skeletal muscles and smooth muscles. It has been found that ATPase activity of myosin correlates proportionally with the maximum velocity of muscle contraction. Even though the original experiments of Barany [25] were conducted to relate the velocity of shortening of unloaded muscle to ATPase activity of myosin, it appears that the ATPase activity of myosin is related to the velocity of shortening with and without load. ATPase activity of myosin is the direct reflection of ATP-binding, hydrolysis and phosphate release in the cycle of attachment and detachment of cross-bridges. Therefore, the rate constants of myosin detachment increased due to isotonic shortening should be proportional to shortening velocity or strain rate, namely for , the biochemical model can be rewritten as
3.5
3.6
3.7
and
3.8where *k*_{8}, *k*_{9} are rate constants. For isometric contraction, i.e. , equations (3.5)–(3.8) are equivalent to the Hai–Murphy model. For contraction with the maximum shortening strain rate, i.e. , detachment rates increase to (*k*_{7}+*k*_{9}) for latch bridges and (*k*_{4}+*k*_{8}) for cross-bridges, respectively. The model described in equations (3.5)–(3.8) implicates that the steady-state level of phosphorylation is mainly governed by *k*_{1} and *k*_{2}, i.e.
3.9For eccentric contraction of smooth muscle cells, i.e. , a more complicated biochemical reaction such as stretch-induced calcium release may be involved [30]. As the focus of this paper is isotonic contraction of smooth muscle cells, for simplicity, the Hai–Murphy model as described in equations (2.1)–(2.4) is used for eccentric contraction. In §5, the above-mentioned active stress–velocity relation and biochemical model are used to simulate shortening-induced deactivation for one-dimensional smooth muscle.

## 4. Constitutive model for two-dimensional cytoskeleton of smooth muscle cells

We now analyse an idealized two-dimensional smooth muscle cell, i.e. stress fibres can form in any direction in two-dimensional space with equal probability.

### (a) The model

Consider a periodic media with *m* stress fibres buried in a matrix with material property equal to passive behaviour of the cytoskeletal network in each smallest periodic unit cell, namely *representative volume element* (RVE). The RVE is a microdisc with out-of-plane thickness *δ* and in-plane radius *R* for two-dimensional problems, as schematically shown in figure 3. For the three-dimensional problem, the RVE is a microsphere with radius *R*, see the electronic supplementary material, appendix C. In this section, we use bold letters to represent vector/tensor. The longitudinal axis of the *k*th, *k*∈[1,*m*], stress fibre is aligned along the radius of the RVE with unit vector in the initial configuration and in the current configuration, where **e**_{j} is the basis vector in Cartesian coordinates with *j*=1 and 2 for two-dimensional space or *j*=1, 2 and 3 for three-dimensional space based on the indicial notation. The time derivative of the macroscopic deformation gradient can be related to the microscopic spatial velocity field **v** and the local coordinate **Y** at initial configuration via volume average
4.1where the superposed dot represents the time derivative and the superposed bar the macroscopic quantities, **F** is the microscopic deformation gradient tensor, *F*_{ij}=∂*y*_{i}/∂*Y* _{j}, where **y** is the local coordinate at the current configuration and *V* is the initial volume of an RVE. By definitions of **F** and , we have
4.2where *λ*^{(k)} denotes the ratio of the current to initial length for the *k*th stress fibre, . The macroscopic nominal stress can be defined as the volume average of the microscopic nominal stress **S**
4.3If both the microscopic stress and velocity field are admissible, then it follows that the macroscopic power density can be related to the volume average of the total power at microscale [31].
4.4Let and denote macroscopic Cauchy stress and rate of deformation, respectively. As mentioned by Nemat-Nasser [31], in the general case, volume average does not hold to relate macroscopic Cauchy stress and rate of deformation to the corresponding quantities at microscale based on the stress measure in equation (4.3). However, and can be calculated as
4.5respectively, where is the microscopic spatial velocity gradient, and . Using compatibility, the true strain rate along the axis of the *k*th stress fibre reads
4.6Hence, the strain and stress within the volume of the *k*th stress fibre, *Ω*^{(k)}, in the current configuration are constant. By equilibrium, the local Cauchy stress, ** σ**, in the RVE reads
4.7where

*σ*

^{(k)}is the Cauchy stress for the

*k*th stress fibre and

*σ*_{p}is the local Cauchy stress in the matrix. The local nominal stress,

**s**, can be calculated, via

**s**=

*J*

**F**

^{−1}

**and**

*σ**s*

^{(k)}=

*σ*

^{(k)}/

*λ*

^{(k)}, as 4.8where

*s*

^{(k)}and

**s**

_{p}=

*J*

**F**

^{−1}

*σ*_{p}is the local nominal stress in the

*k*th stress fibre and matrix, respectively, and . The macroscopic nominal stress can be calculated from equation (4.3) as 4.9where

*c*

^{(k)}is the volume fraction of

*k*th stress fibre in the current configuration. Equation (4.9) satisfies the relation defined in equation (4.4). The macroscopic Cauchy stress can be calculated from equation (4.5) as 4.10For well-developed stress fibres, i.e. , we have 4.11where angle

*θ*is defined in figure 3 and the active contraction stress

*σ*could be related to strain rate via equations (3.1), (3.3) and (3.4).

The passive elasticity, provided mainly by intermediate filaments, nuclei and cell membrane, need to be included in the contractile response of a cell. As shown in equation (4.11), additive decomposition of the active stress and passive stress is assumed as the stress fibres act in parallel with intermediate filaments and the cell membrane. Here, for the sake of simplicity, the passive component of a cell is assumed to behave like a neo-Hookean solid. Following Ogden [32], we consider the multiplicative decomposition
4.12of into a volumetric part and an isochoric part . Let denotes the deviatoric right Cauchy–Green tensor and . A decoupled form of the strain energy density [33] is given by
4.13where *κ* and represent the bulk modulus and the initial elastic modulus of the passive component, respectively. The Cauchy stress in the passive component *σ*_{p} can be calculated via
4.14where *δ*_{ij} is the Kronecker delta and is the deviatoric left Cauchy–Green deformation tensor, . Note that, as the numerical studies presented in Deshpande *et al.* [5], the strains in the cell are relatively small and the neo-Hookean model for the passive elasticity suffices. When warranted, a more refined model for passive elasticity such as Holzapfel *et al.* [34] could be implemented to replace the neo-Hookean model.

The above-described constitutive model has been implemented into commercially available finite-element (FE) software ABAQUS Standard via user-defined subroutine UMAT [35] to solve plane strain/plane stress problems, see the electronic supplementary material, appendix D for details. For interpretation purpose, at each integration point, both the maximal value of attached myosin, i.e. , which occurs at orientation and the averaged value of attached myosin or phosphorylated myosin, i.e.
4.15or
4.16is calculated. The averaged orientation of stress fibres, *θ*_{ave}, is also calculated as
4.17In §6, the UMAT is used to investigate swine carotid media strips under electrical stimulation and contraction of a hollow airway and a hollow arteriole buried in a soft matrix subjected to multiple Ca^{2+} stimulations.

## 5. Simulation for one-dimensional smooth muscle cells: shortening-induced deactivation

One phenomenon has been observed for striated [36,37], heart [38] and smooth muscle [39,40], named shortening-induced deactivation, is that, when a muscle is allowed to shorten during an isometric contraction, the maximum force that redevelops after shortening is smaller than the isometric force at the same length without prior shortening. The underlying mechanism is not fully understood. Possible hypothesis would be that the deactivation is caused by a combination of disorganization of myofilaments, reattachment of cross-bridges and high cytoplasmic calcium concentration after shortening [41]. In this section, we use the one-dimensional smooth muscle model, described in §3, to simulate shortening-induced deactivation for molluscan smooth muscle ABRM [42] and Aplysia I2 smooth muscle [43], respectively. The problems are solved using the fifth-order Runge–Kutta method based ordinary differential equation solver ode45 in MATLAB with initial condition [*M*]=1, [*MP*]=0, [*AMP*]=0, [*AM*]=0 at *t*=0. The shortening velocity is fixed at zero except for the moment that the quick release is imposed. For the sake of simplicity, the passive response of the smooth muscles is ignored in one-dimensional problems.

Figure 4*a* shows the comparison of experimental result for time history of stress, reported by Ruegg [42] on a molluscan smooth muscle ABRM, and the corresponding simulation results. The stress value is normalized by the peak value in the time history of stress, occurred immediately before the quick release. The muscle was subjected to a quick release by 5% of the original length *L*_{o} at *t*=40 s during an isometric contraction, and then kept at the shortened length, see top of figure 4*a*. In the simulation, the rate constants and reference velocity are chosen as ]s^{−1} [44], see in figure 4*a*), *k*_{2}=*k*_{5}=0.1 s^{−1}, *k*_{3}=0.44 s^{−1}, *k*_{4}=0.11 s^{−1}, *k*_{7}=0.006 s^{−1}, *k*_{10}=0.5 s^{−1},*k*_{8}=*k*_{9}=1.7 s^{−1} and *ε*_{o}=0.5 s^{−1}. The time history of normalized isometric stress at the *L*_{o}, predicted by the original Hai–Murphy model or the current model with , is shown in the figure for comparison. As shown for the experimental data, the quick release made normalized stress in the muscle firstly drop to zero and then immediately redevelop from 0.6 at shortened length, which may suggest that the shortening velocity by the quick release achieved the maximal shortening strain rate of the muscle and approximate 40% attached cross-bridges, namely *η*, become detached by the quick release. Comparison between experimental data and predications by the Hai–Murphy model may suggest that the stress redeveloped after quick release is initially less than isometric stress without quick release (40 s<*t*<80 s) and is fully recovered to isometric stress at a later stage (*t*>80 s). The agreement between experiment and simulation by the current model with *k*_{8}=*k*_{9}=1.7 s^{−1} is good. To examine the effects of *k*_{8} and *k*_{9}, figure 4*b* shows predictions of the time histories of normalized stress and cross-bridge phosphorylation, *ξ*, at three selected values of *k*_{8} and *k*_{9}, i.e. *k*_{8}=*k*_{9}=0.5, 1.7, 5.0 s^{−1}. The predictions show that the normalized stress redeveloped immediately after quick release at shortened length is sensitive to the values of *k*_{8} and *k*_{9}: as *k*_{8} and *k*_{9} increase from 0.5 to 5 s^{−1}, the immediately redeveloped normalized stress decreases from 0.8 to 0.18 with *k*_{8}=*k*_{9}=1.7 s^{−1} giving the best fit to experimental data (0.6). The predicted normalized stress for the three selected values of *k*_{8} and *k*_{9} will all be redeveloped to isometric stress at a later stage (*t*>80 s). However, the value of *ξ* is not sensitive the values of *k*_{8} and *k*_{9}.

Figure 5 shows the response of an Aplysia I2 smooth muscle subjected to a quick release with shortening of 2 mm at three velocities, i.e. *V* =−0.5, −1.5 and −2.0 mm s^{−1}, respectively, see top of the figure, during an isometric contraction. The muscle was kept at the shortened length after the quick release. The experimental data for the time history of stress [43] are compared with predictions by the original Hai–Murphy model for isometric contraction and by the current model for isotonic contraction (*v*<0.0 m s^{−1}), as shown in figure 5*a*, in which the stress value is normalized by the peak value in the time history of stress. The rate constants and reference velocity are chosen as s^{−1} [44], see figure 5*b*), *k*_{2}=*k*_{5}=0.39 s^{−1}, *k*_{3}=0.2 s^{−1}, *k*_{4}=0.8 s^{−1}, *k*_{7}=0.11 s^{−1}, *k*_{8}=*k*_{9}=16 s^{−1}, *k*_{10}=0.22 s^{−1} and *ε*_{o}=19.2 s^{−1} in the simulation. Higher shortening velocity in the quick release test results in more reduction in the active stress of the smooth muscle. In all cases, the shortening strain rate induced by the quick release is greater than the maximal shortening strain rate of the muscle as the active stresses during isotonic contraction are all greater than zero. Owing to the low level of *k*_{1} and *k*_{6} at the later stage, the stress redeveloped at shortened length is lower than the predicted isometric stress of the smooth muscle at initial length, as shown in figure 5*a*. The agreement between experimental results and predictions by the current model is reasonably good. The predicted time histories of the attached cross-bridges, *η*, and phosphorylated myosin, *ξ*, are shown in figure 5*b*. Corresponding to quick release-induced reduction in active stress, quick release results in reduction of attached cross-bridges, see in figure 5*b* for *η*: a higher shortening velocity results in more reduction in *η*. However, quick release has negligible influence on phosphorylation of myosin, i.e. *ξ*.

As demonstrated by the examples in this section as well as §3, the Hai–Murphy biochemical model neglects the strain rate effect on the kinetics of cross-bridge interactions with actin filaments in the isotonic contraction of smooth muscles.

## 6. Simulation for two-dimensional smooth muscle cells

In this section, FE package ABAQUS Standard is used to solve boundary value problems based on the constitutive model for a two-dimensional cytoskeletal network, as described in §4. Six-node quadratic plane stress triangle element (CPS6 in ABAQUS notation) and plane strain triangle element (CPE6) are used to analyse plane stress and plane strain problems, respectively, under finite deformation setting, i.e. the effects of geometry changes on momentum balance and rigid body rotations are taken into account. For the sake of simplicity, we ignore the detailed arrangement of smooth muscle cells in a smooth muscle and assume the whole smooth muscle is made of a single continuous smooth muscle cell. The effect of the arrangement of smooth muscle cells will be examined in subsequent papers.

### (a) Swine carotid media strips under electrical stimulation

Singer & Murphy [45] experimentally studied the maximum rate of myosin phosphorylation and stress development in swine carotid media under isometric condition. In their experiment, strips of the swine carotid media, with the dimension of 10×2×0.3 mm, were clamped at both short ends and stimulated using a direct electrical field. The stress along the longitudinal direction at the clamped ends and myosin phosphorylation level were measured for analysis. Relating *η* to the measured stress and *ξ* to the measured myosin phosphorylation level, respectively, Hai & Murphy [12] calibrated the Hai–Murphy model against the experimental data and obtained the rate constants: *k*_{1}, *k*_{6}=0.55 for 5 s^{−1} followed by 0.3 s^{−1}, see top of figure 6*a*, *k*_{2}, *k*_{5}=0.5 s^{−1}, *k*_{3}=0.4 s^{−1}, *k*_{4}=0.1 s^{−1} and *k*_{7}=0.01 s^{−1}. Based on the initial condition [*M*]=1, [*MP*]=0, [*AMP*]=0, [*AM*]=0 at *t*=0, predictions of *η* and *ξ* by the Hai–Murphy model agreed well with the experimental results. However, the effect of transverse deformation (along free ends) in the swine carotid media was not considered in the calculation. In this paper, this effect is examined via a plane stress boundary value problem based on FE simulation. In Hai & Murphy [12], the stress along the longitudinal direction measured at the clamped ends in the experiment [45] was treated as the isometric stress in the media to calibrate the above-mentioned rate constants. However, we will later show that the stress is less than the isometric stress owing to significant rotation of the stress fibres at the clamped ends. In the FE simulation, the same rate constants mentioned above were used and constants *k*_{10} and reference stain rate were obtained via calibration against experimental results [21], i.e. *k*_{10}=0.3 s^{−1} and s^{−1} (figure 2). The passive property of swine carotid media was estimated as *κ*=7840 N m^{−2} and N m^{−2} [12]. The sensitivity of *k*_{8} and *k*_{9} is studied with *k*_{8}=*k*_{9}.

Figure 6*a* shows simulation results for time histories of averaged phosphorylated myosin, i.e. *ξ*_{ave}, and normalized averaged attached myosin, i.e. , at the centre of the strip (point A in figure 6*b*) and normalized longitudinal active stress, , at clamped ends. For swine carotid media, as reported by Rembold & Murphy [20], the maximum tension is 1.8×10^{5} N m^{−2} and . The experimental results for phosphorylation level and longitudinal active stress normalized by at clamped ends are shown in the figure for comparison [12,45]. Three selected values of *k*_{8} and *k*_{9} were used in the FE simulation, namely *k*_{8}=*k*_{9}=0, 0.1 and 1.0, respectively, in which higher value represents greater effect of contraction velocity or strain rate in biochemical reaction, see §3. As shown in the figure, the prediction on phosphorylated myosin, *ξ*_{ave}, is not sensitive to *k*_{8} and *k*_{9}, which is identical to the result for phosphorylated myosin predicted by the original Hai–Murphy model (see figure 2 of Hai & Murphy [12]). Hence, the agreement between *ξ*_{ave} and the measured myosin phosphorylation level is good. The values of both and predicted by FE simulation are sensitive to *k*_{8}, *k*_{9}, i.e. higher values in *k*_{8}, *k*_{9} correspond to lower values in and , in which the prediction for obtained at *k*_{8}, *k*_{9}=0 is identical to the result for attached myosin predicted by the original Hai–Murphy model (see figure 2 of Hai & Murphy [12]). Owing to the existence of a strain gradient in the swine carotid media, the value of is not uniform in the media except for the case with *k*_{8}, *k*_{9}=0. As shown in figure 6*b* for the contours of at *t*=40 s for *k*_{8}, *k*_{9}=0.1 and 1.0, respectively, obtained by FE simulation, the value of increases from the centre to the clamped edges owing to decreased deformation from centre to edges.

As given in equation (2.5), the value of is equivalent to the corresponding isometric stress normalized by . However, as shown in figure 6*a*, the predicted normalized longitudinal active stresses, , at the clamped ends for all the selected values of *k*_{8} and *k*_{9} are less than the corresponding , hence, less than the corresponding normalized isometric stresses. This could be explained as significant rotation of stress fibres at the clamped ends induced by transverse deformation as shown in figure 6*c* for the contour of at *t*=40 s and *k*_{8}, *k*_{9}=0.0, which is identical to the contours of both the orientation of maximum principal stress and averaged orientation of stress fibres, *θ*_{ave}, in the swine carotid media: the direction of stress fibres in the media is aligned with the direction of maximum principal stress, which is rotated significantly at the clamped ends. For brevity, the contours for *k*_{8}, *k*_{9}=0.1 and 1.0 are not shown here, as they follow the same pattern as in figure 6*c*.

The rate constants, *k*_{1},…,*k*_{7}, were obtained by calibrating *η* against the measured longitudinal active stress at the clamped ends and *ξ* against the measured myosin phosphorylation level [12]. As shown in figure 6*a*, the fit between and longitudinal active stress measured at the clamped ends is reasonably good when the values of *k*_{8} and *k*_{9} vary between 0 and 0.1. However, the predicted longitudinal active stress at the clamped end is less than the measured values. The results presented in figure 6 suggest that transverse deformation may have significant influence on the response of the swine carotid media.

### (b) Contraction of an airway and an arteriole subjected to multiple stimulations

Several authors have experimentally examined the contraction of hollow organs subjected to Ca^{2+} stimulation [47,48], which has motivated the FE simulation on the contraction of a hollow airway, with internal radius 100 mm and external radius 140 mm, and a hollow arteriole, with internal radius 40 mm and external radius 66.6 mm, buried in a 400×640 mm passive matrix material, clamped at two opposite edges with remaining edges free, subjected to multiple Ca^{2+} stimulations, as shown in the electronic supplementary material, figure E1*a* for the geometries and boundary condition at the initial configuration. It is assumed that the out-of-plane dimension is much longer than the in-plane ones. Therefore, the problem can be simplified as a plane strain boundary value problem. No attempt was made to calibrate material properties or intracellular Ca^{2+} concentration against existing experimental studies in this paper. The passive matrix material was treated as a neo-Hookean solid with *κ*=7800 N m^{−2} and representing soft tissues such as lung. The material properties of bovine tracheal smooth muscle [12,27] and (figure 2) are used in the simulation for both the airway and arteriole, i.e. *k*_{2},*k*_{5}=0.1 s^{−1}, *k*_{3}=0.44 s^{−1}, *k*_{4}=0.11 s^{−1}, *k*_{7}=0.005 s^{−1}, *k*_{10}=0.31 s^{−1} and . The rate constants of *k*_{1}, *k*_{6} are chosen to represent three consecutive stimulations, namely
6.1as shown in the electronic supplementary material, figure S1E*b*. The passive property of bovine tracheal smooth muscle was assumed to be same as that of swine carotid media, i.e. *κ*=7840 N m^{−2} and . The sensitivity of *k*_{8} and *k*_{9} is studied with *k*_{8}=*k*_{9}.

The time evolutions of relative areas, defined as the area of the hollow organ at the current configuration normalized by that at the initial configuration, of the airway and the arteriole are shown in electronic supplementary material, figure 1E*b* and *c* for selected values of *k*_{8}, *k*_{9}. The contours of the maximum principal true strain at *t*=300 s for selected values of *k*_{8}, *k*_{9} are shown in the electronic supplementary material, figure E2. These figures show (i) lower values of *k*_{8}, *k*_{9} correspond to more reduction of the areas in the hollow organs, (ii) if the same values of *k*_{8}, *k*_{9} were applied to the airway and the arteriole, the relative areal changes of the airway and the arteriole are almost identical, (iii) the deformed shapes of hollow organs may be sensitive to the strain field in the matrix material and (iv) even though the values of *k*_{1}, *k*_{6} decline to zero after each stimulation, these hollow organs still keep in contracted state.

Electronic supplementary material, figure E3 shows the contours of orientation of the maximum principal stress, the orientation of , and the averaged orientation of stress fibres, *θ*_{ave}, in the airway and the arteriole at *t*=300 s and *k*_{8},*k*_{9}=1. For brevity, the contours for other values of *k*_{8}, *k*_{9} are not shown here, as they follow a similar pattern. These results all suggest that the orientation of the maximum principal stress is almost identical to at steady state. However, *θ*_{ave} is not consistent with the orientation of the maximum principal stress.

## 7. Discussion

The model presented here accounts for the nonlinear coupling between intracellular calcium ion signalling, phosphorylation-dependent contractility and the kinetics of myosin phosphorylation and cross-bridge interactions with thin filaments. The critical features of this model are as follows:

(i) It captures the main bio-chemo-mechanical features involved in isotonic contraction of smooth muscle cells, including latch states, phosphorylation-dependent active stress–strain rate relation as well as the strain-rate-dependent kinetics of cross-bridge interactions with thin filaments.

(ii) The cytoskeletal stress fibre network is allowed to evolve based on bio-chemo-mechanical boundary conditions applied to the smooth muscle cells, which enable simulation on dynamic remodelling of cytoskeleton owing to altered boundary conditions.

The critical features enable the model to predict a wide range of experimentally observed phenomena in smooth muscle cells or smooth muscle cells-based tissues. However, a note of caution is appropriate. The model is highly nonlinear and results depend on the choice of parameters. The majority of the rate constants in the current model is inherited from the Hai–Murphy model (*k*_{1},…,*k*_{7}). Systematic studies have been conducted to characterize the effects of these rate constants [12,46] and [20,26,49], which provide a basis for the choice of the rate constants and the interpretation of simulation results for the current model. Further comparisons between simulations and experiments are also needed to justify and validate the microdisc or microsphere RVE-based homogenization approach used for the two-dimensional model or the three-dimensional model.

## 8. Concluding remarks

A group of constitutive equations have been presented for cytoskeletal contractility of idealized one-dimensional smooth muscle cells, which capture the main features of bio-chemo-mechanical responses induced by a rise in the intracellular calcium ion level. The biochemical model is obtained by extending the four-state Hai–Murphy model to isotonic contraction of a smooth muscle and the mechanical model accounts for phosphorylation-dependent active stress–strain rate relation. For one-dimensional numerical examples, the constitutive equations are used to simulate shortening-induced deactivation for smooth muscles. The results obtained by simulation agree well with the experimental measurements reported in the references. The results suggest that the Hai–Murphy biochemical model neglects the strain rate effect on the kinetics of cross-bridge interactions with actin filaments in the isotonic contraction of smooth muscles.

A model is developed to extend the one-dimensional constitutive equations for two-dimensional and three-dimensional cytoskeletal networks of smooth muscle cells. The two-dimensional version of the model has been incorporated, as the user-defined subroutine (UMAT), into the commercial FE package ABAQUS Standard, which is used to investigate swine carotid media strips under electrical stimulation as a plane stress problem. The results obtained by FE simulation are compared with experimental measurements and analytical predictions by the Hai–Murphy model [12]. It is found that the direction of the stress fibres (attached myosin) in the media is aligned with the direction of maximum principal stress, which rotates significantly at the clamped ends. Hence, transverse deformation may have significant influence on the bio-chemo-mechanical response of the swine carotid media. FE simulation is also used to simulate contractions of a hollow airway and a hollow arteriole buried in a soft matrix subjected to multiple Ca^{2+} stimulations as a plane strain problem. The simulation results reveal that the deformed shapes of hollow organs may be sensitive to the strain field in the matrix material and, at steady state, the orientation of the maximal value of attached myosin, , is almost identical to the orientation of the maximum principal stress.

The model is capable of simulating dynamic remodelling of cytoskeletal network in response to altered bio-chemo-mechanical boundary conditions. It could be used to design and interpret appropriate experiments.

## Acknowledgements

The author is grateful for access to the University of Nottingham High Performance Computing Facility.

- Received November 19, 2013.
- Accepted January 3, 2014.

© 2014 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0/, which permits unrestricted use, provided the original author and source are credited.