HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Full Text (PDF) Supplemental Information All Versions of this Article: 288/2/H613    most recent 00425.2004v1 Alert me when this article is cited Alert me if a correction is posted Citation Map Services Email this article to a friend Similar articles in this journal Similar articles in ISI Web of Science Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via ISI Web of Science (2) Citing Articles via Google Scholar Google Scholar Articles by Coutu, P. Articles by Metzger, J. M. Search for Related Content PubMed PubMed Citation Articles by Coutu, P. Articles by Metzger, J. M.

Genetic manipulation of calcium-handling proteins in cardiac myocytes. II. Mathematical modeling studies

Departments of ¹Biomedical Engineering and ²Molecular and Integrative Physiology, University of Michigan, Ann Arbor, Michigan

Submitted 10 May 2004 ; accepted in final form 5 August 2004

	ABSTRACT

TOP ABSTRACT METHODS RESULTS AND DISCUSSION APPENDIX A REFERENCES

 
We developed a mathematical model specific to rat ventricular myocytes that includes electrophysiological representation, ionic homeostasis, force production, and sarcomere movement. We used this model to interpret, analyze, and compare two genetic manipulations that have been shown to increase myocyte relaxation rates, parvalbumin (Parv) de novo expression, and sarco(endo)plasmic reticulum Ca²⁺-ATPase (SERCA2a) overexpression. The model was used to seek mechanistic insights into 1) the relative contribution of two mechanisms by which SERCA2a overexpression modifies Ca²⁺ sequestration, i.e., more pumps and an increase in the SERCA2a-to-phospholamban ratio, 2) the mechanisms behind postrest potentiation and how Parv and SERCA2a influence this response, and 3) why Parv myocytes retain their fast kinetics when endogenous SERCA2a is partially impaired by thapsigargin (a condition used to mimic diastolic dysfunction). The model was also utilized to predict whether Parv metal-binding characteristics might be modified to improve diastolic and systolic functions and whether Parv or SERCA2a might affect diastolic Ca²⁺ levels and myocyte energetics. One outcome of the model was to demonstrate a higher peak and total ATP consumption in SERCA2a myocytes and more even distribution of ATP throughout the cardiac cycle in Parv myocytes. This may have implications for failing hearts that are energetically compromised.

mathematical model; rat; parvalbumin; sarco(endo)plasmic reticulum calcium-adenosinetriphosphatase 2a

In our experimental study (10), we used gene transfer in intact isolated rat adult cardiac myocytes to study, characterize, and compare the effects of SERCA2a overexpression or Parv expression on isolated cardiac myocyte function. We found that when myocytes were subjected to single unloaded contractions, both methods resulted in similar increases in myocyte relaxation rate. However, we established that the response of SERCA2a-expressing myocytes to other conditions, such as -adrenergic stimulation [with isoproterenol (Iso)], diastolic dysfunction [chemically induced with thapsigargin (TG), a SERCA2a inhibitor], or changes in stimulation pacing [during postrest potentiation (PRP) tests], was different from that of Parv-expressing myocytes. The exact mechanisms underlying these differences are unclear. Moreover, other important issues regarding the use of SERCA2a overexpression or Parv expression have not been addressed experimentally. Indeed, issues such as the impact of these two genetic manipulations on energy consumption and Ca²⁺ diastolic level are of critical importance, because energy deficiency has been postulated to contribute to progression of heart failure (13, 39), and elevated diastolic Ca²⁺ levels have been associated with a calmodulin-dependent signaling pathway leading to maladaptive hypertrophy (7, 28).

To lend insight into the mechanisms underlying the SERCA2a overexpression or Parv expression and to address some of the issues mentioned above, we used mathematical modeling. To effectively use data presented in our experimental study (10), we needed a model that was specific to rat ventricular myocytes at physiological temperature that could reproduce sarcomere-shortening and Ca²⁺ fluorescence measurements. To our knowledge, only one model specific to rat ventricular myocytes has been published with a focus on electrophysiological parameters at room temperature, and it does not include sarcomere shortening (27). We therefore created an integrated rat ventricle cardiac myocyte model that contains a detailed mathematical representation of sarcolemmal ionic channels, Ca²⁺ handling (including effects of genetic manipulations), force production (including cooperativity in the force-pCa relation), and unloaded shortening. The model was adapted and integrated in part from component models that have been published previously (24, 29, 37).

With use of this rat integrated myocyte model, it was possible to elucidate the mechanisms underlying some of the results obtained in our experimental study (10). In particular, the model provided insights into 1) how Parv-transduced myocytes retain their already rapid kinetics when endogenous SERCA2a pumps are partially inhibited by TG, 2) whether the effects of SERCA2a overexpression are the same as removal of its natural inhibition by PLB (21), 3) how PRP works in control myocytes, and 4) how Parv and SERCA2a affect PRP. In addition, the model allowed a theoretical investigation of possible genetic modification of Parv's metal-binding site characteristics, such that its accelerating effect on cardiac myocyte relaxation is optimized without compromising the myocyte's systolic function. The model simulations also showed that SERCA2a overexpression might have the negative consequence of increasing energy use compared with control or Parv-expressing myocytes. Finally, the simulation results suggest that wild-type Parv expression might create a small increase in diastolic Ca²⁺ levels because of its slow Ca²⁺ release during the late part of the contraction cycle.

	METHODS

TOP ABSTRACT METHODS RESULTS AND DISCUSSION APPENDIX A REFERENCES

 
Overview

The model was constructed by integrating and adapting previously reported models: 1) ionic channel formulation and Ca²⁺ handling (based on Ref. 37), 2) force development (based on Ref. 29), and 3) sarcomere movement (based on Ref. 24). Each of these components is described in detail. However, first, it seems appropriate to give an overview of how the model is structured. Figure 1 provides an overview of the model. The electrophysiological component of the model is shown in Fig. 1, top. It includes the ionic currents and fluxes that control the membrane potential as well as homeostasis for Na⁺, K⁺, and Ca²⁺ [including sarcoplasmic reticulum (SR) Ca²⁺ handling]. Among the features of this first section of the model is the binding of Ca²⁺ to the low-affinity sites of troponin. The second section of the model represents active force development by the cross bridges (Fig. 1, bottom) and is represented by a four-state model: n0, n1, p0, and p1 (where n represents a nonpermissive cross bridge, p represents a permissive cross bridge, 0 is non-force-generating, and 1 is force-generating). This section of the model is linked to the electrophysiological section by the fact that the fraction of troponin sites occupied by Ca²⁺ controls the transition from the nonpermissive to the permissive state. As in the Negroni-Lascano model (24), the sarcomere is represented by an active element (cross bridges) in parallel with a passive elastic element (cytoskeletal-extracellular matrix). Because these two forces are in parallel, they form the total force (Fig. 1, bottom). Finally, the interaction between sarcomere movement and force production (Fig. 1, bottom) is based on the myofilament sliding theory (18). This section of the model is fully detailed in the original model (24) and is briefly described below.

View larger version (37K): [in this window] [in a new window]   Fig. 1. Block diagram overview of the model. Top: electrophysiological and ionic homeostasis section of the model, consisting of 5 compartments: myoplasm (MYO), subspace (SS), junctional (JSR) and network (NSR) sarcoplasmic reticulum, and extracellular space. Different ionic currents (arrows) and Ca2+ fluxes (arrowheads) are directly linked to equations in APPENDIX A. Arrow direction indicates a positive flux (J) or current (I). JL Trpn, low-affinity site of troponin. Bottom left: 4-state model for active force production by cross bridges. Cross-bridge attachment is permissive (p) or nonpermissive (n) and force-producing (1) or non-force-producing (0). Transitions between nonpermissive and permissive states are controlled by the rate k, which is a function of sarcomere length (L) and the fraction of troponin sites occupied by Ca2+ (TCa). Arrows labeled YVel indicate effect of sarcomere shortening velocity on breaking cross bridges and releasing Ca2+ from troponin. Bottom middle: overall force production model, which is the active force-generating component (bottom left) in parallel with a serial passive elastic component. Total force is the sum of forces produced by the 2 components. Bottom right: mechanisms by which sarcomere movement is handled by the model.

Ionic Channel Formulation

The ionic channel formulation is based on the model of Winslow et al. (37), which was used to represent experiments conducted in canine myocytes. To make the model rat specific, the voltage-activated K⁺ channel representation of Winslow et al. was replaced by a formulation representative of rat cardiac myocytes (16). This novel representation, as well as other minor changes on other currents, is discussed in detail in Tables 1–6      .

The major difference in Ca²⁺ handling and removal between large (e.g., canine) and smaller (e.g., rat) mammalian myocytes occurs mostly in the relative contribution of each protein, rather than in fundamental differences in their basic functionality. In particular, in the rat cardiac myocyte, Ca²⁺ removal is accomplished mostly (97%) by SERCA2a; in the canine cardiac myocyte; however, Ca²⁺ removal is mostly divided between the Na⁺/Ca²⁺ exchanger (NCX) (27%) and SERCA2a (73%) (4). To reflect the rat cardiac myocyte, we modified the canine model of Winslow et al. (37) by reducing NCX density and increasing the number of SERCA2a pumps. We elected to exclude the mitochondrial contribution to Ca²⁺ removal, inasmuch as it is believed to have little effect on a beat-to-beat basis under normal conditions (4). Changes relative to NCX and other Ca²⁺ sarcolemmal currents are discussed in Tables 1–6, whereas changes in SR Ca²⁺ handling, Ca²⁺ buffering, and genetic manipulations are described below.

SR Ca²⁺ Handling

The SERCA2a model proposed by Shannon has been extensively used in the last few years (27, 32–34, 37). However, the values of the SERCA2a parameters vary among these studies. We set the Hill coefficient to 2 and the affinity of the forward mode of the pump to 230 nM to reflect measurements obtained in rodents (5, 32). The other parameters were set to maintain 0.2 mM and 100 nM free Ca²⁺ (4) in the SR and cytosol, respectively, and to match relaxation rates as obtained in our experimental study (10).

The equations for the ryanodine receptor (RyR) remain as reported by Winslow et al. (37), except for the maximal release rate (increased to match higher buffering, see Buffering and Ca²⁺ Indicator) and the off-rate (k) from the closed state (P_c2) to the open state (P_o1; increased to reduce the refractory period). We also included a small leak between the junctional SR and the subspace compartments to reflect data obtained previously (34). This leak remained small compared with the other fluxes (15% of the SERCA2a reverse mode at rest). Finally, we decreased the transfer time constant between the subspace and the cytosol to achieve the fast rising phase of the fluorescence signal obtained in our experimental study.

Buffering and Ca²⁺ Indicator

Cytosolic buffering was set to combine two terms: troponin and other generic buffers. The sum of the two buffers was adjusted such that the buffering affinity and capacity were in the range of data available in the literature (4). Finally, the effect of the Ca²⁺ indicator fura 2-AM was added to the model (binding properties slightly modified from Ref. 3) to compare simulation results with data published in our experimental study (10).

Genetic and Chemical Manipulations

To model the addition of Parv, we used the law of mass action (9). The Parv on- and off-rates were slightly modified from a previous study (9) to obtain a better data match with our experimental study (10). To model the SERCA2a overexpression, we set a multiplying factor (K_SR) for forward and reverse modes. The forward-mode dissociation constant (K_fb) of the pumps under overexpression conditions was also decreased to reflect an increase in the SERCA2a-to-PLB ratio (5, 32). Application of TG was modeled as an allosteric inhibitor of SERCA2a (36). Affinity for TG is very high (K_D < 1 nM) (17); therefore, most of the TG in a solution will diffuse and bind to available SERCA2a sites. Under our experimental conditions (10) (1.0 ml of solution and 20,000 myocytes), 150 nM TG gave a half-maximal effect on Ca²⁺ transient dose-response measurements (unpublished data). Finally, for caffeine tests, the open probability for RyR was set to 1.0 during the drug application.

Force Production Formulation

Negroni and Lascano (24) presented a relatively concise model that includes force development and sarcomere shortening. However, the cooperative nature of the force-pCa relation and the sarcomere length dependence are largely unaccounted for in that model. To incorporate these components, we kept the shortening strategy of that model and combined it with a more complete (still relatively simple) cooperative cross-bridge model of force production (29). The cross-bridge model (model 3 in Ref. 29) is a four-state model in which the transition rate between nonpermissive and permissive cross-bridge states is controlled by the portion of troponin C occupied by Ca²⁺ (Fig. 1, bottom left). Cross-bridge transition rates are also a function of sarcomere length. The original model (29) reflects the strong cooperativity typical of the force-pCa relations observed experimentally (4). The values were adjusted to match the shortening data in our experimental study (10) and other parameters found in the literature, such as the rate of isometric force redevelopment in cardiac muscle (2, 26), twitch force in intact myocytes (38), and the force-pCa curve in intact cardiac muscle (1) (see Fig. B2, E and F, in supplemental data for this article, which may be found at http://ajpheart.physiology.org/cgi/content/full/00425.2004.DC1).

Sarcomere-Shortening Formulation

The Negroni-Lascano model (24) is a simplified version of the Huxley model (18). In the Negroni-Lascano model, the basic mechanical subunit is represented by an active elastic element (the myosin heads) in parallel with a passive elastic element (cytoskeletal proteins plus, in our case, myocyte-laminin attachment resistance). In the Huxley model, a probabilistic distribution of cross-bridge attachment is integrated along the longitudinal axis to evaluate the total force production at a given time. In the Negroni-Lascano model, the population of cross bridges responding to a probabilistic distribution is replaced by a single equivalent cross bridge.

The basic functionality of this model is depicted in Fig. 1 (bottom right). In isometric conditions, the active elastic element is at a stable position of length h_L (Eqs. A125 and A126; with h_L = L – X_L), producing a force proportional to the number of cross bridges attached and to the distance h_L. When the sarcomere shortens, some cross bridges detach and reattach at a different position with a given delay. In the model, the active element compresses to a length h, reducing the amount of force. After a certain delay, the rest of the filament slides back to its resting position, restoring the active element length to h_L. In addition, to match experimental observations, shortening (or relengthening) velocity reduces the number of attached cross bridges and releases Ca²⁺ from troponin C (Fig. 1, bottom left). The passive element force equation, originally represented by a fifth-order polynomial equation, was replaced by three second-order polynomial equations (see APPENDIX A) to simplify the computation of the inverse problem (see Algorithm and Numerical Methods).

Algorithm and Numerical Methods

The program was written in C++ and implemented on a personal computer. The differential equations were solved using a fourth-order Runge-Kutta-Merson algorithm with adaptive time steps (14). The only unusual difficulty arises from solving the force-movement relation. The algorithm is summarized in Fig. 2. During isometric simulations, the length is fixed, and the total force is the sum of the forces produced by the active and passive elements. During isotonic simulations, an external load (F_max) is applied to the myocyte. If the "test" force (passive + active) is smaller than or equal to F_max, the length remains fixed at its original value, and the test force becomes the total force. If the test force (passive + active) is larger than F_max, the sarcomere length is adjusted such that the sum of passive and active forces equals F_max. The inverse quadratic passive force-length relation is solved to find the resulting length. In unloaded shortening conditions, no external load is applied to the myocytes, and F_max = 0.

View larger version (64K): [in this window] [in a new window]   Fig. 2. Algorithm for force(F)-sarcomere length (L) computations. Top: isometric force (constant length). Bottom: isotonic force (constant maximal force).

	RESULTS AND DISCUSSION

TOP ABSTRACT METHODS RESULTS AND DISCUSSION APPENDIX A REFERENCES

To test our intact isolated rat cardiac myocyte mathematical model, we compared and contrasted simulation results obtained under different experimental conditions with previously published data. In particular, we simulated adult rat cardiac myocyte 1) action potential, membrane currents, and Ca²⁺ fluxes during electrical pacing, 2) Ca²⁺ fluxes during caffeine-induced contractions, 3) the force-Ca²⁺ relation in intact myocytes with use of an artificial Ca²⁺ clamp or TG to slow the kinetics in a quasi-steady-state mode, and 4) isotonic force contractions. The results are presented in supplemental data for this article (see above). In general, there was good agreement between the simulation results and literature data.

The model reproduced the Ca²⁺ fluorescence and sarcomere-shortening data from our experimental study (10) for control, Parv-transduced, and SERCA2a-transduced myocytes (Fig. 3). The representative experimental traces were obtained from our experimental study. They were selected (among a set of 30 traces in each group) as the most representative of the average of the different kinetic parameters reported in our experimental study. For the sarcomere length data, the representative traces were rescaled to have a baseline of 1.8 µm and an amplitude that matched the average measurements performed in that study. The parameters of the model were adjusted, each within the limits of related published data, to match the control myocyte data. These parameters remained fixed and were then used for the rest of the simulations presented in this mathematical study (unless specified otherwise). For Parv myocytes, when the on- and off-rates of Ca²⁺ and Mg²⁺ binding to Parv were fixed, only the Parv concentration was allowed to change; 0.062 mM Parv gave the best data match. For SERCA2a overexpression, only the overexpression factor (K_SR) and the forward-mode affinity constant (K_fb) were allowed to vary (for the rationale, see Modeling SERCA2a Overexpression): the best data match was obtained with K_SR = 2.0 and K_fb = 0.00022 mM.

View larger version (28K): [in this window] [in a new window]   Fig. 3. Model validation. A and B: simulation results. C and D: experimental representative traces obtained from our experimental study (10). A and C: normalized traces of Ca2+ fluorescence from Ca2+ indicator fura 2-AM. B and D: sarcomere length traces. There is good agreement between experimental and simulation results for control, parvalbumin-expressing (Parv), and sarco(endo)plasmic reticulum Ca2+-ATPase (SERCA2a)-overexpressing (SERCA2a) myocytes.

The sarcomere-shortening simulations match the experimental data in the control, Parv, and SERCA2a myocyte groups (Fig. 3, B and D). The sarcomere-shortening amplitude was matched in all groups within 0.005 µm. However, the modeling results revealed some limitations. In control myocytes, there were three main differences. The rising phase of sarcomere shortening was delayed by 8 ms in simulations compared with experimental data. Moreover, the experimental data clearly exhibited a biphasic relengthening, whereas in simulation the relengthening phase was clearly monophasic. The largest difference was in Parv myocyte simulations, where the simulation relengthening kinetic properties were slower than the experimental properties.

Aside from the differences in sarcomere shortening, we also sought to understand the differences between expression of Parv and SERCA2a in our experimental study (10) and that in the model. In our experimental study, an experimentally estimated value of 0.138 ± 0.046 mM was obtained for Parv concentration 3 days after gene transfer. This estimated value is somewhat higher than the value used in the model. Several factors could account for these differences, including experimental variations, use of an indirect method for experimentally estimating expression of these proteins, and uncertainties in the model parameters that might influence the value used for Parv concentration.

Similarly, SERCA2a overexpression was estimated to be 46 ± 7% (corresponding to K_SR = 1.46 ± 0.07), which is, in this case, lower than the value used in the model. Again, the model assumptions, such as equivalent forward and reverse rates and a cooperativity factor of 2 for the SERCA2a pumps, might be partly responsible for this discrepancy. Perhaps more importantly, the spatial distribution of the SERCA2a pumps on the junctional SR relative to the RyR Ca²⁺ release sites and whether this distribution is modified during overexpression might not be fully accounted for in the simple common pool model used here.

Mechanistic Insights

Modeling SERCA2a overexpression. When SERCA2a is overexpressed, at least two factors influence the model. First, the pump density is increased. In the model, this is achieved by multiplying the SR uptake flux equation by an overexpression factor, K_SR. Second, in most of the studies where SERCA2a has been overexpressed, no change in PLB has been reported (15, 31). Therefore, the SERCA2a-to-PLB ratio is expected to increase, and with that change, the SERCA2a forward-mode affinity for Ca²⁺ is expected to increase (resulting in a decreased dissociation constant K_fb) (5). As mentioned previously, we included both of these changes to model SERCA2a-overexpressing myocytes (Fig. 3).

In Fig. 4, we examined the individual contribution of these two factors to myocyte performance. In particular, we examined the effect of varying K_SR and K_fb on the Ca²⁺ transient and on the SR free Ca²⁺ content. Decreasing K_fb (increasing affinity) from 230 to 180 nM increased the baseline SR Ca²⁺ content from 0.20 to 0.25 mM (Fig. 4B). The net result was a large increase in the Ca²⁺ transient amplitude (Fig. 4A). When normalized (not shown), the increase in decay rate was small.

View larger version (27K): [in this window] [in a new window]   Fig. 4. Mathematical modeling of SERCA2a overexpression. A and B: effect of modifying Ca2+ sensitivity of forward mode of SERCA2a pumps (Kfb) on free Ca2+ concentration ([Ca2+]) in myoplasm (A) and sarcoplasmic reticulum (SR, B). Decreasing Kfb increased SR Ca2+ content and amplitude of the Ca2+ transient with little effect on kinetics. Kfb = 0.00018, 0.00021, 0.00023, 0.00025, and 0.00028 mM. Thick line represents value used for control myocytes with Kfb = 0.00023. C and D: effect of modifying pump density (KSR) on free Ca2+ in myoplasm (C) and SR (D). Increasing KSR from 1.0 to 3.0 had little effect on Ca2+ transient amplitude but had pronounced effects on decay kinetics. KSR = 1.0, 1.5, 2.0, 2.5, and 3.0. Thick line represents value used for control myocytes with KSR = 1.0.

Two questions arise from these simulations. First, is this uncoupling effect of adding new pumps (accelerating kinetics) vs. increasing the SERCA2a-to-PLB ratio (increasing amplitude) truly reflective of intact cardiac myocytes, or is it just an artifact of the model? Second, is this behavior dependent on experimental conditions such as temperature and animal species? Several experiments affecting the number of SERCA2a pumps and/or the SERCA2a-to-PLB ratio in different species have been published. In rabbit myocytes at 37°C, Chaudhri et al. (6) directly compared SERCA2a overexpression (increased number of pumps and SERCA2a-to-PLB ratio) with results using a PLB antisense (increase in SERCA2a-to-PLB ratio only). Under baseline conditions, they found a very similar increase in sarcomere-shortening amplitude in both groups, especially with stimulation at higher frequency. The effect on sarcomere relaxation was more pronounced in the SERCA2a-overexpressing myocytes than with the PLB antisense, although PLB antisense myocytes still relaxed more rapidly than control myocytes. One of the primary effects of -adrenergic agonists, such as Iso, is initiation of PLB phosphorylation and, therefore, an increase in the functional ratio of SERCA2a to PLB. When 10 nM Iso was applied at 37°C to control rat cardiac myocytes in our experimental study (10), the sarcomere-shortening amplitude was greatly increased (by 107.9%). However, when normalized with respect to their amplitude, the myocytes treated with Iso were only slightly faster than the controls (13.2% difference) (10). Similar results were obtained in another study (11). However, much larger increases in relaxation rates are seen with Iso at lower temperatures (10) or in larger mammals (6). It seems that the results predicted by the model, although more pronounced, qualitatively reflect those documented experimentally in rodent myocytes at physiological temperatures.

Postrest potentiation. Postrest potentiation (PRP) is postulated to reflect SR reuptake of Ca²⁺ after a rest period and the refractory period of the Ca²⁺ release channels (4). To analyze the mechanisms behind PRP, we simulated the same protocol used in our experimental study (10), namely, 10 conditioning pulses (1 s apart), a 120-s rest period, 10 prepulses (1 s apart), a variable rest period, and 10 postpulses (1 s apart). The PRP simulation results are presented in Fig. 5 (control myocytes) and Fig. 6 (Parv and SERCA2a myocytes).

View larger version (31K): [in this window] [in a new window]   Fig. 5. Postrest potentiation (PRP) in control myocytes. Traces, time-dependent signals; symbols, peak values. A: comparison between simulations (traces) and experimental normalized sarcomere-shortening data ( and ) from our experimental study (10). Gray traces and , prepulses 1–10; black traces and , postpulse 1, obtained after rest periods. B: effects of SR Ca2+ reloading with 5-s (solid trace) and 10-s (dotted trace) rest intervals. C and D: effects of modifying ryanodine receptor (RyR) refractory period on ryanodine open probability (C) and normalized sarcomere length (SL, D). Solid traces, effect of prolonged RyR refractory period; and , peak values obtained from original model (same as A); gray traces, prepulses 1–10; black traces, postpulse 1, obtained after rest periods; dotted traces, control myocytes.

View larger version (48K): [in this window] [in a new window]   Fig. 6. PRP in Parv and SERCA2a myocytes. Traces, time-dependent signals; symbols, peak values. A and D: comparison (for Parv and SERCA2a myocytes) between simulations (solid traces) and experimental normalized sarcomere-shortening data ( and ) from our experimental study (10). Gray traces and , prepulses 1–10; black traces and , postpulse 1, obtained after rest periods. Results from control myocyte simulations ( and ) are included for comparison. B: and E: comparison between control (dotted trace) and Parv (B) or SERCA2a (E) simulations (solid trace) on SR reloading with a 10-s rest interval. C: amount of total Ca2+ bound to Parv during a test with a 10-s rest interval. At the end of prepulses, a large amount of Ca2+ is released within 2 s, contributing to faster SR reloading observed in B. F: simulations of PRP test in SERCA2a-overexpressing myocytes, in which a leak between SR and myoplasm has been added. Simulation results show behavior similar to some "atypical" responses in our experimental study.

Next, we studied the effect of Parv on PRP (Fig. 6, A–C). As in Fig. 3, we used 0.062 mM Parv. Figure 6A shows the comparison between data collected in our experimental study and the simulation results. For additional comparisons, we also included the simulation results of control myocytes. The main difference between Parv myocyte simulation and experimental data was in the amplitude of prepulse 1 before normalization. Experimentally, the amplitude of prepulse 1 in Parv myocytes was almost identical to that in control myocytes. In contrast, in the simulations, prepulse 1 in Parv myocytes was only 58% of that in control myocytes. However, once normalized, Parv myocytes exhibited some of the same qualitative behavior that was observed experimentally. Indeed, Parv myocytes showed a greater amplitude deficit between prepulses 1 and 2 than control myocytes, followed by a flatter response. However, most interestingly, Parv myocytes, as seen experimentally, had the strongest early PRP (rest periods of 2–10 s) of all groups. The effect is attributable to a much faster SR Ca²⁺ reloading than in control myocytes (Fig. 6B). This phenomenon can be attributed to the Ca²⁺-buffering capacity of Parv. Indeed, during every cycle, Parv binds and releases a large amount of Ca²⁺. At a stimulation rate of 1 Hz, release of Ca²⁺ from Parv is incomplete between each cycle (Fig. 6C). However, after the final prepulse, all this extra Ca²⁺ is released back into the myoplasm to be rapidly sequestered by the SR, explaining the rapid reloading mentioned above.

We also examined the effect of PRP in SERCA2a-overexpressing myocytes (Fig. 6, D–F). Here, the simulations (Fig. 6D) are in good agreement with the experimental data, showing that SERCA2a-overexpressing myocytes were the least sensitive, among all groups, to change in stimulation rate. To explain this behavior, we plotted the SR Ca²⁺ content (Fig. 6E). The results show that the SERCA2a-overexpressing myocytes exhibited less baseline Ca²⁺ depletion than the control myocytes during the prepulses. This can be explained by the higher rate of sequestration in the SERCA2a myocytes, which allows almost all of the Ca²⁺ to be reloaded between each pulse at a stimulation frequency of 1.0 Hz.

Finally, in our experimental study (10), we showed that some of the SERCA2a-transduced myocytes exhibited a response with increasing sarcomere-shortening amplitude from prepulse 1 to prepulse 10 and were labeled "atypical." These myocytes showed an early strong PRP between rest periods of 2 and 10 s and exhibited negative PRP at longer rest periods. One possible explanation is that, in some myocytes, SERCA2a overexpression favors a leak from the SR. In the main model (Eqs. A89, A110, and A111), a small leak already exists between the junctional SR and the subspace, to reflect a leak flux measured experimentally from the RyR (34). Increasing this leak flux fourfold was not sufficient to reproduce the experimental results; however, increasing it fivefold augmented the Ca²⁺ in the subspace domain sufficiently to create spontaneous contractions (results not shown). As an alternative mechanism, we added another leak between the network SR and the myoplasm in these SERCA2a-transduced myocytes (K_SR = 2.0 and K_fb = 220 nM). The resulting PRP simulation results are presented in Fig. 6F. The simulation results followed the same trend observed experimentally in these atypical SERCA2a myocytes. Interestingly, the same SR-to-myoplasm leak in control myocytes (K_SR = 1.0 and K_fb = 230 nM) also followed the same trend, although with a maximal shortening amplitude of 7.2 nm in control myocytes vs. 37.1 nm in SERCA2a myocytes, which would be difficult to detect experimentally. Whether such a leak would be caused by SERCA2a overexpression or whether SERCA2a overexpression would "rescue" a portion of the population naturally bearing that kind of leak remains to be investigated experimentally.

TG as a chemical tool to mimic cellular diastolic dysfunction. In our experimental study (10), we used TG to mimic impairment in Ca²⁺ removal as seen in heart failure (23). Indeed, TG is a specific high-affinity blocker of SERCA2a (17). We showed that Parv-transduced myocytes, in contrast to control myocytes, were able to maintain their diastolic function in the presence of TG. We also mentioned that use of TG alone to model pure diastolic dysfunction was not ideal, because TG also impaired contractile function. However, use of elevated amounts of external Ca²⁺ helped correct this situation. Here, we will investigate some of these aspects, first in control myocytes (Fig. 7) and then in Parv-transduced myocytes (Fig. 8).

View larger version (31K): [in this window] [in a new window]   Fig. 7. Effects of thapsigargin (TG) in control myocytes. A and B: dose response of TG for normalized Ca2+ fluorescence and normalized sarcomere-shortening traces. Traces, simulation results with increasing TG concentrations; symbols, experimental averages obtained from our experimental study (10) with () or without () 150 nM TG. C: TG dose response of sarcomere length amplitude and relaxation rate index [calculated as 1/(Tp r), where Tp r is time from peak to three-fourths relengthening]. Sarcomere-shortening amplitude was more affected than relaxation kinetics at any TG concentration. D: potential strategies to compensate for loss of contractile amplitude to create a pure diastolic dysfunction model. Three strategies were used in the presence of 150 nM TG: increasing external Ca2+ concentration from 1.8 to 1.99 mM, reducing Na+-Ca2+ function to 72% of its original value, and prolonging the action potential using 0.66 mM 4-aminopyridine (4-AP; 3 quasi-superimposed solid thin traces). All 3 methods produced the similar result of maintaining or restoring systolic function (compared with no special treatment with 150 nM TG) and producing diastolic dysfunction (compared with simulation in the absence of TG).

View larger version (30K): [in this window] [in a new window]   Fig. 8. Effects of TG in Parv myocytes. A and B: dose response of TG for normalized Ca2+ fluorescence and normalized sarcomere-shortening traces. Traces, simulation results with increasing TG concentrations; symbols, averages obtained from our experimental study (10) with or without 150 nM TG. C and D: Ca2+ flux from SERCA2a pumps (C) and by sarcolemmal Ca2+-ATPase and Na+/Ca2+ exchanger combined (D) for control myocytes with and without TG and Parv myocytes with TG. TG affects mostly peak Ca2+ uptake by the SR, whereas sarcolemmal removal processes only increased their transient uptake duration. Note difference in scale between C and D, showing that SR Ca2+ uptake in the rat is the dominant process, even in the presence of 150 nM TG.

P r

Next, we studied the effect of TG on Parv-transduced myocytes. The model was not able to match the fast kinetics observed experimentally in Parv-transduced sarcomere relengthening in the absence of TG (Fig. 3). However, as seen experimentally, Parv-transduced myocytes showed no diastolic dysfunction in sarcomere-shortening data (Fig. 8B) when up to 350 nM TG was applied. Again, as observed experimentally, Ca²⁺ fluorescence in Parv myocytes exhibited diastolic dysfunction in the presence of TG only in the late part of the cycle (Fig. 8A).

From these results, two questions arise: 1) Where does Ca²⁺ go when TG is present? 2) How is Parv able to maintain myocyte relaxation function in the presence of TG? To answer these questions, we tracked the Ca²⁺ sequestration fluxes into the SR (SERCA2a pumps; Fig. 8C) and through the sarcolemma (NCX + sarcolemmal Ca²⁺-ATPase; Fig. 8D). TG at 150 nM reduced the peak SR uptake flux to about one-half of its value in the absence of TG in control myocytes (Fig. 8C). The total Ca²⁺ uptake was also reduced. The differential Ca²⁺ fluxes went through the sarcolemma, not by way of an increase in Ca²⁺ flux peak value but, rather, an increase in the duration of the transient removal (Fig. 8D). When the different Ca²⁺ fluxes were integrated over a 1-s period after stimulation, distribution in the absence of TG was 94.1% for SERCA2a, 4.4% for NCX, and 1.5% for sarcolemmal Ca²⁺-ATPase (see Fig. B2 in supplemental data). In the presence of 150 nM TG, the distribution was 90.6% for SERCA2a, 7.3% for NCX, and 2.1% for sarcolemmal Ca²⁺-ATPase. This change in the relative contribution of each removal mechanism is often seen in heart failure (4). Interestingly, in transduced myocytes, some of the sequestration demand on SERCA2a was eased by Parv during most of the Ca²⁺ transient decay and redistributed later when the demand was less. The total distribution in Parv myocytes in the presence of TG was 91.6% for SERCA2a, 6.2% for NCX, and 2.3% for sarcolemmal Ca²⁺-ATPase.

Model Predictions

Optimizing Parv-binding characteristics. The law of mass action was used to model Parv interaction with Ca²⁺ and Mg²⁺. This results in two first-order differential equations (Eqs. A95 and A96) that are controlled by the four binding rate constants to Ca²⁺ and Mg²⁺: k, k, k, and k. In turn, the concentration of Parv and the fraction of its sites occupied by Ca²⁺ determine its effect on the Ca²⁺ decay rate (Eqs. A93 and A109). In an experimental context, where a given Parv molecule has a predetermined set of binding rate constants, the only controllable variable is the concentration. As predicted with the use of a much simpler model and as observed experimentally in a previous study (9), increasing Parv concentration resulted in an acceleration of the Ca²⁺ transient decay and a decrease in its amplitude (Fig. 9A). However, with the set of parameters used in the present study, the resulting Parv concentration range where attenuation is negligible compared with the increased relaxing effects on the kinetics (Fig. 9D) is rather small compared with that observed experimentally (8). Nonetheless, one can ask the following question: Is it possible, at least theoretically, to modify the binding characteristics of Parv to broaden the concentration range where relaxation rates are increased with little or no effect on sarcomere-shortening and Ca²⁺ transient amplitude?

View larger version (32K): [in this window] [in a new window]   Fig. 9. Mathematical modeling of Parv expression. A: Ca2+ transient traces at different Parv concentrations showing dual effects of amplitude attenuation and accelerated decay. B: amount of Ca2+ binding to Parv during a contraction. Thin solid line, original wild-type model with 0.062 mM Parv (as used in Fig. 3); other lines, modified Parv (off-rates divided by 2.5 at 0.062 and 0.200 mM Parv). C: comparison between original wild-type and modified Parv. Results show that, for the same Ca2+ transient amplitude, modified (or "optimized") Parv has faster decay kinetics than wild-type model. D: analysis of optimal Parv concentration range (i.e., Parv concentration at which amplitude is >90% that of control myocytes), while time from stimulus to one-half decay (Ts d) shows reduction to <90% that of control myocytes. In original wild-type model, this range was very narrow and occurred at relatively low Parv concentrations. Modified (or "optimized") Parv showed a wider optimal range but occurred at higher Parv concentration.

Energetics and diastolic Ca²⁺ levels. Two important issues that have not been addressed experimentally are investigated using the mathematical model. First, in heart failure, energy production is often impaired (13, 39), raising the following important question: What happens to energy utilization when Parv is expressed de novo or when SERCA2a is overexpressed? Second, the diastolic levels of free Ca²⁺ in the myoplasm influence the level of Ca²⁺ bound to calmodulin. In turn, Ca²⁺-calmodulin can initiate a cascade of signaling events that can lead to maladaptive hypertrophy (7, 28). The question is: Does Parv or SERCA2a affect the diastolic Ca²⁺ levels?

The main ATP-dependent proteins involved in the cardiac contraction cycle are the myosin heavy chain, the SERCA2a, the Na⁺-K⁺ pump, and, to a much lesser extent, the sarcolemmal Ca²⁺-ATPase. Unfortunately, because of the formulation of the model, it is not possible to assess ATP consumption by the myosin heavy chain, but it is possible to do so for the three other processes. Indeed, the SERCA2a transports two Ca²⁺ per ATP; the sarcolemmal Ca²⁺-ATPase requires one ATP per Ca²⁺ carried across the membrane, and the Na⁺-K⁺ pump allows three Na⁺ out and two K⁺ in, for each molecule of ATP used (4). Using this information, we traced the combined instantaneous ATP consumption of these three processes, as a function of time during a 0.2-Hz contraction (Fig. 10A). The peak level of ATP utilization was greatest in SERCA2a myocytes, whereas levels for Parv and control myocytes were similar. As explained previously (Fig. 8C), when expressed, Parv dominates Ca²⁺ sequestration during the early part of the Ca²⁺ transient decay, reducing the energy demand during that period and redistributing it later in the cycle, where the demand is less. Baseline ATP use was mostly controlled by the Na⁺-K⁺ pump and was similar in all groups. For that reason, we focused on the transient consumption of ATP by first removing baseline ATP use and then integrating over time to obtain the total amount of ATP used by the three processes (Fig. 10B). Most of the transient ATP use occurs during the first 500 ms (Fig. 10A). When the integration was performed during that period (0–500 ms), SERCA2a myocytes showed the greatest total consumption among the three groups, whereas Parv myocytes used slightly less total ATP than control myocytes. When the integration was performed over the entire cycle (0–5,000 ms), the results were the same, except the Parv level of ATP use reached that of control myocytes. This confirms that Parv myocytes reduce ATP utilization during the transient and redistribute this utilization in the late part of the cycle at a much lower rate.

View larger version (21K): [in this window] [in a new window]   Fig. 10. Energetics and diastolic Ca2+ levels. A: traces of combined instantaneous utilization of ATP by 3 of the major ATP consumers: SERCA2a, Na+-K+ pump, and sarcolemmal Ca2+-ATPase. Results show highest peak consumption of ATP by SERCA2a-overexpressing myocytes. B: total transient ATP utilization during the first 500 ms and during the entire cycle (5,000 ms). Baseline ATP consumption was removed before integration, because it was the same for all myocyte groups. SERCA2a myocytes showed highest total transient ATP consumption. C: effect of Parv and SERCA2a on free Ca2+ diastolic levels at different stimulation frequencies. Diastolic Ca2+ values were measured after 24 beats just before onset of the last stimulation (or at 120,000 ms for 0 Hz).

Model Limitations

One approach to study the performance of the myocardium is to first isolate individual myocytes and then perform electrophysiological, Ca²⁺ transient fluorescence, and/or unloaded sarcomere-shortening measurements. Isolated myocytes can be used to investigate pathological conditions, pharmacological treatments, biological diversity (myocyte lineage, species, gender, and age), environmental conditions, and genetic modifications. In this study, we developed, for the first time to our knowledge, a mathematical model that simultaneously addresses these three modalities (electrophysiology, Ca²⁺ handling, and unloaded sarcomere shortening) in genetically engineered rat ventricular myocytes at physiological temperature.

Although the model was able to replicate several experimental results, such as PRP, inhibition of SERCA2a by TG, and Ca²⁺-handling protein genetic manipulations, it has limitations. In particular, the force-generation/sarcomere-shortening component does not always reveal all the details seen experimentally. For example, the experimental sarcomere-relengthening traces clearly exhibit a first rapid relaxation followed by a somewhat slower return to baseline. Alternatively, the simulation results simply follow a monoexponential time course of relaxation (Fig. 3). Moreover, for equivalent increases in decay rate of Ca²⁺ fluorescence data, the model was not able to reproduce the sarcomere-shortening fast relaxation rates observed in Parv myocytes experimentally. These two points are most likely related and could potentially be attributed to several factors. First, the myofilament compartment of the model is an oversimplification of a more complex system, where millions of independent cross bridges attach and detach, producing force and sarcomere length movement. Moreover, this system is influenced by sarcomere length, availability of ATP, pH, and inorganic phosphate. Second, the lack of force-related data in intact isolated myocytes made the development of that section of the model difficult. Indeed, several parameters had to be based on data from permeabilized myocytes at low temperature or from entire cardiac muscle bundles. Recently, a new technique to measure force in intact cardiac myocytes has been developed using special carbon-fiber force transducers (38). With this technical advancement, it is hoped that tests such as force-pCa measurements [using TG (1); see Fig. B2 in supplemental data], slack tests, isometric and isotonic twitches, and force-frequency relation could provide more insights specific to intact single cardiac myocytes and help enhance the model.

Implications

Along with revealing underlying mechanisms behind the experimental results presented in our experimental study (10), the model simulations led to new insights related to the use of Ca²⁺-handling proteins to correct diastolic dysfunction. From the experimental and mathematical modeling studies presented here, it is apparent that there are different secondary consequences of using Parv or modifying SERCA2a function to increase relaxation rates. Interestingly, the way in which SERCA2a function is modified might also influence the functional outcome of myocyte contractile performance (see Modeling SERCA2a overexpression and Fig. 4). Indeed, Fig. 4 shows that modifying Ca²⁺ sensitivity of SERCA2a (as seen when SERCA2a-to-PLB ratio changes) mostly influences Ca²⁺ transient amplitude in rat myocytes, whereas increasing the number of pumps (without changing the SERCA2a-to-PLB ratio) mostly influences the Ca²⁺ sequestration rate. This result implies that using SERCA2a overexpression, which affects the SERCA2a-to-PLB ratio and the maximal number of pumps, might not be equivalent to only modifying the SERCA2a-to-PLB ratio, as with PLB-knockout and antisense strategies. This may have implications for potential SERCA2a-to-PLB strategies in human heart failure.

-Adrenergic stimulation has three main effects on the heart in vivo: 1) increase in contractility (positive inotropy), 2) increase in relaxation (positive lusitropy), and 3) increase in heart rate (positive chronotropy) (4). In our experimental study (10), we showed that, in SERCA2a-overexpressing myocytes, the -adrenergic response was blunted on stimulation at 0.2 Hz. However, the PRP simulation and experimental results obtained in this study, as well as frequency response obtained in other studies (6, 31), show that SERCA2a overexpression has a positive inotropic effect at increased stimulation pacing compared with control myocytes (see prepulses in Fig. 6D at 1.0 Hz). It is therefore conceivable that SERCA2a overexpression in vivo could appear to have normal behavior in the presence of a -adrenergic agonist, because the positive inotropic response caused by an elevation in heart rate in vivo could possibly compensate for the blunted -adrenergic response.

Results from the PRP simulations also revealed an interesting interplay between Ca²⁺ unloading from Parv and Ca²⁺ reloading in the SR. Indeed, during the first 1 or 2 s after the prepulses, SR Ca²⁺ reloading was much more rapid in Parv-transduced than in control myocytes. Is it possible that, in Parv myocytes, early Ca²⁺ sequestration renders the competition for final removal between NCX and SERCA2a biased toward the latter, allowing less Ca²⁺ to leave the myocyte and more to be taken up by the SR? It seems that this could be the case, because the basal SR Ca²⁺ load is slightly greater in Parv than in control myocytes at the end of prepulse 10 of PRP (Fig. 6B) and during steady-state pacing at 0.2 Hz (basal SR Ca²⁺ = 207 mM for Parv and 200 mM for control). Additionally, TG simulations (Fig. 8, C and D) show that Parv myocytes have an amplitude of SR Ca²⁺ flux similar to that in control myocytes, whereas sarcolemmal Ca²⁺ uptake was lower in myocytes expressing Parv. This may be important, because in heart failure the opposite situation is favored; i.e., the tendency is for Ca²⁺ to leave the myocyte because of a higher NCX-to-SERCA2a sequestration ratio.

In addition, we found through the model that energy consumption was more elevated in SERCA2a-transduced than in control or Parv-transduced myocytes (Fig. 10, A and B). It is estimated that, under normal conditions, the rate of ATP synthesis in cardiac myocytes is 1–2 mM/s (22, 30). If it is assumed that pacing frequency is 1 Hz, the difference observed in our model represents 1% of the total ATP synthesis in normal conditions. However, this difference might become more important in heart failure, where ATP synthesis and transfer efficiency are greatly impaired (35).

In summary, this modeling study and our experimental study (10) provide a direct comparison between two genetic strategies that have been proposed to improve relaxation rate in myocardial tissue: Parv expression and SERCA2a overexpression. Moreover, the model allowed accessibility to variables such as distribution of Ca²⁺ fluxes (SR vs. sarcolemma), instantaneous SR Ca²⁺ load under all conditions, and assessment of energy utilization. Overall, the results showed that these two methods (Parv vs. SERCA2a), or even modification of the SERCA2a-to-PLB ratio using PLB antisense, are not fully equivalent. The best method may be dependent on the cause, the molecular symptoms, and the state of the cardiac disease. Importantly, in vivo experimentation in larger mammalian systems, under normal and diverse pathological conditions, would be required before these results could be extrapolated to possible human therapies.

	APPENDIX A

TOP ABSTRACT METHODS RESULTS AND DISCUSSION APPENDIX A REFERENCES

 
Equations. Na⁺ current (I_Na)

(A1)

(A2)

(A3)

(A4)

(A5)

(A6)

(A7)

For V –40 mV

(A8)

(A9)

(A10)

(A11)

For V < –40 mV

(A12)

(A13)

(A14)

(A15)

Inward rectifier K⁺ current (I_K1)

(A16)

(A17)

(A18)

(A19)

(A20)

Transient outward current (I_to)

(A21)

(A22)

(A23)

(A24)