1 Universidad Simón Bolívar, Departamento de Procesos y Sistemas, Apartado 89000, Valle de Sartenejas, Edo. Miranda 9995, Venezuela.

2 Groupe ESIEE, Laboratoire Algorithmique et Architecture des Systèmes Informatiques (A2SI), Cité Descartes, BP 99 93162, Noisy-Le-Grand, Paris, France.

3 Institut National de la Santé et de la Recherche Médicale (INSERM), U581, Créteil, F-94010 France.

A global model of the cardiac muscle cell is formulated, integrating the electrical and chemical dynamical aspects of the different components in the cytosol, the sarcoplasmic reticulum, and the cellular membrane (L-channels, Ryanodine channels, SERCA-pump, ATP-pump, Na⁺-Ca²⁺-exchanger, myofibrilles). The model is based on the idealization of the cell as a two-interconnected-stirred-tank-reactor system. It produces the expected time responses of the excitation-contraction (E-C) coupling, i.e. the oscillatory-bell-shaped curve in calcium dynamics and the characteristic plateau phase in membrane potential. Validations of the model with experimental data show that it can be adapted to represent cells of different species. Using sensitivity analysis, the phenomenological structure of the model allows the identification of the cellular sub-systems that are related to specific or altered dynamics, and hence the applications of the model in experiment design, drug-effect analysis and cardiac-pathology treatment.

Keywords: Cardiac modeling, Cardiac cell, E-C coupling, Calcium dynamics, Cellular simulation.

En este trabajo, se formula un modelo matemático integral de la célula del músculo cardíaco que toma en cuenta los fenómenos químicos y eléctricos que ocurren en las diversas estructuras que componen los diferentes elementos celulares: el retículo sarcoplásmico (canales de Rianodina, bomba SERCA), la membrana celular (canales L, bomba ATPasa, intercambiador Na⁺-Ca²⁺) y citosol (miofibrillas). El modelo está basado en la interpretación de la célula como dos tanques agitados e interconectados donde ocurren reacciones químicas ("two-interconnected-stirred-tank-reactor system"). Este modelo reproduce la duración del proceso de excitación-contracción ("E-C coupling") esperado, así como las formas características para la curva de la dinámica de calcio (campana asimétrica) y para la curva del potencial de membrana (potencial en forma de meseta). El proceso de validación con datos experimentales mostró que el modelo puede adaptarse para representar el comportamiento de células de diferentes especies de animales. Se realizó un análisis de sensibilidad que mostró que la estructura fenomenológica del modelo permite identificar el comportamiento de los diferentes elementos de la célula y relacionarlos de manera específica con dinámicas alteradas. Esto habilita al modelo para utilizarlo como herramienta en el diseño de protocolos experimentales, análisis de efectos de fármacos y formulación de tratamientos para cardiopatías.

Palabras clave: Modelado cardiaco, Célula cardiaca, Contracción-relajación, Dinámica de calcio, Simulación celular.

 Recibido: octubre de 2007 Recibido en forma final revisado: febrero de 2009

Virtually all forms of heart disease are caused by loss or an inadequate functioning of cardiac muscle cells (cardiomyocytes). Cardiomyocytes are the base of the pumping (contraction-relaxation) movement (Stern, 1992). It is widely accepted that two interdependent aspects modulate the excitation-contraction and relaxation of the cardiac muscle cell (E-C coupling process; Tseg, 1988), namely, the cyclic changes in the cell membrane potential, and the variations in calcium concentration, due to the transport of the calcium and other ions through the cardiomyocyte membranes (Bray et al. 1994; Guyton, 1997).

In order to develop appropriate strategies for pharmacological therapy in heart diseases, the mechanisms underlying the E-C coupling process must be well understood. A large number of "in vitro" and "in vivo" experiments, mostly qualitatively focused, are routinely carried out in biological research to study the role of several sub-cellular elements (Gavaghan et al. 2006). As a consequence, the development of new drugs or treatments consumes large amounts of resources and has not yet produced a complete understanding of the cellular system.

As an alternative, a mathematical model of the cell is an excellent tool to obtain better understanding of the E-C coupling process. It also allows the substitution of an important part of laboratory experiments by "in silico" experiments or simulations. A model can be used to design experimental protocols or test new ones, to reduce experimental variations, to design alternative experiments and to discard "a priori" non-desirable drug effects.

When the model is based on physical and biochemical laws, it also helps us to understand the role of the different cellular elements in alternative dynamic responses of the cell. Furthermore, the model can be adapted (by parametric identification; Söderström & Stoica, 1989) to different kinds of cells using specific sets of experimental data. Special identification techniques allow not only the embedding of experimental data into the model but also the embedding of expert knowledge.

The first models of the cardiac muscle cell were devoted to the description of the membrane potential (electrical models) due to the availability of patch clamp techniques allowing the measurement of ionic currents (Sakmann & Neher, 1995). The development of fluorescents methods helped to obtain specific information about calcium signalling (Bray et al. 1994), and was the reason why calcium homeostasis models (chemical models) were developed afterwards.

However, cardiac models have often been focused on only one aspect: the calcium dynamic or the electrical aspect, in spite of the evident relationship of the two phenomena. Therefore electrical formalism and empirical relations separately are insufficient to represent a complete E-C coupling process.

Among the electrical models, models based on the works of Luo & Rudy (1994) and Noble (2002) can be found, among others (Winslow et al. 1999; Jafri et al. 1998; Fox et al. 2001; Tusscher & Panfilov, 2006; Iribe et al. 2006). They contain a complete description of the electrophysiological aspects of the cell along with some empirical blocks, and the results have been validated using patch clamp techniques. Some of these models are packaged in the commercial tool CardioPrism® for the pharmacological industry (Zenglan, 2001). Generally, they reproduce typical membrane potential, but give a very fast calcium dynamic cycle (300-400 ms) which does not agree with the cardiomyocyte normal behavior.

Calcium dynamics models of the cardiomyocyte have been developed in two ways. One based on a complete description of the chemical activities of the cellular elements, without voltage dependencies (Tang & Othmer, 1994; Hamam et al. 2000; Rocaries et al. 2004) and the other based on space variations in the cardiomyocyte, without the use of detailed chemical descriptions (Michailova et al. 2002). In any case, these models do not provide a prediction of the membrane potential and its interdependence with the calcium homeostasis.

The model presented in this paper contains a complete description of the cardiomyocyte, both calcium dynamics and electrical phenomena. It is based on the analysis of the cardiac muscle cell as a system of two interconnected continuous-stirred tank chemical reactors, and it accounts for the different mass transport mechanisms, the electrical phenomena and the chemical kinetics of the reactions taking place.

Along with the detailed description of the calcium dynamics and electrical phenomena related to the different cellular elements, some empirical structures are introduced in the model where dynamics are not completely understood. These elements will be finally tuned when adjusting the model to expert or experimental data. This "hybrid" formulation allows a certain compromise between precision and complexity of the model, and also helps to adapt it to different conditions.

In the following sections a formulation of the physical model proposed and the corresponding mathematical model of the cardiac muscle cell are presented. Some simulations results allow the comparison of the model with previous ones. Validations with experimental data, after adapting the model to different species, are also included.

In order to formulate the mathematical equations of the model, it is assumed that:

-All chemical reactions occur in liquid phase (Tseng, 1988).

-Both tanks operate at constant volume, being the SR smaller than the cytosol. Volumes have been computed in previous works (Stern, 1992).

-The ionic concentrations are homogeneous throughout the cell, that is, there are not spatial concentration variations due to the very small size of the cardiomyocyte.

-The external membrane (sarcolemma) and the sarcoplasmic reticulum have the same transfer area.

-Extra-cellular concentrations of the different ions are constant (Stern, 1992).

-An external pulse voltage signal triggers the depolarization of the sarcolemma and the calcium homeostasis.

-Calcium concentration and membrane potential are the principal actuators of the E-C coupling process.

The mathematical description of the system is based on the formulation of typical balance equations: mass balances both in the cytosol reactor and in the SR tank and a charge balance in the reactor wall.

Fundamental dynamic equations

Mass balance of calcium ion in the main reactor (cytosol):

Mass balances on other ions (potassium and sodium) are

Mass balance of calcium ion in the SR tank:

not considered because of their faster dynamic response in relation to the calcium, as assumed above.

▪Charge balance in the main reactor wall:

where:

I_stim is the stimulus that triggers the membrane potential (square wave pulse).

The activation of the protein that constitutes a channel, namely the conversion of the chemical reactions that changes the configuration of the protein to allow its opening may depend on the concentration of the present ions and on the membrane potential. Mass balances for the different transition states in the mechanisms of these reactions lead to the computation of the fraction P of open channels.

Ryanodine channels. According to Tang & Othmer (1994), the R-channels present four possible states:

• Open (RC⁺): receptor with Ca2+ only on the positive site.

• Closed (RC⁺C^-): receptor with Ca2+ on both positive and negative sites.

• Refractory (RC^-): receptor with Ca²⁺ on the negative site only.

▪The mass balances related to each transition state are:

▪Finally, the diffusive flow through the R-channels can be computed as:

where:

x₂ is the fraction of open channels.

L-channels. Typically L-channels behavior is modeled using Hodgkin and Huxley theory (1952). Only three states of the channels are assumed, where the fraction d corres

ponds to the open channels among the totality of channels in active state. The transformation rate that describes the behavior of fraction d is a simplified function of two parameters d⁰⁰ and t_d , which are empirically adjusted using data from voltage-clamp experiments (Destexhe & Huguenard, 2000), and is obtained as:

The diffusive flow due to the activity of the L-channels is given by:

▪ Chemical flows:

These flows are a result of the activity of the proteins called SERCA pump, sarcolemma pump and Na⁺-Ca²⁺ exchanger, and are a function of the transport rate of each protein. This rate (reaction rate) is defined by the mechanisms of the chemical reactions involved in the calcium transfer. These processes are different from the transport through the R-channels and the L-channels since they occur against the concentration gradient.

SERCA pump. Using the simplification introduced by Tang & Othmer (1994); Hamam et al. (2000), that proposes a reaction mechanism of the Michaelis-Menten type (Michaelis & Menten, 1913), the reaction rate of the SERCA pump can be obtained as:

Sarcolemma pump. Is similar to the SERCA pump (Maclennan & Kranias, 2003), but its transport velocity is 90% slower (Negretti et al. 1995). Therefore the model by Tang & Othmer (1994) can be used changing the parameters as follows:

where:

q₁: 0.1 times p₁.

q₂: 0.1 times p₂.

Na⁺-Ca²⁺ exchanger: The kinetic mechanisms involved in the calcium flow induced by the Na⁺-Ca²⁺ exchanger enzyme are not well understood, being the model by Michailova et al. (2002) the better alternative:

where:

▪Chemical reactions:

Calcium ions are produced or consumed by chemical reactions in the cytosol due to the activity of the myofibrils and the R-channels.

The dynamical variations in the concentrations of myofibrils M and transition molecule CaM, must also be described in the model. The corresponding mass balances are:

Finally, the total terms of calcium generation and removal in equation (1) are given by the following expressions, where the contributions of the activity of the R-channels have already been added:

▪Ionic currents in the sarcolemma:

The descriptions of the currents generated through the external wall of the cytosol (sarcolemma), namely, the sodium, potassium, calcium and ionic-pumps currents, are taken from the work by Fox et al. (2001), where the mathematical formulations are based on Ohm’s law and Nerst-Plank’s equation. We have used this work because their model improves the description of the electrical currents for the ventricular canine myocyte from the Luo and Rudy model, by means of the application of relatively simple mathematical equations that guarantee the obtaining of numerical rapid and stable solutions, quite unlike the model developed by Jafri et al. (1998).

The equations for the mass balances and main current are finally:

The initial conditions for the solution of the model are summarized in table 2.

The equations of the model are integrated with Simulink^® of Matlab^® (v6.5), using a variable-step numerical solver based on the algorithm by Klopfenstein (1971). This method is relatively fast and accurate, and allows the modification of the relative error tolerance and absolute error tolerance values (Shampine & Reichelt, 1997). The equipment used is a Pentium III 900MHz, 256MB RAM computer. The program takes 5 minutes to reproduce 40 minutes of cell cycles.

The heuristic tuning of the model parameters allows an initial validation of a general model of the cardiac muscle cell, based on expert knowledge or typical experimental values for the dynamics of the output variables (Ca_i²+ and V). Specifically, validations are focused on:

Curb shapes: cyclic, with a plateau phase for V, and bell-shaped for Ca_i²+.

Experimental data for validation is obtained from the literature and also from laboratory experiments.

• Variable Ca_i²⁺ for human cardiomyocytes from Piacentino et al. (2003); and for dog cardiomyocytes from Winslow et al. (1999).

An important amount of experimental data on Ca_i²⁺ in chicken cardiomyocytes has been collected at the Research Laboratory Unit 99 of the INSERM (Paris), based on quantitative calcium fluorescence. Detailed descriptions of the experimental protocols can be found in Hamam et al. (2000).

The model is now validated using different sets of experimental data, corresponding to cardiac muscle cells of different species. Occasionally, modifications of selected parameters are needed to adjust the model.

The following coefficients are changed to obtain this behavior:

i. Period of the input voltage pulse is increased to: α =0.6 s.

ii. Maximum transport rate K of the Na⁺-Ca²⁺ exchanger is increased to K =1.6 μMs-1.

i. Period of the input voltage pulse is increased to: α=1.8 s.

Figure 13 shows stable oscillations for the calcium concentration in the SR, Ca_SR²⁺, with numerical values within the range from 60 μM to 500 μM as expected.

The dynamic behavior of the fraction of open Rianodine channels x₂, presented in Figure 14, is consistent with other results (Tang & Othmer, 1994; Hamam et al. 2000; Rocaries et al. 2004), and also shows values in the correct range (0-1).

A global model of the cardiac muscle cell has been presented in this paper. It provides calcium dynamics and electrical description of the cardiomyocyte, based on its interpretation as a two interconnected chemical reactors system.

Simulation results show that the model predicts the well-known cyclic behavior of the calcium concentration and the membrane potential of a cardiac muscle cell. Not only the shapes of the dynamic responses are accurate, but also the key numerical values are similar to experimental data (Table 5).

The formulation of the physical model as an integrated system of chemical reactors allows relating the calcium dynamics and electrical aspects of the cell dynamics, which has not been completely and coherently done in the past. Effectively, calcium and voltage predictions are more accurate than those obtained by previous models (Table 4).

Analyzing the summary of simulation results presented in Table 5, it can be said that a satisfactory validation of the model is achieved, except in very few situations. When prediction errors go beyond 15%, the possibility of structural errors in the model must be considered. However, since the basic phenomenological structure of the model allows relating parameters with cellular elements (Figure 2), the description of the cellular element involved in a certain type of alteration in the simulations can be improved.

Finally, a collateral result of the extensive validation procedures carried out in this paper, is an expert knowledge of the effects of different alterations in specific model components. Complete results of this sensitivity analysis, allowing

the identification of the cellular elements that determine different kind of dynamic responses, will be presented in a companion paper. It is easy to foresee at this point the applications of the model to support pharmacological research for instance, because it can be used to study the effects of different drugs by determining which alterations in cellular elements produce certain pathological behaviors of the cell.

The work presented in this paper is the result of the French-Venezuelan collaboration research project entitled: "Modeling of the cardiac phenomena: from the cell to the organ". The authors gratefully acknowledge the financial support of ECOS Nord (France), FONACIT, Universidad Simón Bolívar and FUNDAYACUCHO (Venezuela).