Mechanisms of transmurally varying myocyte electromechanics in an integrated computational model.

The mechanical properties of myocardium vary across the transmural aspect of the left ventricular wall. Some of these functional heterogeneities may be related to differences in excitation-contraction coupling characteristics that have been observed in cells isolated from the epicardial, mid-myocardial and endocardial regions of the left ventricle of many species, including canine. Integrative models of coupled myocyte electromechanics are reviewed and used here to investigate sources of heterogeneous electromechanical behaviour in these cells. The simulations (i) illustrate a previously unrecognized role of the transient outward potassium current in mechanical function and (ii) suggest that there may also exist additional heterogeneities affecting crossbridge cycling rates in cells from different transmural regions.

1. Introduction

Although three-dimensional wall strains are highly inhomogeneous across the transmural aspect of the ventricular wall (Waldman ), several studies showed that fibre strain distributions are relatively uniform during contraction (Costa ; Tseng ; Mazhari ) owing to the effects of torsion and myofibre angle distributions (Arts ). However, more recent studies using higher-resolution measurements have revealed significant spatial and temporal differences in myofibre strain in the left ventricular (LV) wall throughout the cardiac cycle (Ashikaga et al. 2004, 2007). Heterogeneities of structure, function and composition have been reported in myocardial tissue and in myocytes isolated from different regions in the LV wall. These heterogeneities are likely to contribute to complex transmural patterns of myofibre strain, though the precise nature of these relationships and their relevance to global heart function are not fully understood.

(a) Transmural heterogeneity

Transmural heterogeneity of electrical activity and Ca2+ handling observed in tissue preparations (Wolk ) can be explained largely by intrinsic differences in individual myocytes isolated from the endocardial, mid-myocardial and epicardial regions of the LV wall (Liu ; Laurita ). Action potential (AP) morphology is different among the three ventricular cell subtypes, with epi- and mid-myocardial cells showing a prominent phase 1 ‘notch’ that is essentially absent in endocardial cells (Liu ). The magnitude of the notch is primarily determined by the transient outward K+ current (Ito1), which is much smaller in endocardial cells (Litovsky & Antzelevitch 1988) as a result of transmural gradients in the expression of K+ channel interacting proteins (Rosati ). AP duration (APD) has a significant transmural variation, with mid-myocardial cells exhibiting the longest APD (Liu ). A transmural gradient in L-type Ca2+ channel current (ICaL) density was observed in one study, with endocardial cells showing higher values than epi- and mid-myocardial cells (Wang & Cohen 2003). This finding remains controversial, however, as several other studies showed no difference in ICaL density among the three cell types (Li ; Bányász ; Cordeiro ). Finally, the expression of sarcoplasmic reticulum Ca2+-ATPase (SERCA) is highest in epicardial cells and most probably explains the faster decay of the Ca2+ transient and increased sarcoplasmic reticulum (SR) load observed in these cells (Laurita ; Cordeiro ).

Fewer studies have focused on transmural variations in the passive and active mechanical properties of isolated myocytes, though some important differences have been documented. Transmural variation of passive myocyte stiffness has been observed in the rat heart and is associated with changes in overall titin expression levels (Cazorla ). The ratio of titin isoform expression is also thought to modulate passive myocyte stiffness (Cazorla ), and the transmural gradients in this ratio have been observed in the canine LV (Bell ). Cazorla identified distinct length–tension relationships between epicardial and endocardial cells isolated from rat and ferret hearts. Epicardial myocytes have been shown to reach peak unloaded shortening more quickly than those from the endocardium in both guinea-pig (Wan ) and dog (Cordeiro ).

Cordeiro have conducted the only experimental study to date examining transmural heterogeneity of ionic currents, Ca2+ handling and mechanical function (electromechanics). They identified several key distinctions in the time course of unloaded shortening among canine LV epi-, mid- and endocardial cells. Endocardial cells displayed a comparatively longer delay in the onset of shortening, while epicardial cells showed the fastest shortening kinetics among the cell subtypes. Significantly, this study did not show any differences in ICaL density or excitation–contraction (EC) coupling gain among the three myocyte subtypes. Furthermore, the amplitudes of the Ca2+ transient and the peak percentage of cell shortening were not seen to be significantly different among epi-, mid- and endocardial cells. The authors concluded that some, but not all, of the differences in shortening time course were related to AP morphology and Ca2+ handling, leaving room for additional, as yet unidentified, sources of mechanical heterogeneity.

Computational models of ventricular myocyte EC coupling have been useful for investigating the mechanisms of functional differences across the LV wall (Viswanathan ; Solovyova ; Saucerman ; Flaim ). With the exception of the work by Solovyova , none of these models includes the myofilament interactions responsible for cell shortening. Several of the heterogeneities mentioned above were mechanistically integrated into a recent model of canine EC coupling (Greenstein ) by our group in the form of three different parameter sets, one for each myocyte subtype (Flaim ). This model successfully reproduced several of the observations of Cordeiro , including differences in AP morphology and Ca2+ handling among cell subtypes.

(b) Models of myocyte electromechanics

Most models of myocyte electromechanics are combinations of previously published models of EC coupling subsystems, coupled via state variables including Ca2+, membrane potential, sarcomere length, pH and ATP. Figure 1 summarizes the mechanisms of coupling and feedback now represented in various system models of ventricular myocytes. While no single published model includes all these mechanisms, this figure demonstrates the feasibility of an integrated model of this type.

Figure 1

Diagram of generalized functional coupling and feedback mechanisms in the models of myocyte electromechanics. Entities that broadly couple cell processes include voltage, Ca2+, length, ATP/ADP and pH. Single arrows represent one-way coupling; double-headed arrows indicate two-way interaction between components. Interactions not shown for clarity include effects of pH on XB attachment, TnC, NCX, RyRs and SERCA; effects of length on RyRs; and effects of ATP on SERCA and NCX. Abbreviations: LCC, L-type Ca2+ channel; RyRs, ryanodine receptors; SERCA, sarcoplasmic reticulum Ca2+-ATPase; PLB, phospholamban; NCX, Na+/Ca2+ exchanger; SACs, stretch-activated channels; TnC, troponin C; XB, crossbridge.

Modelling of ionic currents giving rise to the cardiac AP is the most mature of those subsystems involved in cellular electromechanics (reviewed in detail by Noble & Rudy 2001). Recent developments in the modelling of myocyte electrophysiology (EP) include continuous-time Markov models allowing for more mechanistic representations of ion channel behaviour individually (Irvine ) and in the context of the whole cell (Clancy & Rudy 1999; Flaim ). Detailed understanding of the events underlying Ca2+-induced Ca2+ release (CICR) is relatively recent, and thus models of the same have developed more slowly than those of myocyte EP (see review by Rice & Jafri 2001). The most advanced local control models of CICR (Greenstein & Winslow 2002) are too computationally demanding to be used in multiscale cardiac modelling. More recent models have used time-scale decomposition to capture important properties such as graded Ca2+ release while remaining computationally practicable for integration into larger-scale analyses (Hinch ; Greenstein ). Models of myofilament activation and force production have also been slow to mature owing to gaps that remain in the understanding of actin–myosin interactions and their role in thin filament activation (see review by Rice & de Tombe 2004). Recent myofilament models have addressed these deficiencies to some degree, and reproduce many experimentally observed characteristics of cardiac muscle contraction (Landesberg ; Solovyova ; Niederer ; Rice ).

Integrated models of myocyte electromechanics have been profitably used to investigate many whole-cell phenomena not within the grasp of smaller models. One of the first models to include several EC coupling mechanisms in equal detail was that of Rice , which was used to examine aspects of short-term interval–force relationships in cardiac muscle. Other subsequent models have been used to examine the relationship between force production and pacing interval (Matsuoka ; Iribe ). Iribe used an integrated model of myocyte function to describe the role of Ca2+/calmodulin-dependent kinase II in producing Ca2+ transient alternans during the decay of post-extrasystolic potentiation.

Models of myocyte electromechanics have also been extensively used to investigate the dynamic relationship that exists between the length or tension of myocytes and their electrical activity. Acute alteration of the AP by changes in myocyte length (mechanoelectric feedback) has been modelled in a number of studies (Kohl ; Kohl & Sachs 2001; Nickerson ; Solovyova ). Others have used integrated models to evaluate proposed mechanisms underlying the slow force response to stretch (Bluhm ; Niederer & Smith 2007).

Other models have added components such as mitochondrial metabolism and Ca2+ uptake (Cortassa ) and pH regulation (Crampin & Smith 2006) to fully coupled myocyte models. These modifications enabled analysis of metabolic control in the heart under normal conditions (Cortassa ) and of the mechanisms by which acidosis affects myocyte contractility during ischaemia (Crampin & Smith 2006).

Here, we describe our use of a coupled computational model of electromechanics to investigate mechanisms underlying transmural variations in myocyte function in the canine LV. We coupled recent models of canine EC coupling (Flaim ) and myofilament force production (Rice ) in order to test the hypothesis that known heterogeneities in AP morphology and Ca2+ transients are sufficient to explain experimentally observed differences in the time courses of unloaded cell shortening. The analysis suggests an important role for transmural variations in Ito1 in determining regional shortening and also indicates that variations in crossbridge function may be present.

2. Methods

(a) Model of myocyte electromechanics

We started with the endocardial canine myocyte model formulated by Flaim and coupled it to the recent model of myofilament force production and cell mechanics of Rice . Forward coupling between EP/Ca2+ handling and myofilaments was accomplished by driving the myofilament model with the variable representing [Ca2+]i. Reverse coupling was achieved by replacing the static Ca2+ buffering by TnC low-affinity regulatory sites used in the model of Flaim with the dynamic flux of Ca2+ on and off the myofilaments described by the Rice model. The modified differential equation describing [Ca2+]i as well as other details of model coupling can be found in the electronic supplementary material.

The complete electromechanics model was programmed in Matlab and run for more than 50 beats to achieve steady-state conditions for each parameter set used. A basic cycle length (BCL) of 2000 ms was used for all simulations to allow valid comparison with the experiments of Cordeiro .

AP-clamp protocols were simulated by prescribing the membrane voltage with an appropriate time-varying function rather than allowing it to be determined by membrane currents. The epicardial and endocardial AP waveforms used in AP-clamp simulations were digitized from representative traces reported by Cordeiro .

(b) Modification of published model parameters

All adjustments to model parameters from those originally published are presented in the electronic supplementary material. These adjustments were intentionally kept to a minimum, and were only made to correct the largest differences between experimental data and model behaviour. EP/Ca2+ handling parameter sets for epi-, mid- and endocardial cells set by Flaim were retained, along with those described in the original canine EC coupling model of Greenstein . The lone exception was an increase in the number of Ca2+ release units from 50 000 to 75 000 to increase the amplitude of the Ca2+ transient at a BCL of 2000 ms to more typically observed levels. This same parameter value has been used by others (Niederer & Smith 2007). A few minor changes to the original myofilament model parameters (Rice ) were made to increase the stability of unloaded shortening responses in the coupled electromechanics model. Alterations to parameters governing myofilament Ca2+ sensitivity and cooperativity were also made in order to give a time course and magnitude of shortening that agreed qualitatively with experiments (Cordeiro ). These and other adjustments to myofilament model parameters described throughout the remainder of this paper are listed in detail in the electronic supplementary material.

(c) Modelling of the exogenous Ca2+ buffer fluo-3

We represented the buffering of intracellular Ca2+ by fluo-3 dye using a standard differential equation for simple buffering (equation (S2) in the electronic supplementary material). Ca2+ transients measured experimentally reflect the time course of Ca2+ bound to the indicator molecule ([Ca2+]fluo) and are often calibrated to correspond to the concentration of cytosolic Ca2+. To mimic this measurement, the model reports the time course of [Ca2+]fluo scaled to vary between the minimum and maximum values of [Ca2+]i. We have called this quantity [Ca2+]i/app, the apparent cytosolic Ca2+ concentration. This modification enables direct comparison of model output with experimental measurements.

In all simulations where fluo-3 was present, the concentration of this indicator in the cytosol was assumed to be 15 μM, corresponding to the concentration of fluo-3 used in bathing solutions by Cordeiro . New steady-state conditions for the coupled model had to be found when fluo-3 was included, as the buffering action of this indicator was substantial enough to alter SR load and other model values. This was accomplished as before by running the model for more than 50 beats in the presence of 15 μM fluo-3.

(d) Simulation of shortening under clamped [Ca2+] and Ca2+ flux conditions

The myofilament model of Rice was used independently of the coupled electromechanics model to study Ca2+ contraction dynamics under various conditions. Two Ca2+ input modes were used in these simulations. The first allows myocyte contraction and shortening to be driven with a prescribed Ca2+ transient, which automatically sources as much Ca2+ flux as needed to maintain the command [Ca2+]i signal. This mode was used to study the shortening response to specific Ca2+ transients.

A second type of Ca2+ clamp was created whereby a fixed Ca2+ flux, Jclamp, was used to drive shortening. This approach is similar to the one originally used to estimate Ca2+ flux across the SR (Sipido & Wier 1991). A differential equation describing the time course of [Ca2+]i was added to the original myofilament model equations in which d[Ca2+]i/dt is set equal to the sum of different Ca2+ fluxes,where Jfluo, JTnC and JSERCA are [Ca2+]i and time-dependent Ca2+ fluxes describing buffering by fluo-3, buffering by TnC and uptake by SERCA, respectively. The fluxes representing the buffering of Ca2+ by fluo-3 and TnC are calculated as described above and in the original work of Rice , while uptake of SERCA is modelled as in the EP model with the assumption that SR content is fixed (see electronic supplementary material). This ‘flux-clamp’ mode allows the Ca2+ transient driving contraction to be flexibly determined under changing conditions such as the presence and absence of fluo-3.

3. Results

(a) Model validation

The fully coupled electromechanical model with endocardial myocyte parameters produced a time course of EC coupling events comparable with those reported experimentally (figure 2a,b). The model displayed unloaded cell shortening characteristics of latency to onset, time to peak and peak percentage of shortening that agreed with mean experimental values to within the reported standard deviations (Cordeiro ). The model showed qualitative agreement with other experimentally reported values. Model time from peak of unloaded shortening to 90 per cent relaxation (RT90) was 222 ms, compared with a mean value of 265 ms (Cordeiro ). Modelled Ca2+ transients (time course of Ca2+ bound to fluo-3 and scaled to peak free cytosolic Ca2+) showed a latency to the onset of 16 ms (versus 17.4 ms experimental mean), time to peak of 161 ms (versus 193 ms), RT90 of 489 ms (versus 566 ms) and peak cytosolic Ca2+ concentration of 580 nM (versus 663 nM). The APD of the endocardial model was 326 ms, compared with a mean value of 266 ms reported in experiments at the same BCL of 2000 ms (Cordeiro ).

Figure 2

Comparison of simulated and experimentally measured canine ventricular myocyte electromechanics. (a) Simulated model and (b) experimentally measured time courses of endocardial myocyte membrane potential (solid curves), [Ca2+]i/app (dashed curves, see §2 for definition) and unloaded cell shortening (dotted curves). [Ca2+]i/app and shortening traces have been scaled and are displayed in arbitrary units to facilitate comparison. (c) Comparison of model (triangles) and experimentally reported (circles) transmural trends in EC coupling. (All experimental traces and values are taken from Cordeiro .) (i) APD (ms); (ii) Sh latency (ms); (iii) Sh TTP (ms); (iv) RT90 (ms); (v) CaT TTP (ms); (vi) CaT RT90 (ms). Abbreviations: Sh, shortening; CaT, Ca2+ transient; TTP, time to peak; RT90, time from peak to 90 per cent relaxation.

(b) Transmural heterogeneity in a model of myocyte electromechanics

Parameter sets describing heterogeneous properties of membrane currents and Ca2+ handling for epicardial and mid-myocardial cells were proposed by Flaim . We performed simulations of unloaded shortening using these parameter sets in the electromechanics model to evaluate the potential for these proposed differences to explain distinctions in mechanical function reported by Cordeiro .

The transmural trend in APD achieved in the original model (Flaim ) was retained in the electromechanics model and agrees with Cordeiro and other reports (Liu ). Trends in cell mechanical function between the three cell types agreed qualitatively with those reported experimentally, including (i) longer latency to the onset of contraction in endocardial cells, (ii) shorter time to peak shortening in epicardial cells and (iii) longer RT90 in mid-myocardial cells (figure 2c). Peak percentage cell shortening matched measurements within the reported standard deviation for all three cell types (not shown), indicating a reasonable relation between peak systolic [Ca2+]i and myofilament Ca2+ sensitivity.

Conversely, the model failed to reproduce the significantly reduced RT90 observed experimentally in epicardial cells compared with the other cell types (figure 2c). The epicardial and mid-myocardial parameter sets both yielded Ca2+ transients with a rapid time to peak that was not fully consistent with measured representative traces or mean values (Cordeiro ). The mid-myocardial cell model also produced a Ca2+ transient RT90 of approximately 200 ms longer than that of the endocardial model, while this difference was not seen in experiments.

(c) AP-clamp simulations

Since the model showed a strong dependence of SR Ca2+ release on the AP waveform in epicardial and mid-myocardial cells, the influence of AP morphology on Ca2+ release and shortening was examined in the absence of other differences (such as SR load) using AP-clamp conditions. Initial conditions for all model state variables were taken from the mid-myocardial model run to steady state, after which a single beat was elicited in response to a representative endocardial or epicardial AP digitized from the literature (Cordeiro ).

L-type Ca2+ channel currents (ICaL), Ca2+ transients and cell shortening in response to epicardial and endocardial AP waveforms showed large differences in amplitude and time course (figure 3). The epicardial AP generated a peak current five times greater than that of the endocardial AP. By contrast, the endocardial AP produced a marked prolongation of the secondary phase of ICaL, though the respective currents ultimately terminated at about the same time, around 230 ms.

Figure 3

Influence of AP morphology on (a) ICaL, (b) Ca2+ transient and (c) unloaded shortening. Results shown are simulated responses to typical epi- or endocardial APs (inset to (a), digitized from Cordeiro ) in the electromechanics model. Solid and shaded curves correspond to epi- and endocardial AP-clamp, respectively. Dashed curve in (c) shows unloaded shortening by the model in response to a clamped Ca2+ transient shown in (b) as the dashed curve (see text for details).

Differences in ICaL had a clear effect on downstream EC coupling events. The peak of the Ca2+ transient under epicardial AP-clamp arrived approximately 120 ms faster and was approximately 40 per cent larger than that of the endocardial AP (figure 3b), corresponding to the greater amount of triggering Ca2+ in the initial phase of ICaL. The endocardial AP-elicited Ca2+ transient reached a peak magnitude above 50 per cent of that of the other transient, despite the fivefold smaller initial trigger current through a longer secondary phase of ICaL (figure 3a, shaded curve). The onset of unloaded shortening was plainly different between the two cases, with shortening in response to the endocardial AP starting approximately 30 ms after that of the epicardial AP (figure 3c). Peak unloaded cell shortening in response to the epicardial AP was reached approximately 100 ms faster than that observed under endocardial AP, and was approximately 30 per cent greater in magnitude.

To understand better the relationship between Ca2+ transient time course and unloaded shortening dynamics, we simulated an unloaded shortening event under a specified or ‘clamped’ Ca2+ transient corresponding to a scaled version of the epicardial AP-elicited transient (figure 3b, dashed curve). The transient was scaled such that its peak corresponded in magnitude with that of the endocardial transient in an attempt to isolate the effects of transient shape from those of absolute magnitude. Unloaded shortening during Ca2+ clamp (figure 3c, dashed curve) displayed shortening latency, time to peak and time from peak to 90 per cent relaxation all well within 10 per cent of those same values in epicardial AP-driven shortening (figure 3c, solid curve), showing that the time course of the Ca2+ transient was the primary determinant of shortening dynamics on this scale. Surprisingly, the contraction driven by the scaled epicardial Ca2+ transient resulted in a 30 per cent reduction in peak shortening compared with that from the endocardial AP-clamp, in spite of sharing the same absolute magnitude of peak cytosolic Ca2+.

(d) Influence of an exogenous buffer on Ca2+ transients and unloaded shortening

Discrepancies between model-predicted and experimentally observed Ca2+ transients for mid-myocardial and epicardial cells limited the ability of the electromechanics model to discern potential sources of transmural heterogeneity downstream of Ca2+ handling. We therefore simulated unloaded shortening driven directly by representative Ca2+ transients measured experimentally in the three cell types.

Before doing this, we added a simple representation of fluo-3 buffering and kinetics to the electromechanics model (as described in §2) to assess its effects on the Ca2+ transient and shortening. The concentration of fluo-3 was set to 15 μM, consistent with the concentration used in experiments (Cordeiro ), and the model was run to steady state. We then compared [Ca2+]i/app with [Ca2+]i from the model containing no fluo-3 (figure 4a). Addition of the fluo-3 had negligible impact on the time to peak Ca2+, but lowered the magnitude of the peak by approximately 20 per cent through buffering of free Ca2+. Time from peak to 90 per cent decay of [Ca2+]i/app was increased by 40 per cent over the same value calculated from [Ca2+]i. This prolongation of the Ca2+ transient is caused by competition for free Ca2+ between the high-affinity indicator and the usual Ca2+ uptake pathways, but also possesses an apparent component that reflects the dissociation kinetics of the dye. Overall, this simulation suggests that the apparent Ca2+ transient as assessed by indicator fluorescence is blunted at its peak and decays more slowly in comparison with the actual transient in the absence of this exogenous buffer.

Figure 4

Effect of fluo-3 on the apparent Ca2+ transient and potential errors in the prediction of shortening in response to [Ca2+]i/app. (a) [Ca2+]i (solid curve) taken from the mid-myocardial cell model with no fluo-3 present plotted with [Ca2+]i/app (dashed curve) from the same model in the presence of 15 μM fluo-3. (b) Simulated cell shortening in response to the respective Ca2+ signals shown in (a), e.g. solid curve is shortening in response to the solid curve in (a) and so on. Shaded curves in (a) and (b) correspond to flux-clamp validation as described in the text. (c) This is similar to (a) but now compares a calibrated fluo-3 fluorescence signal (dashed curve, digitized from Cordeiro ) with [Ca2+]i predicted via the flux-clamp protocol (solid curve). (d) Shortening in response to the respective Ca2+ inputs shown in (c). Tick marks on curves in (b) and (d) indicate peak shortening and 90 per cent relaxation.

To assess the type and magnitude of errors that could be introduced by driving cell shortening directly with measured indicator fluorescence, [Ca2+]i/app (figure 4a, dashed curve) was used to drive shortening in the myofilament model (figure 4b, dashed curve). Peak shortening was reduced by approximately 10 per cent, while the time to peak shortening and RT90 were increased by 32 and 15 per cent over original values, respectively. These differences indicated the need for a method to estimate the true time course of intracellular Ca2+ from fluorescence measurements.

(e) Flux-clamp validation

The electromechanics model was used to generate a test case for the predictive function of the flux-clamp input mode. An input flux, Jclamp, was determined from the electromechanics model run to steady state in the presence of 15 μM fluo-3, and then used to drive the flux-clamp myofilament model in the absence of fluo-3. The resulting Ca2+ transient and unloaded shortening (figure 4a,b, shaded curves) were compared with traces from the base model. The predicted Ca2+ transient was close in magnitude and shape to that of the original simulation, and unloaded shortening time course showed only slight differences (less than 8% change from baseline for magnitude of shortening, time to peak and RT90).

(f) Shortening responses to realistic Ca2+ transients

We used the flux-clamp protocol to estimate the true Ca2+ transient from a representative time course of indicator fluorescence from a mid-myocardial cell (Cordeiro ). Differences similar to those seen in simulations using the fully coupled model (figure 4a) can be observed between experimentally measured [Ca2+]i/app and [Ca2+]i estimated using flux-clamp (figure 4c). These included increased magnitude and faster relaxation for the estimated transient. The estimated shortening in the absence of fluo-3 generated with the flux-clamped model (figure 4d, solid curve) was compared with shortening driven directly with the fluorescence data (figure 4d, dashed curve) using the same myofilament model parameters in both cases. The contraction driven by fluorescence displayed differences in magnitude of shortening (−15%), time to peak (+7.6%) and RT90 (+28.5%) compared with flux-clamp shortening.

The parameters of the myofilament model were adjusted until time to peak shortening and shortening magnitude under mid-myocardial cell flux-clamp were within the standard deviations of measured mean values (Cordeiro ). The adjusted parameter set corresponded to a steady-state force–Ca2+ relationship possessing a half-activation value (Ca50) and Hill coefficient (nH) of approximately 1.0 μM and 3.5, respectively. Latency to the onset of shortening using these parameters was 12 ms slower than the experimentally reported mean; however, this may be attributed to noise in the original fluorescent signal or error in digitizing the trace, as even small differences in the initial rising phase of the Ca2+ transient have a strong effect on shortening latency. Additionally, RT90 was 54 ms shorter than the experimentally reported mean of 312 ms. This is most probably due to a secondary, slower time constant of relaxation seen in the representative shortening traces for mid- and endocardial cells (cf. figure 1 in Cordeiro ), which the model did not reproduce. Retaining the myofilament model parameter set used to fit reported mid-myocardial cell shortening, the flux-clamp protocol was used to determine driving fluxes matching representative epicardial and endocardial fluorescence traces, and unloaded shortening in the absence of indicator was predicted as previously (figure 5a). The myofilament model parameters adjusted to match the mid-myocardial shortening response also resulted in good reproduction of reported values for endocardial cells (figure 5b). Shortening RT90 for the mid- and endocardial cell responses showed the same systematic deviation from reported mean values and displayed the correct transmural trend (figure 5b). Peak cell shortening under endocardial flux-clamp was within the experimentally reported variation for this value (not shown).

Figure 5

Simulated unloaded shortening for epi-, mid- and endocardial cells in response to measured Ca2+ transients corrected via flux-clamp protocol. (a) Time course comparison of [Ca2+]i (shaded curves) and resulting cell shortening (solid curves) in each of the three cell types: (i) epi-, (ii) mid- and (iii) endocardial. The dotted curve shows shortening after an increase in crossbridge cycling rate consistent with the expression of the V1 myosin isoform (see text for details). (b) Comparison of simulated (i) shortening TTP (triangles) and (ii) RT90 (triangles) with mean values reported by Cordeiro et al. (2004; circles) for the three cell types. The dashed lines leading to triangles were determined from the dotted curve in (a) and once again represent the effect of increased crossbridge cycling rates.

Shortening simulated using the base myofilament model parameter set in response to epicardial flux-clamp did not show good agreement with experimental values of time to peak and RT90 (figure 5b). Time to peak shortening and RT90 were both slower than experimental means by 38 and 101 ms, respectively. These results were strongly suggestive of an underlying difference in myofilament properties between epicardial cells and the other two cell subtypes.

(g) Testing hypothesized myofilament heterogeneities

We considered and quantitatively tested several different hypotheses regarding changes to the myofilaments that might reconcile the discrepancies in Ca2+ contraction dynamics between the model and experimental measurements in epicardial cells. As many of the proposed changes affected myofilament Ca2+ sensitivity, the flux-clamp protocol had to be repeated for each one in order to fit a new driving flux. Figure 6 shows the difference between model and experimental values of shortening time to peak and RT90 for each hypothesis tested in comparison with the original parameter set (BASE in figure 6). We tested the potential effects of transmural variation in myosin isoform expression (suggested by Cordeiro ) by changing the crossbridge kinetic rates in the myofilament model to those described by Rice as corresponding to the faster V1 myosin isoenzyme (V1 in figure 6). Passive cell stiffness was increased (TITIN in figure 6), consistent with increased expression of the stiffer N2B titin isoform observed in the canine epicardium (Bell ). We hypothesized that phosphorylation of myosin-binding protein C, which is thought to increase crossbridge attachment rate (Stelzer ), could be altered in epicardial cells (MyBPC in figure 6). We also tested the sensitivity of the model to increased myofilament Ca2+ sensitivity (SENS), increased myofilament cooperativity (COOP) and a combination of increased titin-based stiffness with steeper cooperative activation (TITIN+COOP). The myofilament model parameter sets used to represent each of these hypotheses are shown in the electronic supplementary material.

Figure 6

Testing of hypothesized mechanisms of altered Ca2+ contraction coupling in epicardial myocytes. Description of the hypothesis represented by each abbreviation is presented in the text. Other abbreviations are the same as given in the legend of figure 2. Black bars, Sh TTP; grey bars, Sh RT90.

Of all the hypotheses tested in the model, none was able to account completely for the differences between the epicardial cell model and experimentally reported values (figure 6). The V1 and COOP parameter sets reduced RT90 below that shown by the BASE model, but not significantly. Overall, the relaxation of shortening in the model was not sensitive to any of the parameters changed in this process. Time to peak shortening was reduced to some extent by all but the SENS parameter set, which increased time to peak. Time to peak was most sensitive to the V1 parameter set, which brought this value to within the experimentally observed variation (figure 5a, dotted curve; figure 5b, dashed line and triangle).

4. Discussion

The present work suggests that much of the mechanical heterogeneity observed by Cordeiro in canine myocytes may be explained by previously identified and modelled differences in EC coupling and clarifies the mechanisms by which this control is likely to operate. We have presented evidence that, where heterogeneity of Ca2+ transients does not completely explain differences in shortening, a specific alteration in myofilament protein expression reconciles experimental and modelling results.

(a) Modulation of contractile dynamics by early repolarization and Ito1

Cordeiro reported significantly prolonged time between stimulus and the onset of unloaded cell shortening in endocardial cells in comparison with mid- and epicardial cells. In the coupled model, this was the result of a large difference in the slope of the rising phase of the Ca2+ transient. Early repolarization of the AP was seen as the primary modulator of this aspect of Ca2+ morphology (figure 3). As Ito1 magnitude determines distinctions in early repolarization (Litovsky & Antzelevitch 1988), we conclude that it ultimately controls contraction latency to a large extent.

The process by which Ito1 changes the time course of contraction is mechanistically explained by the electromechanics model. The slowed rate of early repolarization under endocardial AP-clamp allows less extracellular Ca2+ into the dyadic subspace, and results in only a small amount of Ca2+ release through the RyRs. However, a lower concentration of Ca2+ in the dyadic subspace causes Ca2+-dependent inactivation of L-type Ca2+ channels to proceed at a slower rate, prolonging the secondary phase of ICaL relative to that seen under epicardial AP-clamp. The sustained nature of this phase causes enough SR Ca2+ release to generate a cytosolic transient that is delayed but still of significant amplitude. A similar relationship between the rate of early repolarization and Ca2+ transient shape has been seen experimentally (Sah ).

Given the steep force–Ca2+ relationship of cardiac muscle, it seems counterintuitive that the magnitudes of shortening between the two AP-clamp cases would be closer than those of their respective Ca2+ transients. Simulations using realistic myofilament cooperativity reveal that the extent of unloaded shortening is strongly affected by the shape of the Ca2+ transient, not just the absolute magnitude (figure 3b,c). Cooperative recruitment of force slows the rate of force development at submaximal levels of myofilament activation (Campbell 1997), meaning that rapidly peaking Ca2+ transients allow less time for the system to reach maximal force than their slower counterparts. Steep cooperativity of myofilament force production may therefore function as a mechanism to preserve contractile function in endocardial cells in the absence of Ito1. Similarly, Ca2+ transient morphology may explain why endocardial cells showed slightly larger peak unloaded shortening than the other cell types in spite of no significant difference in Ca2+ transient amplitude (Cordeiro ).

(b) Transmural heterogeneity of myofilament protein expression

Flux-clamp simulations suggested that mid- and epicardial cells could not share the same Ca2+ contraction dynamics. The increased SR Ca2+ content observed in epicardial cells (Cordeiro ), probably due to increased SERCA expression (Laurita ), may be responsible for faster time to peak of the Ca2+ transient compared with mid-myocardial cells. Increased SERCA also has the effect of speeding decay and (to a lesser extent) the time to peak of the Ca2+ transient. Yet the combined effects of Ito1, increased SR content and increased SR reuptake of Ca2+, manifest in the experimentally measured epicardial Ca2+ transient, appear to be insufficient to account for the observed increase in contraction speed.

We showed using the model that increasing crossbridge kinetics in a manner consistent with increased expression of the faster V1 isoform of myosin produces a relationship between time course of free Ca2+ and time to peak shortening, which agreed well with experiments. This change alone was the most effective compared with several alternative hypotheses and suggests a transmural gradient in the expression ratio of V1 to V3 myosin, with V1 expression highest at the epicardium. Although no study of which we are aware has examined transmural gradients of myosin isoform expression in canine, others have been performed in rat (Carnes ), in rabbit (Litten ) and most recently in pig (Stelzer ). These studies universally report higher expression of V1 myosin on the epicardium than endocardium, supporting our hypothesis.

We were unable to find a parameter set for the myofilament model of Rice that would completely describe the epicardial myocyte Ca2+ contraction dynamics reported experimentally (Cordeiro ). While speeding the crossbridge cycling rate allowed the model to correctly predict time to peak shortening, relaxation still proceeded at too slow a rate. Indeed, none of the parameter modifications we attempted was able to account for the faster relaxation rate reported experimentally for epicardial cells. This is not a general failure of the model, as the relaxation–Ca2+ relationship is apparently well reproduced in mid- and endocardial cells. Instead, this points to the absence of some mechanism not included in the model whereby faster crossbridge cycling influences deactivation of the thin filament.

(c) Insights from tight coupling between EP and myofilament models

Modelling of cardiac EC coupling usually stops short of including contraction. Models and experiments have already examined the relationship between early repolarization (primarily via Ito1) and the Ca2+ transient (Sah ; Greenstein ). These studies have emphasized the adverse effects of reduced Ito1 on contraction seen in the diseased myocardium. However, it is clear that heterogeneous Ito1 expression is present in the normal heart with little effect on the strength of contraction (Cordeiro ). The electromechanics model suggests that several mechanisms working together can compensate for lower initial ICaL in endocardial cells to preserve contractile function. Coupling of the EP and myofilament models allowed us to elucidate this subtle behaviour.

The majority of models of myocyte EC coupling assume a constant affinity of TnC for Ca2+, while experiments suggest that thin filament activation and crossbridge binding alter TnC buffering capacity (Davis ). Owing to the steep cooperative regulation of crossbridge binding by Ca2+, the buffering capacity of the myofilament system may also be viewed as cooperative, undergoing large time- and Ca2+-dependent shifts within a single cardiac cycle. In the coupled electromechanics model, this dynamic buffering shortened APD relative to the original EC coupling model with static TnC buffering by reducing the buffering capacity of the myofilaments during the diastolic interval and allowing more rapid decline of [Ca2+]i via SR uptake and extrusion through other mechanisms. Dynamic buffering of Ca2+ by the myofilaments played an especially important role in the flux-clamp protocol we had designed to estimate the Ca2+ transient after removal of the fluorescent indicator. In the absence of the exogenous buffer, the myofilament buffering capacity adjusted itself in a complex nonlinear manner, resulting in changes to Ca2+ transient morphology.

Most models of myofilament force production are tested using a clamped Ca2+ transient based on the time course of fluorescent indicators. We used the electromechanics model to show that such indicators can cause large errors in Ca2+ transient shape. Thus aspects of the EP/Ca2+ handling model allowed us to make greater quantitative use of the myofilament model. In summary, including tightly coupled models of all three subsystems improved the quality of all models involved.

(d) Limitations

The models presented here have several limitations. Each of the published models that form the components of the present model contains limitations of its own (Flaim ; Greenstein ; Rice ). The most important of these is the strong cycle-length dependence of Ca2+ transient time to peak in the Greenstein model (cf. fig. 6b in Greenstein ), especially when using the parameter sets for mid- and epicardial cells presented by Flaim . This dependence caused the electromechanics model to produce qualitatively different Ca2+ transients unlike those reported experimentally when run at a matching cycle length of 2000 ms. This discrepancy did not affect the use of the model in elucidating the modulation of contraction time course by Ito1, but did limit the ability of the model to identify potential heterogeneities of myofilament properties. This shortcoming was effectively overcome via the flux-clamp simulations, which allowed the prediction of shortening directly from experimentally reported Ca2+ transients.

We did not fully explore the effects of increased SERCA and higher SR Ca2+ load reported in canine epicardial cells (Laurita ; Cordeiro ). The model of Flaim exhibits a higher SR Ca2+ load in epicardial cells relative to the other types at a cycle length of 1000 ms, which is consistent with experimental results. However, the electromechanics model at slower cycle lengths produced only a moderately elevated SR load, once again indicating an incorrect frequency dependence of Ca2+ handling. The result achieved in the model of Flaim came by doubling the activity of SERCA relative to the other cell types to reflect the experimental reports of much higher SERCA expression in epicardial cells. When this hypothesis was implemented in similar fashion in the flux-clamp protocol, it resulted in unreasonable responses after withdrawal of fluo-3. As a result, we left SERCA parameters unchanged during epicardial cell flux-clamp and allowed any additional uptake of Ca2+ by SERCA present in these cells to be implicitly incorporated into the driving flux during the fitting process.

The flux-clamp protocol involves several assumptions. The cytosolic concentration of fluo-3 is assumed to be equal to that reported for the bathing solution, and may not reflect the actual value. Furthermore, Ca2+ binding kinetics for fluo-3 were taken from experiments performed in vitro (Naraghi 1997) and may not represent kinetics within an intact myocyte. Flux-clamp was developed on the assumption that most cytosolic Ca2+ fluxes would not be strongly affected by withdrawing fluo-3 and that the Ca2+ dependence of myofilament buffering and SERCA activity were adequately described by their respective model representations. This series of assumptions seems reasonable, however, as, through the use of the flux-clamp protocol, a set of myofilament parameters was determined that was able to describe the Ca2+ contraction dynamics of both mid- and endocardial cells. Furthermore, this parameter set displayed a steady-state force–Ca2+ relationship similar to those reported experimentally in isolated canine myocytes (Wolff ).

5. Conclusion

We have used an integrative model of ventricular myocyte function that includes electrophysiology, Ca2+ handling and myofilament dynamics to explore potential mechanisms by which myocytes within the same heart modulate their EC coupling characteristics. Our work identifies Ito1 and crossbridge cycling kinetics as primary mechanisms underlying transmural variation of myocyte contractile function.