## Abstract

Isolated cardiac myocytes exhibit spontaneous patterns of rhythmic contraction, driven by intracellular calcium waves. In order to study the coupling between spatio-temporal calcium dynamics and cell contraction in large deformation regimes, a new strain-energy function, describing the influence of sarcomere length on the calcium-dependent generation of active intracellular stresses, is proposed. This strain-energy function includes anisotropic passive and active contributions that were first validated separately from experimental stress–strain curves and stress–sarcomere length curves, respectively. An extended validation of this formulation was then conducted by considering this strain-energy function as the core of an integrated mechano-chemical three-dimensional model of cardiac myocyte contraction, where autocatalytic intracellular calcium dynamics were described by a representative two-variable model able to generate realistic intracellular calcium waves similar to those observed experimentally. Finite-element simulations of the three-dimensional cell model, conducted for different intracellular locations of triggering calcium sparks, explained very satisfactorily, both qualitatively and quantitatively, the contraction patterns of cardiac myocytes observed by time-lapse videomicroscopy. This integrative approach of the mechano-chemical couplings driving cardiac myocyte contraction provides a comprehensive framework for analysing active stress regulation and associated mechano-transduction processes that contribute to the efficiency of cardiac cell contractility in both physiological and pathological contexts.

## 1. Introduction

Rhythmic cardiac myocyte contractions are controlled by intracellular Ca^{2+} oscillations, which are triggered by membrane depolarization under normal conditions (Kort *et al.* 1985; Stern *et al.* 1988; Viatchenko-Karpinski *et al.* 1999; Vinogradova *et al.* 2005). This depolarization induces trans-sarcolemmal Ca^{2+} influx that generates synchronous Ca^{2+} release from the sarcoplasmic reticulum (SR) to bring about a uniform rise in cytosolic Ca^{2+} (calcium-induced calcium release (CICR) mechanism; Fabiato 1983). This autocatalytic process not only triggers intracellular calcium release, but also the propagation of calcium waves that modulate the overall cell mechanical contraction. A large number of experimental studies attempted to clarify this excitation–contraction coupling, and especially the characteristics of Ca^{2+} oscillations and wave propagations using single-cell preparations (Capogrossi & Lakatta 1985; Sasse *et al*. 2007; Rapila *et al.* 2008). Such experimental results stimulated different theoretical studies, which provide a mechanistic explanation for these spatio-temporal phenomena, as well as realistic simulations of calcium intracellular dynamics within simple or more-refined cellular geometries (Means *et al.* 2006; Dupont *et al.* 2007; Solovey *et al.* 2008). However, an essential requirement for understanding the regulation of cardiac cell contractility is the ability to couple mechanical and biochemical processes within integrative mechano-chemical models of the cell, which will offer a comprehensive view of rather large sets of experimental data, while providing a reliable basis for extensive simulations predicting the outcome of any model parameter on cardiac cell contractility.

Thus, even if it only represents one part of the global electro-mechano-chemical regulation of heart dynamics (Campbell *et al.* 2008), modelling cardiac myocyte mechano-chemical couplings is of particular importance because the analysis of cardiac cell response to mechanical stress may shed light on various pathologies of the rhythmic contractility of the cardiac tissue. One particular aim is to integrate comprehensive experimental data for cardiac myocyte structure and interactions between mechanical and biochemical processes within appropriate constitutive models, which may provide a relevant basis for finite-element analysis computation and prediction of cell contraction in varying conditions. Recent work has analysed the theoretical properties of isotropic continuous contracting excitable media driven by travelling nonlinear waves of excitation at a tissue level (Panfilov *et al.* 2005, 2007), but only a small number of models of the contractile cardiac cell have been proposed in a continuum mechanics framework (Okada *et al.* 2005; Tracqui *et al.* 2008). Indeed, most of the developed modelling approaches usually consider mechanical couplings as a superimposition of different subsets of phenomenological relationships (Stuyvers *et al.* 2002; Niederer *et al.* 2006), instead of considering more fundamental constitutive laws that may be applicable to real three-dimensional cell geometries.

Therefore, this paper proposes a global description of cardiac cell contractile behaviour, which integrates in a single strain-energy function both passive and active components of cell biomechanical response. Indeed, such formulation ensures a natural way of dealing with the complexity of continuum media deformations without introducing virtual intermediate configurations (Stalhand *et al.* 2007). Thus, we have developed an integrated and constitutive model of adult cardiac myocyte contraction in which nonlinear cell elasticity properties and sarcomere-length dependence of contractile intracellular stresses have been gathered within a global strain-energy function that accounts for both contractile and cytoskeleton anisotropies. Interestingly, this approach allows us to provide a representative mechano-chemical model of the cardiac myocyte that is able to account for the main features of the three-dimensional cell contraction. In our model, mechano-chemical couplings were defined by using the simplified two-variable model of Goldbeter *et al.* (1990), which was considered as the biochemical kernel generating the different spatio-temporal patterns of intracellular calcium dynamics. The contribution of cardiac myocyte elasticity to the strain-energy function was derived from the experimental stress–strain curves obtained from passively stretched rat cardiac myocytes (Cazorla *et al.* 2003). Active mechanical properties have been originally introduced in the strain-energy function by considering stress–sarcomere length curves recorded during the isometric contraction of rat cardiac myocytes immersed in solutions buffered with different Ca^{2+} concentrations (Weiwad *et al.* 2000).

We will first show that the active stress tensor derived from the proposed energy function fits very well the successive sets of experimental data obtained by Weiwad *et al.* (2000) when increasing extracellular Ca^{2+} concentrations. Considering then a simplified but realistic cylindrical cardiac myocyte geometry, the strain-energy function was introduced as the core of a three-dimensional finite-element model simulating the progressive contraction of the cardiac myocyte from part to part in response to a twitch stress propagation originating locally either from the cell extremities or close to the cell centre. The successful comparison of such simulated behaviours with the contraction patterns exhibited by spontaneously contracting adult rat cardiac myocytes provides a satisfactory validation of the proposed mechano-chemical model. More fundamentally, the formalization of the explicit coupling between biochemical and biomechanical intracellular processes also suggests that the proposed hyperelastic description of the contractile cell behaviour can be used for describing stress and strain distributions at the sarcomere level under a wide range of physiological and pathological conditions affecting the efficiency of cardiac myocytes contractility, and thus the heart performance.

## 2. Model formulation

### (a) Strain-energy function modelling cardiac myocyte elasticity

In the light of the nonlinear stress–strain relationships revealed by uniaxial stretch test measurements (Cazorla *et al.* 2003), the cardiac myocyte could be considered as a hyperelastic, yet homogeneous, continuum. However, micromanipulations using magnetic tweezers and performed on adult mouse cardiac myocytes show that the specific alignment of the sarcomeres is associated with a significant mechanical anisotropy of the cell (Lammerding *et al.* 2003). With this technique, Lammerding *et al.* (2003) reported that the cardiac cell cytoskeleton appears to be almost twice as stiff in the sarcomere direction (principal cell axis) than in the transverse direction.

Among the various anisotropic strain-energy functions that can describe such a mechanical medium property (Holzapfel 2001), and assuming a quasi-incompressibility of the cell, we retained the following three-parameter strain-energy function *W*_{p}, previously suggested for modelling the passive myocardium mechanical behaviour by several groups (Humphrey *et al.* 1990*a*,*b*; Sacks 2003), which reads
2.1

In equation (2.1), *a*_{1}, *a*_{2} and *a*_{3} are three material constants, *I*_{1} is the first invariant of the right Cauchy–Green strain tensor **C** (*I*_{1}=Trace(**C**)) and the pseudo-invariant *I*_{4} accounts for the anisotropic behaviour of the cell (Holzapfel 2001). The pseudo-invariant *I*_{4} is given by *I*_{4}=**f**_{0} **C** **f**_{0}, where **f**_{0} is the unit vector in the sarcomeres direction at zero stress state in the reference configuration.

### (b) Strain-energy function for cardiac myocyte Ca^{2+}-dependent contraction

#### (i) Active energy function and derivation of active stress

In order to consider Ca^{2+}-dependent active contraction in any deformed configuration of the cell, we introduced the active strain-energy function *W*_{act}(*I*_{4},Ca) given by (Lin & Yin 1998; Usyk *et al.* 2001; Bourdarias *et al.* 2003)
2.2
with *I*_{4}=*λ*^{2}, where *λ* defines the sarcomere elongation along the deformed fibre direction specified by the unit vector **f**_{s}. Hence, the anisotropic contractile Cauchy stress tensor *τ*_{act} is defined by
2.3
where ⊗ denotes the dyadic or tensor product. The amplitude *T*_{act}(*λ*,Ca) of the active stress is given by
2.4
where *J*=det(**F**) is the determinant of the deformation gradient tensor **F** characterizing the medium elastic deformation from the reference to the deformed configuration. Substituting the expression of *W*_{act}(*I*_{4},Ca) in equation (2.2) gives
2.5
Assuming cardiac myocyte incompressibility (*J*=1), the amplitude of the contractile stress tensor is thus related to the sarcomere elongation *λ* and Ca^{2+} concentrations by the relationship
2.6

#### (ii) Modelling sarcomere-length-dependent activation

It is known from experimental data that cardiac myocyte contraction depends on the local cytosolic Ca^{2+} concentrations, which regulate the conformational change of myosin heads binding to actin filaments, and thus the mechanical tension produced by pulling myosin and actin filaments in opposite directions. In addition, cardiac myocyte exhibits a feedback regulation process that adapts the tension-generating capacity of the cell to variations of the sarcomere length SL from its resting length SL_{0}. This sarcomere-length dependence of force development, which underlies the relationship commonly referred to as the Frank–Starling mechanism of heart contraction adaptation (Levy *et al.* 2005), is characterized by asymmetric bell-shaped curves when the contractile stress is plotted as a function of sarcomere length for given Ca^{2+} concentrations: first, the stress increases almost linearly with SL, then decreases abruptly to zero (Weiwad *et al.* 2000).

In order to propose an original formulation of the active strain-energy function *W*_{act}(*λ*,Ca) that includes these regulatory processes, we first start from a decoupled formulation of the amplitude *T*_{act}(*λ*,Ca) as a product of two independent functions of the cytosolic Ca^{2+} concentrations and sarcomere elongation *λ*=SL/SL_{0}, respectively,
2.7
The second term on the right-hand side corresponds to the usual description of sigmoid-like experimental force Ca^{2+} curves by Hill’s function, with half-maximal activation concentration Ca_{50} and Hill’s exponent *n**H* (Weiwad *et al.* 2000). The possible dependency of Ca_{50} and *n**H* with respect to *λ* has been neglected in a first approximation.

Second, we make explicit the function in equation (2.7) by introducing an original function that specifies the nonlinear variation of the maximal active-stress with sarcomere length as
2.8
where *a* (in Pa) and *b* are positive numbers, while *λ*_{R}=SL_{R}/SL_{0} corresponds to the sarcomere extension at the maximal activation state that occurs for some sarcomere length SL_{R}. These parameter values will be identified from experimental data, as further reported in §4. Considering equations (2.6) and (2.7), it is then possible to propose a general formulation of the active strain-energy function *W*_{act}(*λ*,Ca).

When *λ*=1, *C*_{act}(1,Ca) has a finite value and *W*_{act}(*λ*,Ca) is equal to zero. When *λ*≠1, substituting the expression of *T*_{act}(*λ*,Ca) into equation (2.6) and solving the corresponding differential equation for *C*_{act}(*λ*,Ca) gives, after some algebra,
2.9

The energy function *W*_{act}(*λ*,Ca) is then given by substituting the above expression in equation (2.2). The continuity of the energy function *W*_{act}(*λ*,Ca) for the value *λ*=1 is insured by choosing the integration constant *k* such that *W*_{act}(1,Ca)=0. Thus, the active energy function we propose reads
2.10

Finally, the Cauchy stress tensor ** τ** within the cardiac myocyte can be derived directly from the global nonlinear strain-energy function

*W*(

*λ*,Ca)=

*W*

_{p}(

*λ*)+

*W*

_{act}(

*λ*,Ca) gathering passive and active contributions, such as 2.11 where

*p*is the Lagrangian multiplier resulting from the considered incompressibility of the cardiac myocyte, while

**I**,

**F**and

**E**are the identity, deformation gradient and Green–Lagrange strain tensors, respectively.

### (c) Modelling intracellular calcium waves

The experimentally observed spontaneous and rhythmic contraction of an isolated cardiac myocyte is induced by the well-documented autocatalytic CICR process, which is triggered by the gating of ryanodine receptors, and which insures a massive release of calcium from the SR (Fabiato & Fabiato 1975). Cell contraction then proceeds, followed by a relaxation phase insured by Ca^{2+} pumping back into the SR, which makes this calcium pool available for the next cellular contraction.

We considered the simplified two-variable model of Goldbeter *et al.* (1990) and the subsequent coupling with diffusion processes (Dupont & Goldbeter 1992, 1994) as a biochemical core for the generation of these self-sustained intracellular calcium oscillations and the associated propagating calcium waves. The model accounts for the ability of excitable cells to generate dynamical instabilities and different types of propagating waves, and especially solitary waves (Sneyd *et al.* 1993).

Briefly, this biochemical model considers as variables the cytosolic *Z*(**x**,*t*) and sarcoplasmic *Y* (**x**,*t*) calcium concentrations, which vary simultaneously according to the balance between different calcium fluxes involving: (i) Ca^{2+} pumped into (flux *ν*_{2}) and released (flux *ν*_{3}) from the SR and (ii) flux *ν*_{1}*β* of calcium into the cytosol, basal leak *k*_{f}*Y* of Ca^{2+} from the SR and inward *ν*_{0} and outward *k**Z* calcium fluxes taking place through specific calcium channels at the cell membrane.

In addition to these biochemical processes, anisotropic diffusion of cytosolic Ca^{2+} (Engel *et al.* 1994) has been introduced on the basis of the measurements of Subramanian *et al.* (2001), who reported that Ca^{2+} diffusion is twice as fast along the sarcomere orientation (principal cell axis) than in the transverse direction. Accordingly, we considered a diagonal diffusion tensor **D**=diag(*D*_{xx},*D*_{yy},*D*_{zz}), where the ratios *D*_{yy}/*D*_{xx} and *D*_{zz}/*D*_{xx} between the diagonal tensor components are equal to 0.5, with the non-diagonal components being zero. In the deformed configuration domain, the submodel of calcium dynamics we considered is thus given by the following system of nonlinear partial differential equations:
2.12
and
2.13
where **v** is the spatial velocity of convective intracellular calcium fluxes generated by cell deformations. The calcium fluxes *ν*_{2} and *ν*_{3} are modelled by the following enzymatic Michaelis–Menten reaction terms (Goldbeter *et al.* 1990):
2.14
in which integers *n*, *m* and *p* represent the cooperativity degrees of the biochemical processes, while the positive constants *V*_{M2} and *V*_{M3} define the maximum reaction rates and *K*_{2}, *K*_{3} and *K*_{4} are affinity constants (see table 1).

### (d) Finite-element simulations

The preceeding equations were numerically solved in a three-dimensional domain defined from adult rat cardiac myocyte morphology (figure 1). The real rod-like three-dimensional cell morphology has been approximated by a cylinder with a typical cell length of approximately 110 *μ*m and a cell mean cross-sectional diameter of approximately 27 *μ*m. This simplified three-dimensional geometry provides a relevant compromise between reasonable computation time and realistic simulations of single-cell contraction initiated by calcium spark originating from any part of the cell. In particular, cell-bending motions initiated by intracellular torques can be simulated. In addition, for computational reasons, the cell nucleus has not been introduced as a supplementary subdomain in our present computations. Assuming that inertial effects and volumic forces are negligible, the mechanical balance equation we solved on the so-defined three-dimensional deformed configuration domain is
2.15
where the expression of the Cauchy stress tensor ** τ**(

*λ*,

*Z*) is given by equation (2.11), with Ca

^{2+}concentrations corresponding to the variable

*Z*(

**x**,

*t*). Thus, the value of the intracellular Cauchy stress is obtained by solving equations (2.11)–(2.13) simultaneously in order to compute the spatio-temporal evolution of

*W*

_{act}(

*λ*,

*Z*). In the reference configuration, the cardiac myocyte is assumed to be in a stress-free and strain-free states. Initial conditions for calcium concentrations in the resting cell were defined by a cytosolic Ca

^{2+}value of 0.1

*μ*M for

*Z*(

**x**,

*t*

_{0}), with an associated value of 1.6

*μ*M for initial sarcoplasmic concentrations

*Y*(

**x**,

*t*

_{0}).

As boundary conditions, we considered a cell stuck on its left (or bottom) side on a plastic dish, and thus zero displacement conditions were imposed on this boundary. Stress-free conditions were prescribed on all other cell boundaries, together with zero-flux (Neumann) conditions for variables *Z*(**x**,*t*) and *Y* (**x**,*t*). The generation of calcium sparks in different locations of the cardiac myocyte was simulated by applying a local increase of parameter *β* from its initial zero value. The so-defined global nonlinear partial differential system was numerically solved by a finite-element method (Comsol Multiphysics software) in a three-dimensional domain meshed by 446 nodes and 1599 triangle elements.

## 3. Time-lapse videomicroscopy of cardiac myocytes contraction

To get a realistic quantitative comparison with real cardiac myocyte contraction dynamics, time-lapse videomicroscopy observations were conducted on cardiac myocytes kindly provided by the Laboratoire de Bioénergétique Fondamentale et Appliquée (J. Olivares, Université Joseph Fourier, Grenoble). Cardiomyocyte contractions were observed at room temperature. The spontaneous contractions of the isolated rat cardiac myocytes have been recorded using an imaging workstation composed of an inverted microscope (Zeiss Axiovert 135) equipped with a Phase 1 Achrostigmat 5X and 10X objectives, automated shutters (Uniblitz) and a charge-coupled device camera (CoolSNAP ES, Roper Scientific) monitored via a computer by image acquisition and analysis software (Meta Vue, Roper Scientific) given in the electronic supplementary material, S1.

## 4. Results

### (a) Identification of elasticity constants for adult cardiac myocyte

According to our working hypothesis, the passive cell elasticity modulus can be derived from the three material constants *a*_{1}, *a*_{2} and *a*_{3} previously introduced in the strain-energy function defined by equation (2.1). This was done through a data-fitting analysis of previously published experimental data on a single cardiac myocyte stretching along the sarcomere orientation and where a highly nonlinear relationship between the sarcomere extension and the imposed external tension has been reported (Cazorla *et al.* 2003). However, the lack of data regarding the cell response in a direction transverse to the sarcomere orientation may hinder the determination of the coefficients *a*_{i}. We therefore took advantage of the results previously obtained by Lammerding *et al.* (2003) on the cardiac myocyte probed by magnetic beads, who estimated that Young’s modulus of the cardiac myocyte in the sarcomere alignment direction is twice that measured in the transverse direction. Considering that this anisotropy ratio holds whatever the sarcomere elongation is defines a constrained parameter space within which we identified optimal values of *a*_{1}, *a*_{2} and *a*_{3} through a least mean-squared error data-fitting procedure. Moreover, restrictions on rheological parameters, such as Baker–Ericksen’s inequality discussed in detail in Humphrey *et al.* (1990*a*), have been incorporated as additional constraints for the optimization process.

As a result, a very satisfactory fit of the experimental data was obtained with parameter values *a*_{1}=1.631 kPa, *a*_{2}=6.949 kPa and *a*_{3}=0.280 kPa (figure 2), which then makes possible an extended description of the passive nonlinear anisotropic strain–stress response of the cardiac myocyte in large deformation domains. Furthermore, it is then possible to infer the variations of Young’s moduli values with the extension ratio when considering the cell response to uniaxial traction tests performed both along and across the sarcomere orientation.

#### (i) Identification of an active strain-energy function including sarcomere-length-dependent activation

Different studies have been conducted to determine the full range of SL-active tension curves for cardiac myocytes, especially up to the point where active tension development was zero. We concentrated here on the data reported by Weiwad *et al.* (2000) in experiments conducted on skinned rat cardiac myocytes.

In these experiments, cardiac myocytes, activated with increasing values of extracellular Ca^{2+} concentrations, exhibited a typical sarcomere-length dependence on force development with a steep force rise for SL values between 2.2 and 2.3 *μ*m. At *p*Ca=4.9 (), full activation of the cell contractile apparatus occurred, whereas no activation took place at *p*Ca=7.0.

In order to get an active strain-energy function *W*_{act}(*λ*,Ca) that reproduces this adaptive contraction mechanism, we fit the overall set of data points reported by Weiwad *et al.* (2000) for *p*Ca ranging from 4.9 to 5.7 with the expression of active stress *T*_{act}(*λ*,Ca) proposed in equations (2.7) and (2.8). The least-squared fits to the experimental SL-active tension data points at different Ca^{2+} activation levels are shown in figure 3, for a resting sarcomere length SL_{0}=1.9 *μ*m and a Hill exponent *n**H*=2.6, in agreement with the experimental data (Weiwad *et al.* 2000).

Let us remark that we did not analyse the three sets of data points separately. Instead, the set of fitting curves presented in figure 3 correspond to results obtained when minimizing the least-squared error over all the three sets of data points, with the same set of unknown parameters *a*, *b* and *λ*_{R}. With this procedure, a quite satisfactory fit was obtained with the parameter values *a*=55.33 kPa, *b*=0.0371 and *λ*_{R}=1.179, corresponding to SL_{R}=2.24 *μ*m. The value Ca_{50} was left as a free parameter to check the validity of the fitting procedure. One indeed obtained an identified optimum value Ca_{50}=3 *μ*M, which is quite close to the value Ca_{50}=3.27 *μ*M identified by Weiwad *et al.* (2000) from experimental Ca^{2+}-active tension curves. As expected, a slightly better fit was obtained if each experimental curve is analysed separately.

Owing to this very satisfactory description of the experimental SL-tension curves by the proposed relationship (2.8), we then plotted in figure 4 the active strain-energy function *W*_{act}(*λ*,Ca) for *p*Ca values equal to 4.9, 5.46 and 5.7. This function is nonlinear, tending to plateau for sarcomere elongation above 1.25.

### (b) Finite-element simulations of different rhythmic contraction patterns

The relevance of the proposed active strain-energy function was further investigated by considering its coupling with nonlinear calcium dynamics within a three-dimensional model able to generate spontaneous and self-sustained rhythmic contractions of isolated cardiac myocytes. For this purpose, we considered a parametrization of the intracellular calcium fluxes that ensures the existence of solitary calcium waves when localized perturbations of intracellular Ca^{2+} concentrations are imposed in different intracellular locations. Such perturbations correspond to the experimentally observed calcium sparks attributed to the local and concerted opening of a small number of ryanodine receptors (Sheehan *et al.* 2006; Cheng & Lederer 2008). Kinetic parameter values used in the simulations (table 1) have been taken in the range of values reported by Dupont & Goldbeter (1992, 1994) and chosen in order to get cell contraction features that correspond to the measurements extracted from time-lapse videomicroscopy experiments (Pustoc’h *et al.* 2005). Typically, contraction periodicity is approximately 10–18 s, with a contraction/relaxation phase lasting approximately 0.9–1.5 s and producing a typical cell shortening of 4–10% of cell length.

In the light of experimental observations, different contraction patterns have been simulated (see electronic supplementary material, S2). First, when a localized Ca^{2+} spark originated from the left-hand side of the cell, it triggers a Ca^{2+} wave that spreads from left to right (figure 5). The cell contraction proceeds according to two intricate phases mixing contraction/relaxation on the left side and contraction/relaxation of the right side of the cell (figure 5*a,c*). Notably, asymmetric propagation induces cell bending, a behaviour consistent with our experimental observations (figure 1, times *t*_{2} and *t*_{6}). Figure 6 illustrates the temporal evolution of principal strain fields in longitudinal and radial directions. It reveals the coupling between chemical and mechanical cell dynamics, keeping the cell under stress during all the propagation of the calcium wave. Alternatively, when the Ca^{2+} spark site was initiated in the centre of the cell, then two contraction waves propagate at equal velocity in opposite directions (figure 7*a*(i–iii)).

We additionally simulated the quasi-simultaneous occurrence of multiple calcium sparks in different regions of the cell. When two sparks originate from opposite sides of the cell, they trigger two independent Ca^{2+} waves that propagate towards the centre where they collide before disappearing (figure 7*b*(i–iii)). This behaviour agrees quite well with reported experimental data on colliding Ca^{2+} waves (Zimmermann & Walz 1997). From a mechanical point of view, the associated strain waves follow the same dynamics, with the collision of Ca^{2+} waves inducing a complete relaxation of the cardiac myocyte.

## 5. Discussion

A number of advances in experimental techniques, including the extended use of Ca^{2+} probes coupled with confocal or total internal reflection fluorescence microscopy (Blatter *et al.* 2003; Bai *et al.* 2008; Alonso *et al.* 2009), together with adapted digital image processing (Kamgoue *et al.* 2009) have promoted an accurate analysis of intracellular Ca^{2+} wave dynamics and cardiac cell contractions. Along with these developments, simulation studies have provided significant contributions to the understanding of the complexity of the excitation–contraction process that drives and regulates the mechanical behaviour of cardiac myocytes.

In this context, the present work proposes an integrated and comprehensive mechano-chemical core for modelling self-sustained contractions of cardiac myocytes that is based on an original formulation of a strain-energy function that explicits the dependence of the contractile efficiency with regard to cell elasticity, sarcomere length and cytosolic Ca^{2+} concentrations. This approach may be viewed as an alternative to the theoretical framework recently proposed by Stalhand *et al.* (2007) for modelling smooth muscle contraction. Indeed, Stalhand *et al.* (2007) modelled active contraction by introducing a virtual deformed configuration in which only active elements contract to an extension ratio *λ*_{a}, without deforming the passive structural components. The compatibility of the overall medium deformation, which is lost in this splitting process, is then recovered by further elongating elastically contractile elements up to the actual size of the deformed configuration. The contraction dynamics is additionally driven by a differential equation that prescribes the evolution of the time derivative of *λ*_{a}. The modelling approach we proposed here does not introduce such intermediate configuration and directly derives the cell contraction from the variation with time of the strain-energy function in response to an underlying spatio-temporal variation of calcium concentrations.

Using, as a biochemical core, a simple autocatalytic model in which the diffusion of Ca^{2+} released by Ca^{2+} sparks induces a regenerative Ca^{2+} release from the adjacent SR, we tested this novel strain-energy function by inducing distinct contraction patterns of the three-dimensional cardiac myocyte model. The diverse contractile patterns we simulated reproduced quite satisfactorily, both qualitatively and quantitatively, the spontaneous periodic contraction of single adult rat cardiac myocytes observed with time-lapse videomicroscopy. Notably, our simulations highlight the influence of the calcium-triggering regions on the cell contraction efficiency, with cell contraction amplitude being larger when the cell twitch is initiated in the middle part of the cell rather than at cell extremities. Such results extend previous studies analysing the influence of sparks localization on wave propagation in purely biochemical models (Spiro & Othmer 1999).

Although regulation of intracellular calcium pools by inositol 1,4,5-trisphosphate receptors and by mitochondrial Ca^{2+} buffering, as well as sarcomere-length dependence of the Hill coefficient and Ca_{50} value, was ignored, our model could successfully reproduce the monodirectional or bidirectional propagation and disappearance of strain waves, together with relevant changes in cell length. In addition, we did not consider the perturbation of intracellular Ca^{2+} dynamics induced by the cell nucleus (Alonso *et al.* 2006) and its role in the generation of spiralling Ca^{2+} waves (Dupont *et al.* 1996; Ishida *et al.* 1999). Such more-complex propagation dynamics would give a more perturbed cell contraction pattern, associated with a weaker efficiency of cell contraction.

Thus, the active strain-energy function presented here has been shown to provide a fair description of active stress regulation at the sarcomere level. It may thus serve as a basis for further theoretical analysis of modifications of intracellular calcium kinetics and cardiac cell rheological properties that may enhance the efficiency of cardiac myocyte contractility or that lead to abnormal cell contraction, with direct consequences on cardiac tissue contractility and cardiac rhythm (Stuyvers *et al.* 2000). Even if all the complexity of the electro-mechanical regulatory feedbacks (Campbell *et al.* 2008) is not considered here, the present model proposes an integrative and a rather simplified description of the cardiac myocyte contraction that may help to focus, without excessive complexity, on open questions related to the modulation of cardiac myocyte rythmic beating by the micro-environment stiffness (Engler *et al.* 2008; Tracqui *et al.* 2008) and more generally, to the non-homogeneous development of intracellular stress as a regulator of the mechano-transduction pathways involved in the control of the cardiac myocyte contraction.

## Acknowledgements

We thank Dr J. Olivares (LBFA-UJF, Grenoble, France) for providing isolated cardiomyocytes. This work is supported by a grant from IXXI Institut Rhône-Alpin des Systèmes Complexes.

## Footnotes

One contribution of 17 to a Theme Issue ‘From biological and clinical experiments to mathematical models’.

- © 2009 The Royal Society