Introduction

The purpose of the plan described in this contribution is the simulation and subsequently the review of the dynamic processes running on subcellular level in the electromechanical heart contractile system. These processes can be seen on the very simplified scheme of the electromechanical heart contractile system (see fig. 1). In the first part of our research project we have focused on the formulation of the mathematical model of the biochemic-mechanical subsystem (which is described in this paper).

Figure 1. Electromechanical heart contractile system - simplified scheme

Basic information that enables us to build a suitable model is the description of the biochemical bindings in a sarcomere [1,2,14]. This description is divided into two main blocks of the biochemical reactions. The first block of the reactions, placed in the Table 1, describes the bindings in the contractile system (i.e. among Actin (A), Myosin (M) and structures of Adenosinetriphosphate (ATP)) The second block, stated in the Table 2, describes the regulation system, thus relations between intracellular calcium Ca²⁺_i and troponin C (TnC).

TABLE 1. The actomyosin system

TABLE 2. The regulation system

By expressing of the biochemical reaction's kinetic behaviour, the system of ten ordinary non-linear 1^st order differential equations can be derived with constant coefficients - rate constants, denoted as "k", of biochemical reactions mentioned above and which represent also the parameters of the mathematical model (see Table 3).

TABLE 3. Mathematical Model

The model, presented in Matlab, allows precarious definition of a set of experiments differing from each other by values of observed parameters, initial conditions, etc. Detailed investigation of its dynamic properties has shown that the model is stable, stiff with plenty of linear loops. This fact is characteristic for a specific class of mathematical models of biological systems. The full-scale description of model's dynamic properties can be found in [3] and [4].

Dynamic changes of intracellular calcium concentraction [Ca²⁺]_i, which are resultants of the outward action potential course and relating ion currents, were in this model due to simplification substituted by the time course of an ideal rate change of a concentration of the intracellular calcium released out of the JSR (Junctional Sarcoplasmic Reticulum) d[Ca²⁺_F]/dt (see Figure 2) - the only one input signal of the model, whose mathematical characteristics (e.g. amplitude) were chosen to fit the steady state value of [Ca²⁺]_i and maximal value of [Ca²⁺]_i (before maximal contraction) during simulation according to present physiological experiments in vitro [5, 6]. Here should be denoted, that the contributions of intracellular calcium ions non-accumulated in the sarcoplasmic reticulum and therefore passing through the cell membrane directly into sarcomere and other phenomenons like mitochodria are in this model presently neglected.

Figure 2. The input signal curve representing the ideal time course of the rate of the calcium ions release out of JSR (here is, naturally, the amount of calcium ions released out of the JSR [Ca²⁺_F] identical to the amount of calcium ions entering the intracellular space [Ca²⁺_i])

In the Figure 2 are the points of simulation start, stop, model's end of the transient response and the value ensuring the right heart beat frequency (70 1/min), marked by means of indexes START, STOP, S and b. The puls width (100 ms) represents the duration of calcium release out of the JSR.

The model's basic ability to reach a significant conformance with physiological experimental measurements was qualitatively verified by comparing the simulation of the [Ca²⁺]_i course (see Figure 3) with the same one, but measured in vitro (see Figure 4) on a rabbit's heart cell [7]. Quantitative verification can only be performed constraintly with knowledge of great variability of physiological measurement results especially dependent on the method applied.

Figure 3. Transient curve of the intracellular calcium [Ca²⁺_i]- simulation result (t_S shows the point of the model's end of transient response)

Figure 4. Transient curve of intracellular calcium [Ca²⁺_i] measured in vitro The calcium steady state value is not in this figure denoted

Methods

Following dynamics depends on the quantity of biochemical parameters some of which have physiological and medical meaning and thus offer great number of degrees of freedom on selection the simulation experiments. Variability of simulation experiments with respect to mentioned biochemical parameters embodies a number of variations following from the definition of the parameter's vector (1):

P = (kan, kab, kDPA, k –DPA, k –DAP, k –AD, kAT, k –AT, k –TA, kTA, kAH, k –AH, kH, k –H, kr)

(1)

where the rate constant k_an represents the rate of association of Ca²⁺ ions on the specific asociation sites on troponin C (TnC), k_ab represents the rate of dissociation of Ca²⁺ ions from the specific binding sites on TnC, k_H represents the rate of rending of ATP (adenosintriphosphate) to ADP (adenosindiphosphate) and a phosphor molecule, k _–H represents the rate of reverse reaction, k_DPA represents the rate of forming the crossbridge in the weak binding form, k _–DPA represents the rate of reverse reaction, k _–DAP represents the rate of forming the crossbridge in the strong binding form, k _–AD represents the rate of rending of ADP molecule from the myosin's head, k_AT represents the rate of transformation of crossbridge from strong to weak form after binding the ATP molecule, k _–AT represents the rate of reverse reaction, k_TA represents the rate of binding of myosin's head with ATP molecule on the actin molecule, k _–TA represents the rate of reverse reaction, k_AH represents the rate of hydrolysis of ATP, k _–AH represents the rate of reverse reaction, k_r represents the rate of active transport of Ca²⁺ ions form sarcomere into sarcoplasmic reticulum (into NSR - Network Sarcoplasmic Reticulum [8]).

Standard physiological reality obtained by in vitro experiments reflects the element's values in vector (1) marked as P_ref [9, 10, 11, 12, 13]. Naturally values of these elements in P_ref were not measured during a single physiological experiment and thus show some variance with respect to definition of laboratory conditions (t, pH) and to animal class (heart myocyte of rabbit, guinea pig, cattle, etc.). With respect to difficulty of laboratory experiments and thus availability of corresponding measuring results we have to admit this non-uniformity in the definition of the vector P_ref.

Simulated changes in dynamics of phenomenons running in the electromechanical heart contractile system can be achieved by apriori changes of positions of the P point (defined by vector P) in 15 dimensional parametric space. A well known problem of every mathematical model simulating a biological system is the problem of ranges of physiological interpretation, it means definition of suitable set of positions of P regarding experimentally verified position of P_ref. While mathematics allows to simulate with geat ranges of parameter's values and is only restricted by model's dynamic properties (e.g. numerical unstability), borders of physiological reality are the subset of mathematical borders and can only be discovered as a result of experimental measurements in vitro, which are not in many cases available. From these reasons the variability of parameters k_r, k_an, k_ab, which influence significantly [Ca²⁺]_i and [TnC] and they are succefully physiological interpretable, was investigated. By means of this selection the number of degrees of freedom of the observed dynamics was also indispensably reduced. A great problem before own simulation is setting all initial conditions in the model because most of them are not known. We found experimentally, that after implementing the P_ref into model and after setting those well-known initial conditions, the model reaches the steady state in time no longer than 2 seconds. Therefore the input signal (defined by rectangular shape with amplitude of 480 [uM/s] and pulse width of 100 msec) starts at 2 seconds after begin of simulation (Figure 2). All simulations presented in this contribution were performed under the same conditions: begin of simulation = 0 seconds, total time = 4 seconds, sample period = 0.00001 seconds, integration method Gear, tolerance = 0.0001 seconds. The model was simulated on a two PENTIUM III processors workstation based on Windows NT and Matlab 5.

Results

Variability of k_r. The value of this parameter represents the rate of active transport [Ca²⁺]_i from a sarcomere back into the NSR. In our model we use the linear description of calcium pump (2):

(2)

The figures 5 and 6 show the maximal value of [Ca²⁺]_i and minimal value of [TnC] (both concentrations are in M) as a function of k_r, that is changing from 40 to 400 [1/s] with reference value of 200 [1/s] according to [6]. One can see, that with increasing rate of calcium pump, concentraction amplitude of free [Ca²⁺]_i is decreasing and minimal value of concentration of free (i.e. not occupied) binding sites on TnC is increasing. This result corresponds with reality because quicker outflow from sarcomere (increase of k_r) results in rapid decrease of maximal concentration of free intracellular calcium ions and on the contrary increase of concentration of free Ca-specific binding sites on TnC. The figures 7 and 8 demonstrate the dependence of concentration of free intracellular calcium ions and free Ca-specific binding sites on TnC (both in M) in the steady state on the parameter k_r. Hence follows that for reference value of k_r is steady state for [Ca²⁺]_i equals to 0,1 [uM] [5]. Although the time dependence of [TnC] is not still experimentally verified, one can read from the graph the steady state value for [TnC] (approximately 6^.78.10^-5 [M]). Picture 9 shows the dependence of [Ca²⁺]_i on the k_r and on the time. Here again has to be denoted that the efflux of the calcium ions out of the JSR is substituted by the input signal (see Figure 2) and the contributions of other calcium processes (Na-Ca exchanger, mitochondria) are not presently in the model included.

Figure 5. Dependence of the intracellular free calcium ions on the k_r parameter values

Figure 6. Dependence of the free Ca-specific binding sites on the TnC on the k_r parameter values

Figure 7. Dependence of the free intracellular calcium in the steady state on the k_r parameter values

Figure 8. Dependence of the Ca-specific binding sites on the TnC in the steady state on the k_r parameter values

Figure 9 - Dependence of the free intracellular calcium ions on the parameter k_r and on the time

Variability of k_ab. The value of this parameter reflects the rate of calcium ions dissociation from Ca-specific binding sites on TnC, which is described by the biochemical reaction (3) and the differential equation (4).

(3)

(4)

The figures 10 and 11 show the maximal value of [Ca²⁺]_i and minimal value of [TnC] (both concentrations in M) as a function of k_ab, that is changing from 20 to 200 [1/s] with reference value of 120 [1/s] according to [9]. From the figure 10 results that the change of k_ab initiates non-linear change of [Ca²⁺]_i so that increasing value of k_ab results in decreasing rate of [Ca²⁺]_i growth. Pictures 12 and 13 demonstrate the dependence of concentration of free intracellular calcium ions and free Ca-specific binding sites on TnC (both in M) in the steady state on the parameter k_ab. From the picture 12 results that real physiological value of k_ab should be higher than 60 [1/s], because in this range does not influence the steady state value of intracellular calcium concentration. Picture 14 shows the dependence of [Ca²⁺]_i on k_ab and on time.

Figure 10. Dependence of the intracellular free calcium ions on the kab parameter values

Figure 11. Dependence of the free Ca-specific binding sites on the TnC on the kab parameter values

Figure 12. Dependence of the free intracellular calcium in the steady state on the kab parameter values

Figure 13. Dependence of the Ca-specific binding sites on the TnC in the steady state on the kab parameter values

Figure 14. Dependence of the free intracellular calcium ions on the parameter kab and on the time

Variability of k_an. The value of this parameter reflects the rate of calcium ions association on Ca-specific binding sites on TnC, which is described by he biochemical reaction (5) and the differential equation (6).

(5)

(6)

The figures 15 and 16 show the maximal value of [Ca²⁺]_i and the minimal value of [TnC] (both concentrations in M) as a function of k_an, that is changing from 1^.10¹⁰ to 1^.10¹¹ [M-1s-1] with reference value of 4^.10¹⁰ [M-1s-1] according to [9]. Figure 15 shows that rate of decrease max[Ca²⁺]_i is decreasing with the increasing value of parameter k_an. Figures 17 and 18 demonstrate the dependence of concentration of free intracellular calcium ions and free Ca-specific binding sites on TnC (both in M) in the steady state on the parameter k_an. From the picture 17 results that real physiological value of k_an should be lower than 6^.10⁷ [M-1s-1], because in this range does not influence the steady state value of intracellular calcium concentration. Figure 18 gives us a relatively surprising result because with increasing value of k_an concentration of free Ca-specific binding sites on TnC is approximately linearly decreasing. Figure 19 shows the dependence of [Ca²⁺]_i on the k_an and on the time.

Figure 15. Dependence of the intracellular free calcium ions on the k_an parameter values

Figure 16. Dependence of the free Ca-specific binding sites on the TnC on the k_an parameter values

Figure 17. Dependence of the free intracellular calcium in the steady state on the k_an parameter values

Figure 18. Dependence of the Ca-specific binding sites on the TnC in the steady state on the k_an parameter values

Figure 19. Dependence of the free intracellular calcium ions on the parameter k_an and on the time

Conclusions

The results of simulations show that several of investigated dynamic processes have monoexponential character and offer ability of compartment formulation with regard to weak non-linearities. Presented changes of [Ca²⁺]_i and [TnC] also verify rightness of used method. It is clear that other parameters can be changed or their combinations can be studied and thus other phenomenons could be investigated, but presently they can not be verified in vitro. Performed simulation experiments show the dynamic behaviour of the observed part of the heart contractile system outside the area hitherto recognized by means of physiological in vitro experiments, i.e. for values of the vector P # P_ref.

In the next stage of the research project we are about to substitute the model's input signal by enhanced description of intracellular calcium ions dynamics in cardiomyocyte and then perform coupling with the action potential.

 

Acknowledgement

This contribution is the part of a research project supported by Grant's Agency of Czech republic under No. 106/99/0975.

 

References