MCA (metabolic control analysis) was originally developed to deal with steady-state systems. In the present theoretical study, the control analysis is applied to the cyclic quasi-steady-state system of ion transport in cardiac myocytes. It is demonstrated that the metabolic control of particular components (channels, exchangers, pumps) of the system over such quasi-steady-state variables as action potential amplitude, action potential duration, area under the Ca2+ peak and average fluxes through particular channels during one oscillation period can be defined and calculated. It is shown that the control over particular variables in the analysed, periodical system is distributed among many (potentially all) components of the system. Nevertheless, some components seem to exert much more control than other components, and different variables are controlled to the greatest extent by different channels. Finally, it is hypothesized that the Na+ and K+ transport system exerts a significant control over the Ca2+ transport system, but not vice versa.
MCA (metabolic control analysis) [1–3] appears to be an immensely useful quantitative tool for analysing the dynamic behaviour of biochemical systems. It has been used in a great number of experimental and theoretical studies concerning the control of metabolic pathways (see [4–9] for a few examples). In its original form MCA concerns steady-state system properties. Its most fundamental idea is based on the ratio of the relative change (dY/Y) in some variable Y caused by a small relative change (dX/X) in some parameter/variable X (eqn 1):
In principle these changes are infinitesimal but in practice sufficiently small finite changes can be considered. Several coefficients have been defined within MCA. Perhaps the most important and most frequently used is the FCC (flux control coefficient) expressed as a relative change in the flux J through a given system divided by a relative change in a given enzyme concentration/activity Ei that caused the change in the flux (eqn 2):
The FCC for a given enzyme is a global parameter: its value depends on the kinetic properties of the whole system and not only on the properties of the enzyme. The values of the FCC can be different in different steady states. Usually, in linear pathways, the value of the FCC for a given enzyme is between zero (no control over flux) and one (entire control over flux). It has been demonstrated that in linear metabolic pathways the FCC values for all components of the system sum up to unity (this is the so-called summation property; eqn 3):
However, this property does not apply to branched pathways. MCA demonstrates that the control over the flux is distributed among many different steps and therefore there are no rate-limiting steps. FCCs may be attributed not only to (single) enzymes but also to metabolic blocks (groups of enzymes), carriers (e.g. the ATP/ADP carrier) and physical processes (e.g. proton leak through the inner mitochondrial membrane).
Another important coefficient is the CCC (concentration control coefficient) defined as the ratio of the relative change in a given metabolite concentration Mj to the relative change in a given enzyme concentration/activity Ei that caused the change in the metabolite concentration (eqn 4):
The values of CCC can be either positive or negative. In a given metabolic system the sum of all the CCC values equals zero (this is the so-called connectivity property; eqn 5):
Yet another coefficient defined within MCA is the elasticity coefficient expressing the sensitivity of a given enzymatic reaction rate vk to the concentration of a given metabolite Mj (eqn 6):
Generally the FCC values of particular components in some systems are inversely proportional to the values of their elasticity coefficients. Some other coefficients, for instance the response coefficient, have been defined within MCA. However, they are not directly related to the present theoretical study.
As mentioned above, MCA is traditionally applied to steady-states. It has also been shown that the MCA framework is not straightforwardly applicable to transient and oscillatory systems [13–16]. However, several theoretical approaches have been defined in order to tackle problems arising in the analysis of such systems [13–17]. It has been demonstrated that in quasi-steady-state systems with a constant oscillation period (e.g. imposed through external stimulation) it is possible to define variables whose unique relative changes can be calculated and analysed under the framework of MCA [15,16]. In such systems, variable values return periodically to the initial value, average fluxes during one cycle can be determined and such variables as the amplitude, duration and area under the oscillation (pulse) of some metabolite concentration (or thermodynamic force) can be defined.
A good example of such a cyclic quasi-steady-state system is the ion circulation and action potential generation in isolated cardiac myocytes. In this system heart cells are periodically stimulated electrically with a constant stimulation interval. This stimulation transiently increases some ion channel conductivities, generates an ion movement across the cell membrane and SR (sarcoplasmic reticulum) membrane, and thus induces changes in the membrane potential. Several kinetic models of this system have been developed [10–12] that are generally similar, although some details differ. The elements of the system taken into account explicitly in the model developed by Matsuoka et al. , used in the present study, are presented in Figure 1. The system is composed of several pumps and ion channels that conduct Na+, K+ and/or Ca2+ currents.
Scheme of the ion transport system in cardiac myocytes
The most important events occurring in the system during one oscillation period are as follows. External electrical stimulation causes a very high, but very short activation of a Na+ channel, that leads to an intensive short-term inward Na+ current, INa. This is related to a sudden increase in the membrane potential from about −85 mV (resting potential) to over +40 mV (action potential). As a result, the K1 (inward rectifier) channel closes, allowing a plateau phase in the membrane potential to generate. An increase of the membrane potential above −30 mV activates the CaL (L-type Ca2+-dependent) channel and causes an influx of Ca2+ ions from the extracellular space to the cytosol (ICaL). The increased Ca2+ concentration promotes CICR (calcium-induced calcium release) from the SR. This flux occurs through RyR (ryanodine) channels (IRyR current) with its amplitude being about two orders of magnitude higher than that of the ICaL; nevertheless RyR channels are quickly deactivated and the ratio of total Ca2+ entry through RyR and CaL channels is 15:1. The high membrane potential also causes a slow, time-dependent activation of K+ [Ks (delayed rectifier, slow component)+Kr (delayed rectifier, rapid component)] channels and an increase in K+ ion efflux from the cell (IKs+IKr). As a result the membrane potential decreases, the K1 channel is activated and the decrease in the membrane potential is further accelerated. Finally, the NaK (sarcolemmal Na+/K+-ATPase) pump and the NaCa (sarcolemmal Na+/Ca2+) exchanger bring the system back to the initial resting state. The changes in the membrane potential result from the ion (charge) movement across the cellular membrane and their magnitude is dependent on the capacitance of the membrane. Then 400 ms after a previous external stimulation the next stimulation is applied and the whole cycle is repeated. Figure 2 shows the simulated time course of the membrane potential and cytosolic Ca2+ concentration during one cycle. Of course, the calcium pulse is crucial for the stimulation of the myocyte contraction (driven by actomyosin-ATPase) and at least some components of the ATP-producing system. The above description of the behaviour of the system is a generalized one and a detailed description can be found in [10–12].
Simulated changes in transmembrane potential (
Vm) and cytosolic calcium concentration during a single oscillation period of 400 ms
In the present theoretical study the metabolic control exerted by particular components of the ion transport system in cardiac myocytes (channels, pumps, exchangers) over the amplitude and duration of the action potential, the area under the Ca2+ peak and average fluxes through particular components is analysed. It was found that, (i) the control is distributed between different steps; (ii) some steps tend to have more control than other steps; and (iii) different steps have the greatest control over different variables. Generally, it is demonstrated that MCA is very useful in describing not only steady-states, but also cyclic quasi-steady-states in biochemical/biophysical systems.
In the computer simulations performed in the present study the model of ion transport in isolated cardiac myocytes developed previously by Matsuoka et al.  was used without any modifications. Within this model the dependences of the activities of particular elements of the system (channels, pumps, exchangers) on ion concentrations, membrane potential and time are expressed in the form of kinetic equations. The changes in ion concentration and membrane potential over time are described as a set of ordinary differential equations. The authors of the model validated it successfully by a comparison of computer simulations with experimental data concerning the action potential time course, internal Ca2+ concentration transients and sarcomere length changes occurring during the contraction event. In a theoretical way they studied the effect of large-scale changes in some selected channel activities on the shape of the action potential and the magnitude of some selected currents; however, they did not perform the control analysis of the system and did not deal with the variables analysed in the present study.
The discussed model was implemented by us in the FORTRAN programming language. In order to check that no errors were made during the implementation several system properties, including those presented in Figure 2, were simulated and compared with the simulations presented in the original paper . The complete description of the model is also given in .
The control coefficients of particular steps Xi (channels, pumps, exchangers) over different variables Yj (amplitude and duration of the action potential, area under the Ca2+ peak, average fluxes through particular channels) were calculated using the following equation:
In subsequent computer simulations the original rate constants of particular steps were increased by a relative factor of 0.01 (by 1%) and the relative changes in the variable Yj between the original and the new quasi-steady-state were recorded. In a subsequent series of variable estimations, Y represented (i) action potential amplitude (the maximal difference, in mV, between the peak action potential and the resting potential: see Figure 2); (ii) action potential duration (apd50; width of the action potential peak, in ms, at 50% of the maximal difference between the action potential and the resting potential: see Figure 2) representing the actual time span of electrical excitation during each cycle; (iii) area under the Ca2+ pulse (the integral of Ca2+ concentration over time after subtracting the resting Ca2+ concentration) generally representing the potential magnitude of influence on Ca2+ sensitive subsystems (see Figure 2); and (iv) average fluxes (the net amounts of translocated charges related to some particular ion divided by the oscillation period, in pA) through selected elements of the system during the whole pulse period (400 ms), which quantitatively represent the contribution of particular elements to the overall ion-cycling activity of the system. The first three of the above mentioned variables are often easily assessed during a typical cell-behaviour recording experiment.
The ion circulation system in cardiac myocytes presented in Figure 1 is an extremely ‘branched’ system in the sense that there are many independent ways in which particular ions can be transported across the plasma membrane and SR membrane. Additionally, apart from two exchangers, namely the NaK pump and NaCa exchanger, the transport of Na+, K+ and Ca2+ ions is independent. Therefore different channels only affect each other via a common ion transported and, indirectly, via the changes in the membrane potential related to the ion redistribution across the plasma membrane. For this reason it is very important to know the values of the ionic currents through particular channels. Figure 3 presents the simulated average currents (or fluxes; in pA) during one cycle (400 ms). The ion flow through particular channels is either exclusively or predominantly unidirectional during the cycle (results not shown, see also ), and therefore the net ion movement is (almost) identical with the total ion movement. The greatest currents are conducted by the SR pump and channels: SRU (Ca2+-ATPase SR pump), SRL (SR Ca2+ leakage), RyR and SRT (transcompartmental SR Ca2+ transport). From among the plasma membrane channels, the greatest K+ currents are conducted by NaK and K1 channels, the greatest Na+ currents are conducted by NaK and bNSC (background non-selective cation) channels, whereas the greatest Ca2+ currents are conducted by NaCa and CaL channels. One can expect that these channels play an important role in determining the dynamic properties of the system.
Average currents (net amounts of translocated charges related to some particular ion divided by the duration of the oscillation period, in pA) through particular channels during one cycle period of 400 ms
Figure 4 presents the simulated control coefficients of particular elements of the system over the action potential amplitude. An increase in the activity of the NaK pump increases this amplitude, while an activation of the Na channel, the bNSC channel and the CaL channel moderately decreases the amplitude. The activities of other channels have a much smaller impact on the magnitude of the action potential amplitude. Some channels, namely Ks, Kr, CaT (T-type Ca2+), Cab (sarcolemmal Ca2+ leakage), Ito (transient outward current), Kpl [non-specific, voltage-dependent outward (plateau)], l(Ca) (Ca2+-activated background cation) and KATP (ATP-sensitive K+ channel), have apparently no control on the action potential amplitude.
Simulated control coefficients of particular system elements over the action potential amplitude
Figure 5 shows the control exerted by particular channels on the duration of the action potential (peak width at 50% of the action potential amplitude). The pattern of control is significantly different here from that in the case of the action potential amplitude. bNSC, K1 and Na channels have a strong negative control, whereas the NaK pump exerts a moderate positive control over the duration of the action potential. Moreover, the absolute values of control coefficients are generally greater here than in the previous case (the action potential duration is more sensitive to the channel activities than the action potential magnitude).
Simulated control coefficients of particular system elements over the action potential duration at 50% height
The simulated control coefficients of particular channels over the area under the Ca2+ peak are presented in Figure 6. Again, the control is distributed among several steps, from among which the NaK pump exerts a strong negative control, the bNSC exchanger and CaL channel exert a moderate positive control, while the control kept by other channels is significantly smaller. Again, several channels, namely Ks, Kr, Kpl, l(Ca), Ito, KATP, Cab and CaT seem to essentially show no control.
Simulated control coefficients of particular system elements over the area under the Ca2+ transient
It also has to be noted that these control coefficient values are much greater than in the previous two cases, meaning that the behaviour of the Ca2+ handling subsystem is much more sensitive to the changes in the activities of the ion transporting elements than the action potential generating subsystem.
Figure 7 shows the simulated FCCs of particular elements of the system quantifying the control over average ion currents exerted by the ion transporters that have the most significant control over the action potential amplitude and duration, and over the calcium pulse (Na, NaK, K1 and bNSC) and, not surprisingly, also have highest values of ion fluxes. The control over the ion fluxes through these channels is distributed mostly among just these channels (Na, NaK, K1 and bNSC). The only major exception is the relatively high control of the CaL channel over the flux through the NaK channel.
Simulated FCCs of particular channels over average currents (net amounts of translocated charges related to some particular ion divided by the oscillation period) through four selected channels (Na, NaK, K1, bNSC)
The present paper analyses the control exerted by particular components (channels, pumps, exchangers) of the cyclic quasi-steady-state ion transport system in isolated cardiac myocytes over several variables including the amplitude and duration of the action potential, the area under the Ca2+ transient and average fluxes through particular channels during one oscillation period. It is demonstrated that, first, it is possible to determine the pattern of metabolic control not only in a steady-state system, but also in a quasi-steady-state system with periodic changes of fluxes, metabolite concentration and thermodynamic forces. Secondly, the control is distributed among many (essentially all) components of the system, as in the majority of the systems analysed using MCA. Thirdly, different variables are controlled to a different extent by different components of the system. K+, Na+ and Ca2+ are transported by several routes: channels, pumps and exchangers. However, as shown in Figure 3, the average fluxes during one cycle can be very different for a given ion. In some cases the transport of two ions is coupled and has a fixed stoichiometry (NaK pump, NaCa exchanger), and in some cases an independent (without a fixed stoichiometry) flow of two or three ions through one channel occurs [CaL, bNSC, l(Ca), Ito, Ks, Na]. Of course, one can suppose that the channels with the highest average fluxes play the most important role in the system, although the time variations of the fluxes during one oscillation period are also of a great importance, especially for the shape of the action potential and Ca2+ transient.
The amplitude of the action potential is controlled to the greatest extent by the NaK pump (Figure 4), the element of the system that conducts the highest K+ and Na+ fluxes. Also bNSC and Na channels which conduct the second and third greatest Na+ currents respectively, exert a significant control over the action potential amplitude. Interestingly, the action potential amplitude is mainly modified by changing the resting potential value and not the peak value of the action potential. This is why the outward Na+ currents (hyperpolarizing currents) have a positive control over the amplitude, whereas the inward Na+ currents have a negative control over it.
The negative control over the action potential duration is mostly exerted by the bNSC, K1 and Na channels and the positive control by the NaK channel (Figure 5). All of these channels conduct significant K+ and/or Na+ fluxes. A negative control seems to be associated with the Na+ inward and K+ outward fluxes, whereas a positive control is associated with the opposite fluxes.
The area under the Ca2+ peak is significantly increased by an increased activity of the bNSC, CaL and Na channels and by a decreased activity of the NaK and NaCa exchangers. The CaL channel conducts an influx of Ca2+ ions from the extracellular space to the cytosol, while the NaCa exchanger conducts the opposite Ca2+ current, and therefore the effect of their activity on the calcium peak size is straightforward. The bNSC and Na channels transport Na+ ions to the cytosol, while the NaK pump transports these ions outside the cell. Because the outward Ca2+ transport by the NaCa exchanger is coupled with a Na+ influx to the cytosol, a high external and low internal Na+ concentration promotes the outward Ca2+ transport and therefore diminishes the Ca2+ peak. Interestingly, although the SR-related Ca2+ channels conduct the highest fluxes (see Figure 3), their control over the area under the Ca2+ peak is very moderate. This is mostly caused by the fact that the release of Ca2+ from the SR to the cytosol through RyR and SRL and the ATP-driven Ca2+ uptake by SRU are significantly controlled by the cytosolic Ca2+ itself; as discussed above, high sensitivities to metabolite concentrations (high elasticity coefficients) are associated with a low control over system variables.
The activities of four channels, namely Na, NaK, K1 and bNSC, to a significant extent control the values of the variables analysed above (the action potential amplitude, the action potential duration and the area under the calcium peak). Therefore it is interesting to test what controls the fluxes only through these channels. The simulated FCCs of all components of the system over the currents through the Na, NaK, K1 and bNSC channels, presented in Figure 7, demonstrate that the control is distributed mostly among these channels themselves. They are elements of the Na+ and K+ transport system, but not the Ca2+ transport system. This suggests that the Na+ and K+ transport system is kinetically superior in relation to the Ca2+ transport system: the former controls the latter via the membrane potential and Na+ gradient across the plasma membrane (driving the NaCa exchanger). In fact, the system composed of only these four channels is able to generate the shape of the membrane potential pulse similar to the original one (shown in Figure 2) although, of course, it cannot generate the original Ca2+ pulse (simulations not shown).
The simplistic formulation of the adenosine phosphate conversion used in the model makes the concentrations of ATP, ADP and Pi almost constant during the simulated oscillation cycle. However, in reality, some important changes in these concentrations may take place. Therefore a more sophisticated ATP production block might also contribute to the distribution of the control over the system variables.
Eight channels out of eighteen, namely Ks, Kr, CaT, Cab, Ito, Kpl, l(Ca) and KATP, seem to have essentially no control over the analysed system variables (the amplitude and duration of the action potential, the area under the Ca2+ peak and the average fluxes through particular channels). This fact seems to be strictly associated with the very low currents conducted by these channels (see Figure 3). In fact, in the simulation in which these channels were eliminated from the system (their activity was assumed to be zero), the general behaviour of the system remained essentially unchanged (simulations not shown). Therefore a question arises as to why these channels did appear during biological evolution and still exist in cardiac myocytes. It is possible that they play some unknown function and/or are activated by some unknown factors in conditions of stress, for instance during an energetic crisis and/or hypoxia. For example, the current conveyed by the KATP channel only becomes important under pronounced cytosolic ATP depletion to the level of 0.5 mM; in such conditions the channel activity increases 100-fold when compared with its basal activity at a concentration of typically 5 mM ATP. It is also interesting that some of these channels (namely Kr, CaT, Cab, KATP and Kpl) tend to control only their own activities, with values of very close to one. This is contrary to the self control coefficient values for the significant channels, for which the change in their respective activities are compensated, thus leading to only a minor change in fluxes (values equal to or smaller than 0.2).
The ion transport system analysed in the present study is in many aspects very different from a ‘typical’ metabolic biochemical pathway (leaving apart the fact that it is a cyclic, quasi-steady-state system). First, there are essentially no (or, at best, very few) linear sequential reaction chains. Secondly, the system is extremely ‘branched’, many different ways of the flow of a given ion through the plasma membrane or SR membrane are present. Thirdly, there is no chemical conversion of one metabolite into another (if ATP hydrolysis by NaK and SRU is ignored), but only a purely physical flow of different ions from one compartment into another. This system architecture has some important dynamic consequences. First, while the channels conducting a common ion can ‘communicate’ through the ion concentration and the membrane electrical potential, the channels not conducting common ions can exert an influence on each other only via this potential (in a sense, the membrane potential is a ‘common metabolite’ for all channels in a given membrane). Secondly, the summation property cannot be applied to the system, the FCCs of all components over particular (average) fluxes do not sum up to unity. Thirdly, the currents of positively charged ions in one direction, force the currents of positively charged ions in the opposite direction (this is why the system is cyclic and not unidirectional; however, one must bear in mind that the cyclic ion movement is driven by an unidirectional ATP hydrolysis by NaK and SRU).
The presented coefficients of the control over the action potential amplitude, the action potential duration and the area under the Ca2+ peak exhibit some resemblance to the CCCs: for instance, they can be either positive or negative and are related to metabolite concentrations and/or thermodynamic potentials. On the other hand, they refer to quasi-steady-state properties of the system and do not fulfill the connectivity property. Therefore they are certainly not identical with the CCCs.
One must also bear in mind that the evaluated coefficients represent only a single point in the parameter (activities) space. Thus a significant change in some activities will usually redefine the values of different control coefficients. This means, that under non-standard conditions the control over different parameters may be shifted to other components of the system. Such a situation may take place, for example, when some genes coding particular ion channels are knocked out or overexpressed, or when a permanent inhibition of particular enzymes occurs.
In the present paper we demonstrate only a few system variables selected for control analysis; in our simulations, however, several other variables were considered, for instance the area under the action potential, the Ca2+ peak amplitude or ATP usage by NaK and SRU, and many more (results not shown). In fact, any conceivable, arbitrary variable describing a particular property of the quasi-steady-state system in mind can be used as a target of the control by the dynamic elements of the entity. Of course, the selected variables should hold a potentially useful biological and/or physical meaning. The purpose of the present study was to analyse only some representative examples. In the present study, we show the plausibility and usefulness of the evaluation of control coefficients over arbitrarily defined parameters of systems oscillating with a constant period, thus enabling a quantitative analysis of interdependencies between the elements of the entity. However, our approach does not deal with changes in the frequencies of oscillations; in fact, the oscillation period is defined as a constant in the employed model.
Of course, in the heart in vivo the cardiac cells are not excited by an artificial electrical stimulation, but by neural signals coming from the SA (sinoatrial) node. The frequency of this excitation may vary depending on conditions, for instance on the presence of different hormones or neural signals coming from the brain. However, the discussed excitation is still external and hierarchically superior in relation to the ion circulation system taken into account within the model (additionally, the model does not involve any feedback signals from cardiac myocytes to the pacemaker cells in the SA node). Therefore the model cannot be used to study the control of different system elements over the frequency of the excitation and contraction events.
The knowledge of the control of particular elements of the system over different system variables gives a better insight into the functioning of the system and therefore is important for the pure science. It also potentially has very important practical applications , possibly in formulating treatment for heart diseases, for instance those related to the heart-beating cycle. The model predictions can potentially be tested by analysing, in an experimental way, the effect of different mutations leading to small changes in particular channel activities or concentrations on different system variables. They may also be validated by using selective inhibitors of particular ion transporters in experiments recording selected system variables.
In conclusion, in the present study the control analysis of chosen variables, namely the action potential amplitude, the action potential duration, the area under the Ca2+ peak and the average fluxes through particular ion channels in the cyclic ion circulation system in isolated cardiac myocytes was performed. It has been shown that: (i) the control is distributed among different components of the system (channels, exchangers, pumps); (ii) different variables are controlled by different components to a different extent; (iii) the Na+ and K+ transport system exerts a significant control over the Ca2+ transport system, but not vice versa. Overall, it has been demonstrated that an approach based on MCA can be very useful in analysing not only steady-state, but also cyclic quasi-steady-state systems.
background non-selective cation channel
sarcolemmal Ca2+ leakage
L-type Ca2+ dependent channel
T-type Ca2+ channel
concentration control coefficient
flux control coefficient
transient outward current channel
total current (flux) through system element X (in pA)
inward rectifier channel
ATP-sensitive K+ channel
non-specific, voltage-dependent outward (plateau) channel
delayed rectifier channel, rapid component
delayed rectifier channel, slow component
Ca2+ activated background cation channel
metabolic control analysis
sarcolemmal Na+/Ca2+ exchanger
ryanodine Ca2+ channel
SR Ca2+ leakage
transcompartmental SR Ca2+ transport
Ca2+-ATPase SR pump