Journal of Pharmacological and Toxicological Methods 81 (2016) 233–239

Contents lists available at ScienceDirect

Journal of Pharmacological and Toxicological Methods
journal homepage: www.elsevier.com/locate/jpharmtox

A temperature-dependent in silico model of the humanether-à-go-go-related (hERG) gene channel
Zhihua Li ⁎, Sara Dutta, Jiansong Sheng, Phu N. Tran, Wendy Wu, Thomas Colatsky
Division of Applied Regulatory Science, Ofﬁce of Clinical Pharmacology, Ofﬁce of Translational Sciences, Center for Drug Evaluation and Research, Food and Drug Administration, United Statesa r t i c l e

i n f oArticle history:
Received 29 January 2016
Received in revised form 20 April 2016
Accepted 9 May 2016
Available online 11 May 2016
Keywords:
CiPA
hERG
Kv11.1
Markov model
Methods
Temperaturea b s t r a c t
Introduction: Current regulatory guidelines for assessing the risk of QT prolongation include in vitro assays
assessing drug effects on the human ether-à-go-go-related (hERG; also known as Kv11.1) channel expressed incell lines. These assays are typically conducted at room temperature to promote the ease and stability of recording hERG currents. However, the new Comprehensive in vitro Proarrhythmia Assay (CiPA) paradigm proposes touse an in silico model of the human ventricular myocyte to assess risk, requiring as input hERG channel pharmacology data obtained at physiological temperatures. To accommodate current industry safety pharmacology
practices for measuring hERG channel activity, an in silico model of hERG channel that allows for the extrapolation of hERG assay data across different temperatures is desired. Because temperature may have an effect on bothchannel gating and drug binding rate, such models may need to have two components: a base model dealing
with temperature-dependent gating changes without drug, and a pharmacodynamic component simulating
temperature-dependent drug binding kinetics. As a ﬁrst step, a base mode that can capture temperature effects
on hERG channel gating without drug is needed.
Methods and results: To meet this need for a temperature-dependent base model, a Markov model of the hERG
channel with state transition rates explicitly dependent on temperature was developed and calibrated using
data from a variety of published experiments conducted over a range of temperatures. The model was able to reproduce observed temperature-dependent changes in key channel gating properties and also to predict the results obtained in independent sets of new experiments.
Discussion: This new temperature-sensitive model of hERG gating represents an attempt to improve thepredictivity of safety pharmacology testing by enabling the translation of room temperature hERG assay data
to more physiological conditions. With further development, this model can be incorporated into the CiPA paradigm and also be used as a tool for developing insights into the thermodynamics of hERG channel gating mechanisms and the temperature-dependence of hERG channel block by drugs.
Published by Elsevier Inc.

1. Introduction
The human ether-à-go-go-related gene (hERG, also known as
Kv11.1) encodes the alpha pore-forming subunit of the rapidly activating component of the cardiac delayed rectiﬁer current (IKr)(Sanguinetti, Jiang, Curran, & Keating, 1995). IKr is critical for cardiac action potential (AP) repolarization(Jurkiewicz & Sanguinetti, 1993) andsuppressing IKr by drugs can be associated with electrocardiographic
QT prolongation and an elevated risk of serious proarrhythmia, speciﬁcally Torsade de Pointes (TdP) (Malik & Camm, 2001). Because of this,current regulatory guidelines consider hERG channel block a potential
safety concern and recommend the evaluation of drug effects on hERG
channel as part of the preclinical development program (“ICH S7B:
⁎ Corresponding author at: WO Bldg 64 Rm 2084, 10903 New Hampshire Ave, Silver
Spring, MD 20993, United States.
E-mail address: Zhihua.li@fda.hhs.gov (Z. Li).

http://dx.doi.org/10.1016/j.vascn.2016.05.005
1056-8719/Published by Elsevier Inc.

Note for Guidance on the Nonclinical Evaluation of the potential for Delayed Ventricular Repolarization (QT Interval Prolongation) by Human
Pharmaceuticals”). The evaluation of hERG channel block is typicallycarried out early in drug discovery to determine the concentration atwhich hERG channel activity is reduced by 50% (IC50) (Redfern et al.,2003). The methodology used to conduct hERG assays varies widely
across laboratories, as do the results obtained, with voltage protocol
and temperature being among the top sources of variability (Kirsch etal., 2004).
Some effort has been made to address this issue of methodologicaldifferences. For example, computational models have been developed
to extract kinetic information from dynamic drug-hERG channel interaction data that are intrinsic to the drug/channel interaction and less dependent a priori on the protocol used (Di Veroli, Davies, Zhang,
Abi-Gerges, & Boyett, 2013, 2014). However, unless explicitly introducing temperature-dependent parameters, these models cannot accountfor temperature-dependent hERG channel activity changes and234

Z. Li et al. / Journal of Pharmacological and Toxicological Methods 81 (2016) 233–239

therefore are limited in their ability to permit an extrapolation of the results across temperatures.
Because temperature may have an effect on both channel gating anddrug binding rate, such models may need to have two components: abase model dealing with temperature-dependent gating changes without drug, and a pharmacodynamic component simulating temperaturedependent drug binding kinetics. Despite extensive studies on temperature effects on hERG channel gating (Vandenberg et al., 2006; Zhou etal., 1998) and early attempts to build a temperature-dependent hERG
channel model (Fink, Noble, Virag, Varro, & Giles, 2008), a base modelthat can quantitatively reproduce these effects across different labs
and independently predict new experimental data across different temperatures is currently unavailable.
Here we developed a Markov model that can recapitulate macroscopic hERG channel gating behavior for a temperature range of 20 °C
to 37 °C. Based on literature data, this model was calibrated using extensive experimental results obtained at different temperatures, and wasvalidated using independent sets of experiments. Such a model would
improve on existing hERG models by yielding more insight into the molecular basis of temperature-dependent channel activities, and also
could be used as the base model to study temperature-dependent
drug binding kinetics.
2. Methods
2.1. Model development
Several Markov models for hERG channel have been described in theliterature with varying structure and number of states (Clancy & Rudy,
2001; Di Veroli et al., 2013; Fink et al., 2008; Kiehn, Lacerda, & Brown,1999). We based our model on the Di Veroli et al. (Di Veroli et al.,2013) model because this is the only model designed to study drug
binding kinetics. The Markov model uses distinct states to represent different hERG channel conformations and the states transitions aremodeled as ordinary differential equations. All state transition rates
are deﬁned as:
R ¼ A  eBV  qðT−20Þ=10
where R is the state transition rate, V is the membrane potential (electrical ﬁeld across the channel), A and B reﬂect the energy barrier heightin the absence and presence of electrical ﬁeld respectively (Balser,
Roden, & Bennett, 1990), T is the temperature, and q is the commonlyused temperature extrapolating Q10 value deﬁned as the change inrate for each 10 °C change in temperature (Hille, 2001). We chose Q10
over more complex thermodynamic formulations (Fink et al., 2008;
Irvine, Jafri, & Winslow, 1999) to represent temperature sensitivity because it simpliﬁes model parameterization and interpretation. A similar
strategy was used by others in modeling cardiac sodium channel gating(Clancy & Rudy, 1999). Microscopic reversibility was ensured by ﬁxing
the products of forward and reverse transition rates of closed loops(Hille, 2001). The model was implemented in R (http://www.rproject.org), and the equations are numerically solved by the R package
deSolve (http://desolve.r-forge.r-project.org). Model ﬁtting was doneon a computer cluster hosted by High Performance Computing Lab of
FDA using GAPSM (Genetic Algorithm-based Parameterization for Systems Modeling), a genetic algorithm-based ﬁtting method we developed in the process of building other large scale systems biology/
pharmacology models (Li, Zhou, Lu, & Colatsky, 2014). The differential
equations and parameter values are described in the Online
supplement.

2013; Vandenberg et al., 2006; Zhou et al., 1998). The current value
can be approximated using the model's open state as follows: I =G * O * (V − EK), where I is the current, G is the conductance, O is theopen state of the model, V is the voltage and EK is the potassium reversal
potential. For normalized current and conductance this equation can be
simpliﬁed to O/Omax, where Omax is the maximum open state value.
2.2.1. Fitting the model to hERG channel biophysical data at 20 °C from Di
Veroli et al. (2013) study
2.2.1.1. Deactivation time course protocol. Membrane potential was heldat −80 mV for 2 s and then a step to +60 mV was applied for 4 s before
stepping to Vtest (−30, −45, −60, −75 and −90 mV) for 3 s. The ﬁrst
1000 ms upon stepping to Vtest were plotted normalized to the peakcurrent value.
2.2.1.2. Activation time course protocol. An envelope of tails protocol wasused to elucidate activation kinetics. Membrane potential was held at− 100 mV for 3 s and then a step to Vtest (+ 60, + 40, + 20, + 10,+0 mV) for t1 was applied. Another step to −100 mV for 3 s was applied followed by another step to Vtest for t2 N t1. The process was repeated eight times for different durations. The activation time course wasobtained by taking peak values upon each repolarization to −100 mV
and normalizing to the value obtained for Vtest = +60 mV.
2.2.1.3. Activation current voltage protocol. Membrane potential was heldat − 80 mV for 2 s and then a step to Vtest was applied for 4 s before
stepping to − 80 mV for 2 s. The protocol was applied 10 times with
Vtest ranging from −70 mV to +70 mV. Peak values upon repolarizationto − 80 mV were used to build the activation curve by normalizing to
the value obtained with Vtest = +70 mV.
2.2.1.4. Inactivation curve protocol. Membrane potential was held at−80 mV for 2 s and then a step to +60 mV was applied for 4 s. Anotherstep to a voltage of Vtest (9 values between −100 mV and +60 mV) wasapplied for 3 s afterwards. Peak values during Vtest were used to build
the inactivation curve by normalizing to the value obtained with
Vtest = −100 mV.
2.2.1.5. Inactivation rates protocol. For Vtest ranging from − 100 mV to−40 mV, the same protocol as the inactivation curve protocol (deﬁned
above) was used. The initial current shape (until it reaches its maximum) upon repolarization to Vtest was used to generate reactivation
rates by ﬁtting a single exponential function and extracting the time
constant (the exponential recovery equation is described below). For
Vtest ranging from −30 mV to +50 mV, membrane potential was heldat − 80 mV for 2 s and then a step to + 60 mV was applied for 5 s.
Then a brief step to −100 mV lasting 6 ms was applied before stepping
to Vtest for 4 s. Finally, another step to −80 mV for 4 s was added afterward. The initial part of the current shape obtained upon depolarization
to Vtest (the ﬁrst 100 ms) was used to derive inactivation rates by ﬁtting
a single exponential function and extracting the time constant (the exponential decay equation is described below).
2.2.2. Fitting the model to hERG channel biophysical data at various temperatures from Vandenberg et al. (2006) study
2.2.2.1. Activation curve protocol. Membrane potential was initially set at−80 mV, and then it was stepped up to Vtest (ranging from −80 mV to+40 mV) for 30 s before stepping down to −120 mV. Peak values during Vtest were used to build the inactivation curve by normalizing to thevalue obtained with Vtest = +40 mV.

2.2. Model simulations
A series of voltage protocols were applied to the model in order tocalibrate and validate it against hERG biophysical data (Di Veroli et al.,2.2.2.2. Inactivation curve protocol. Membrane potential was initially set
at −80 mV, and then it was stepped up to +40 mV for 1 s to fully activate and inactivate the channels, before stepping to Vtest (ranging from

Z. Li et al. / Journal of Pharmacological and Toxicological Methods 81 (2016) 233–239−150 mV to +60 mV). The peak amplitude of the tail current at eachpotential, after recovery from inactivation, was calculated and corrected
for channel deactivation (Vandenberg et al., 2006). For Vtest b −20 mV,
an exponential curve was ﬁtted to the decay part and extrapolated thepeak value. Data was then normalized to Vtest = −150 mV.
2.2.2.3. Activation time course protocol. Membrane potential was
depolarized to 0 mV for various durations (between 100 and3500 ms) before stepping down to − 50 mV. Peak values upon
repolarizing to − 50 mV were normalized to the step with maximum
duration 3500 ms.
2.2.2.4. Protocol to calculate time constants of inactivation, deactivation
and recovery. A triple pulse protocol was applied, consisting of a depolarization step to +40 mV for 1 s, a repolarization step to −80 mV for
20 ms, followed by a voltage step to +40 mV for 100 ms and ﬁnally astep to −120 mV. To calculate the recovery constant, a single exponential (as described below) was ﬁtted to the ﬁrst part of the rising phase
after stepping to −120 mV. The deactivation time constant was calculated by ﬁtting the decay part after the maximum is reached when
stepping to − 120 mV with a single exponential. For the inactivation
time constant, we ﬁt a single exponential function to the second+40 mV step.
2.2.3. Predicting independent hERG biophysical data at 35 °C from Zhou etal. (1998)
2.2.3.1. Activation current voltage protocol. Membrane potential was heldat − 80 mV and then a depolarization step to Vtest (ranging from− 60 mV to +50 mV) was applied for 4 s, after which the membranepotential was held at − 50 mV for 5 s. The peak tail current upon
stepping to −50 mV was normalized to Vtest = +50 mV.
2.2.3.2. Activation time course protocol. Membrane potential was held at− 80 mV, followed by a depolarizing step to 0 mV for different times(between 10 and 450 ms) and a repolarizing step to − 60 mV. Peaktail current upon repolarization to −60 mV was normalized to the longest step time.
2.2.3.3. Deactivation time course protocol. From − 80 mV, the potential
was stepped to +60 mV and then to −70 mV. The decay of the current
upon stepping to −70 mV was recorded and normalized to the maximum value reached.
2.2.3.4. Protocol to calculate recovery from inactivation time constants.
Membrane potential was held at − 80 mV, before stepping to+ 60 mV for 200 ms and then to Vtest (ranging from − 100 mV to− 20 mV) for 100 ms. Recovery time constants were calculated byﬁtting a single exponential function (as deﬁned below) upon stepping
to Vtest.
2.2.3.5. Protocol to calculate inactivation time constants. Membrane potential was held at − 80 mV, followed by a depolarizing step to+ 60 mV for 200 ms, then to − 100 mV for 2 ms and ﬁnally to Vtest(ranging between −20 mV and +60 mV). The inactivation time constant was calculated by ﬁtting a single exponential function upon
stepping to Vtest.
2.2.4. Predicting the three-step protocol at room temperature
Membrane potential was held at − 80 mV, followed by adepolarizing step to + 20 mV for 1 s. The membrane potential wasthen stepped down to − 40 mV for 500 ms and ﬁnally to − 120 mV
for 500 ms before stepping back to −80 mV. The current was approximated using the equation described in Section 2.2 (I = G * O * (V − EK))and normalized to the absolute maximum value reached. The potassium

235

Koreversal potential (Ek) was calculated as ðRTF Þ  logð Ki Þ , where R is thegas constant (8314.0 J/kmol/K), T is the absolute temperature, F is the
Faraday constant (96,485.0 coul/mol), Ko is the extracellular potassium
concentration (5 mM), and Ki is the intracellular potassium concentration (140 mM).

2.2.5. Equations
The equations described here were used to calculate recovery, inactivation and deactivation time constants, as well as voltage dependence(Vhalf) of I-V relationship as in Milnes et al. (Milnes, Witchel, Leaney,
Leishman, & Hancox, 2010). Recovery time constants were calculated
by ﬁtting to the following equation:

O ¼ Span  1−e−kt þ O min ;where O is the open state, t is time, Span is the difference between theminimum fraction of open channels, Omin, and the extrapolated plateau
value, k is the rate-constant of the exponential recovery and the time
constant can be calculated as 1 k . Inactivation and deactivation time
constants were calculated by ﬁtting to the following equation:
 t
O ¼ A  e−τ þ C;
where A is the total fraction of open channels, t is time, τ is the time constant of decay and C is the residual fraction of open channels. The inactivation I-V relationship curves were ﬁtted by a modiﬁed Boltzmann
equation as follows:I¼ 1−1= 1 þ eðVhalf−VmÞ=S
Imax
where I is the current corresponding to membrane potential Vm, Imax isthe maximal current observed, Vhalf is the half inactivation voltage, andS is the slope factor describing the I-V relationship.
2.3. Electrophysiology
For hERG current recording, a HEK293 cell line stably expressinghERG1a protein (generously supplied by Dr. Gail Roberton, Universityof Wisconsin-Madison) was cultured in MEM supplemented with
non-essential amino acids solution, 10% FBS and 10 μM G418 (GIBCO,
Grand Island, NY). During electrophysiological recording, cells wereconstantly perfused with extracellular solution containing (in mM):
130 NaCl, 5 KCl, 10 HEPES, 1 MgCl2, 1 CaCl2, 12.5 Dextrose (pH 7.4, adjusted with NaOH). The internal pipette solution contains (in mM):
120 K-Gluconate, 20 KCl, 10 HEPES, 5 EGTA, 1.5 MgATP (pH 7.3, adjusted with KOH). MultiClamp 700B, Axon Digidata 1550A and pClamp software (Molecular Devices, Sunnyvale, CA) were used for whole-cellrecording and data acquisition. All voltage commands were correctedwith a − 15 mV liquid junction potential resulting from using theabove internal solution.
3. Results
As shown in Fig. 1, the hERG model has six states corresponding toopen, closed, or inactivated channel conformations, with each statetransition governed by membrane voltage and temperature. Each transition is governed by a different set of free parameters, in contrast to theoriginal Di Veroli hERG model in which many state transitions were assumed to have the same parameters. By forcing some state transition
parameters to be identical, the number of free parameters is reduced
but the model becomes less ﬂexible and robust. We loosened the constraints on the state transition parameters by allowing them to be independent from each other in order to build a robust model that accounts

236

Z. Li et al. / Journal of Pharmacological and Toxicological Methods 81 (2016) 233–239

Fig. 1. Structure of the hERG model. C1 and C2 are closed states, O is open state, and thecorresponding inactivated states are IC1, IC2, and IO. The transition between adjacentstates is a ﬁrst order reaction dependent on membrane voltage, temperature, with three
free parameters A, B, and q (see Methods). Each state transition has a different set offree parameters, which are distinguished between each other by numeric sufﬁxes.

for the breadth of published data across different temperatures andstudies.
The new hERG model was ﬁtted simultaneously to experimental
data from two comprehensive studies to estimate free parameters.
The Di Veroli study used CHO cells to study the biophysical propertiesof hERG channel at 20 °C on an automated patch clamp system (DiVeroli et al., 2013), while the Vandenberg study used the same cellline and a different set of protocols to probe channel gating processes
at various temperatures using a manual patch clamp system(Vandenberg et al., 2006). The ﬁtting of our model to the Di Veroli
dataset at 20 °C and the Vandenberg study at multiple temperaturesare shown in Figs. 2 and 3 respectively. As shown in Fig. 2, while somedata points were ﬁtted well, others show less good a ﬁt compared tothe original Di Veroli model (Di Veroli et al., 2013), especially the inactivation I-V curve (Fig. 2 D) and recovery time constants (ﬁrst half of Fig.
2 E). Careful examination suggests that this is due to the variability between the room temperature data from Di Veroli and Vandenberg studies. For instance, the Di Veroli inactivation I-V curve (circles in Fig. 2 D)
would result in a half inactivation voltage (Vhalf) of −20 mV, while theroom temperature inactivation I-V curve in Vandenberg study (nonﬁlled circles in Fig. 3 B) gave a Vhalf of −88 mV. Because the inactivation curve in Fig. 3 B was corrected for deactivation (Vandenberg etal., 2006) while Fig. 2 D was not, we also calculated a Vhalf using noncorrected inactivation I-V curve based on Fig. 5C of the original Vandenberg paper (Vandenberg et al., 2006) to make it more comparable to Di
Veroli's value in Fig. 2 D. This resulted in an inactivation Vhalf of− 73 mV, still much more negative than Di Veroli's value (− 20 mV).
This discrepancy is probably due to lab-to-lab and system-to-system(manual vs automated) variability. Trying to ﬁt values from both dataset, the model computed an inactivation I-V relationship (solid line in
Fig. 2 D) with a Vhalf of − 56 mV to achieve an overall best ﬁt to all
the data, which also leads to a seemingly bad ﬁtting on Fig. 2 D. Interestingly, this calculated Vhalf is close to at least two other studies reporting
room temperature inactivation Vhalf values (Lu, Mahaut-Smith, Huang,& Vandenberg, 2003; B. D. Walker et al., 2000), suggesting our model
was able to “average out” the occasional biased values found in any particular study. A similar pattern was found for the room temperature

Fig. 2. Fitting to hERG biophysical data at 20 °C. Experimental data, taken from Di Veroli et al. (2013), are shown as symbols, and simulated data are shown as solid lines. A, fraction of open
channels over time during deactivation. B, fraction of open channels over time during activation. C, activation curve. D, inactivation curve. E, recovery and inactivation rate. The top panels
show the voltage protocols used. Note that points in panel E were calculated using two separate protocols: recovery time constants (−100 to −40 mV) using the protocol 2D while
inactivation time constants (−30 to +50 mV) using the protocol in 2E.

Z. Li et al. / Journal of Pharmacological and Toxicological Methods 81 (2016) 233–239

237

Fig. 3. Fitting to hERG channel experiments at various temperatures. Experimental data, taken from Vandenberg et al. (Vandenberg et al., 2006), are shown as symbols and simulated data
are shown as solid lines. A, activation curve. B, inactivation curve. C, activation time course. D, time constants of various processes at room or physiological temperatures; Measured areexperimentally determined data from Vandenberg paper while Fitted are simulated using our model. The insets show the voltage command protocol used. Note that for A and B, theoriginal Vandenberg study contains 14 °C data, which are not used for model ﬁtting here because the main purpose is to extrapolate between room temperature and physiologicaltemperature.

recovery data. Di Veroli data (circles in Fig. 2 E) suggested a time constant (τ) of ~1.9 ms at −100 mV. In contrast, Vandenberg study reported a room temperature recovery τ of 3.8 ms at −120 mV (Fig. 3 D), andjudging by the voltage dependency of recovery speed, one would expect
an even larger τ at −100 mV. Even though Di Veroli study temperature(20 °C) is lower than that of Vandenberg study (24 °C), this small temperature difference cannot explain the N 2 fold difference in recovery
speed at the same voltage (−100 mv). Again, the model found a middle
ground solution (5.8 ms at −100 mV) that minimized the overall error,
but at the same time created a seemingly bad ﬁtting (Fig. 2 E). In terms
of ﬁtting to multiple temperature data from Vandenberg study, ourmodel correctly reproduced the left shift of the activation curve andthe right shift of the inactivation curve as temperature increases (Fig.
3 A and B). The signiﬁcant increase in activation speed after temperature increases was quantitatively captured by the model (Fig. 3 C),
and so were the temperature sensitivities of deactivation, inactivation,and recovery from inactivation (Fig. 3 D). Even though a few data points
were not ﬁtted perfectly, the overall ﬁtting to the Vandenberg study in
Fig. 3 is reasonable, suggesting our model can reproduce temperaturedependent changes in various hERG channel gating processes.
The data presented above were used to parameterize the model (theﬁtted parameters can be found in the Online supplement). While it isencouraging to see that a single set of parameters can reproduce the results of an ensemble of different voltage protocols across different temperatures and labs, it is more important to determine whether this same
parameter set can predict the results of an independent set of experimental data not used in the ﬁtting process. For the purpose of validation,
we chose a third set of data using hERG-transfected HEK cells recordedat near physiological temperature (Zhou et al., 1998). As shown in Fig. 4,

the current model is able to predict quantitatively most of the experimental results in this independent data set, including the activation IV relationship (Fig. 4 A), activation/deactivation time course (Fig. 4 B
and C), as well as time constants for recovery (Fig. 4 D) and inactivation(Fig. 4 E). In addition, to test if our model can also independently predict
new experimental results from room temperature, we challenged themodel with a new three-step protocol (Fig. 5) and then compared thepredicted traces to those experimentally obtained using HEK cells atroom temperature in our laboratories. As shown in Fig. 5, the model
was able to predict the amplitude and shape of the currents across thethree steps. Both the outward (− 40 mV step) and the inward(−120 mV step) tail currents were predicted very well. The computed
pre-pulse currents (+20 mV step) are a little low compared to experimental values, but are within the standard deviation of the data.
Taken together, this suggests that the hERG model is able to predict
the experimental results from protocols and cell lines not used in theﬁtting process. This feature is very useful if the model were to be used
to study drug-channel interactions, where new voltage protocols will
be designed to study different aspects of drug binding/unbinding
kinetics.
4. Discussion
Here we presented a computational hERG model for which the channel gating parameters are explicitly dependent on temperature. To ourbest knowledge, this is the ﬁrst hERG model that can quantitatively reproduce and independently predict the temperature sensitivities ofhERG channel gating properties. Incorporating temperature into themodel allows it to predict how temperature changes affect gating

238

Z. Li et al. / Journal of Pharmacological and Toxicological Methods 81 (2016) 233–239

Fig. 4. Prediction of independent experiments at 35 °C in HEK cells. The hERG model was used to predict the outcome of experiments from an independent study (Zhou et al., 1998). Solidlines are model prediction while symbols are experimental data. A. Activation I-V relationship. B. Activation time course. C. Deactivation time course. D. Recovery from inactivation time
constants. E. Inactivation time constants. The insets show the voltage command protocol used.

processes under various voltage protocols. After supplementing thebase model with a pharmacodynamics component describing temperature-dependent drug binding kinetics, the full model can potentiallyallow for the extrapolation of hERG assay data obtained at room temperature to models of the human cardiac ventricular myocyte at 37 °C.
The Markovian nature of this model makes it straightforward to simulate state-dependent blocking and dissociation, an important featureof many hERG-blocking drugs (Stork et al., 2007; B.D. Walker et al.,Fig. 5. Prediction of independent experiments at room temperature in HEK cells. A new
three-step voltage protocol was applied to HEK cells at room temperature (24 °C) andthe results (symbols) are compared to independent prediction of the model (lines).
Error bars indicating standard deviation are also shown (N = 3). The inset shows thevoltage command protocol used. The three sections of symbols correspond to the prepulse (+20 mV), the ﬁrst test pulse (−40 mV), and the second test pulse (−120 mV)
respectively.

1999). Thus the current model provides an initial framework for studying temperature-dependent hERG channel gating changes and the kinetics of drug block.
As with any computational approach to simulate the real world behaviors, it is critical to bear in mind the strength and weakness of thecurrent hERG model. First of all, because the goal is to build a “generic”
hERG model that can reproduce the widest range of data possible, themodel was calibrated and validated using data across different laboratories, different recording systems (manual and automated patch system),and different cell lines (CHO and HEK cells). While this strategy allows
for broader application of the model, it also introduces signiﬁcant potential inter-laboratory and inter-platform variability into the data. As
can be seen from Figs. 2 and 3, sometimes the ﬁtting to a particular
dataset was not optimal because the best ﬁt to that dataset was
sacriﬁced to achieve a better overall ﬁtting to the whole set. Thus applications in a particular experimental setting may require further adjustment of the parameters. Secondly, for a model with such complexity
and a training data set with such variability, inevitably some parameters
are not well constrained (Irvine et al., 1999). As a result, the set of parameters presented here is not the only possible solution during the parameterization process. However, the fact that this parameter set can
predict the results of independent experiments suggests that thismodel captures some underlying biological mechanisms from the training data rather than “memorizes” or overﬁts the data. Lastly, this is abase model with only channel gating processes, and drug binding/unbinding processes will need to be added later for each speciﬁc drug. Currently for most compounds it is unknown how much an effecttemperature has on drug binding kinetics. It has been suggested thatthe drug binding processes might be less temperature sensitive (DiVeroli et al., 2013) and modeling the temperature effects on gating

Z. Li et al. / Journal of Pharmacological and Toxicological Methods 81 (2016) 233–239

processes, but not drug binding processes, may be enough to explain the
drug-channel interactions at different temperatures (Vandenberg et al.,2006). However, some drugs/toxins do show temperature sensitivity intheir binding/unbinding rates (Hill, Sunde, Campbell, & Vandenberg,2007) and in those cases, drug-speciﬁc temperature parameters(Q10s) will have to be introduced into the drug binding component
and estimated by doing channel blocking experiments at multiple temperatures. Currently a series of compounds are being tested at bothroom and physiological temperatures to add a pharmacodynamic component to the base model and to evaluate the necessity of modeling the
drug binding processes in a temperature-dependent manner.
In summary, here we presented a computational hERG model thatcan quantitatively reproduce and predict temperature-induced changes
in major gating processes. This model improves on existing hERG Markov models by incorporating temperature dependence of state transitions, and could be used to probe dynamic drug-channel interactions
in a temperature dependent manner. Currently the model is being further developed (e.g. to make its gating behavior more closely resemble
human native IKr currents in vivo) and integrated into the consensus
human ventricular cardiomyocyte model (O′Hara Rudy, or ORd)
(O'Hara, Virag, Varro, & Rudy, 2011) to extract dynamic drug-channelinteraction for the compounds on the CiPA training list (Fermini et al.,2016)
Conﬂict of interest statement
The authors declared no conﬂict of interest.
Sources of funding
This project was funded by FDA Critical Path Award. The study sponsor(s) played no role in the study design; in the collection, analysis, andinterpretation of data; in the writing of the report; and in the decision tosubmit the paper for publication.
Disclaimer
This report is not an ofﬁcial U.S. Food and Drug Administration guidance or policy statement. No ofﬁcial support or endorsement by the U.S.
Food and Drug Administration is intended or should be inferred.
Appendix A. Supplementary data
Supplementary data to this article can be found online at http://dx.
doi.org/10.1016/j.vascn.2016.05.005.
References
Balser, J. R., Roden, D. M., & Bennett, P. B. (1990). Global parameter optimization for cardiac potassium channel gating models. Biophysical Journal, 57, 433–444.
Clancy, C. E., & Rudy, Y. (1999). Linking a genetic defect to its cellular phenotype in a cardiac arrhythmia. Nature, 400, 566–569.
Clancy, C. E., & Rudy, Y. (2001). Cellular consequences of HERG mutations in the long QT
syndrome: Precursors to sudden cardiac death. Cardiovascular Research, 50, 301–313.
Di Veroli, G. Y., Davies, M. R., Zhang, H., Abi-Gerges, N., & Boyett, M. R. (2013). Highthroughput screening of drug-binding dynamics to HERG improves early drug safetyassessment. American Journal of Physiology - Heart and Circulatory Physiology, 304,
H104–H117.

239

Di Veroli, G. Y., Davies, M. R., Zhang, H., Abi-Gerges, N., & Boyett, M. R. (2014). hERG inhibitors with similar potency but different binding kinetics do not pose the sameproarrhythmic risk: Implications for drug safety assessment. Journal of
Cardiovascular Electrophysiology, 25, 197–207.
Fermini, B., Hancox, J. C., Abi-Gerges, N., Bridgland-Taylor, M., Chaudhary, K. W., Colatsky,T., ... Vandenberg, J. I. (2016). A new perspective in the ﬁeld of cardiac safety testing
through the comprehensive in vitro proarrhythmia assay paradigm. Journal of
Biomolecular Screening, 21, 1–11.
Fink, M., Noble, D., Virag, L., Varro, A., & Giles, W. R. (2008). Contributions of HERG K+
current to repolarization of the human ventricular action potential. Progress inBiophysics & Molecular biology, 96, 357–376.
Hill, A. P., Sunde, M., Campbell, T. J., & Vandenberg, J. I. (2007). Mechanism of block of thehERG K+ channel by the scorpion toxin CnErg1. Biophysical Journal, 92, 3915–3929.
Hille, B. (2001). Ion channels of excitable membranes (3rd ed.).
ICH S7B: Note for guidance on the nonclinical evaluation of the potential for delayed ventricular repolarization (QT interval prolongation) by human pharmaceuticals. In. London: The international conference on harmonisation of technical requirements forregistration of pharmaceuticals for human use (pp. 3–9).
Irvine, L. A., Jafri, M. S., & Winslow, R. L. (1999). Cardiac sodium channel Markov modelwith temperature dependence and recovery from inactivation. Biophysical Journal,
76, 1868–1885.
Jurkiewicz, N. K., & Sanguinetti, M. C. (1993). Rate-dependent prolongation of cardiac action potentials by a methanesulfonanilide class III antiarrhythmic agent. Speciﬁc
block of rapidly activating delayed rectiﬁer K+ current by dofetilide. Circulation
Research, 72, 75–83.
Kiehn, J., Lacerda, A. E., & Brown, A. M. (1999). Pathways of HERG inactivation. American
Journal of Physiology, 277, H199–H210.
Kirsch, G. E., Trepakova, E. S., Brimecombe, J. C., Sidach, S. S., Erickson, H. D., Kochan, M. C.,... Brown, A. M. (2004). Variability in the measurement of hERG potassium channel
inhibition: Effects of temperature and stimulus pattern. Journal of Pharmacologicaland Toxicological Methods, 50, 93–101.
Li, Z., Zhou, H., Lu, Y., & Colatsky, T. (2014). A critical role for immune system response inmediating anti-inﬂuenza drug synergies assessed by mechanistic modeling. CPT:
Pharmacometrics & Systems Pharmacology, 3, e135.
Lu, Y., Mahaut-Smith, M. P., Huang, C. L., & Vandenberg, J. I. (2003). Mutant MiRP1 subunits modulate HERG K+ channel gating: A mechanism for pro-arrhythmia in long
QT syndrome type 6. Journal of Physiology, 551, 253–262.
Malik, M., & Camm, A. J. (2001). Evaluation of drug-induced QT interval prolongation: Implications for drug approval and labelling. Drug Safety, 24, 323–351.
Milnes, J. T., Witchel, H. J., Leaney, J. L., Leishman, D. J., & Hancox, J. C. (2010). Investigating
dynamic protocol-dependence of hERG potassium channel inhibition at 37 degrees C:
Cisapride versus dofetilide. Journal of Pharmacological and Toxicological Methods, 61,
178–191.
O'Hara, T., Virag, L., Varro, A., & Rudy, Y. (2011). Simulation of the undiseased human cardiac ventricular action potential: Model formulation and experimental validation.
PLoS Computational Biology, 7, e1002061.
Redfern, W. S., Carlsson, L., Davis, A. S., Lynch, W. G., MacKenzie, I., Palethorpe, S., ...
Hammond, T. G. (2003). Relationships between preclinical cardiac electrophysiology,
clinical QT interval prolongation and torsade de pointes for a broad range of drugs:
Evidence for a provisional safety margin in drug development. Cardiovascular
Research, 58, 32–45.
Sanguinetti, M. C., Jiang, C., Curran, M. E., & Keating, M. T. (1995). A mechanistic link between an inherited and an acquired cardiac arrhythmia: HERG encodes the IKr potassium channel. Cell, 81, 299–307.
Stork, D., Timin, E. N., Berjukow, S., Huber, C., Hohaus, A., Auer, M., & Hering, S. (2007).
State dependent dissociation of HERG channel inhibitors. British Journal of
Pharmacology, 151, 1368–1376.
Vandenberg, J. I., Varghese, A., Lu, Y., Bursill, J. A., Mahaut-Smith, M. P., & Huang, C. L.
(2006). Temperature dependence of human ether-a-go-go-related gene K + currents. American Journal of Physiology - Cell Physiology, 291, C165–C175.
Walker, B. D., Singleton, C. B., Bursill, J. A., Wyse, K. R., Valenzuela, S. M., Qiu, M. R., ...
Campbell, T. J. (1999). Inhibition of the human ether-a-go-go-related gene (HERG)
potassium channel by cisapride: Afﬁnity for open and inactivated states. British
Journal of Pharmacology, 128, 444–450.
Walker, B. D., Singleton, C. B., Tie, H., Bursill, J. A., Wyse, K. R., Valenzuela, S. M., ...
Campbell, T. J. (2000). Comparative effects of azimilide and ambasilide on thehuman ether-a-go-go-related gene (HERG) potassium channel. Cardiovascular
Research, 48, 44–58.
Zhou, Z., Gong, Q., Ye, B., Fan, Z., Makielski, J. C., Robertson, G. A., & January, C. T. (1998).
Properties of HERG channels stably expressed in HEK 293 cells studied at physiological temperature. Biophysical Journal, 74, 230–241.