Sinusoidal voltage protocols for rapid characterisation of ion channel kinetics

Key points Ion current kinetics are commonly represented by current–voltage relationships, time constant–voltage relationships and subsequently mathematical models fitted to these. These experiments take substantial time, which means they are rarely performed in the same cell. Rather than traditional square‐wave voltage clamps, we fitted a model to the current evoked by a novel sum‐of‐sinusoids voltage clamp that was only 8 s long. Short protocols that can be performed multiple times within a single cell will offer many new opportunities to measure how ion current kinetics are affected by changing conditions. The new model predicts the current under traditional square‐wave protocols well, with better predictions of underlying currents than literature models. The current under a novel physiologically relevant series of action potential clamps is predicted extremely well. The short sinusoidal protocols allow a model to be fully fitted to individual cells, allowing us to examine cell–cell variability in current kinetics for the first time. Abstract Understanding the roles of ion currents is crucial to predict the action of pharmaceuticals and mutations in different scenarios, and thereby to guide clinical interventions in the heart, brain and other electrophysiological systems. Our ability to predict how ion currents contribute to cellular electrophysiology is in turn critically dependent on our characterisation of ion channel kinetics – the voltage‐dependent rates of transition between open, closed and inactivated channel states. We present a new method for rapidly exploring and characterising ion channel kinetics, applying it to the hERG potassium channel as an example, with the aim of generating a quantitatively predictive representation of the ion current. We fitted a mathematical model to currents evoked by a novel 8 second sinusoidal voltage clamp in CHO cells overexpressing hERG1a. The model was then used to predict over 5 minutes of recordings in the same cell in response to further protocols: a series of traditional square step voltage clamps, and also a novel voltage clamp comprising a collection of physiologically relevant action potentials. We demonstrate that we can make predictive cell‐specific models that outperform the use of averaged data from a number of different cells, and thereby examine which changes in gating are responsible for cell–cell variability in current kinetics. Our technique allows rapid collection of consistent and high quality data, from single cells, and produces more predictive mathematical ion channel models than traditional approaches.

1 Figure 1 of the main text features simulations from 29 literature hERG or I Kr models. In Table A1 2 we list these models, give references, and show the seven different structures that they feature in 3 Figure Figure A1: Different mathematical model structures for the literature models listed in Table A1. The model we use in the main text takes structure D as shown in Figure 4B. Note that, depending on their parameterisations, models D, E, F and G could satisfy independent gating (Hodgkin-Huxley assumptions) and then the states annotated as 'I' could also be considered 'IO' -but as drawn here this is not a requirement, and in neither case does current flow, so we removed the 'O' from these for simplicity.
There are two models in this table that are not in Figure   indicates the total number of free parameters (the number given includes a G Kr parameter for the conductance). The Model Type is 'HH' for Hodgkin-Huxley models and 'MM' for Markov Models, or a hybrid of the two (MM/HH, which generally means a MM with some symmetry in transition rates). Temp. represents temperature and is given as 'Room' for room temperature or 'PT' for physiological temperature. Structure respresents the model structure as shown iin Figure A1 We see that models are calibrated to experimental datasets from different conditionsspecies, cell types and temperatures. These conditions may not be the same temperature, cell type or even species as the models are intended to represent (see Niederer et al., 2009), and so the final column indicates the species of the cell model in which the hERG model was used.  Kiehn et al. (1999) model are defined at specific voltages, so for this model there are 8 parameters (and 1 conductance parameter) for each voltage at which the model is defined. each repeat. From an initial holding potential of −80 mV, this protocol comprised a 5 s step to 46 10 mV followed by a 5 s step to −50 mV before returning again to a holding potential of −80 mV. This protocol is depicted in Figure B2. We repeated this protocol while dofetilide was added (see 48 Figure 2) and the current traces recorded from this protocol were used to assess when a steady level 49 of dofetilide block had been reached. 50 Protocols 1,2 -Activation Kinetics 51 After the initial period at holding potential incorporating the −120 mV leak step, a step to V step 1 52 followed and was held at that voltage for T step ms, before a step to −120 mV for 2.5 s, before 53 returning to holding potential of −80 mV for 1 second. The protocol was repeated 6 times with a 54 different T step on each repeat. T step took the values of 3, 10, 30, 100, 300 and 1000 ms.

55
For Protocol 1, V step 1 is 0 mV. This protocol is depicted in Figure B3.

56
For Protocol 2, V step 1 is +40 mV. This protocol is depicted in Figure B4.    Protocol 3 -Steady-State Activation 60 From the initial period at holding potential incorporating the −120 mV leak step, a step to V step 61 was applied for 5 seconds, followed by a 1 s step to −40 mV, before a 500 ms step to −120 mV, and 62 then returning back to holding potential for one second. This process was repeated 7 times with a 63 different V step on each repeat. V step ranged from −60 mV to +60 mV in 20 mV increments. This 64 protocol is depicted in Figure 5A (left column).

66
From the initial period at holding potential incorporating the −120 mV leak step, a step to 50 mV 67 for 600 ms, and a step to −90 mV for 60 ms, followed by a step to V step for 150 ms, before a 500 ms 68 step to −120 mV, and a 1 s step back to holding potential of −80 mV; This was repeated 16 times 69 with a different V step on each repeat. V step ranged from −100 mV to 50 mV in 10 mV increments.

70
This protocol is depicted in Figure 5A (middle column).

71
From the initial period at holding potential incorporating the −120 mV leak step, a step to 50 mV for 2 s was applied, followed by a step to V step for 6 s, before a 500 ms step to −120 mV, and then 74 returning back to holding potential for one second. This process was repeated 9 times with a 75 different V step on each repeat. V step ranged from −120 mV to −40 mV in 10 mV increments. This 76 protocol is depicted in Figure 5A (right column).

80
The full protocol is comprised of 250 ms at holding potential of −80 mV, followed by a 50 ms 'leak 81 detection' step to −120 mV, and then 200 ms back at −80 mV. This was followed by an 'activation' 82 step of 1 s step to 40 mV; a 'closing' 500 ms step to −120 mV; and a return to −80 mV for 1 second.

83
The 3.5 s sinusoidal portion of the protocol then followed (the form of which is described below), 84 before a 'closing' 500 ms step to −120 mV, and a return to −80 mV for 1 s.

85
The sinusoidal portion of the protocol takes the form of a sum of three sine waves as discussed in the main text Equation (1), and repeated here for convenience: where A 1 = 54 mV, A 2 = 26 mV, A 3 = 10 mV, ω 1 = 0.007 ms −1 , ω 2 = 0.037 ms −1 and ω 3 = 86 0.19 ms −1 , and t is time measured in milliseconds.

87
The protocol was initially designed with just the −120 mV leak step and not the additional 88 'activation' steps to 40 mV and −120 mV (which were included after preliminary experiments as 89 described in Section B2.2) and so the sine wave was shifted using t 0 = 2500 ms to begin at the same 90 phase after we incorporated the additional steps. This offset is not expected to be important to 91 include, but was included here for clarity for anyone attempting to reproduce our study.

93
All of the protocols described in this section were adjusted on the amplifier to account for the liquid 94 junction potential which was calculated to be 4.1 mV from the ionic composition of our physiological 95 solutions which are described in Section B1.2. The liquid junction potential was calculated using 96 the junction potential calculator in the pClamp software.

98
To show the effect of 0.3 µM dofetilide subtraction and to demonstrate the lack of endogenous 99 currents remaining following this step, we compare the current traces before and after dofetilide 100 subtraction in Figure B5. We show currents from Cell #5 which features in the figures of the main 101 text.

102
There was little contribution from endogenous currents in these cells. The application of 0.3 µM 103 dofetilide eliminates almost all current, indicating that both the vehicle recording and 'dofetilide 104 subtracted' currents are almost entirely due to hERG. While levels of endogenous currents and the 105 impact of leak subtraction on the current traces may vary from cell-to-cell, overall we found only 106 small endogenous currents were observed in the cells.  Figure B5: Raw recordings, recordings with leak subtraction, and 0.3 µM dofetilide controls for Cell #5. Top row: input voltage trace from the sine wave protocol (same on both sides for comparison with traces below). In A-D we show the current in response to this voltage protocol in four situations -each panel shows the same four traces, with a different one highlighted in red. In A we highlight the raw recording of whole current; in B the whole current recording after leak subtraction, in C the raw recording after the addition of dofetilide; and in D the dofetilide recording after leak subtraction.
An example of where dofetilide subtraction may play a more important role is shown in Fig-108 ure B6, where the hERG current is much lower and so any contribution from endogenous currents  Figure B6: Raw recordings, recordings with leak subtraction, and 0.3 µM dofetilide controls for Cell #6. Top row: input voltage trace from the sine wave protocol (same on both sides for comparison with traces below). In A-D we show the current in response to this voltage protocol in four situations -each panel shows the same four traces, with a different one highlighted in red. In A we highlight the raw recording of whole current; in B the whole current recording after leak subtraction, in C the raw recording after the addition of dofetilide; and in D the dofetilide recording after leak subtraction.
V step section of the deactivation protocol. Again, we used the faster time constant as the recovery 128 from inactivation time constant and the slower time constant as that for deactivation.

129
To produce the peak current-voltage relationship for the steady state activation protocol for the 130 simulated data traces we wrote MatLab code to identify the peak current in the region between 131 5.6292 and 5.7292 seconds on each sweep of the protocol, which corresponds to the current response 132 just after the 5 second V step when the voltage is stepped to −40 mV. We then normalised the peak 133 current data to the maximum overall peak identified in this region to produce the current-voltage 134 relationship curve. For the simulated data we wrote a MatLab script (included in code download) 135 to identify the peak-current voltage relationship for this protocol but for the experimental traces 136 we verified these peak points manually to avoid incorrect peaks being identified due to noise or 137 capacitive effects. We also identified the peak currents in the currents evoked by the activation 138 kinetics protocols manually for the same reason. In the activation kinetics protocol we identified 139 the peak currents during the V step for each interval of T step duration. For an observed experimental recording which we will denote y, we can infer the probability of dif- When the prior distribution is assumed to be uniform (as it is in this study), we can make inferences 160 based on just the likelihood, as the prior P (θ) is either constant or zero. If a proposed parameter 161 is outside our chosen prior then likelihood is 0 and we simply record that this parameter set has a 162 likelihood of 0 and propose another parameter set. 163 We assume that the errors at each time point are independent and so the conditional probability We assume that the experimental noise is independently and normally distributed with a mean 167 of zero and variance of σ 2 . The likelihood is then expressed as In our case f t (θ) is the predicted current at each time point given the parameters, this is given 169 by equation (7) after solving the model system (equations (3)- (6)). Calculating equation (B.6) 170 requires the evaluation of the product of many numbers less than 1, so it is more numerically 171 convenient to calculate the log-likelihood instead. As our aim is to identify parameter sets θ which 172 maximise the likelihood in equation (B.6), maximising the likelihood is equivalent to maximising 173 the log-likelihood: In practice, the sums over time in equation (B.7) are formulated so that we exclude time points 175 from regions where the data are affected by capacitive spikes. To be precise, we exclude 5 ms 176 intervals following step-changes in the imposed voltage clamp. In the sine wave protocol (Pr7) these step-changes occur at 0.25 seconds, 0.3 seconds, 0.5 seconds, 1.5 seconds, 2 seconds, 3 seconds, 6.5 seconds and 7 seconds (spikes are seen in experimental recordings at these times in Figure 3).

B2.2 Conductance estimation to inform the prior 180
Preliminary work revealed that using sine wave protocols alone often allowed kinetic parameters in 181 the hERG model to be recovered, but there was potential for identifiability problems (or at least 182 we encountered difficulties in finding a global optimum due to a rugged likelihood surface) when 183 simultaneously fitting the conductance parameter and transition rate parameters P 1 to P 8 (although 184 previous work suggests all parameters are theoretically identifiable (Walch and Eisenberg, 2016)).

185
To add extra information on conductance, we incorporated a voltage-step to +40 mV followed by a  2012)). The conductance we estimated was used as a lower bound for the 193 prior distribution of the conductance, as we describe below. 195 In this section we describe our prior assumptions on the values that each model parameter can take.

196
The prior for the conductance G Kr is assumed to be independent of the kinetic parameters, and to 197 take a uniform distribution. As discussed above, the lower bound is formed by estimating a lower where V is voltage and A and B are model parameters (P 1 to P 8 for k 1 to k 4 , as shown in Figure 4).

202
For parameters of the form A we assumed that the prior distribution is uniform between 10 −7 203 and 1000 ms −1 , again to cover (and extend beyond) the full physiological range expected with hERG 204 channel gating. 205 We assume that the prior distributions for B parameters are uniform between 10 −7 and 0.4 206 mV −1 . The lower bound for this parameter was selected as the voltage-dependence becomes prac-207 tically redundant when B becomes small: when B = 10 −7 the value of exp(BV ) will change by less 208 than 0.0015% across the voltages we reach in this study. The upper value is beyond the physiolog-209 ically expected range. 210 We also impose a prior on the maximum rate of transition k between any states (maximum 211 across the full voltage range in the protocol (that is from −120 to 58.25 mV)). If the maximum 212 rate k is greater than 1000 ms −1 , or less than 1.67 × 10 −5 ms −1 , the pair of parameter values that 213 give rise to this are assigned prior probability zero (strictly, this is equivalent to defining 2D prior fr/~hansen/cmaes.m. We imposed bounds based on the prior as we describe above in Section B2.3.

235
We run the CMA-ES algorithm from different starting points and continue to do so until we 236 identify the same region of parameter space for optimal parameter sets for each experimental data 237 trace when starting from many different starting points. In this way, we can be confident that we 238 identify the same region of high likelihood consistently (not simply the first local minimum that is 239 found), and we have more confidence that this corresponds to the globally optimal likelihood.

240
These initial starting points for the CMA-ES algorithm are sampled from within the prior 241 defined for each parameter, described in section B2.3. To sample from the prior we simply select 242 the voltage-dependent transition rate parameters (of the form B described above) uniformly from 243 the defined range. The same approach is used to sample the conductance parameter.

244
For the parameters of the form A above we sample starting points in a logarithmic fashion across 245 the range of the uniform prior. This approach helps to restrict the initial guesses of parameters to 246 the region of measurable time scales we imposed by defining the maximum and minimum ranges on 247 the overall transition rate, as described above. We also run a small selection of starting points with  (1996)).

275
In order to verify that there was sufficient information within the sinusoidal protocol to parameterise 277 our model we performed a synthetic data study. The aim in such a study is to ascertain whether we 278 can recover the parameters used in the simulation from a simulated data trace (with added noise 279 in this case).

C1 Producing synthetic data 281
In order to produce synthetic data we simulate with some fixed known parameter set. We performed 282 this with both our best initial parameter set estimate (those parameters in literature HH models), 283 and also from the parameters we obtained after fitting to the experimental data trace as we present

292
C2 Inferring parameters from synthetic data 293 We then attempt to infer parameters from this synthetic data trace, using the CMA-ES algorithm 294 followed by MCMC as described in Section 2.5. In Figure C7 we present probability density distri-295 butions obtained when using both synthetic and experimental traces. We are able to recover the 296 original parameters underlying the synthetic trace with high accuracy.

297
The synthetic data study provides us with confidence in the suitability of our protocol for 298 accurately identifying parameters of the model presented in Figure 4B in the main text, and also 299 that the parameter inference protocol(s) we are using are suitable for the task. We believe such an 300 approach should always be used to test whether there is sufficient information in the experimental 301 data being proposed for calibration of a mathematical model. The test should be performed twice: 302 before conducting the experiment (with the pre-existing best guess at the parameters); and also after 303 conducting the experiment (with the new maximum posterior density estimate of the parameters 304 -as we illustrate in Figure C7). 305 Figure C7: Probability density distributions for each parameter estimates from fitting to both experimental data (red) and simulated data (blue). Crosses indicate the parameter set with the maximum posterior density.

306
In Tables D2-D10 we Figure E8 shows the summary curves for Pr1 (voltage clamp shown in Figure B3) 333 and Pr2 (voltage clamp shown in Figure B4).

334
Traditionally these peak current curves would be plotted by normalizing to the peak current  Figure E8: Predictions of peak current-voltage relationship derived from experiment and model predictions in response to; A) Activation Kinetics Pr1, B) Activation Kinetics Pr2, with comparison of our model prediction with predictions from existing literature models. Currents have been normalized to the peak current in the initial deactivation step in the sine wave protocol (around 1.6 seconds) as we do not expect the channel to be fully open at the longest T step in these activation kinetics protocols.
In Figures E9 and E10   In addition to creating cell specific models as described so far we also created an averaged model 344 by first normalizing each experimental trace to one reference trace (so that each trace was given 345 equal weight in the averaging regardless of the conductance of the channel) and then summing and 346 averaging the current value at each time point along the protocol.

347
The parameter values obtained when calibrating each cell-specific and averaged model are shown 348 in Table F11. These values correspond to the parameter sets with maximum posterior density 349 identified in the MCMC chain. The full posterior density distributions for each parameter for each 350 of the 9 cells are shown in Figure F11. Table F11: Table of parameter values at the maximum posterior density for each cell-specific model, and the model fitted to averaged data (N.B. not the average of the cell-specific parameters). Here the model parameter numberings correspond to those detailed in Figure 4B, and G Kr represents the conductance value fitted for each model. *Note that the conductance fitted for the 'Averaged' model reflects mainly the conductance for the reference experimental trace (used for scaling all other traces before averaging), and should not be considered the 'average' conductance, hence its omission from Figure 7A.  Table F12 with a comparison between the experimental result and the 357 average model predictions with the cell-specific predictions. 358 We note that we have ordered the cells in this table (as in Figure 7) according to the percentage 359 change in leak resistance between performing the vehicle and dofetilide repeats of the sine wave 360 voltage protocol used to construct the model. This ordering acts as an estimated ranking for the 361 quality of each recording. The benefit of a cell-specific approach occurs when using the highest 362 quality data for both model construction and validation. We should note that even though in cells 363 #4 and #6 the average model provides the better prediction of the steady-state activation peak 364 current-voltage relationship than the cell-specific model, the cell-specific models are still providing 365 very good predictions in these cases, it is just that the experimental behavior is more like the average 366 model behavior for these cells. We also note that for eight out of the nine cells, the cell-specific 367 model provides a better prediction of the current response to the action potential protocol than the 368 average model, however, in the case where the cell-specific model is worse the difference is only a 369 small amount.  Figure F11: Distributions for each parameter for each of the 9 cell-specific models and the averaged data model. To aid comparison these are all histograms with 100 bars (plotting probability distributions here leads to very different maxima, obscuring the spread information), and so the y-axis is in arbitrary units related to the number of samples. We see that the parameter values tend to be given distinct distributions and so we would consider most of them to be 'significantly different', indicating that the variation we see in Figure 7 is due to cell-cell variability in the recordings rather than noise or unidentifiability in our parameter estimates.  Table F12: Table showing the error measure defined by equation F.10 between cell-specific or average models and the experimental current recording for fit (sine wave Pr7) and predictions with validation protocols (all other columns). Cells are ordered in ascending order according to the percentage change in leak resistance R leak . Here the color scale is set so that within each pair of columns T represents lowest error and T represents the highest error for each protocol/pair of columns. Note that the cells with larger currents will show larger errors, but the left column cell-specific predictions tend to perform better than the average model, particularly for cells where the average model gives a relatively large error. For predictions of the action potential protocol currents, Table F12 demonstrates that the cell-375 specific modeling approach yields predictions that are very close to or better than the average model.

376
Additionally, for the predictions of the steady-state activation protocol the cell-specific approach 377 generally yields very good and more accurate (for 4/5) predictions of validation data when the 378 highest quality data is used (cells #1-5). This benefit is absent when lower quality experimental 379 data is used where the average model provides very similar, but slightly better, predictions (cells 380 #6-9).

381
We also compare cell-specific and average predictions for each of the 9 cells for the deactivation, 382 recovery from inactivation and instantaneous inactivation time constants as were shown for one 383 cell in Figure 5. We show this comparison for each cell in Figure F12 and F13 for 8/9 cells and in 384 Figure F14 for all cells. Cell #6 was omitted in the first two plots because this cell had a particularly 385 low current and it was difficult to accurately fit exponential curves to the experimental data for this  Table F12. Voltage (mV) Cell 9 Figure F13: Cell-specific model predictions of time constant/voltage relationships for recovery from inactivation in Pr5. Each plot represents a different cell, with cell-specific model prediction depicted by the bold line, and the dashed line showing the cell's experimental data. Black lines on each plot represents the average model prediction. Cells are ordered as in Table F12.  Table F12.
to better quality data.

G Synthetic data study for an I Ks model 395
In order to demonstrate the wider applicability of our approach we also show its suitability for 396 parameterising a model of a different ion current -the slowly activating delayed rectifier potassium 397 channel (I Ks ). We tested the approach using a synthetic data study similar to that described in 398 Section C2 for the model in the main text. The I Ks model we used was taken from the Ten , (G.12) In Table G13 we compare the parameters we used to generate the synthetic data and the 403 parameters we identified as producing the best fit to these synthetic data (which included added 404 noise). We see that the maximum likelihood parameters are very similar to those from which 405 the synthetic data were produced, demonstrating the theoretical capability of this approach to 406 parameterise a model of I Ks . The practicality of this approach with real data is still to be explored.  (1997), as shown in Figure H15B. This model has 14 kinetic parameters, 411 rather than the 8 in our Hodgkin-Huxley style model, we maintain the non-voltage-dependent tran-412 sition rates between states C 2 and C 1 (states and parameters labelled as depicted in Figure H15B)  Figure H15A shows the fit to the sinusoidal protocol Pr6 for Cell #5 data 415 (analogous to Figure 4 of the main text). 416 We see a very good fit to the training protocol with almost all parameters having narrow posterior 417 distributions, as we did with the Hodgkin-Huxley (HH) model. In Figure H16 we show the predictions of the parameterised five state model for the action 432 potential protocol Pr6, akin to that shown in Figure 6 in the main text. We see that the five state 433 model is quite predictive for this protocol, and has much less error than with its original parameters 434 (shown in purple in panel E and quantified as 0.0764 nA in Table D6  is slightly more than the model in the main text (0.0475 nA vs. 0.0471 nA respectively).

440
The larger structure also makes slightly larger errors when predicting the standard voltage-441 step protocols, shown in Figure H17 below. We draw attention to the activation, inactivation and 442 recovery from inactivation summary curves, all of which are slightly worse under the five-state model 443 than under the HH model (see Figure 5 of the main text).

444
Overall, despite an improvement in model fit (comparing Figure Figure H16: Validation prediction for the five-state model -the current in response to the action potential protocol, for comparison with Figure 6 of the main text.