Repositorio Institucional de Documentos

An Improved Human Ventricular Cell Model for Investigation of Cardiac Arrhythmias under Hyperkalemic Conditions Author: Jesús Carro Fernández Directors: Esther Pueyo and José Félix Rodríguez Master's Degree in Biomedical Engineering December, 2010 Academic Year: 2010-2011 Communications Technology Group Group of Structural Mechanics and Material Modeling An Improved Human Ventricular Cell Model for Investigation of Cardiac Arrhythmias under Hyperkalemic Conditions Author: Jesús Carro Fernández Directors: Esther Pueyo and José Félix Rodríguez Master's Degree in Biomedical Engineering December, 2010 Academic Year: 2010-2011 Communications Technology Group Group of Structural Mechanics and Material Modeling To Débora. To my family. An Improved Human Ventricular Cell Model for Investigation of Cardiac Arrhythmias under Hyperkalemic Conditions Abstract The use of experiments for studying cardiac arrhythmias or the effect of drugs on cardiac electrophysiology is mostly limited to measurements obtained from electrograms (EGMs, measured on the heart surface) or, more often, electrocardiograms (ECGs, measured on the body surface). Despite the fact that many diagnostic and therapeutical decisions rely only upon interpretation of ECG patterns, the cellular and subcellular mechanisms underlying pathophysiological ECG changes remain mostly unclear. Among the different approaches aimed to connect the ECG with its underlying basis, multi-scale computational modeling of the heart arises as a powerful tool to understand cardiac functioning from the ionic to the whole organ level. With the increase in computational resources available to the scientic community, mathematical modeling and simulation of heart's electrical activity is becoming a fundamental tool to understand cardiac behavior. In this study several modications were introduced to a recently proposed action potential (AP) cell model so as to render it suitable for the study of ventricular arrhythmias. These modications were based on new experimental data and in the results of several cellular arrhythmic risk biomarkers reported in the literature. Five stimulation protocols were applied to the original and improved models of isolated cell, and a number of cellular arrhythmic risk biomarkers were computed. The stimulation protocol included a steady-state protocol, abrupt changes in cycle length (CL) protocol, S1S2 and dynamic restitution protocols, and concentration rate dependence protocol. In addition, the behavior of the proposed model under hyperkalemic conditions was simulated in a one dimensional ber by increasing the extracellular [K+], measuring the AP duration (APD), conduction velocity (CV) and effective refractory period (ERP) after steady-state conditions had been reached. Our modications led to: a) further improved AP triangulation (78.1ms); b) APD rate adaptation curves characterized by fast and slow time constants within physiological ranges (10.1sand 105.9s); c) maximum S1S2 restitution slope in accordance with experimental data (SS1S2 = 1.0). Underhyperkalemia, ourresultsshowedthatAPDprogressivelydecreasedwiththelevelofhyperkalemia, while ERP increased after a threshold in the extracellular [K+]was reached ([K+]o= 6 mM). Conduction velocity decreased with hyperkalemia and the conduction was blocked above [K+]o= 10.4mM. Above [K+]o= 9.8mM, alternans appeared in the APD. These results suggest that the longer ERP values and the conductionblock above[K+]o= 10.4mM foundinthe centralzoneofacutelyischemictissueas compared to the normal zone could create areas of block that could set a substrate for reentrant arrhythmias. 7 Mejora de un Modelo Celular de Ventrículo Humano para la Investigación de Arritmias Cardiacas bajo Condiciones de Hiperpotasemia Resumen La utilización de experimentos para el estudiode arritmias cardiacas o el efecto de medicamentos en la electrosiología cardiaca está limitada principalmente a medidas obtenidas de los electrogramas (EGMs, medidos en la supercie del corazón) o, más habitualmente, electrocardiogramas (ECGs, medidos en la superce del cuerpo). A pesar de que muchos diagnósticos y decisiones terapéuticas están fundadas solamente en la interpretación de los patrones de ECG, los mecanismos celulares y subcelulares subyacentes a cambios siopatológicos en el ECG, continúan no quedando claros en la mayoría de casos. Entre los diferentes acercamientos para conectar el ECG con los mecanismos subyacentes, el modelado computacional multiescala del corazón se plantea como una herramiente potente para entender el funcionamiento cardiaco, desde el nivel iónico al órgano completo. Con el incremento de los recursos computacionales disponibles para la comunidad cientíca, el modelado matemático y la simulación de la actividad eléctrica del corazón se está convirtiendo en una herramienta fundamental para comprender el comportamiento cardiaco. En este estudio, varias modicaciones fueron introducidas a un modelo de potencial de acción (AP) propuesto recientemente, para hacerlo apto para el estudio de arritmias ventriculares. Estas modicaciones se basaron en nuevos datos experimentales y en el resultado de varios biomarcadores arrítmicos celulares presentados en la literatura. Se aplicaron cinco protocolos de estimulación al modelo origninal y a la versión mejorada para calcular los diferentes biomarcadores de riesgo arrítmico. Los protocolos de estimulación incluían un protocolo de estado estático, un protocolo de cambios abruptos en el período de estimulación (CL), los protocolos de restitución S1S2 y dinámico, y un protocolo de depencencia de las concentraciones a la frecuencia cardiaca. Además, se simuló el comportamiento del modelo propuesto bajo condiciones de hiperpotasemia en una bra unidimensional mediante el incremento de la concentración extracelular de K+, midiendo la duracióndelAP(APD),lavelocidaddecondución(CV)yelperíodorefractarioefectivotrasalcanzarcondiciones de estado estable. Estas modicaciones llevaron a: a) mejorar aún más la triangulación del APD (78,1ms); b) Caracterización del APD mediante una constante de tiempo rápida y una lenta dentro del rango sológico (10,1sy 105,9s); c) máxima pendiente de la curva de restitución S1S2 en concordancia con los datos experimentales (SS1S2 = 1,0). Bajo condiciones hiperpotasémicas, los resulatos mostraron que el APD decrece progresivamente con el nivel de hiperpotasemia, mientras que el ERP aumenta tras un umbral en la concentración extracelular de K+([K+]o= 6 mM). La CV disminuyó con hiperpotasemia y la conducción se bloqueó por encima de [K+]o= 10,4mM. Por encima de [K+]o= 9,8mM, aparecieron alternantes en el APD. Estos resultados sugieren que los altos valores de ERP y el bloqueo de conducción sobre [K+]o= 10,4mM encontrado en la zona central de tegido isquémico agudo, puede crear un área de bloqueo que sea el substrato para un arrhitmia reentrante. 9 List of Tables 4.1 Biomarkers of arrhythmic risk for the TP06, GPB and CRLP human models. . . . . . . . . . . . 36 4.2 Percentages of variation in the APD caused by blocking different currents. . . . . . . . . . . . 37 17 Document Chapter 1 Introduction Cardiovascular disease (CVD) is the main cause of death in industrialised countries. More than 2 million people in the European Union (EU) and nearly 4.5 million people in all Europe die from CVD every year (European Heart Network 09), a number which is well above that of people dying from cancer, dementia and HIV/AIDS all together (World Health Organization 09). CVD is thus responsible for 42% of total deaths in the EU (38% of female, and 45% of male deaths). Spain, despite being one of the countries with the lowest rate of mortality due to CVD within Europe, is showing a tendency to a worst perspective, with a 2.1% increment per year in the number of such deaths. 1.1 Cardiac Arrhythmias Among CVD, cardiac arrhythmias represent a common cause of mortality. Cardiac arrhythmias are conditions in which there is a failure in the normal rhythm of the heartbeat. Arrhythmias can be caused either by an abnormal formation of the electrical impulse, by abnormalities in the conduction of that impulse throughout the heart, or by both. There are different types of cardiac arrhythmias. The so-called re-entrant arrhythmias occur when a propagating impulse fails to die out after normal activation of the heart and persists to re-excite other regions that have already recovered excitability. Some arrhythmias of this type are dangerous and may have fatal consequences. Frequently occurring re-entrant arrhythmias in the ventricles include monomorphic and polymorphic ventricular tachycardia (VT) and ventricular brillation (VF). Monomorphic VT and polymorphic VT (among which Torsades de Pointes) both lead to increased ventricular excitation and contraction rates and a decrease of cardiac output. In addition, both may destabilize into VF. VF leads to a further increase in excitation rate and a loss of coherence of ventricular contraction, which results in almost zero cardiac output. This will cause death to occur within minutes unless normal sinus rhythm is restored. In Europe sudden cardiac death (SCD) accounts for over 300,000 deaths per year, of which 75-80% are due to VF. VT and VF thus belong to the most dangerous cardiac arrhythmias. 1.2 Ischemia Ischemia isa pathological conditionin whichthe deliveryof substrates, mainlyoxygen, tothe myocardium is reduced. This causes a progressive deterioration of the electrical activity in the injured region. Three are the main components of the cellular alterations under ischemia: hyperkalemia, acidosis and hypoxia (Dodama et al., 1984; Gilmour and Zipes, 1980; Kagiyama et al., 1982; Kishida et al., 1979; Senges et al., 1979; Vleugels et al., 1980). Under hyperkalemia, potassium starts to pour out of ischemic myocardial cells causing a depolarization of themyocardiumatrest, slowing conduction and decreasing excitability. Under these conditions, recovery of excitability is known to outlast repolarization, a phenomenon termed post-repolarization refractoriness. 21 22 Jesús Carro Fernández Hypoxia causes a reduction in the ATP/ADP ratio than can slow down calcium pumps and exchangers (Tian and Ingwall, 1196), increasing intracellular sodium and calcium concentrations as well as activating ATP sensitive potasium channels (IKATP channels) (Carmeliet, 1999). Acidosis reduces the availability of sodium and calcium channels, reducing excitability and slowing conduction. Theseabnormalitiescauseinhomogeneities ofrestingpotentialanddispersion inrecoveryofexcitability in the tissue. In addition, these inhomogeneities set the stage for unidirectional block that can become the substrate for cardiac arrhythmias. 1.3 New Methods for Diagnosis and Treatment Due to the high incidence of CVD in the European population, and among this of cardiac arrhythmias, the associated cost with it turns out to be huge. It is estimated that CVD costs the EU economy €192,000 million a year, and €7,000 million a year to the Spanish Government, thus becoming the main cause of disease burden (European Heart Network 09). Based on these numbers it is clear that industrialised countries, and in particular Europe and Spain, need to face the challenge of appropriately managing prevention, diagnosis and treatment of cardiac diseases but at affordable cost. It is in this context where computational tools may provide relevant capabilities to improve quality and efficiency of healthcare services. Mathematical modeling and simulation of heart’s electrical activity (so-called cardiac electrophysiology) and signal processing of bioelectrical signals provide an ideal framework where to join the information from clinical and experimental studies with the understanding of the mechanisms underlying them so as to produce tools able to test different hypothesis and predict potential abnormalities in heart's behavior. In a relatively near future it is not unlikely that those tools might even become to start being used in clinical practice as complementary instruments to help in the prevention of cardiac diseases and the improvement of their diagnosis and therapy. 1.4 Computational Modeling and Simulation At present, clinical data is mostly limited to measurements obtained from electrograms (EGMs, measured on the heart surface) or, more often, electrocardiograms (ECGs, measured on the body surface). The main reason behind the extensive use of the ECG signal in the clinics can be found in the fact that, in addition to being a non-invasive procedure, is simple and of low-cost. Despite the fact that many diagnostic and therapeutical decisions rely only upon interpretation of ECG patterns, the cellular and subcellular mechanisms underlying pathophysiological ECG changes remain mostly unclear. Among the different approaches aimed to connect the ECG with its underlying basis, multi-scale computational modeling of the heart arises as a powerful tool to understand cardiac functioning from the ionic to the whole organ level. There are a number of important reasons that support the use of cardiac models to explain ECG-derived outcomes. First, modeling represents an alternative to experimental investigation, which is pretty much limited by ethical considerations, particularly when dealing with human hearts. Animal experimentation on other species, such as guinea pigs, rabbits or pigs, turns out to be more feasible, but corresponding results are not always easily transferable to humans, due to differences in key aspects such as heart size, heart rate, or duration and shape of the action potential (AP, representative of the electrical activity of a cell) (ten Tusscher et al., 2004). Another advantage of modeling is that it offers the possibility to test hypothesis or investigate conditions that are not easily reproducible in experiments, as it is the case with the study of three-dimensional phenomena, where experimental recording techniques offer limited possibilities due to poor spatial resolution, poor depth resolution, or applicability only to wedge preparations (ten Tusscher, 2004). Finally, computer models are of great importance in elds such as drug testing, where the use of simulations can help to detect adverse properties of compounds before their commercialization, and thus reduce the cost of bringing a new drug to the market. Computational modeling of the heart has experienced a great development during the last ve decades. In 1960 Noble formulated the rst electrophysiological model for cardiac myocytes (Noble, 1960), and since UniversidadZaragoza An improved human ventricular cell model 23 then more complex models of cellular electrical activity have been proposed for different species and regions of the heart (Courtemanche et al., 1998; Grandi et al., 2010; Luo and Rudy, 1991, 1994; Maleckar et al., 2009; Noble et al., 1998; ten Tusscher and Panlov, 2006). Additional modeling of the communication betweenneighbouring cellsmakesit possibletoextendthosecell modelsandconvert themintotissue oreven whole organ models. The electrical properties of the myocardium are generally described by the bidomain equations, a set of coupled parabolic and elliptic partial differential equations (PDEs) that represents the tissue as two separate, distinct continua - one intracellular and the other extracellular. The intracellular and the extracellular media are connected via the cell membrane, and thus the two PDEs are coupled at each point in space through a set of non-linear ordinary differential equations (ODEs), which describe the ionic transport across the cell membrane. In some cases, the monodomain representation of cardiac activity is used, which involves solving a single parabolic PDE, by assuming either that the extracellular potentials are negligible, or thattheanisotropyratios are equalinthe intracellularand theextracellulardomains (Bernabeu et al., 2009). 1.5 Human Ventricular Cell Models In the last years some mathematical models of electrical and ionic homeostasis in human ventricular myocytes have been proposed. One of the most extensively used models is the one proposed by ten Tusscher and Panlov (ten Tusscher and Panlov, 2006) (TP06), which is an improved version of the model described in (ten Tusscher et al., 2004) (TNNP04), where the calcium dynamics, the slow delayed rectier current (IKs) and the L-type calcium current (ICaL) were reformulated. The major advantage of the TP06 model is that it accurately reproduces the behavior of the AP duration restitution curve (see Section 2 for more information about this). Other relevant model for human ventricular cells was described in the paper by Iyer et al. (Iyer et al., 2004). In this modelthe mainionic currents are describedusing Markovianchains, thereasonforwhich the model becomes very complex, having more than 60 state variables. Recently, a new model of human ventricular AP has been proposed by Grandi et al. (Grandi et al., 2010) (GPB), which improves the response to frequency changes and has a better performance against current blocks with respect to the TP06 model. 1.6 This Work Despite the fact that the human ventricular AP models currently available in the literature are able to reproduce certain electrophysiological properties, they fail in representing other relevant aspects, as will be later described. In order to obtain a model suitable for the study of cardiac arrhythmias, in this study we have worked on building a new human ventricular AP model based on the analysis of well-known arrhythmic risk biomarkers. For that, we have used the GPB model as a starting point, and we have modied it by accounting for recent experimental measurements of potassium currents (Fink et al., 2008) and by introducing fast and slow voltage-dependent inactivation gates for the L-type calcium current (ten Tusscher and Panlov, 2006). All the modications were made with the aim of improving the model performance on several arrhythmic risk biomarkers (see Chapter 2). The proposed model, herein denoted by CRLP, has been tested against the GPB and TP06 models by applying four stimulation protocols and computing twelve cellular arrhythmic risk biomarkers following the methodology proposed in (Romero et al., 2009). The results show that the introduced modications have brought most biomarkers into the physiological range, and considerably improved others with respect to the GPB and TP06 models of AP. The proposed model has also been tested under hyperkalemic conditions to evaluate the feasibility of using this model for studying the acute regional ischemic heart. This document is structured in four chapters, each one with the following contents: •Chapter 2 introduces the importance of using arhythmic risk biomarkers to develop models and describes each of the arrhythmic risk biomarkers used in the evaluation of the model and how they are computed. BiomedicalEngineering 24 Jesús Carro Fernández •Chapter 3 establishes the basis for the development of the new CRLP model and describes the modi- cations introduced to the GPB model to render it more suitable for the study of cardiac arrhythmias. •Chapter 4 describes the simulation tests conducted to evaluate model's performance and presents the corresponding results. •Chapter 5 presents the discussion and the conclusions of the study. In addition, other contents can be found in the appendices to the main document: •Appendix A: Description of the CRLP human ventricular cell model. •Appendix B: Contributions of the study to international conferences. •Appendix C: Compilation of the acronyms used in the document. UniversidadZaragoza Chapter 2 Biomarkers of arrhythmic risk Arrhythmic risk biomakers are commonly used to predict drug cardioxicity in the development process of a new pharmacological compound. In a study developed by Romero et al. (Romero et al., 2009) an evaluation of how the principal ionic currents of a human ventricular cell model inuence different biomarkers was presented. Despite these evidences, the traditional way of developing a cellular model is based on isolated measurements on cells (i.e., individual ionic currents or intracellular concentrations). The exclusive use of current measurements leads to many uncertainties, as it only provides information about the behavior of each current, but not about the contribution of each current to the global cellular behavior. The main goal of this study is to propose a new methodology for the development of ionic cell models that joins the two types of approaches described above for model development. Specically, in this study biomarkers of arrhythmic risk were used to drive the modications of the model. While the formulation of individual currents was based on specic measurements, model parameters were adjusted based on different biomarkers that must remain within a physiological range. This allowed the proposed model to provide results in better concordance with experimental data. In the following sections, the main cellular arrhythmic risk biomarkers are presented, and the protocols used to measured them are described. 2.1 Selected Biomarkers In this study, a number of arrhythmic risk biomarkers proposed in the literature were selected to be used for model evaluation purposes, following the methodology described in (Romero et al., 2009): •APD: The AP duration (APD) is considered the main preclinical biomarker of drug cardiotoxicity becasuse of its link to long QT syndrome (Hondeghem et al., 2001; Volders et al., 2000). Unless otherwise specied, APD is used to denote AP duration at 90% repolarization. •Triangulation: This biomarker represents the shape of the end of the AP, calculated as the difference between the APD at 50% and 90% of the repolarization (Figure 2.1.a). Low triangulation values are indicative of square APs, while high values indicate triangular APs. It has been proposed as a predictor of serious proarrhythmia in (Hondeghem et al., 2001). •Systolic and diastolic [Ca2+]ilevels : The calcium dynamics at different frequencies has been reported as a risk biomarker in (Bers and Despa, 2006; Pieske et al., 1999). These biomarkers evaluate the diastolic (calcium level at rest) and the systolic (calcium level at peak) values for different frequencies (Figure 2.1.b). •Time constants of the APD adaptation to abrupt changes in cycle length (CL): The adaptation to abrupt changes in CL was proposed as an arrhythmic risk biomarker in (Pueyo et al., 2004). In (Pueyo et al., 2010) the APD adaptation to abrupt changes in CL was tted by two exponentials with time constants τfast and τslow (Figure 2.1.c). These time constants were reported as cellular arrhythmic risk biomarkers. 25 32 Jesús Carro Fernández f2ss =0.67 1 + eVm+35 7 + 0.33 αf2= 300 ·e −(Vm+25)2 170 βf2=31 1 + e25−Vm 10 γf2=16 1 + eVm+30 10 τf2=αf2+βf2+γf2 df2 dt=f2ss −f2 τf2 Additional modications to the L-type calcium current include redenition of the time constant τdof the activation gate d, which was replaced with the description given in (ten Tusscher and Panlov, 2006). The equations for the activation gate dof the L-type calcium current in the CRLP model are: dss =1 1 + e −(Vm+5) 6 αd=1.4 1 + e −35−Vm 13 + 0.25 βd=1.4 1 + eVm+5 5 γd=1 1 + e50−Vm 20 τd=αd·βd+γd dd dt=dss −d τd The calcium-dependent inactivation gates (fCaBjand fCaBsl ) denitions were taken as in the GPB model: dfCaBj dt= 1.7·Caj·(1−fCaBj)−11.9·10−3·fCaBj dfCaBsl dt= 1.7·Casl ·(1−fCaBsl )−11.9·10−3·fCaBsl The relative permeabilities of the L-type calcium channels were adjusted to: pCa = 1.9887 ·10−4cm/s pK= 5.4675 ·10−8cm/s pNa = 3.0375 ·10−9cm/s UniversidadZaragoza An improved human ventricular cell model 33 a) −100 −50 0 50 0 500 1000 Vm(mV ) (ms) τfCRR τf2CRR τfGR exp inact exp rec slow exp rec fast b) −40 −20 0 20 40 60 −10 −5 0 Vm(mV ) ICaL(pA/pF) GR CRR exp c) −140 −120 −100 −80 −60 −40 −20 −10 0 Vm(mV ) IK1(pA/pF) exp GR CRR d) −100 −80 −60 −40 −1 0 1 Vm(mV ) IK1(pA/pF) Figure 3.2. Current characteristics in simulations and experiments: a) ICaL inactivation time constants as function of Vm; b) ICaL current density versus voltage; c) maximum IK1 versus voltage and experimental data. [K+]i= 138 mM and [K+]o= 4 mM; d) Detail of c. Blue dots and markers: experimental data. ¯ ICaj=pCa ·Vm·F rdy ·FoRT ·(Caj·e2·Vm·F oRT −Cao) e2·Vm·F oRT −1 ¯ ICasl =pCa ·Vm·F rdy ·FoRT ·(Casl ·e2·Vm·F oRT −Cao) e2·Vm·F oRT −1 ¯ INaj=pNa ·Vm·F rdy ·FoRT ·(Naj·eVm·F oRT −Nao) eVm·F oRT −1 ¯ INasl =pNa ·Vm·F rdy ·FoRT ·(Nasl ·eVm·F oRT −Nao) eVm·F oRT −1 ¯ IK=pK·Vm·Frdy ·FoRT ·(Ki·eVm·F oRT −Ko) eVm·F oRT −1 The results obtained with the new denition of the ICaL current were compared against experimental data and against the results derived from the original denition of ICaL in the GPB model. These results are presented in Figure 3.2.b. Experimental data is from (Li et al., 1999; Magyar et al., 2000; Pelzmann et al., 1998). 3.3.2 Inward Rectier K Current (IK1) The IK1 current was readjusted based on experimental data from (Fink et al., 2008). In the CRLP model, the IK1 current is formulated as follows: aK1 =4.094 1.0 + e0.1217·(Vm−EK−49.934) bK1 =15.720 ·e0.0674·(Vm−EK−3.257) +e0.0618·(Vm−EK−594.31) 1.0 + e−0.1629·(Vm−EK+14.207) BiomedicalEngineering 34 Jesús Carro Fernández K1ss =aK1 aK1 +bK1 IK1 = 0.5715 ·√[K+]o 5.4·K1ss ·(Vm−EK) where [K+]odenotes the extracellular K+concentration, and EKis the K+reversal potential. Figure 3.2.c shows the modied IK1 compared against the experimental data from (Fink et al., 2008) and IK1 from the GPB model. Note that Figure 3.2.c is a detail of Figure 3.2.c. 3.3.3 Na/K Pump Current (INa,K ) In order to improve the APD adaptation response to abrupt changes in CL, maximal INa,K conductance was reduced by 45% with respect to the GPB model. 3.3.4 [K+]iand GNa A more physiological value for [K+]iof 138 mM was used in our model. Considering this value for [K+]i and in order to get physiological values for the maximal upstroke velocity in tissue, dV /dt, the maximum conductance of the sodium current, GNa, was reduced to 18.86 mS/µF . UniversidadZaragoza Chapter 4 Results 4.1 Conducted Simulations Our proposed CRLP model was tested against the arrhythmic risk biomarkers described in chapter 2 and the results were compared against those obtained using the previously published TP06 and GPB human ventricular cell models. Also, in order to test the model in the same situations described in (Grandi et al., 2010), current blocks were simulated as reported in (Grandi et al., 2010). Finally, the behavior of the model under hyperkalemicconditionswassimulatedin aone-dimensional ber, and conduction velocity (CV),APD and effective refractory period (ERP) were measured. Cells were stimulated with square transmembrane current pulses twice the diastolic threshold and 1-ms duration. Forward Euler integration with a time step ∆t= 0.002 ms was used to integrate the system of differential equations governing the cellular electrical behavior. The Rush and Larsen integration scheme (Rush and Larsen, 1978) was used to integrate the Hodgkin-Huxley type equations for the gating variables of the various time-dependent currents (m,h, and jfor INa;xtof,xtos,ytofand ytosfor Ito;xkr for IKr ;xks for IKs and d,fand f2for ICaL). In the hyperkalemia study, a 4cm long multicellular ber was used. For the simulations, the operator splitting techique was implemented for solving the reaction-diffusion equation that describes the conduction along a cable. The same integration scheme used in isolated cells was used in this case. Implicit Euler witha timestep∆t= 0.002 ms andamesh size∆x= 0.1mm wasused tointegratethe parabolicequation. To obtain a maximum planar CV of 74 cm/s(Taggart et al., 2000), a cellular conduction of 0.0013 cm2/ms was required. Five positions within the ber were selected for further analysis in each simulation. These positions were located at 1.5, 1.75, 2, 2.25 and 2.5 from the stimulation end of the cable. The APD and the CV were computed as their mean values over these positions. 4.2 Results for Control Conditions 4.2.1 Biomarkers of arrhythmic risk Table4.1 showsthe computedbiomarkersforthe modiedCRLP model. Resultsare comparedagainstthose obtained with previous models (TP06 and GPB) and against a variety of physiological data described in (Beuckelmann et al., 1992; Drouin et al., 1995; Franz et al., 1988; Li et al., 1998, 1999; Nash et al., 2006; Pieske et al., 1999, 2002; Schmidt et al., 1998). These results indicate that the GPB model performs better than the TP06 model with respect to the following biomarkers: triangulation, systolic and diastolic [Ca2+]iat 0.5 Hz and rate dependence of steadystate [Na+]iand [Ca2+]i. However, according to other analyzed biomarkers, the GPB model renders worse results than the TP06 model. Specically, in (Franz et al., 1988) two phases of APD rate adaptation were 35 36 Jesús Carro Fernández Biomarker TP06 GPB CRLP Physiolo. APD90 301.2 285.0 305.6 271-366 Triangulation 28.4 51.5 78 44-112 Sys [Ca2+]i1 Hz 0.886 0.383 0.602 1.59-2.01 Sys [Ca2+]i0.5 Hz 0.199 0.345 0.523 0.71-1.68 Dia [Ca2+]i1 HZ 0.104 0.089 0.097 0.20-0.33 Dia [Ca2+]i0.5 Hz 0.068 0.085 0.091 0.14-0.32 SmaxS1S2 1.3 0.2 1.0 0.79-4.25 SmaxDYN 1.0 1.1 0.9 --- τfast 13.3 --- 10.2 --- τslow 124.8 56.3 105.6 70-110 Max. sys. [Ca2+]i1157 178 147 130-170 Max. sys. [Na+]i217 132 134 145 Table 4.1. Biomarkers of arrhythmic risk for the TP06, GPB and CRLP human models. Green color indicates within physiological range. Blue color indicates out of physiological range but better than previous models. Red color indicates out of physiological range. reported: a rst phase with a short time constant and a second phase with a longer time constant. The APD rate response in the GPB model does not present this rst phase. Also, the S1S2 restitution curves of the TP06 model (Figure 4.1.c) are in better agreement with experimental data. For the CRLP model, most of the twelve analyzedbiomarkers are within physiological range or otherwise closer to the physiological range as compared to GPB and TP06 models, except for systolic and diastolic [Ca2+]ilevels at 1 Hz, where TP06 outperforms CRLP. 4.2.2 Blocking Currents We have also studied the behavior of the model under total and partial block of potassium currents, comparing the results against those reported in (Jost et al., 2008) from experiments in human ventricle: •IKs:To simulate the effect of 1µM of HMR-1556, IKs was completely blocked. •IKr :To simulate the effect of 50 nM of dofetilide, IKr was completely blocked. •IK1:To simulate the effect of 10 µM of BaCl2,IK1 was reduced by 50%. Table 4.2 summarizes the effect of current block on the different models. Results are given as percentage of variation of the APD value in control conditions. In those cases where more than one current is blocked, the variation refers to the APD obtained when the rst current is blocked. Experimental values have been taken from (Jost et al., 2008). In Figure 4.2 the efects of the different blocks for the original and the modied models are represented. 4.3 Results for Hyperkalemic Conditions Forthesimulationofhyperkalemiaa1-dimensionalber wasused,asdescribed insection4.1. Hyperkalemic conditions were simulated by increasing the extracellular K+concentration from 4.0 mM to 15 mM (depending on the degree of simulated hyperkalemia). The stimulation protocol consisted of a train of ten basic stimulations (S1) with a cyclic length of 1,000 ms, followed by an extrasystole stimulus (S2) delivered at different times taken in 0.1 millisecond steps in order to determine ERP. CV and APD were measured after the last S1 stimulus. UniversidadZaragoza An improved human ventricular cell model 37 a) 0 200 400 600 −80 −60 −40 −20 0 20 40 Time (ms) Vm(mV) Steady-State AP Shape b) 0 5 10 15 20 260 280 300 APD Rate Adaptation Time (min) APD90 (ms) c) 0 200 400 600 800 200 250 300 S1S2 DI (ms) APD90 (ms) GR TTP06 CRR d) 0 200 400 600 800 200 250 300 Dynamic DI (ms) APD90 (ms) e) 0.5 1 1.5 2 2.5 3 200 400 600 800 1000 1200 [Ca2+ ]iRate Dependence Frequency (Hz) Systolic [Ca2+]i(%) f) 0.5 1 1.5 2 2.5 3 100 150 200 [Na+]iRate Dependence Frequency (Hz) Systolic [Na+]i(%) Figure 4.1. Biomarkers of arrhythmic risk from TP06, GPB and CRLP models: a) AP for CL = 1000 ms; b) APD rate adaptation to abrupt changes in CL (From 1000 ms to 600 ms); c) S1S2 restitution protocol curve; d) Dynamic restitution protocol curve; e) rate dependence of steady-state [Na+]i; f) Rate dependence of steady-state [Ca2+]i Current Ref. TT GPB CRLP Exper. IKr 0% Control 74.9 0.7 0.6 <2.8 IKr 0% Control 15.9 18.6 14.1 44±4 IKs 0% IKr B. --- 1.4 0.8 9 IKr 0% IK1 50% Control 3.6 11.5 14.6 4.8±1.5 IKr 0% IKr B. 8.9 15.6 21.1 33 IK1 50% Table 4.2. Percentages of variation in the APD caused by blocking different currents. Green color indicates within physiological range. Blue color indicates out of physiological range but better than previous models. Red color indicate out of physiological range. BiomedicalEngineering 38 Jesús Carro Fernández GPB Model a) 0 200 400 600 −80 −60 −40 −20 0 20 40 Time (ms) Vm(mV) Control IKr IKs IKr &IKs b) 0 200 400 600 −80 −60 −40 −20 0 20 40 Time (ms) Vm(mV) Control IKr IK1 IKr &IK1 CRLP Model c) 0 200 400 600 −80 −60 −40 −20 0 20 40 Time (ms) Vm(mV) Control IKr IKs IKr &IKs d) 0 200 400 600 −80 −60 −40 −20 0 20 40 Time (ms) Vm(mV) Control IKr IK1 IKr &IK1 Figure 4.2. Blocking currents for the GPB and CRLP models: a) Effect of completely blocking IKr and IKs currents in the GPB model, b) Effect of completely blocking IKr and partialy blocking IK1 in the GPB model, b) Effect of completely blocking IKr and IKs currents in the CRLP model, b) Effect of completely blocking IKr and partially blocking IK1 in the CRLP model. UniversidadZaragoza An improved human ventricular cell model 39 a) 4 6 8 10 20 40 60 80 [K+]0(mM) C.V. (cm/s) b) 6000 6500 7000 7500 8000 8500 −60 −40 −20 0 time (ms) V (mV) c) 4 6 8 10 0 500 1000 [K+]o(mM) (ms) ERP APD90 Figure 4.3. Behavior of the model under hyperkalemic conditions in 1-D ber: a) conduction velocity versus [K+]o; b) example of alternans for [K+]o= 10.0mM, each color represents a different position in the ber; c) dependence of ERP and APD with the [K+]o. The basic stimuli consisted of rectangular pulses of 3 ms duration and 1.5 times the diastolic threshold. This threshold was determined as follows. For a given [K+]ovalue, the model was rst stabilized until when the product of the gates h·jreached 99% of the steady-state value hss ·jss. Once the model had reached steady-state conditions, an stimulation current was applied at the left end of the cable. Diastolic threshold wasdenedas theminimumstimulationcurrentrequiredin orderforanactionpotentialtopropagatealong the cable. 4.3.1 Conduction Velocity Figure 4.3.a shows CV obtained for different values of [K+]o. In our proposed CRLP model, as in TP06, there is no supernormal conduction (Shaw and Rudy, 1997) (Kagiyama et al., 1982). Above 10.4 mM there is no conduction and over 9.8 mM alternans appeared (Figure 4.3.b). 4.3.2 Action Potential Duration and Effective Refractory Period Figure 4.3.c shows APD and ERP measurements for different values of [K+]o. APD decreases when the extracellular K+concentration increases. For high [K+]values, the APD curve is divided into two branches, which is explained by the ocurrence of alternans at those high concentrations. The ERP has three phases. In the rst one, ERP decreases with APD. After this, ERP begins to increase (the excitability of the cells decreases due to the increase in [K+]o). Finally, the last phase is caused by the ocurrence of alternans. Depending on whether the beat that is measured is an odd or an even beat, ERP takes one or another value. BiomedicalEngineering Chapter 5 Discussion and conclusions In this study a new human ventricular cell model has been developed, which is suitable for investigation of cardiac arrhythmias. Taking the recently proposed GPB model as a starting point, several modications have been introduced by accounting for recent experimental measurements of potassium currents, and by reformulating the L-type calcium current so as to introduce fast and slow voltage-dependent inactivation. All the modications have been made with the aim of maintaining the advantages that the GPB model presents over previous models, such as the TP06 model, while improving the performance on other electrophysiological aspects. The novelty of this study is that the development of our proposed CRLP model is not only based on new data from experiments of ionic currents. More importantly, the development is also based on the results of several arrhythmic risk biomarkers reported in the literature. For the proposed CRLP model, the introduction of fast and slow ICaL inactivation has led to a more physiological rate adaptation responseas compared to the original GPB model. The other modications have rendered most biomarkers to be within the physiological range or otherwise closer to the physiological range as compared to the GPB and TP06 models, except for systolic and diastolic [Ca2+]ilevels at 1 Hz where TP06 outperforms CRLP. These modications, however, have altered the behavior of the model against current blocks minimally, thus providing results in very good agreement with available experimental data. The results of the analysis conducted in this study points out the importance that calcium dynamics have over different biomarkers. This suggests the need for continuing with the development of more reliable calcium dynamic models that allow improving the performance of whole cell AP models. The hyperkalemia simulations showed that APD progressively decreased with the level of hyperkalemia, while ERP increased after a threshold in the extracellular [K+]was reached ([K+]o= 6 mM). Conduction velocity decreased with hyperkalemia and the conduction was blocked above [K+]o= 10.4mM. Additionally, alternans appeared in the APD above [K+]o= 9.8mM. These results suggest that the longer ERP values and the conduction blocking above [K+]o= 10.4mM found in the central zone of acutely ischemic tissue as compared to the normal zone could create areas of block that could set the stage for reentrant arrhythmias. Future develpments of the model will include a complete ischemic model with ATP sensitive potassium current (IKATP ) and simulations of acidosis. The simulations with the model will be made, not only in a 1-D ber but in 2-D tissue and in a realistic human ventricular geometry. 41 Clinical Practice Clinical Practice The clinical practice of the study was conducted at Hospital Clínico de Valencia. It was divided into two parts. The rst part was carried out in collaboration with the group directed by Dr. Ricardo Ruiz Granell in the Cardiostimulation Unit. The second part was conducted in the research group led by Dr. Francisco Javier Chorro about the investigation of different pathologies using rabbit hearts. In the rst part of the clinical practice with the group of Dr. Ricardo Ruiz, several interventions were witnessed. Some cardioversions were performed by team doctors, where different types of arrhythmias were terminated by delivering a synchronized shock. Likewise, different checks and follow-ups of pacemakers and debrillators were witnessed. Also, an intervention in a patient with history of chronic ischemic arrhythmia took place during the days of the internship. In that intervention an ablation of the part of the heart causing reentrant circuits was performed. For this, rst of all, a 3D map of the inside of the heart and its electrical activity was measured. An arrhythmia was induced with an extra-stimulus test to check its focus and, nally, all the identied zone was burned in order to prevent the spread of late fronts. During the days of the clinical practice, different conversations on the applicability of the investigation with computational models took place between the student and the doctors. In this conversations the doctors proposed the application of our model to the investigation of arrhythmias caused by genetic diseases like Brugada syndrome, where some channels present a pathological behavior and are common in the clinical routine. Other pathologies like arrhythmias caused by acute ichemia were very rare to nd in the clinic because usually the patient does not survive to reach the hospital and in case it happens, the procedure is clear, to deliver an electroshock. During the second part of the clinical practice with the group led by Dr. Francisco Javier Chorro, some experiments in rabbit hearts were witnessed. This group is currently working in the study of the differences between the heart response in rabbits trained and not trained. Likewise, they are working on comparing the electrophysiology of the heart in control conditions and when it is subjected to stretch. In this studies they are using an increasing stimulation frequency to generate arrhythmias and the extra-stimulus test to study CV and ERP. With a sensor they perform an electrical mapping in the surface of the heart and they process the measures to create activation maps. For all these experiments they use a Langendorff system that allows them to mantain the heart alive after having been extracted from the body. With this clinical practice it was posible to see a real application of the same procedures used in the simulations conducted in this study. 51 Appendices Appendix A Human Ventricular Cell Model A.1 Model Parameters A.1.1 Physical Constants Name Value Units R8314 J/(mol ·K) Frdy 96485 C/mol T310 K FoRT F rdy/(R·T)mV −1 Cmem 1.381 ·10−10 F A.1.2 Enviromental Parameters Name Value Units Lengthcell 100 ·10−5dm Radiuscell 10.25 ·10−5dm Vcell π·Radius2 cell ·Lengthcell L Vmyo 0.65 ·Vcell L Vsr 0.035 ·Vcell L Vsl 0.02 ·Vcell L Vjunc 5.39 ·10−4·Vcell L JCajunc,sl 8.2413 ·10−13 L/ms JCasl,myo 3.7243 ·10−12 L/ms JNajunc,sl 1.8313 ·10−14 L/ms JNasl,myo 1.6386 ·10−12 L/ms A.1.3 Fractional Currents Name Value Units Fjunc 0.11 [−] Fsl 1−Fjunc [−] FjuncCaL0.9 [−] FslCaL1−FjuncCaL[−] 55 56 Jesús Carro Fernández A.1.4 Ion Concentrations Name Value Units Ki138 mM Ko5.4mM Cli15 mM Clo150 mM Cao1.8mM Mgi1mM Nao140 mM A.1.5 Sodium Transport Name Value Units GNa 18.86 mS/µF GNaB 0.597 ·10−3mS/µF ¯ INaK 0.99 pA/pF KmKo 1.5mM KmNaip 11 mM A.1.6 Potassium Currents Name Value Units pNaK 1.833 ·10−2[−] GK1 5.7153 ·10−1mS/µF GKp 2·10−3mS/µF GKr 3.5·10−2mS/µF GKsjunc 3.5·10−3mS/µF GKssl 3.5·10−3mS/µF Gtof,EPI 1.144 ·10−1mS/µF Gtos,EPI 1.56 ·10−2mS/µF Gtof,ENDO 1.404 ·10−3mS/µF Gtos,ENDO 3.7596 ·10−2mS/µF A.1.7 Chlorine currents Name Value Units GClCa 0.054813 mS/µF GClB 9·10−3mS/µF KdClCa 100 ·10−3mM UniversidadZaragoza An improved human ventricular cell model 57 A.1.8 Calcium Transport Name Value Units pCa 1.9887 ·10−4L/(F·ms) pK5.4675 ·10−8L/(F·ms) pNa 3.0375 ·10−9L/(F·ms) ¯ INCX 4.5pA/pF KmCai3.59 ·10−3mM KmCao1.3mM KmNai12.29 mM KmNao87.5mM ksat 0.32 [−] nu 0.27 [−] Kdact 0.15 ·10−3mM ¯ IPMCA 0.0673 pA/pF KmpCa 0.5·10−3mM GCaB 5.513 ·10−4mS/µF A.1.9 SR Calcium Fluxes Name Value Units VmaxSRCaP 5.3114 ·10−3mM/ms ec50SR 0.45 mM hillSRCaP 1.787 [−] kiCa 0.5mM−1·mS−1 kim0.005 mS−1 koCa 10 mM−2·ms−1 kom0.06 mS−1 ks 25 mS−1 Kmf0.246 ·10−3mM Kmr1.7mM MaxSR 15 [−] MinSR 1 [−] BiomedicalEngineering 64 Jesús Carro Fernández Currents ¯ ICaj=pCa ·Vm·F rdy ·FoRT ·(Caj·e2·Vm·F oRT −Cao) e2·Vm·F oRT −1 ¯ ICasl =pCa ·Vm·F rdy ·FoRT ·(Casl ·e2·Vm·F oRT −Cao) e2·Vm·F oRT −1 ¯ INaj=pNa ·Vm·F rdy ·FoRT ·(Naj·eVm·F oRT −Nao) eVm·F oRT −1 ¯ INasl =pNa ·Vm·F rdy ·FoRT ·(Nasl ·eVm·F oRT −Nao) eVm·F oRT −1 ¯ IK=pK·Vm·Frdy ·FoRT ·(Ki·eVm·F oRT −Ko) eVm·F oRT −1 ICajunc =FjuncCaL·¯ ICaj·d·f·f2·(1−fCaBj) ICasl =FslCaL·¯ ICasl ·d·f·f2·(1−fCaBsl ) ICaNajunc =FjuncCaL·¯ INaj·d·f·f2·(1−fCaBj) ICaNasl =FslCaL·¯ INasl ·d·f·f2·(1−fCaBsl ) ICaK=¯ IK·d·f·f2·(FjuncCaL·(1−fCaBj)+FslCaL·(1−fCaBsl )) ICa =ICajunc +ICasl ICaNa =ICaNajunc +ICaNasl ICaL =ICa +ICaK+ICaNa UniversidadZaragoza An improved human ventricular cell model 65 A.2.13 Incx: Na-Ca Exchanger Current Kajunc =1 1 + (Kdact Caj)2 Kasl =1 1 + (Kdact Casl )2 s1junc =enu·Vm·F oRT ·Na3 j·Cao s1sl =enu·Vm·F oRT ·Na3 sl ·Cao s2junc =e(nu−1)·Vm·F oRT ·Na3 o·Caj s2sl =e(nu−1)·Vm·F oRT ·Na3 o·Casl s3junc =KmCai·Na3 o·(1 + (Naj KmNai)3)+Km3 Nao·Caj·(1 + Caj KmCai) +KmCao·Na3 j+Na3 j·Cao+Na3 o·Caj s3sl =KmCai·Na3 o·(1 + (Nasl KmNai)3)+Km3 Nao·Casl ·(1 + Casl KmCai) +KmCao·Na3 sl +Na3 sl ·Cao+Na3 o·Casl Incxjunc =Fjunc ·¯ INCX ·Kajunc ·(s1junc −s2junc) s3junc ·(1 + ksat ·e(nu−1)·Vm·F oRT ) Incxsl =Fsl ·¯ INCX ·Kasl ·(s1sl −s2sl) s3sl ·(1 + ksat ·e(nu−1)·Vm·F oRT ) Incx =Incxjunc +Incxsl A.2.14 IpCa: Sarcolemmal Calcium Pump Current IpCajunc =Fjunc ·¯ IPMCA ·Ca1.6 j Km1.6 pCa +Ca1.6 j IpCasl =Fsl ·¯ IPMCA ·Ca1.6 sl Km1.6 pCa +Ca1.6 sl IpCa =IpCajunc +IpCasl A.2.15 ICaBk: Background Calcium Current ICaBkjunc =Fjunc ·GCaB ·(Vm−ECajunc ) ICaBksl =Fsl ·GCaB ·(Vm−ECasl ) ICaBk =ICaBkjunc +ICaBksl BiomedicalEngineering 66 Jesús Carro Fernández A.2.16 SR Fluxes: CalciumRelease, SRCalciumPump, SR Calcium Leak kCaSR =MaxSR −MaxSR −MinSR 1 + (ec50SR CaSR )2.5 koSRCa =koCa kCaSR kiSRCa =kiCa ·kCaSR RI = 1 −RyRr −RyRo −RyRi dRyRr dt=kim·RI −kiSRCa ·Caj·RyRr −(koSRCa ·Ca2 j·RyRr −kom·RyRo) dRyRo dt=koSRCa ·Ca2 j·RyRr −kom·RyRo −(kiSRCa ·Caj·RyRo −kim·RyRi) dRyRi dt=kiSRCa ·Caj·RyRo −kim·RyRi −(kom·RyRi −koSRCa ·Ca2 j·RI) JSRCarel =ks ·RyRo ·(CaSR −Caj) JserCa = VmaxSRCaP ·((Cai Kmf)hillSRCaP −(CaSR Kmr)hillSRCaP ) 1 + (Cai Kmf)hillSRCaP +(CaSR Kmr)hillSRCaP JSRleak = 5.348 ·10−6·(CaSR −Caj) A.2.17 Ion Homeostasis Sodium Buffers dNaBj dt=konNa ·Naj·(BmaxNaj−NaBj)−koffNa ·NaBj dNaBsl dt=konNa ·Nasl ·(BmaxNasl −NaBsl )−koffNa ·NaBsl UniversidadZaragoza An improved human ventricular cell model 67 Cytosolic Calcium Buffers dTnCl dt=konT nCl·Cai·(BmaxT nClow −T nCl)−koffT nCl·T nCl dTnChc dt=konT nChCa ·Cai·(BmaxT nChigh −TnChc−TnChm)−koffT nChCa ·TnChc dTnChm dt=konT nChMg ·Mgi·(BmaxT nChigh −TnChc−T nChm)−koffT nChMg ·TnChm dCaM dt=konCaM ·Cai·(BmaxCaM −CaM)−koffCaM ·CaM dMyoc dt=konmyoCa ·Cai·(Bmaxmyosin −Myoc−Myom)−koffmyoCa ·Myoc dMyom dt=konmyoMg ·Mgi·(Bmaxmyosin −Myoc−Myom)−koffmyoMg ·Myom dSRB dt=konSR ·Cai·(BmaxSR −SRB)−koffSR ·SRB JCaBcytosol =dT nCl dt+dTnChc dt+dTnChm dt+dCaM dt+dMyoc dt+dMyom dt+dSRB dt Junctional and SL Calcium Buffers dSLLj dt=konsll·Caj·(BmaxSLlowj −SLLj)−koffsll·SLLj dSLLsl dt=konsll·Casl ·(BmaxSLlowsl −SLLsl)−koffsll·SLLsl dSLHj dt=konslh·Caj·(BmaxSLhighj −SLHj)−koffslh·SLHj dSLHsl dt=konslh·Casl ·(BmaxSLhighsl −SLHsl)−koffslh·SLHsl JCaBjunction =dSLLj dt+dSLHj dt JCaBsk =dSLLsl dt+dSLHsl dt Sodium Concentrations dNaj dt= −INatotjunc ·Cmem Vjunc ·F rdy +JNajunc,sl Vjunc ·(Nasl −Naj)− dNaBj dt dNasl dt=−INatotsl ·Cmem Vsl ·Frdy +JNajunc,sl Vsl ·(Naj−Nasl) + JNasl,myo Vsl ·(Nai−Nasl)−dNaBsl dt dNai dt=JNasl,myo Vmyo ·(Nasl −Nai) BiomedicalEngineering 68 Jesús Carro Fernández Calcium Concentrations dCsqnb dt=koncsqn ·CaSR ·(Bmaxcsqn −Csqnb)−koffcsqn ·Csqnb dCaj dt= −ICatotjunc ·Cmem Vjunc ·2·F rdy +JCajunc,sl Vjunc ·(Casl −Caj)−JCaBjunction +JSRCarel ·Vsr Vjunc +JSRleak ·Vmyo Vjunc dCasl dt=−ICatotsl ·Cmem Vsl ·2·Frdy +JCajunc,sl Vsl ·(Caj−Casl) + JCasl,myo Vsl ·(Cai−Casl)−JCaBsk dCai dt=−JserCa ·Vsr Vmyo −JCaBcytosol +JCasl,myo Vmyo ·(Casl −Cai) dCaSR dt=JserCa −(JSRleak ·Vmyo Vsr +JSRCarel )−dCsqnb dt A.2.18 Membrane Potential ICatotjunc =ICajunc +ICaBkjunc +IpCajunc −2·Incxjunc ICatotsl =ICasl +ICaBksl +IpCasl −2·Incxsl ICatot =ICatotjunc +ICatotsl IKtot =Ito +IKr +IKs +IK1 −2·INa,K +ICaK+IKp INatotjunc =INajunc +INaBkjunc + 3 ·Incxjunc + 3 ·INa,Kjunc +ICaNajunc INatotsl =INasl +INaBksl + 3 ·Incxsl + 3 ·INa,Ksl +ICaNasl INatot =INatotjunc +INatotsl ICltot =IClCa +IClBk Iion =INatot +ICltot +ICatot +IKtot dVm dt=−(Iion −Istim) UniversidadZaragoza An improved human ventricular cell model 69 A.3 Initial Conditions Name EPI ENDO Units CaSR 6.093596 ·10−16.138856 ·10−1mM Cai9.658067 ·10−59.719632 ·10−5mM Caj2.038197 ·10−42.048633 ·10−4mM Casl 1.184305 ·10−41.188246 ·10−4mM Csqnb1.258048 1.262853 mM CaM 3.267494 ·10−43.288063 ·10−4mM Myoc2.520383 ·10−32.522168 ·10−3mM Myom1.369529 ·10−11.369514 ·10−1mM SRB 2.373753 ·10−32.38683 ·10−3mM TnChc1.225914 ·10−11.225802 ·10−1mM TnChm8.12201 ·10−38.128604 ·10−3mM TnCl9.757237 ·10−39.811535 ·10−3mM h7.13497 ·10−17.126555 ·10−1[−] j7.128671 ·10−17.119893 ·10−1[−] m2.163678 ·10−32.176608 ·10−3[−] SLHj8.053908 ·10−28.078504 ·10−2mM SLHsl 1.235381 ·10−11.238366 ·10−1mM SLLj8.563314 ·10−38.606485 ·10−3mM SLLsl 1.097424 ·10−21.101044 ·10−2mM d1.871177 ·10−61.879996 ·10−6[−] f9.804391 ·10−19.789409 ·10−1[−] f29.99401 ·10−19.993986 ·10−1[−] fCaBj2.847118 ·10−22.861794 ·10−2[−] fCaBsl 1.692189 ·10−21.69833 ·10−2[−] xKr 1.516232 ·10−21.896559 ·10−2[−] RyRi 1.411382 ·10−71.43831 ·10−7[−] RyRo 1.126209 ·10−61.149876 ·10−6[−] RyRr 8.886338 ·10−18.888214 ·10−1[−] xKs 3.549354 ·10−33.55636 ·10−3[−] NaBj3.796195 3.785209 mM NaBsl 8.283308 ·10−18.259271 ·10−1mM Nai1.007825 ·1011.001989 ·101mM Naj1.007931 ·1011.00211 ·101mM Nasl 1.00781 ·1011.001974 ·101mM xtof3.584625 ·10−43.592405 ·10−4[−] xtos3.584727 ·10−43.592503 ·10−4[−] ytof9.999976 ·10−19.999976 ·10−1[−] ytos8.087629 ·10−18.161309 ·10−1[−] Vm−84.13368 −84.10546 mV BiomedicalEngineering Appendix B International Conferences 71 72 Jesús Carro Fernández B.1 ComputinginCardiology 2010, Belfast(United Kingdom) Analysis and Improvement of a Human Ventricular Cell Model for Investigation of Cardiac Arrhythmias J Carro1,2, JF Rodr´ ıguez1, P Laguna1,2, E Pueyo1,2 1Instituto de Investigaci´ on en Ingenier´ ıa de Arag´ on (I3A), Universidad de Zaragoza, Spain 2CIBER de Bioingenier´ ıa, Biomateriales y Nanomedicina (CIBER-BBN), Spain Abstract The use of experiments for studying cardiac arrhythmias, the effect of drugs, or pathologies on cardiac electrophysiology is very limited. This has made mathematical modeling and simulation of heart’s electrical activity a fundamental tool to understand cardiac behavior. In this study several modifications were introduced to a recently proposed human ventricular cell model. Four stimulation protocols were applied to the original and improved models of isolated cell, and a number of cellular arrhythmic risk biomarkers were computed: steady-state action potential (AP) and [Ca2+]transient properties, AP duration (APD) restitution curves, APD adaptation to abrupt changes in heart rate, and intracellular [Ca2+]and [Na+] rate dependence. Our modifications led to: a) further improved AP triangulation (78.1ms); b) APD rate adaptation curves characterized by fast and slow time constants within physiological ranges (10.1sand 105.9s); c) maximum S1S2 restitution slope in accordance with experimental data (SS1S2 = 1.0). 1. Introduction Ventricular arrhythmias can have their origin in diseases, genetic disorders, drug cardiotoxicity, and a number of other causes. Research in this field has identified a number of potential biomarkers of arrhythmic risk related to cellular electrophysiological properties. However, performing experimental and clinical studies involving human hearts to quantify these biomarkers is very difficult. On the other hand, animal hearts used for experimental studies may differ significantly from human hearts. In addition, ventricular cardiac arrhythmias are three-dimensional phenomena whereas experimental observations are still largely constrained to surface recordings. Under this circumstances, mathematical models of myocardial cells and computer simulations of cardiac activity in the human heart can overcome some of these problems. Recently, a new model has been proposed by Grandi et al. [1] (GR) that improves the response to frequency changes and has a better performance against current blocks with respect to the model proposed by ten Tusscher and Panfilov [2] (TTP06). This article introduces several modifications to the GR model by accounting for recent experimental measurements of potassium currents [3] and the introduction of fast and slow inactivation of L-type calcium current [2]. The performance of the proposed model, herein denoted by CRR, has been tested against the GR and TTP06 models by applying four stimulation protocols and computing twelve cellular arrhythmic risk biomarkers following the methodology proposed in [4]. The results show that the introduced modifications have brought most biomarkers into the physiological range, and considerably improved others with respect to the GR and TTP06 models of AP. 2. Methods 2.1. Biomarkers of arrhythmic risk Four stimulation protocols have been applied and twelve cellular biomarkers of arrhythmic risk have been computed as in [4]: Steady-state AP and intracellular [Ca2+]concentration properties. Systolic and diastolic [Ca2+]ilevels under steady-state pacing at frequencies of 0.5 and 1 Hz and APD and AP Triangulation at 1 Hz have been calculated. These biomarkers have been proposed as arrhythmic risk biomarkers in the literature [5–7]. APD restitution (APDR) curves. APDR curves have been obtained using the S1S2 and the dynamic restitution protocols as in [4]. In both cases, the maximum slope, suggested as a risk marker in the literature [8,9], has been computed. APD rate adaptation to abrupt changes in cycle length (CL). APD rate adaptation dynamics have been proposed as a clinical marker in [10]. Dynamics have been fitted to two exponentials with time constants τfast and τslow [11]. Rate dependence of steady-state [Na+]iand [Ca2+]i. UniversidadZaragoza An improved human ventricular cell model 73 a) −100 −50 0 50 0 500 1000 Vm(mV ) (ms) τfCRR τf2CRR τfGR exp inact exp rec slow exp rec fast b) −40 −20 0 20 40 60 −10 −5 0 Vm(mV ) ICaL(pA/pF) GR CRR exp Figure 1. ICaLcharacteristics in simulations and experiments: a) Inactivation time constants; b) current density versus voltage. Markers: experimental data. The importance of [Na+]iand [Ca2+]idynamics in arrhythmogenesis has been reported [12, 13]. The systolic values of both concentrations have been measured at different frequencies (0.25, 0.5, 1, 1.5, 2, 2.5 and 3 Hz) and normalized to the level of the minimum frequency. The maximum value has been computed and used as biomarker. 2.2. Modifications of the model Two ionic currents have been reformulated and some model parameters have been redefined in the GR model: L-Type Calcium Current (ICaL). The voltage-dependent inactivation gate fhas been replaced with the product of a fast, f, and a slow, f2, inactivation gates as in [2]. Based on available experimental data, the value of the opening rate αf2, associated with f2, has been modified to: αf2= 300 ·e− (Vm+25)2 170 , where Vmdenotes the transmembrane potential. Figure 1a depicts τfand τf2along with experimental data [2].The time constant (τd) of the activation gate dhas been replaced with the description given in [2]. The relative permeabilities of the L-type calcium channels have been adjusted to: pCa = 1.458 ·10−4cm/s, pK= 7.290 ·10−8cm/s,pNa = 4.050 ·10−9cm/s. This new definition of ICaLis compared against experimental data and the original definition of the GR model in Figure 1b. Experimental data is from [14–16]. Inward Rectifier K Current (IK1). This current has been re-adjusted based on experimental data from [3]. The readjusted current is formulated as follows: VEk=Vm−Ek aK1 =4.094 1.0 + e0.1217·(VEk −49.934) bK1 =15.720 ·e0.0674·(VEk −3.257)+e0.0618·(VEk −594.31) 1.0 + e−0.1629·(VEk+14.207) K1ss =aK1 aK1 +bK1 a) −140 −120 −100 −80 −60 −40 −20 −10 0 Vm(mV ) IK1(pA/pF) exp GR CRR b) −100 −80 −60 −40 −1 0 1 Vm(mV ) IK1(pA/pF) Figure 2. Maximum IK1 versus voltage and experimental data. [K+]i= 138 mM and [K+]o= 4 mM. Blue dots: experimental data. IK1 = 0.5715 ·r[K+]o 5.4·K1ss ·VEk where [K+]odenotes the extracellular K+concentration, and EKis the K+reversal potential. Figure 2 shows the modified IK1 compared against the experimental data from [3] and IK1 from the GR model. Note that Figure 2b is a zoom from Figure 2a. INa,K: Na/K Pump Current. Maximal INa,K conductance has been reduced by 45%. [K+]iand GNa. A more physiological value of [K+]i (138 mM) has been used in the model. In this regard, in order to get physiological values for the maximal upstroke velocity, dV/dt, the maximum conductance of the sodium current, GNa, has been reduced to 18.86 mS/µF . 2.3. Blocking Currents We have also studied the behavior of the model under total and partial block of potassium currents, comparing the results against those reported in [17]: IKs: To simulate the effect of 1µM of HMR-1556, IKs has been completely blocked. IKr: To simulate the effect of 50 nM of dofetilide, IKr has been completely blocked. IK1: To simulate the effect of 10 µM of BaCl2,IK1 has been blocked at 50%. 2.4. Numerical implementation and simulation Model differential equations were implemented in Fortran. Cells were stimulated with square transmembrane current pulses twice the diastolic threshold at CL = 1000 ms and 1-ms duration. Forward Euler integration with a time step ∆t= 0.002 ms was used to integrate the system of differential equations governing the cellular electrical behavior. The Rush and Larsen integration scheme was used to integrate the Hodgkin-Huxley type equations for the gating variables of the various time-dependent currents. BiomedicalEngineering Communications Technology Group Group of Structural Mechanics and Material Modeling