Mathematical models of canine right and left atria cardiomyocytes

Mathematical models of canine right and left atria cardiomyocytes 1 first-author †E-mail:xialing@zju.edu.cn †E-mail:xialing@zju.edu.cn xialing@zju.edu.cn 1 1 1 2 3 Key Laboratory of Biomedical Engineering of Ministry of Education, Department of Biomedical Engineering, Zhejiang University, Hangzhou 310027, China Key Laboratory of Biomedical Engineering of Ministry of Education, Department of Biomedical Engineering, Zhejiang University, Hangzhou 310027, China Center of Arrhythmia Diagnosis and Treatment, Cardiovascular Institute & Fuwai Hospital, Chinese Academy of Medical Sciences & Peking Union Medical College, Beijing 100037, China Center of Arrhythmia Diagnosis and Treatment, Cardiovascular Institute & Fuwai Hospital, Chinese Academy of Medical Sciences & Peking Union Medical College, Beijing 100037, China Biological Physics Group, School of Physics and Astronomy, The University of Manchester, Manchester M13 9PL, UK Biological Physics Group, School of Physics and Astronomy, The University of Manchester, Manchester M13 9PL, UK The aim of this study is to build two mathematical models of canine ionic currents specific to right atria and left atria. The canine left atria mathematical model was firstly modified from the Ramirez-Nattel-Courtemanche (RNC) model using the recently available experimental data of ionic currents and was further developed based on our own experimental data. A model of right atria was then built by considering the differences between right atria and left atria. The two developed models well reproduced the experimental data on action potential morphology, the rate dependence, and action potential duration restitution. They are useful for investigating the mechanisms underlying the heterogeneity of canine regional action potentials and would help the simulation of whole heart excitation propagation and cardiac arrhythmia in the near future. Atrial fibrillation Canine atria Mathematical model Action potential duration restitution Introduction 1 Many heart diseases, such as myocardial infarction, heart failure, and ventricular fibrillation, are fatal. Research on understanding the mechanisms behind becomes significantly important. As an important complementarity of clinical research, different computer models from cellular to whole heart level and simulations of electrophysiology of cardiac myocyte were introduced to better understand the cellular pathological mechanism and sequentially to make predictions of potential clinical changes due to variations at subcellular level. After Hodgkin-Huxley equations were described (Hodgkin and Huxley, 1952 1952 ), some mathematical models of diverse cells have been proposed. The first cardiac cell model was built by Noble ( 1960 1960 ) for Purkinje fibre cells. He used the total current flowing through ion channels in the membrane to determine the transmembrane potentials. Since then, models have been built for ventricular, atrial, and sinoatrial node action potentials in different species. Due to the severity of ventricular arrhythmia compared to atrial arrhythmia, modeling of ventricular myocyte was firstly investigated. Beeler and Reuter ( 1977 1977 ) built the first action potential model of mammalian ventricular cell (B-R model). Their model solved a differential equation for each current and transmembrane potential. Based on B-R model, Luo and Rudy ( 1991 1991 ; 1994a 1994a ; 1994b 1994b ) and Rudy and Luo ( 1993 1993 ) developed a new model of ventricular cell and made progress in providing framework for modeling other types of ventricular cells. Recently, a mathematical model of action potential (AP) of human ventricular myocytes was published by ten Tusscher et al. ( 2004 2004 ). There are also mathematical models of AP of rabbit sinoatrial node cells published by Demir et al. ( 1994 1994 ) and Zhang et al. ( 2001 2001 ). Compared with the modeling of ventricular myocytes, there is relatively less concern about the modeling of atrial myocytes partially because of its structural complexity. Recently, with growing interest in atrial diseases, such as atrial fibrillation (AF), more and more research studies have been conducted on the atrium. Consequently, some mathematical models of atrial myocyte based on different species, such as rabbit, canine, and human, have been published (Lindblad et al., 1996 1996 ; Courtemanche et al., 1998 1998 ; Nygren et al., 1998 1998 ; Ramirez et al., 2000 2000 ; Kneller et al., 2002 2002 ). These published models would help us better understand the physiological and pathological mechanisms of atrial region and are of great value for clinical research. However, due to the limited experimental condition and complicated structure of the atrium, it is difficult to obtain accurate experimental data of ionic currents; thus, the difficulties in developing these atrial models should be resolved. After introducing a model of human atrial cell (Courtemanche et al., 1998 1998 ), Ramirez et al. ( 2000 2000 ) built a model of canine atrial cell, named as the Ramirez-Nattel-Courtemanche (RNC) model. The RNC model is the first mathematical model of the canine atrial AP, which provides a useful tool to study the mechanisms underlying electrical activity of canine atrium. However, some results produced by the RNC model were different from other published data. For example, its resting membrane potential (~−83 mV) is more negative than that published by Li et al. ( 2001 2001 ) and Ehrlich et al. ( 2003 2003 ) (~−73 mV), and the kinetics of some ionic currents of the RNC model is also not the same as these newly available data. Furthermore, the RNC model does not account for the heterogeneous electrophysiology of diverse cardiac cells that have different AP shapes, AP durations, sets of ionic currents, etc. It has been recognized that the AP heterogeneity is a major determinant of atrial reentrant arrhythmias (Feng et al., 1998 1998 ); thus, it is necessary to build mathematical models specific to different regions of canine atrium, while the RNC model was only specific to the right atrial regions. The aim of this study is to build two mathematical models of canine ionic currents specific to both the right atrium (RA) and left atrium (LA). The canine LA mathematical model was firstly modified from the RNC model and refined by some recently available experimental data of ionic currents and our own experimental ionic current data. A model of the RA was then built by considering the differences between the RA and the LA. Finally, the agreement between the model simulation results and the data from different experiments was evaluated, and the results simulated with the two mathematical models were compared. Materials and methods 1 Animals 2 Eight adult mongrels (1–2 years old) of either sex weighing 18–26 kg were used. These dogs were subjected to 12-lead electrocardiogram to confirm sinus rhythm. Echocardiograms were performed to exclude the structural heart disease. The LA size and the left ventricular ejection fraction (LVEF) were also measured by echocardiogram. The animal experiments were reviewed and approved by the local authorities and were performed in the Centre of Arrhythmia Diagnosis and Treatment, Fuwai Hospital, Beijing, China. Cellular electrophysiography 2 Single atrial myocyte was dispersed from the right atrial appendage (RAA) and left atrial appendage (LAA) using a previously described method (Dun et al., 2003 2003 ). Cells in a 1-ml bath were continuously superfused (2–3 ml/min) with normal Tyrode’s solution containing 137 mmol/L NaCl, 5.4 mmol/L KCl, 1.0 mmol/L MgCl 2 , 2.0 mmol/L CaCl 2 , 0.33 mmol/L NaH 2 PO 4 , 10 mmol/L hydroxyethyl piperazine ethanesulfonic acid (HEPES), and 10 mmol/L glucose (pH 7.4, NaOH). The solution was bubbled with 100% O 2 . The whole-cell patch-clamp technique was used to record ionic currents in the voltage-clamp mode (Axonpatch 200B, Axon Instruments, Foster City, CA, USA). Data were sampled at 13.3 kHz with an A/D converter (Digidata 1200, Axon Instruments), and the recordings were filtered with a low-pass corner frequency of 2 kHz. Patch pipette resistances ranged from 1.0 to 2.0 MΩ, when filled with an internal solution. Junction potentials were zeroed before the pipette touched the cell. The internal solution for the L-type Ca 2+ current ( I CaL ) recording contained (mmol/L): CsCl, 120; MgCl 2 , 1.0; Mg-adenosine triphosphate (MgATP), 5.0; 1,2-bis( o -aminophenoxy) ethane- N , N , N ′, N ′-tetraacetic acid (BAPTA), 10; HEPES, 10; and tetraethylammonium chloride (TEA-Cl), 10 (pH 7.30, CsOH). The external solution for I CaL recording contained (mmol/L): choline-Cl, 137; CaCl 2 , 2.0; MgCl 2 , 1.0; HEPES, 5; glucose, 10; CsCl, 4.6; TEA-Cl, 10; and 4-AP, 5 (pH 7.30, CsOH). In order to record the I CaL , currents were elicited by 240-ms voltage stepping to the true potential difference across the membrane ( V t ) −60 and +60 mV from a holding potential of −70 mV at 0.2 Hz. The internal solution for the transient outward K + current ( I to ) recording contained (mmol/L): KCl, 20; K-aspartate, 120; MgCl 2 , 1.0; ethyleneglycol bis(2-aminoethyl ether)tetraacetic acid (EGTA), 10; Na 2 -phosphocreatine, 2.0; K 2 ATP, 4.0; and HEPES, 5 (pH 7.30, KOH). The external solution for I to recording contained (mmol/L): choline-Cl, 137; KCl, 5.4; MgCl 2 , 1.0; NaH 2 PO 4 , 0.33; HEPES, 10; and glucose, 10 (pH 7.30, NaOH). Currents were elicited by 250-ms voltage stepping to V t −40 and +60 mV from a holding potential of −90 mV at 0.2 Hz. Data acquisition started 10 min after membrane rupture. Voltages were not corrected for liquid junction potentials between the bath and pipette solutions. I CaL and I to magnitudes were normalized by each cell’s membrane capacitance (pF) and expressed as current density (pA/pF). Statistical analysis 2 Statistical analysis was performed by SPSS 10.5 software (IBM, USA). All values were expressed as the mean±standard deviation (SD). Comparisons among the groups were made using the one-way analysis of variance (ANOVA) and the Student-Neuman-Keuls procedure. A P value less than 0.05 was considered statistically significant. Model description 2 As described by the RNC model, the cell membrane, as a capacitor, is connected in parallel with variable resistances and batteries representing the different ionic currents and pumps (Xia et al., 2006 2006 ). In this work, the following differential equation was used to describe the electrophysiological behavior of a myocyte cell: $$\frac{{{\text{d}}V}}{{{\text{d}}t}} = - ({I_{{\text{ion}}}} + {I_{{\text{stim}}}})/{C_{\text{m}}}$$ , where V is transmembrane voltage, t is time, I ion is the total transmembrane ionic current, I stim is the stimulus current applied externally, and C m is the cell membrane capacitance. Furthermore, I ion can be calculated by the following equation: $${I_{{\text{ion}}}} = {I_{{\text{Na}}}} + {I_{{\text{K1}}}} + {I_{{\text{to}}}} + {I_{{\text{Kur,d}}}} + {I_{{\text{Kr}}}} + {I_{{\text{Ks}}}} + {I_{{\text{CaL}}}} + {I_{{\text{Cl,Ca}}}} + {I_{{\text{p,Ca}}}} + {I_{{\text{NaCa}}}} + {I_{{\text{NaK}}}} + {I_{{\text{b,Na}}}} + {I_{{\text{b,Ca}}}}$$ , where the definitions of these currents are similar to those in the RNC model (Ramirez et al., 2000 2000 ): I Na is the fast Na + current, I K1 is the inward rectifier K + current, I to is the transient outward K + current, I Kur,d is the dog ultra-rapid delayed rectifier K + current, I Kr and I Ks are the rapid and slow delayed rectifier K + current components, respectively, I CaL is the L-type Ca 2+ current, I Cl,Ca is the Ca 2+ -activated Cl − current, I p,Ca is sarcoplasmic Ca 2+ pump current, I NaCa is the Na + /Ca 2+ exchanger current, I NaK is the Na + -K + pump current, and I b,Na and I b,Ca are the background Na + and Ca 2+ currents, respectively. The general approach to model the electrical AP of canine atrial myocytes is similar to that of the RNC model (Ramirez et al., 2000 2000 ); however, our model is based on more recently available experimental data of some major ionic currents from canine atrium, including I K1 , I Kr , I Ks , I to , I Kur,d , and I CaL . A detailed description of these changes and the experimental data used in this work is given in the Section 2.5. Physical units used in this work are the same as the RNC model, and the computer software encoding the model is written in Fortan 95 and Matlab 7.0 with the use of double-precision arithmetic. The modified equations are given in the Appendix. Model modification 2 A number of articles have been published regarding electrophysiological properties of canine atrial myocytes (Feng et al., 1998 1998 ; Wang et al., 1999 1999 ; Li et al., 2000 2000 ; 2001 2001 ; Dun et al., 2003 2003 ; Ehrlich et al., 2004 2004 ; Cha et al., 2005 2005 ). Certain ionic currents in the membrane were investigated with their corresponding experimental data published. Considering these valuable results and using the experimental data measured directly from isolated canine atrial cells, we reformulated some equations describing the ionic dynamic of the channels across the cell membrane with the use of Origin Pro 7.5, while the other equations remained unchanged. Table 1 1 summarizes the major ionic currents used in our modified model and RNC model. The details for each current were discussed below. The data sources were also listed. For the optimization of the model equations to fit the data, Boltzman equations of Origin 7.5 (OriginLab, Northampton, MA, USA) were used and asymptotic-symmetry-based method was applied. The results of RNC model and our modified model were finally plotted together for comparison. 2 Current 2 Conductance G 2 Data source Our model RNC model I K1 G K1 =0.090 nS/nF (LA/RA) G K1 =0.15 nS/nF Our experiment; Li et al., 2000 2000 ; 2001 2001 ; Ehrlich et al., 2003 2003 ; 2004 2004 I Kr G Kr =0.065 nS/nF (LA), G Kr =0.052 nS/nF (RA) G Kr =0.06984 nS/nF Our experiment; Feng et al., 1998 1998 ; Li et al., 2000 2000 ; 2001 2001 ; Ehrlich et al., 2003 2003 ; Cha et al., 2005 2005 I Ks G Ks =0.055 nS/nF (LA/RA) G Ks =0.0561 nS/nF Our experiment; Cha et al., 2004 2004 ; 2005 2005 ; Ehrlich et al., 2003 2003 ; Li et al., 2000 2000 ; 2001 2001 ; Feng et al., 1998 1998 I to G to =0.145 nS/nF (LA/RA) G to =0.19824 nS/nF Our experiment; Cha et al., 2004 2004 ; 2005 2005 ; Ehrlich et al., 2003 2003 ; Li et al., 2000 2000 ; 2001 2001 ; Feng et al., 1998 1998 ; Yue et al., 1997 1997 I Kur,d We formed a new equation for the G Kur,d based on recently available literature and our results. Ramirez et al., 2000 2000 Our experiment; Li et al., 2001 2001 ; Ehrlich et al., 2003 2003 ; Cha et al., 2005 2005 I CaL G CaL =0.27 nS/nF (LA/RA) G CaL =0.24 nS/nF Our experiment; Cha et al., 2004 2004 ; 2005 2005 ; Ehrlich et al., 2003 2003 ; Li et al., 2000 2000 ; Yue et al., 1997 1997 Left atrium 3 1. Inward rectifier K + current ( I K1 ) 4 I K1 plays a significant role in maintaining the resting membrane potential and in shaping the late repolarization phase of the AP (Zobel et al., 2003 2003 ). I K1 can be calculated by the following formulation: $${I_{{\text{K1}}}} = \frac{{{G_{{\text{K1}}}} \cdot ({V_{\text{m}}} - {E_{\text{K}}} - 10)}}{{0.5 + \exp \left( {0.07({V_{\text{m}}} + 75)} \right)}}$$ , where G K1 is the maximal conductance, V m is the transmembrane potential, and E K is equilibrium potential for K + . In the RNC model, the equation of I K1 was similar to that of its human AP model, but with some maximal conductance adjusted. Fig. 1 1 summaries the published experimental data of I K1 . The I - V relationship from our modified LA model was also plotted in Fig. 1 1 with its peak point marked. When compared to the RNC model, the peak point of our data is about 10 mV positive; however, our results match better with the experimental data (Li et al., 2000 2000 ; 2001 2001 ; Ehrlich et al., 2003 2003 ; 2004 2004 ). In addition, the experimental data show a lower peak-amplitude in I - V curve, and a decreased maximal conductance of I K1 (Table 1 1 ) was used in our work. Corresponding experimental values are also shown. Error bar is the scale of experimental data from our group and other investigators (Table 1 1 ), and filled squares are values used in our model to fit the I - V relationship 2. Rapid delayed rectifier K + current ( I Kr ) 4 I Kr plays an important role in repolarization, which can be given by the following equation: $${I_{{\text{Kr}}}} = {G_{{\text{Kr}}}} \cdot {x_{\text{r}}} \cdot {p_{\text{i}}} \cdot ({V_{\text{m}}} - {E_{\text{K}}})$$ , where x r is the activation gate and p i is the inactivation gate. The experimental data of voltage-dependent activation steady states for this current are rarely available, and there is no significant difference between the values used in the RNC model (Li et al., 2000 2000 ) and our own experimental data. Also there is an obvious difference between RA and LA. These experimental data were summarized from several independent researchers (Feng et al., 1998 1998 ; Li et al., 2000 2000 ; 2001 2001 ; Ehrlich et al., 2003 2003 ; Cha et al., 2005 2005 ). In this paper, the experimental data of peak of tail currents were used to get the peak I - V relationship for I Kr . To fit I Kr to the experimental I - V curves without significant changes in steady state curves, we decreased G Kr in Eq. ( 4 4 ) (Table 1 1 ). The simulated current traces using a voltage-clamp method are shown in Fig. 2a 2a , and the peak I - V relationships from both simulation and experiments are shown in Fig. 2b 2b . The experimental values of RA and LA were also presented. The error bars indicate their variability obtained from our group and other investigators (Table 1 1 ) 3. Slow delayed rectifier K + current ( I Ks ) 4 I Ks is a very slowly activating current, which does not show appreciable inactivation. It can be described as: $${I_{{\text{Ks}}}} = {G_{{\text{Ks}}}} \cdot x_{\text{s}}^2 \cdot ({V_{\text{m}}} - {E_{\text{K}}})$$ , where x s is the activation variable, but no inactivation variable is defined here. The result of previous investigation (Li et al., 2001 2001 ) suggests that there is no significant difference of I Ks in dynamics of activation gate and the current density between LA and RA, so the same maximal current densities of I Ks were used for both LA and RA. Fig. 3 3 shows the simulated current traces with LA model by using voltage-clamp experiments (Fig. 3a 3a ) and the peak I - V relationships from both simulation and experiments (Fig. 3b 3b ). Similar to I Kr , there is no peak current for I Ks , and the peak of tail currents from the experimental data was used to get the peak I - V relationship for I Ks . The experimental data were also summarized from the studies conducted by several investigators (Table 1 1 ). G Ks of 0.055 nS/nF was selected for both RA and LA to fit the experimentally obtained I - V relationships. The corresponding I - V relationship from the RNC model is also shown, which is slightly higher than that from our LA model. The experimental values of LA were also presented. The error bars indicate the variability of the experimental data obtained by our group and other investigators (Table 1 1 ) 4. Transient outward K + current ( I to ) 4 Based on the measured data obtained in isolated canine atrial myocytes, formulation of I to is given below: $${I_{{\text{to}}}} = {G_{{\text{to}}}} \cdot {o_{\text{a}}} \cdot {o_{\text{i}}} \cdot ({V_{\text{m}}} - {E_{\text{K}}})$$ , where o a and o i are voltage-dependent activation and inactivation gates. The steady-state activation curve is fitted to data of canine atrial myocytes from our experiments and other investigators (Table 1 1 ). Fig. 4a 4a summarized the experimental data with their variability given. Because the activation between LA and RA cells was almost the same, the same formulation was therefore used. The steady-state activation curves used in our LA model and the RNC model are also shown in Fig. 4a 4a . The activation time constants (Fig. 4b 4b ) were derived from the RNC model due to its similarity. Fig. 4c 4c also shows the steady-state inactivation curves of our LA model and the RNC model. Fig. 4d 4d shows the inactivation time constants. Fig. 4e 4e shows the simulation result of voltage-clamp experiments of LA model. Fig. 4f 4f is of peak I - V relationships of model simulations from the RNC model and our LA model, with experimental data shown together for comparison. Note, in Fig. 4f 4f , only the I - V relationship for LA is shown. Since experimental data for I - V relationships show no significant difference between RA and LA, we therefore used the same G to values to fit the experimental results (Table 1 1 ). The experimental data for the relevant parameters are also included. They were summarised from our own experiments and other published studies (Table 1 1 ) 5. Ultra-rapid delayed rectifier K + current ( I Kur,d ) 4 I Kur,d was reported later than I Kr and I Ks . Because of diverse experimental specimen, different terms were once applied, such as steady state current ( I ss ) and sustained potassium current ( I sus ) (Wang et al., 1993 1993 ; Dun et al., 2003 2003 ). I Kur,d activates much earlier and more rapidly than both I Kr and I Ks , and it may be of special importance in early plateau stage where I Kr and I Ks are not yet activated. The following formulation was used in our study: $${I_{{\text{Kur,d}}}} = {G_{{\text{Kur,d}}}} \cdot {u_{\text{a}}} \cdot {u_{\text{i}}} \cdot ({V_{\text{m}}} - {E_{\text{K}}})$$ , $${G_{{\text{Kur,d}}}} = 0.006 + \frac{{0.0325}}{{1 + \exp \left( {({V_{\text{m}}} + 20)/( - 12)} \right)}}$$ , where G Kur,d is the maximal conductance, which is voltage-dependent, u a is the activation variable, and u i is the inactivation variable. Again, no significant differences in kinetic or voltage-dependent properties between LA and RA were found in our experiments or by other studies (Li et al., 2001 2001 ; Ehrlich et al., 2003 2003 ; Cha et al., 2005 2005 ). Therefore, we only adjusted the equations based on the RNC model to fit the experimental data. However, according to the experimental data, the maximal current density is smaller, compared to the RNC model. In order to fit to the I - V relationship, the expression of G Kur,d was reformatted (Eq. ( 8 8 )). Fig. 5b 5b plotted the peak I - V curves of our LA model and the RNC model together with experimental data that were derived from Li et al. ( 2000 2000 ; 2001 2001 ) and Feng et al. ( 1998 1998 ). Our experimental results agree with the published data, while the result of RNC model shows a big difference. The experimental data summarized from Li et al. ( 2000 2000 , 2001 2001 ) and Feng et al. ( 1998 1998 ) are also shown as error bar scales 6. L-type Ca 2+ current ( I CaL ) 4 I CaL is known as the initiator of excitation-contraction coupling which is the most important mechanism in heart cells. I CaL also plays an essential role in shaping the AP morphology, especially in the plateau phase. Our formulation of I CaL is given by $${I_{{\text{CaL}}}} = {G_{{\text{CaL}}}} \cdot d \cdot f \cdot {f_{{\text{Ca}}}} \cdot ({V_{\text{m}}} - 60)$$ , where f and f Ca are voltage- and Ca 2+ -dependent inactivation gates and d is the voltage-dependent activation gate. Both the steady-state activation curve and inactivation curve are based on steady-state data from canine atrial myocytes reported by Cha et al. ( 2004 2004 ; 2005 2005 ), Ehrlich et al. ( 2003 2003 ), Li et al. ( 2000 2000 ; 2001 2001 ), and Feng et al. ( 1998 1998 ). Simulated results of our LA model of steady-state activation and inactivation curves are shown in Figs. 6a and 6c 6a and 6c . The corresponding results of the RNC model and experimental data were also included for comparison. Because the experimental data of activation time constants are rarely available for this current, the curve derived from the RNC model was used, as shown in Fig. 6b 6b . Next, the summary of these experimental data of inactivation time constants (Cha et al., 2004 2004 ; 2005 2005 ; Ehrlich et al., 2003 2003 ; Li et al., 2000 2000 ; Yue et al., 1997 1997 ) shows that RA cells have a larger inactivation time constant than LA cells in the scale from −10 to 50 mV (Fig. 6d 6d ); we therefore used two distinct expressions to model these two kinds of cells. Fig. 6e 6e shows I CaL of LA elicited by a simulated pulse protocol. The peak I - V relationships simulated with our LA model are shown in Fig. 6f 6f , which agree well with published experimental data. We also realized that these experimental data varied in a much larger scale than other currents, so median values were taken and assumed as typical data. For a comparison, the peak I - V relationship of RNC model was also included. From the Fig. 6f 6f , it can be seen that the I - V curve from the RNC model is about 10 mV negatively shifted. The experimental data for the relevant parameters are also included. They were summarised from our own experiments and other published studies (Table 1 1 ) Right atrium 3 No significant difference was found in most of the major ionic currents between LA and RA cells, except for the I Kr . Findings of Li et al. ( 2001 2001 ) suggest that the density of I Kr was greater in the LA than in the RA, while I Ks , I Kur,d , I K1 , I CaL , and I to were all comparable in the LA and RA. There were no significant differences in kinetic or voltage-dependent properties of currents in LA versus RA (Li et al., 2001 2001 ). Our experimental results are in good agreement with the published findings. Therefore, the RA model was constructed by decreasing the maximal conductance of I Kr to fit the amplitude of the experimental values (Fig. 2 2 ), while keeping the other conditions the same with LA model. Numerical algorithms 2 All simulations were written in Compaq Visual Fortran 95 and Matlab 7.0 and ran on a PC Intel Pentium4 2.80-GHz CPU with 2.00-GB Ram. Time derivatives were integrated with a double-precision Livermore solver for ordinary differential equations (DLSODE) by employing a maximal time step of 0.05 ms. Results 1 Canine AP model 2 To get the AP morphology under stable condition, the AP was recorded after the 10th stimulus (the stimulus is with a frequency of 1 Hz, a duration of 2 ms, and an amplitude of 2 nA). Fig. 7 7 shows the simulated AP, calcium transient, and some major ionic currents of RA model. Repeated simulation results in a stable state with intracellular ionic concentrations of Na + (11.7 mmol/L), K + (138.4 mmol/L), and Ca 2+ (0.0001 mmol/L), and a stable resting potential of −73.2 mV during diastoles of cardiac cycles, which is in good agreement with our experimental data and other published studies (Feng et al., 1998 1998 ; Li et al., 2001 2001 ). As shown in Fig. 7a 7a , the AP of LA model had shorter action potential duration (APD) than RA, while the resting potential and the AP overshot were similar between RA and LA. These results were fit to the experimental data reported by other investigators as mentioned previously (Li et al., 2001 2001 ; Ehrlich et al., 2003 2003 ; Cha et al., 2005 2005 ). The shorter APD in LA model was caused by the larger I Kr in LA cells, which is an important current in repolarization of AP. The repolarization was faster (or a shorter APD) in LA than in RA, leading to more stable reentry circuits with a shorter intrinsic period in the LA, which implies the potential role that LA plays in the initiation and maintenance of AF. The differences of I K1 and I Ks between LA and RA models are consequences of the different repolarization speeds. No significant difference was seen in other major currents, as shown in Fig. 7 7 . Rate dependence of APD 2 The rate dependence of APD and refractoriness to changes in pacing rate play an essential role in excitation and propagation, particularly in the reentrant arrhythmias. Simulation results are shown in Table 2 2 , which were generated by using RA model at various pacing intervals. Typically, APDs to 90% and 50% repolarization (APD 90 and APD 50 ) were used to show the change of APD. It can be seen that there is no significant increase in APD when the simulation period is larger than 1000 ms. However, as the simulation period decreases, the APD obviously shortens. In the case of pacing interval at 200 ms, AP overshoot (APO) and maximal AP upstroke velocity ( v max ) were significantly smaller than those in other cases, which may be caused by the incompleteness of the recoveries of ionic channels, especially for sodium current, due to a very short simulation compared to the AP excitation. These results suggest that shortening of APD to a certain threshold will cause significant variation in electrophysiological properties of cardiac cells, and therefore will result in serious heart problems. APD 90 : APD to 90% repolarization; APD 50 : APD to 50% repolarization; APO: AP overshoot; APA: AP amplitude; RMP: resting membrane potential; v max : maximal AP upstroke velocity Pacing interval (ms) APD 90 (ms) APD 50 (ms) APO (mV) APA (mV) RMP (mV) v max (V/s) 200 136.9 50.6 41.7 113 −71.3 283 300 150.0 63.7 44.8 117 −72.4 301 500 163.2 85.6 44.9 117 −72.7 316 1000 169.9 99.3 44.9 118 −73.1 318 3000 170.3 101.9 44.7 118 −73.4 318 Our simulation results were compared with the mean values of APD 90 in vitro of canine right atrial cells obtained by Li et al. ( 2000 2000 ; 2001 2001 ), as shown in Fig. 8 8 . It can be seen that at a lower frequency our simulation results agree well with the experimental data. However, their difference increases as the pacing frequency increases, likely due to the imperfect calcium dynamic in our models. In the present study, we reformulated new equations for several ionic currents (e.g., I K1 , G Kur,d ) but directly adopted the calcium handling from RNC model, which may cause the incompatibility or instability of the whole AP model. Weak calcium dynamic results in a weak response of cardiac cells to stimulus pacing, so refining our models by improving the modeling of calcium handling needs to be considered in the future. APD restitution 2 The restitution of action potential duration (APDR) has recently been demonstrated in both experimental and mathematical studies to be an important determinant of cardiac arrhythmias (Karma, 1994 1994 ; Koller et al., 1998 1998 ; Qu et al., 1999 1999 ; Garfinkel et al., 2000 2000 ). It is hence an essential AP property and we also simulated the APDR curves with a S1-S2 protocol here. The S1-S2 restitution protocol consists of 10 S1 stimuli applied at a frequency of 1 Hz and a strength of two times the threshold value (2000 pA), followed by a S2 extra-stimulus delivered at some diastolic interval (DI) after the AP generated by the last S1 stimulus (ten Tusscher et al., 2004 2004 ). During the simulation, the DI was decreased from 500 ms to 10 ms, and the APD generated by S2 stimulus was recorded simultaneously. APDR curve was then obtained by plotting APD against DI. Fig. 9 9 shows the APDR curves for both the RA and LA cells. It can be seen that the APD of RA model is longer than that of LA model, which matches the results shown in Fig. 7a 7a . For the slopes of the APDR curve, no significant difference between RA and LA was found. Discussion 1 Mathematical models serve as an important complement to experimental work in attempts to reveal the ionic mechanisms underlying the AP and other electrophysiological phenomena, with which experimental and clinical data can be integrated and reviewed as a holistic system (Nygren et al., 1998 1998 ). Once a model has been developed and validated, it can be used to predict the response of the cell to certain drugs and different experimental conditions (Nygren et al., 1998 1998 ). Furthermore, a specific model for a distinct type of cell is necessary since the heterogeneity of electrophysiological property plays an important role in cardiac arrhythmias (Feng et al., 1998 1998 ). In the present study, based on the RNC model (Ramirez et al., 2000 2000 ), which is only specific to canine right atrial region, we firstly introduced a model of left atrial myocytes, using data from our own experiment and other previously published data. To our knowledge, no mathematical model was built specifically for the left atrial cells. Our LA model reproduced a variety of behaviors that have been observed in experimental research. For example, the resting potential for LA cells is around −73 mV, which is in good accordance with published data (Feng et al., 1998 1998 ; Li et al., 2001 2001 ). The LA model also shows good AP morphology (Fig. 7a 7a ) and rate dependence, which is also in good agreement with published experimental values (Ehrlich et al., 2004 2004 ; Cha et al., 2005 2005 ). Next, by considering the differences between the RA and LA, a mathematical model of the RA with some appropriate modifications was developed. Note that, the general structure of our RA model follows the RNC model. For example, being unaware of detailed underlying mechanism, we derived the calcium dynamic in our models directly from the RNC model. Second, some membrane currents that are in absence of experimental data were simply derived from the RNC model, e.g., I NaK and I NaCa . Furthermore, I Na is so well-modeled that we just adopted the equations in the RNC model with the purpose to keep the whole model stable. However, there are some important differences between the RNC model and our RA model. The latest experimental data were used in our RA model, which is more convincing. The major ionic currents of the RNC model were refined to fit to our measurement results of canine atrial myocytes and/or recent data published. These updates resulted in several differences between the RNC model and our RA model: (1) the more positive reverse potential of I - V relationship of I K1 , which leads to a more positive resting membrane potential (~−72 mV) and agrees better with the available experimental data (Li et al., 2000 2000 ; 2001 2001 ; Ehrlich et al., 2003 2003 ; 2004 2004 ); and (2) smaller I to and I Kur,d , while larger I Kr and I CaL , which lead to a slightly wider plateau and longer APD. Limitations 2 The summarized experimental data from different published studies were used to get a set of equations that describe the change of ionic currents and concentrations along with time. However, limitations are inevitable due to the following reasons: (1) there are insufficient experimental data for some ionic currents, or even the underlying mechanism is still under investigation; (2) experimental data are largely variable due to individual heterogeneity of experimental material; (3) the experimental conditions are very difficult to control; and (4) the various protocols used in experiments (e.g., cell isolation) are likely to result in different experimental findings. The I Kur,d is the specific current in atrial myocytes, which is very important in the repolarization of the AP. However, there are no data obtained in the measurement on I Kur,d inactivation time constant because of its non-inactivating characteristics. The formulation for inactivation time constant used in the RNC model follows its human AP model of atrium. Although existence of difference between canine atrial cells and human atrial cells is obvious, it is one practical way to build the canine AP model by using equations in human AP model. In this study, we just followed the RNC model; that is, we adopted their equation for inactivation time constant. Unfortunately, there are currently no data available for I CaL activation time constant specific to canine atria, so the formulation used in the RNC model had to be adopted. To our knowledge, the precise mechanism of calcium-dependent inactivation of I CaL is still unclear, so the mathematical expression of calcium-dependent inactivation was described with a function of intracellular calcium (see in Appendix). Actually, calcium handling is one of the most important mechanisms in the cardiac cells, due to its function in excitation-contraction coupling (ECC), which translates electrical activity into mechanical activity. But it is also one of the most difficult issues during the modeling of AP. Obviously, dysfunction of calcium handling causes abnormality in AP, such as delayed afterdepolarization (DAD). Thus, further experimental and theoretical studies are needed to better model the mechanism of calcium handling. Furthermore, the experiments performed by different researchers are generally performed in different conditions, for example, different temperatures. The nature of myocytes may also be affected by different isolation methods and pharmacological agents. So the averaged experimental data from different published studies were used in our model (Cha et al., 2004 2004 ; 2005 2005 ; Ehrlich et al., 2003 2003 ; Li et al., 2000 2000 ; 2001 2001 ; Feng et al., 1998 1998 ; Yue et al., 1997 1997 ) and also compared with our own data and modeled results. Potential significance 2 It is known that AF is commonly-observed arrhythmia in clinical practice, and that AP heterogeneity plays an essential role in the generation and maintenance of atrial reentrant arrhythmias (Yue et al., 1997 1997 ; Fareh et al., 1998 1998 ; Ramirez et al., 2000 2000 ). Therefore, it is of considerable value of our work to build models specific to different regions in canine atria. We proposed two mathematical models for RA and LA cells based on ionic current data measured in canine atrial cells by ourselves and some recent published ionic current data. Our models reproduced a number of experimental observations, such as I - V curves of ionic currents and AP morphologies, and provided a valuable tool for analyzing the different roles of ionic currents in canine atrial APs. The repolarization was faster (or a shorter APD) in LA than in RA, leading to more stable reentry circuits with a shorter intrinsic period in the LA, which implies the potential role that LA plays in the initiation and maintenance of AF. In the future, with the combination of the anatomic models of the whole heart, our models can be used in simulation of excitation and propagation at 3-dimensional (3D) level in a more precise approach. These works will certainly give potential insights into more accurate mechanisms underlying cardiac arrhythmias and even assist with diagnosis and treatment in daily clinic practice. References GW Beeler , H Reuter . Reconstruction of the action potential of ventricular myocardial fibres . The Journal of Physiology , 1977 . 268 ( 1 ): 177 - 210 . 874889 . Reconstruction of the action potential of ventricular myocardial fibres TJ Cha , JR Ehrlich , L Zhang , S Nattel . Atrial ionic remodeling induced by atrial tachycardia in the presence of congestive heart failure . Circulation , 2004 . 110 ( 12 ): 1520 - 1526 . 15381657 DOI: 10.1161/01.CIR.0000142052.03565.87 http://doi.org/10.1161/01.CIR.0000142052.03565.87 . Atrial ionic remodeling induced by atrial tachycardia in the presence of congestive heart failure TJ Cha , JR Ehrlich , L Zhang , D Chartier , TK Leung , S Nattel . Atrial tachycardia remodeling of pulmonary vein cardiomyocytes: comparison with left atrium and potential relation to arrhythmogenesis . Circulation , 2005 . 111 ( 6 ): 728 - 735 . 15699259 DOI: 10.1161/01.CIR.0000155240.05251.D0 http://doi.org/10.1161/01.CIR.0000155240.05251.D0 . Atrial tachycardia remodeling of pulmonary vein cardiomyocytes: comparison with left atrium and potential relation to arrhythmogenesis M Courtemanche , RJ Ramirez , S Nattel . Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model . American Journal of Physiology , 1998 . 275 ( 5 ): H301 - H321 . 9688927 . Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model SS Demir , JW Clark , CR Murphey , WR Giles . A mathematical model of a rabbit sinoatrial node cell . American Journal of Physiology , 1994 . 266 ( 3 ): C832 - C852 . 8166247 . A mathematical model of a rabbit sinoatrial node cell W Dun , T Yagi , MR Rosen , PA Boyden . Calcium and potassium currents in cells from adult and aged canine right atria . Cardiovascular Research , 2003 . 58 ( 3 ): 526 - 534 . 12798425 DOI: 10.1016/S0008-6363(03)00288-8 http://doi.org/10.1016/S0008-6363(03)00288-8 . Calcium and potassium currents in cells from adult and aged canine right atria JR Ehrlich , TJ Cha , L Zhang , D Chartier , P Melnyk , H Hohnloser , S Nattel . Cellular electrophysiology of canine pulmonary vein cardiomyocytes: action potential and ionic current properties . The Journal of Physiology , 2003 . 551 ( 3 ): 801 - 813 . 12847206 DOI: 10.1113/jphysiol.2003.046417 http://doi.org/10.1113/jphysiol.2003.046417 . Cellular electrophysiology of canine pulmonary vein cardiomyocytes: action potential and ionic current properties JR Ehrlich , TJ Cha , L Zhang , D Chartier , L Villeneuve , TE Hebert , S Nattel . Characterization of a hyperpolarization-activated time-dependent potassium current in canine cardiomyocytes from pulmonary vein myocardial sleeves and left atrium . The Journal of Physiology , 2004 . 557 ( 2 ): 583 - 597 . 15020696 DOI: 10.1113/jphysiol.2004.061119 http://doi.org/10.1113/jphysiol.2004.061119 . Characterization of a hyperpolarization-activated time-dependent potassium current in canine cardiomyocytes from pulmonary vein myocardial sleeves and left atrium S Fareh , C Villemaire , S Nattel . Importance of refractoriness heterogeneity in the enhanced vulnerability to atrial fibrillation induction caused by tachycardia-induced atrial electrical remodeling . Circulation , 1998 . 98 ( 20 ): 2202 - 2209 . 9815876 . Importance of refractoriness heterogeneity in the enhanced vulnerability to atrial fibrillation induction caused by tachycardia-induced atrial electrical remodeling J Feng , L Yue , Z Wang , S Nattel . Ionic mechanisms of regional action potential heterogeneity in the canine right atrium . Circulation Research , 1998 . 83 ( 5 ): 541 - 551 . 9734477 . Ionic mechanisms of regional action potential heterogeneity in the canine right atrium A Garfinkel , YH Kim , O Voroshilovsky , Z Qu , JR Kil , MH Lee , HS Karagueuzian , JN Weiss , PS Chen . Preventing ventricular fibrillation by flattening cardiac restitution . Proceedings of the National Academy of Sciences , 2000 . 97 ( 11 ): 6061 - 6066 . DOI: 10.1073/pnas.090492697 http://doi.org/10.1073/pnas.090492697 . Preventing ventricular fibrillation by flattening cardiac restitution AL Hodgkin , AF Huxley . A quantitative description of membrane current and its application to conduction and excitation in nerve . The Journal of Physiology , 1952 . 117 ( 4 ): 500 - 544 . 12991237 . A quantitative description of membrane current and its application to conduction and excitation in nerve A Karma . Electrical alternans and spiral wave breakup in cardiac tissue . Chaos , 1994 . 4 ( 3 ): 461 - 472 . 12780121 DOI: 10.1063/1.166024 http://doi.org/10.1063/1.166024 . Electrical alternans and spiral wave breakup in cardiac tissue J Kneller , RJ Ramirez , D Chertier , M Courtemanche , S Nattel . Time-dependent transients in an ionically based mathematical model of the canine atrial action potential . American Journal of Physiology , 2002 . 282 ( 4 ): H1437 - H1451 . 11893581 . Time-dependent transients in an ionically based mathematical model of the canine atrial action potential ML Koller , ML Riccio , RF Gilmour . Dynamic restitution of action potential duration during electrical alternans and ventricular fibrillation . American Journal of Physiology , 1998 . 275 ( 5 ): H1635 - H1642 . 9815071 . Dynamic restitution of action potential duration during electrical alternans and ventricular fibrillation D Li , P Melnyk , J Feng , Z Wang , K Petrecca , A Shrier , S Nattel . Effects of experimental heart failure on atrial cellular and ionic electrophysiology . Circulation , 2000 . 101 ( 22 ): 2631 - 2638 . 10840016 . Effects of experimental heart failure on atrial cellular and ionic electrophysiology D Li , L Zhang , J Kneller , S Nattel . Potential ionic mechanism for repolarization differences between canine right and left atrium . Circulation Research , 2001 . 88 ( 11 ): 1168 - 1175 . 11397783 DOI: 10.1161/hh1101.091266 http://doi.org/10.1161/hh1101.091266 . Potential ionic mechanism for repolarization differences between canine right and left atrium DS Lindblad , CR Murphey , JW Clark , WR Giles . A model of the action potential and underlying membrane currents in a rabbit atrial cell . American Journal of Physiology , 1996 . 271 ( 4 ): H1666 - H1696 . 8897964 . A model of the action potential and underlying membrane currents in a rabbit atrial cell CH Luo , Y Rudy . A model of the ventricular cardiac action potential: depolarization, repolarization, and their interaction . Circulation Research , 1991 . 68 ( 6 ): 1501 - 1526 . 1709839 . A model of the ventricular cardiac action potential: depolarization, repolarization, and their interaction CH Luo , Y Rudy . A dynamic model of the cardiac ventricular action potential. I. Simulations of ionic currents and concentration changes . Circulation Research , 1994 . 74 ( 6 ): 1071 - 1096 . 7514509 . A dynamic model of the cardiac ventricular action potential. I. Simulations of ionic currents and concentration changes CH Luo , Y Rudy . A dynamic model of the cardiac ventricular action potential. II. Afterdepolarizations, triggered activity, and potentiation . Circulation Research , 1994 . 74 ( 6 ): 1097 - 1113 . 7514510 . A dynamic model of the cardiac ventricular action potential. II. Afterdepolarizations, triggered activity, and potentiation D Noble . Cardiac action and pacemaker potentials based on the Hodgkin-Huxley equations . Nature , 1960 . 188 ( 4749 ): 495 - 497 . DOI: 10.1038/188495b0 http://doi.org/10.1038/188495b0 . Cardiac action and pacemaker potentials based on the Hodgkin-Huxley equations A Nygren , C Fiset , L Firek , JW Clark , DS Lindblad , RB Clark , WR Giles . Mathematical model of an adult human atrial cell: the role of K + currents in repolarization . Circulation Research , 1998 . 82 ( 1 ): 63 - 81 . 9440706 . Mathematical model of an adult human atrial cell: the role of K+ currents in repolarization Z Qu , JN Weiss , A Garfinkel . Cardiac electrical restitution properties and stability of reentrant spiral waves: a simulation study . American Journal of Physiology , 1999 . 276 ( 1 ): H269 - H283 . 9887041 . Cardiac electrical restitution properties and stability of reentrant spiral waves: a simulation study RJ Ramirez , S Nattel , M Courtemanche . Mathematical analysis of canine atrial action potentials: rate, regional factors, and electrical remodeling . American Journal of Physiology , 2000 . 279 ( 4 ): H1767 - H1785 . 11009464 . Mathematical analysis of canine atrial action potentials: rate, regional factors, and electrical remodeling Y Rudy , CH Luo . Cellular responses to electrical stimulation: a study using a model of the ventricular cardiac action potential . Advances in Experimental Medicine and Biology , 1993 . 346 ( 1 ): 79 - 90 . 8184783 . Cellular responses to electrical stimulation: a study using a model of the ventricular cardiac action potential KHWJ ten Tusscher , D Noble , PJ Noble , AV Panfilov . A model for human ventricular tissue . American Journal of Physiology , 2004 . 286 ( 4 ): H1573 - H1589 . 14656705 DOI: 10.1152/ajpheart.00794.2003 http://doi.org/10.1152/ajpheart.00794.2003 . A model for human ventricular tissue Z Wang , B Fermini , S Nattel . Sustained depolarization-induced outward current in human atrial myocytes. Evidence for a novel delayed rectifier K + current similar to Kv1.5 cloned channel currents . Circulation Research , 1993 . 73 ( 6 ): 1061 - 1076 . 8222078 . Sustained depolarization-induced outward current in human atrial myocytes. Evidence for a novel delayed rectifier K+ current similar to Kv1.5 cloned channel currents Z Wang , J Feng , H Shi , A Pond , JM Nerbonne , S Nattel . Potential molecular basis of different physiological properties of the transient outward K + current in rabbit and human atrial myocytes . Circulation Research , 1999 . 84 ( 5 ): 551 - 561 . 10082477 . Potential molecular basis of different physiological properties of the transient outward K+ current in rabbit and human atrial myocytes L Xia , Y Zhang , H Zhang , Q Wei , F Liu , S Crozier . Simulation of Brugada syndrome using cellular and three-dimensional whole-heart modeling approaches . Physiological Measurement , 2006 . 27 ( 11 ): 1125 - 1142 . 17028406 DOI: 10.1088/0967-3334/27/11/006 http://doi.org/10.1088/0967-3334/27/11/006 . Simulation of Brugada syndrome using cellular and three-dimensional whole-heart modeling approaches L Yue , J Feng , R Gaspo , GR Li , Z Wang , S Nattel . Ionic remodeling underlying action potential changes in a canine model of atrial fibrillation . Circulation Research , 1997 . 81 ( 4 ): 512 - 525 . 9314832 . Ionic remodeling underlying action potential changes in a canine model of atrial fibrillation H Zhang , AV Holden , MR Boyet . Gradient model versus mosaic model of the sinoatrial node . Circulation , 2001 . 103 584 - 588 . 11157726 . Gradient model versus mosaic model of the sinoatrial node C Zobel , HC Cho , T Nguyen , R Pekhletski , RJ Diaz , GJ Wilson , PH Backx . Molecular dissection of the inward rectifier potassium current (IK1) in rabbit cardiomyocytes: evidence for heteromeric co-assembly of Kir2.1 and Kir2.2 . The Journal of Physiology , 2003 . 550 ( 2 ): 365 - 372 . 12794173 DOI: 10.1113/jphysiol.2002.036400 http://doi.org/10.1113/jphysiol.2002.036400 . Molecular dissection of the inward rectifier potassium current (IK1) in rabbit cardiomyocytes: evidence for heteromeric co-assembly of Kir2.1 and Kir2.2 The authors are very grateful to Prof. Alan MURRAY and Dr. Ding-chang ZHENG of Newcastle University, UK for their valuable suggestions and English language assistance in preparing the manuscript. 10.1631/jzus.B0900346 Project supported by the National Basic Research Program (973) of China (No. 2007CB512100), the National High-Tech R & D Program (863) of China (No. 2006AA02Z307), and the National Natural Science Foundation of China (No. 30570484) 2009-10-30 2010-05-11 2010-03-03 2022-01-05 Copyright © Zhejiang University and Springer-Verlag Berlin Heidelberg 2010 2010 Journal of Zhejiang University-SCIENCE B (Biomedicine and Biotechnology) 06 11