Sinusoidal voltage protocols for rapid characterization of ion channel kinetics

Rationale 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 topredict how ion currents contribute to cellular electrophysiology is in turn critically dependent on our characterization of ion channel kinetics — the voltage-dependent rates of transition between open, closed and inactivated channel states. Objective The hERG potassium channel plays a fundamental role in controlling electrical activity in the heart and many other tissues. There is a well-established link between hERG mutations, or block by pharmaceuticals, and increased arrhythmic risk. We present a new method for rapidly exploring and characterizing ion channel kinetics, applying it to the hERG channel as an example, with the aim of generating a quantitatively predictive representation of the ion current. Methods & Results We fit a mathematical model to currents evoked by a novel 8 secand sinusoidal voltage clamp in CHO cells over-expressing hERG1a. The model is 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 comprised of 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. Conclusions 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. The approach will be widely applicable to other voltage-gatedion currents both in the heart and other electrophysiological systems. Subject Terms Ion Channels/Membrane Transport; Electrophysiology; Computational Biology. Non-standard Abbreviations and Acronyms CHO — Chinese Hamster Ovary [cells]. GKr — maximal conductance of IKr. HEK — Human Embryonic Kidney [cells]. hERG — human Ether-a-go-go Related Gene. IKr — rapid delayed rectifying potassium current, carried by the Kv11.1 ion channel whose primary subunit is encoded by hERG.

1 Introduction pathological settings. The action potentials were shifted slightly so that their resting potentials were exactly −80 mV (see the Supplementary Code for full details and code to reproduce this protocol).

128
The protocol used to characterise the current and train the model is a voltage clamp comprised of 129 simple steps and a main sinusoidal section that is in the form of a sum of three sine waves of different 130 amplitudes and frequencies, designed to rapidly explore hERG channel kinetics. The underlying 131 rationale is to force the protocol to 'sweep' both the time and voltage dependence of the current 132 gating over physiological voltage ranges. 133 The start of the protocol takes the form of a leak step followed by a simple activation step which is 134 similar to Protocol 0. This activation step was included to improve the identifiability of the maximal  The main sinusoidal portion of the protocol takes the form of a sum of three sine waves as shown 140 in Equation (1): 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 = 142 0.19 ms −1 , and t is time measured in milliseconds.

143
In terms of frequencies, existing models and I Kr recordings include characteristic timescales of 144 order 10 ms to one second 26,27 . Therefore we designed the sinusoidal protocol's three frequencies 145 to probe channel kinetics across all these orders of magnitude (10 ms, 100 ms and 1 s timescales). 146 We selected frequencies that were co-prime rather than exactly multiples of ten: ω 1 to ω 3 are 147 ordered slow to fast and correspond approximately to sine waves of period 900, 170 and 33 ms, 148 respectively. The aim was that the three distinct frequencies should not become 'in phase': the 149 protocol never repeats patterns that the cell has experienced before (ensuring new information is 150 supplied throughout). The offset t 0 is 2500 ms as explained in Online Supplement B1.4. If one was 151 to study other ion channels, these frequencies may need adjustment to examine relevant timescales.

152
To decide the amplitudes, the oscillations are centered around −30 mV so that a physiological 153 range is explored (−120 < V < 60 mV). The amplitudes of the sine waves were selected to keep 154 the protocol within this range (A 1 + A 2 + A 3 = 90 mV) and to ensure that A 1 > A 2 > A 3 so that 155 the fastest timescale had the smallest oscillations (to avoid the faster gating processes masking the 156 voltage-dependence of slower ones).

168
We used the leak-step from −80 mV to −120 mV in order to leak-correct the experimental data, 169 according to: We identified the most appropriate R leak value to minimize the difference between the mean current 171 value during the leak step (to −120 mV) compared to the mean value at a holding potential of 172 −80 mV, whilst ensuring that the trace was not over-corrected (which would result in negative 173 currents during the initial stages of activation). 174 We manually selected leak resistances to correct the current evoked by the sinusoidal protocol 175 in both vehicle and dofetilide conditions. We then applied this leak resistance to the remaining 176 protocols performed in the same condition on each cell. In preliminary work, we observed that our sinusoidal protocols could elicit endogenous voltage-183 dependent background currents within expression-system cells. We observed that the levels of 184 endogenous currents the protocols elicited varied from cell to cell. These currents could adversely 185 affect the predictive ability of the resulting mathematical models, as the fitting process attempted 186 to create a model that described both the endogenous and I Kr components of the recorded currents.

187
To overcome this technical issue we made a number of alterations to our pilot experiments. 188 Firstly, we constrained the design of the sinusoidal protocol, as discussed above, so that only 189 voltages within a physiological range of −120 mV to +60 mV were explored, as endogenous currents 190 were much more prominent at voltages above +60 mV that we explored in pilot studies. 191 Secondly, we changed to use CHO cells in this study, rather than the HEK cells we used in pilot 192 studies, as CHO cells generally had lower endogenous currents. Whilst our model is equivalent to a two gate Hodgkin-Huxley formulation 20 , we use a Markov model description in practice (simply to generalize the computational code for other model structures).
The system of ordinary differential equations underlying the mathematical model structure shown in Figure 3B is then: where 211 [IC] = 1 − (C + O + I) .
The eight parameters P 1 to P 8 determine the rates k 1 to k 4 according to the exponential relationships 212 shown in Figure 3B. The current, I Kr , is modeled with a standard Ohmic expression: where G Kr is the maximal conductance, E K is the Nernst potential for potassium ions, and O is the 214 open probability, given by the solution to the system of equations above. E K is not inferred, but 215 is calculated directly from the ratio of ion concentrations on each side of the cell membrane using 216 the Nernst equation:

Note on Normalization
Where existing literature model simulations were plotted alongside experimental traces, or one 240 experimental trace was compared with another, we first had to normalize to account for differences 241 in conductance values. This was achieved by selecting a scaling factor for the conductance value for We calibrate a mathematical model using only data recorded under the sinusoidal protocol (Pr7).

252
The Hodgkin-Huxley 20 -style structure of the model we use, and its corresponding model parameters, 253 can be seen in Figure 3B. We independently fitted this model to each of the experimental current 254 traces shown in Figure 2. For each cell, we obtain a probability distribution of estimates for each 255 parameter that captures any observational uncertainty in the parameter values 31, 32 .

256
The result of the fitting procedure for one cell is shown in Figure 3. The parameter set with 257 maximum posterior density is shown in Figure 3A, demonstrating an excellent fit between experi-258 mental and simulated data. The resulting posterior probability density for the parameters obtained 259 from this Bayesian inference approach is projected across each parameter in Figure 3C. We also 260 tested that our approach is theoretically appropriate for inferring all parameters by using synthetic 261 data studies, as described in Online Supplement C. The plausible parameter space is very narrow:

262
if multiple parameter set samples are taken from the distribution shown in Figure 3C, the resulting 263 simulated current traces are indistinguishable to the eye. To quantify this, taking 1000 samples we 264 find that the 95% credible intervals for the simulated currents were always within at most either 265 3.47% or, in absolute terms, 0.0043 nA of the simulated current given by the maximum posterior 266 density parameter set.

267
The results we present in Figure 3 are from a single cell with a good quality recording and a high 268 signal:noise ratio (this choice of cell, and other cells' predictions, are discussed later). We fit models 269 on a cell-specific basis, and then also use averaged experimental data to create a single 'averaged' 270 model as described in Online Supplement F. We will compare these approaches below. We provide 271 all parameter values with the maximum posterior density for all models in Supplementary Table F11 using the experimental procedure shown in Figure 1B.

279
To make the predictions for Protocols Pr1-6 we performed simulations using the parameter 280 set with the maximum posterior density in the fit to the sinusoidal protocol (Pr7). As with the 281 calibration protocol, all the predictions we will discuss below are indistinguishable by eye from the 282 result of taking multiple samples from the distributions in Figure 3C and plotting a prediction for 283 each of these parameter sets. We also compare the predictions from our new model with those from 284 a sample of widely-used literature models 26,34-37 .

285
In Figure 4, we show traditional voltage step protocols, experimental recordings and the sim-286 ulated predictions from the model. We also show some of the most commonly plotted summary 287 curves for experimental data under these protocols, together with predicted summary curves from 288 both our model and a range of existing literature models (methods used to derive summary plots 289 are given in Methods B1.7, results for Pr1&2 in Online Supplement E). We can predict a wide 290 range of current behavior in response to the standard voltage-step protocols, without having used 291 any of this information to fit the model.

292
There are a number of points to draw attention to in Figure 4.  Figure 5 shows the model prediction of the currents invoked in response to the physiologically-306 inspired action potential protocol Pr6, compared with the experimental recording (as shown in 307 Figure 1B we used the first repeat of Pr6 for validation purposes, and the second as a quality 308 control measure). Replicating behavior under action potentials is perhaps the most important 309 requirement for a hERG channel model for use in physiological or pharmacological studies. The 310 model is able to predict the response to all of the complex action potential protocol extremely well, 311 and much better than existing models (even though we have scaled literature models' maximal 312 conductances (G Kr ) to fit this trace as well as possible in Figure 5). 313 We provide a quantitative comparison of predicted current traces for our model and each of In Figure 6A we present the maximum posterior density parameter values when repeating the above 322 approach using data from nine different cells. The clustered parameter values demonstrate that 323 parameters derived from different cells take similar values, giving us confidence that the procedure is 324 reproducible and biophysically meaningful. There is more cell-to-cell variability in some parameters 325 than others, which may be related to variability in the underlying physiological processes that 326 they represent; supporting the value, and perhaps necessity, of a cell-specific approach. We also 327 acknowledge that some parameters may be more or less sensitive to variability in experimental conditions such as temperature, residual background/endogenous currents, and imperfect dofetilide 329 and/or leak subtraction. 330 We order the cells in Figure 6 based on the lowest difference in leak resistance between the 331 vehicle and dofetilide recordings of Pr7. This ordering gives a measure of recording stability, and is 332 intended to be a surrogate for data quality. The cell presented above, in Figures 3-5, corresponds to 333 Cell #5 of 9 under this ranking, so we obtain very good predictions even with our 'median' quality 334 data. We show cell-specific predictions of the current-voltage relationship for the peak steady-state 335 activation current for each cell-specific model in Figure 6B. While we focused on Cell #5 in the 336 results section, Cells #1-4 also produce excellent cell-specific predictions (similar comparisons for 337 other summary plots are in Supplementary Figures F12-F14). 338 We also investigated the benefit of a cell-specific approach by building a model using averaged 339 experimental data from all nine cells. We describe this approach in Online Supplement F, and 340 summarize the results in Supplementary Table F12. Generally, for the cells with the highest data 341 quality (Cells #1-5) the cell-specific models provide better predictions than the average model,

342
as we see for Pr3 when comparing colored cell-specific predictions and experiment with the black 343 line for the average model in Figure 6B. The same trend holds for the action potential protocol  We saw that our model is able to replicate the experimental training data very well ( Figure 3).

361
This is often the point at which literature approaches stop, and conclude that a mathematical 362 model is a good representation of ion channel kinetics (something that is also true more generally  The extremely good prediction from all our cell-specific models of the response to the complex 374 action potential protocol is particularly remarkable ( Figure 5). Cell-to-cell variability in ion chan-375 nel kinetics was captured by fitting different underlying kinetic parameters. These parameter sets 376 were shown to have modest variation, and this variation in kinetics was quantitatively predictive 377 of variation observed in independent validation experiments ( Figure 6). Cell-specific predictions 378 were particularly strong when using the highest quality data, highlighting the necessary data qual-

383
There are still some aspects of the experimental behavior that are not replicated by our model.

384
In particular, there is only one time constant of deactivation, and low voltage-dependence in the   Figure A1 for the range of possibilities).

426
In summary, we have demonstrated significant advantages in our cell-specific mathematical The authors declare that the research was conducted in the absence of any commercial or financial 461 relationships that could be construed as a potential conflict of interest.        Table A1 in Online Supplement A). All models have been simulated with their original published parameters, with a temperature of 21.5℃ and our experimental concentrations used for calculating the reversal potential. B: top -a schematic representation of the experimental procedure used for this study over time (not to scale). A simple activation step protocol is repeated in the sections marked 'Pr0', before moving on to the highlighted section (enlarged below) where data used in the study were recorded. The recording protocols 'Pr1-7' are performed twice, once before dofetilide addition, and once after, with the hERG current isolated by subtraction. For full details of the protocols please refer to Methods B1.4. Simulations of expected behavior in response to this protocol from existing I Kr and hERG models, normalized by scaling the conductance value for each model to minimize the absolute difference between each trace and a reference trace. For calculation of the reversal potential, a temperature of 21.5°C was used to match the mean experimental conditions. Bottom row: Raw data (following leak and dofetilide subtraction) from experimental repeats at room temperature from 9 cells. Experimental traces have been scaled, to remove the effect of different maximal conductances, by a factor chosen to minimize the absolute differences between each trace and a reference experimental trace (that with the peak current during the sinusoidal portion of Pr7). B: an enlargement of the highlighted sections of panel A. Whilst there is some variation between cells in the experimental results, they are much more consistent than the predictions from the different models. Note that the deactivation time constant we plot here is a weighted tau, described in Methods B1.7.
Note that some literature model predictions are missing from the summary plots as we were either unable to fit exponential curves to 'flat' simulation output reliably; or the exponential decay occurred in the opposite direction to experimental traces, and we considered the comparison unwarranted. Conductance, G Kr , is scaled for each of the literature models to give the least square difference between their prediction and these experimental data, i.e. we display a best-case scaling for each of these models. A quantification of the error in our model prediction versus these literature models is given in Supplementary Table D6: the performance shown in panels D and E holds for the whole trace, so the mean error in predicted current across the whole protocol is between 69% and 264% larger for the literature models' predictions than for our sine-wave fitted model. Voltage (mV) Cell 9 P 1 P 2 P 3 P 4 P 5 P 6 P 7 P 8 G Kr  Figure 1A of the main text features simulations from 29 literature hERG or I Kr models. In Table A1 605 we list these models, give references, and show the seven different structures that they feature in 606 Figure A1.

607
There are two models in this table that are not in Figure 1A:  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 3B. 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.

B Additional Methods
This section contains further description of the methods that were used, with a particular focus on 612 details of the Bayesian Inference scheme in Section B2. These sections do not feature in the main 613 Methods section due to space constraints.  The protocols were first created as text files and then converted to .abf stimulus files to make 633 corresponding .pro protocol files in the pClamp 10 software. A CV 203BU amplifier headstage 634 and a Digidata 1440A were used. A Sutter MP225 micromanipulator was used for positioning of 635 the microelectrode. The current signal was sampled at a rate of 10 kHz. 75-80% series resistance 636 compensation was applied and data were 5 kHz low pass Bessel filtered by the hardware. No 637 software filtering was applied. Whole-cell capacitance compensation was applied electronically.

638
Leak subtraction was applied offline by using a 50 ms leak step to allow correction. To make a series 639 of successive recordings using different protocols on the same cell, the pClamp "Sequencing Keys" 640 tool was utilized, with a .sks file detailing the sequence the protocols should be performed in.

642
Here we list the details of the voltage-clamp protocols that were not included in the main text. The 643 protocols can also be found encoded in the software, available to download as described at the end 644 of the main text. repeat. From an initial holding potential of −80 mV, this protocol comprised a 5 s step to 10 mV followed by a 5 s step to −50 mV before returning again to a holding potential of −80 mV. This 650 protocol is depicted in Supplementary Figure B2. We repeated this protocol while dofetilide was 651 added (see Figure 1B) and the current traces recorded from this protocol were used to assess when 652 a steady level of dofetilide block had been reached.

654
After the initial period at holding potential incorporating the −120 mV leak step, a step to V step 1 655 followed and was held at that voltage for T step ms, before a step to −120 mV for 2.5 s, before 656 returning to holding potential of −80 mV for 1 second. The protocol was repeated 6 times with a 657 different T step on each repeat. T step took the values of 3, 10, 30, 100, 300 and 1000 ms.

663
From the initial period at holding potential incorporating the −120 mV leak step, a step to V step 664 was applied for 5 seconds, followed by a 1 s step to −40 mV, before a 500 ms step to −120 mV, and 665 then returning back to holding potential for one second. This process was repeated 7 times with a 666 different V step on each repeat. V step ranged from −60 mV to +60 mV in 20 mV increments. This 667 protocol is depicted in Figure 4A (left column).

669
From the initial period at holding potential incorporating the −120 mV leak step, a step to 50 mV 670 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 671 step to −120 mV, and a 1 s step back to holding potential of −80 mV; This was repeated 16 times 672 with a different V step on each repeat. V step ranged from −100 mV to 50 mV in 10 mV increments.

Protocol 5 -Deactivation
From the initial period at holding potential incorporating the −120 mV leak step, a step to 50 mV 676 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 677 returning back to holding potential for one second. This process was repeated 9 times with a 678 different V step on each repeat. V step ranged from −120 mV to −40 mV in 10 mV increments. This 679 protocol is depicted in Figure 4A (right column). The full protocol is comprised of 250 ms at holding potential of −80 mV, followed by a 50 ms 'leak 684 detection' step to −120 mV, and then 200 ms back at −80 mV. This was followed by an 'activation' 685 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.

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

688
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 = 689 0.19 ms −1 , and t is time measured in milliseconds.

690
The protocol was initially designed with just the −120 mV leak step and not the additional 691 'activation' steps to 40 mV and −120 mV (which were included after preliminary experiments as 692 described in Section B2.2) and so the sine wave was shifted using t 0 = 2500 ms to begin at the same 693 phase after we incorporated the additional steps. This offset is not expected to be important to 694 include, but was included here for clarity for anyone attempting to reproduce our study.    To produce the peak current-voltage relationship for the steady state activation protocol for the 733 simulated data traces we wrote MatLab code to identify the peak current in the region between 734 5.6292 and 5.7292 seconds on each sweep of the protocol, which corresponds to the current response 735 just after the 5 second V step when the voltage is stepped to −40 mV. We then normalized the peak 736 current data to the maximum overall peak identified in this region to produce the current-voltage 737 relationship curve. For the simulated data we wrote a MatLab script (included in code download) 738 to identify the peak-current voltage relationship for this protocol but for the experimental traces 739 we verified these peak points manually to avoid incorrect peaks being identified due to noise or 740 capacitive effects. We also identified the peak currents in the currents evoked by the activation 741 kinetics protocols manually for the same reason. In the activation kinetics protocol we identified 742 the peak currents during the V step for each interval of T step duration. We assume that the experimental noise is independently and normally distributed with a mean 770 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 772 by equation (7) In practice, the sums over time in equation (B.7) are formulated so that we exclude time points 778 from regions where the data are affected by capacitive spikes. To be precise, we exclude 5 ms 779 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 781 seconds and 7 seconds (spikes are seen in experimental recordings at these times in Figure 2).

B2.2 Conductance estimation to inform the prior 783
Preliminary work revealed that using sine wave protocols alone often allowed kinetic parameters 784 in the hERG model to be recovered, but there was potential for identifiability problems (or at 785 least we encountered difficulties in finding a global optimum due to a rugged likelihood surface) 786 when simultaneously fitting the conductance parameter and transition rate parameters P 1 to P 8 787 (although previous work suggests all parameters are theoretically identifiable 32 ). To add extra 788 information on conductance, we incorporated a voltage-step to +40 mV followed by a step down to 789 −120 mV, as described in the definition of the sinusoidal protocol above. The aim being to provoke In this section we describe our prior assumptions on the values that each model parameter can take.

799
The prior for the conductance G Kr is assumed to be independent of the kinetic parameters, and to 800 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 3).

805
For parameters of the form A we assumed that the prior distribution is uniform between 10 −7 806 and 1000 ms −1 , again to cover (and extend beyond) the full physiological range expected with hERG 807 channel gating. 808 We assume that the prior distributions for B parameters are uniform between 10 −7 and 0.4 809 mV −1 . The lower bound for this parameter was selected as the voltage-dependence becomes prac-810 tically redundant when B becomes small: when B = 10 −7 the value of exp(BV ) will change by less 811 than 0.0015% across the voltages we reach in this study. The upper value is beyond the physiolog-812 ically expected range. 813 We also impose a prior on the maximum rate of transition k between any states (maximum 814 across the full voltage range in the protocol (that is from −120 to 58.25 mV)). If the maximum 815 rate k is greater than 1000 ms −1 , or less than 1.67 × 10 −5 ms −1 , the pair of parameter values that 816 give rise to this are assigned prior probability zero (strictly, this is equivalent to defining 2D prior In order to verify that there was sufficient information within the sinusoidal protocol to parameterize 880 our model we performed a synthetic data study. The aim in such a study is to ascertain whether we 881 can recover the parameters used in the simulation from a simulated data trace (with added noise 882 in this case).

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

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

900
The synthetic data study provides us with confidence in the suitability of our protocol for 901 accurately identifying parameters of the model presented in Figure 3B in the main text, and also 902 that the parameter inference protocol(s) we are using are suitable for the task. We believe such an 903 approach should always be used to test whether there is sufficient information in the experimental 904 data being proposed for calibration of a mathematical model. The test should be performed twice: 905 before conducting the experiment (with the pre-existing best guess at the parameters); and also after 906 conducting the experiment (with the new maximum posterior density estimate of the parameters 907 -as we illustrate in Figure C7). 908 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.

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

937
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

F Comparing Cell-Specific with Average Model
In addition to creating cell specific models as described so far we also created an averaged model 947 by first normalizing each experimental trace to one reference trace (so that each trace was given 948 equal weight in the averaging regardless of the conductance of the channel) and then summing and 949 averaging the current value at each time point along the protocol.

950
The parameter values obtained when calibrating each cell-specific and averaged model are shown 951 in Table F11. These values correspond to the parameter sets with maximum posterior density 952 identified in the MCMC chain. The full posterior density distributions for each parameter for each 953 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 3B, 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 6A.  Table F12 with a comparison between the experimental result and the 960 average model predictions with the cell-specific predictions. 961 We note that we have ordered the cells in this table (as in Figure 6) according to the percentage 962 change in leak resistance between performing the vehicle and dofetilide repeats of the sine wave 963 voltage protocol used to construct the model. This ordering acts as an estimated ranking for the 964 quality of each recording. The benefit of a cell-specific approach occurs when using the highest 965 quality data for both model construction and validation. We should note that even though in cells 966 #4 and #6 the average model provides the better prediction of the steady-state activation peak 967 current-voltage relationship than the cell-specific model, the cell-specific models are still providing 968 very good predictions in these cases, it is just that the experimental behavior is more like the average 969 model behavior for these cells. We also note that for eight out of the nine cells, the cell-specific  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 6 is due to cell-cell variability in the recordings rather than noise or unidentifiability in our parameter estimates.
We use a measure of  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-978 specific modeling approach yields predictions that are very close to or better than the average model.

979
Additionally, for the predictions of the steady-state activation protocol the cell-specific approach 980 generally yields very good and more accurate (for 4/5) predictions of validation data when the 981 highest quality data is used (cells #1-5). This benefit is absent when lower quality experimental 982 data is used where the average model provides very similar, but slightly better, predictions (cells 983 #6-9). 984 We also compare cell-specific and average predictions for each of the 9 cells for the deactivation, 985 recovery from inactivation and instantaneous inactivation time constants as were shown for one 986 cell in Figure 4. We show this comparison for each cell in Figure F12 and F13 for 8/9 cells and in 987 Figure F14 for all cells. Cell #6 was omitted in the first two plots because this cell had a particularly 988 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.

997
G Wider applicability: synthetic data study for an I Ks model 998 In order to demonstrate the wider applicability of our approach we also show its suitability for 999 parameterising a model of a different ion current -the slowly activating delayed rectifier potassium 1000 channel (I Ks ). We tested the approach using a synthetic data study similar to that described 1001 in Section C2 for the model in the main text. The I Ks model we used was taken from the ten α Xs = 1 1 + P 6 exp(P 7 V ) , (G.13) β Xs = 1 P 5 1 + P 3 exp(−P 4 V ) , (G.14) τ Xs = α Xs β Xs , (G.15) dX s dt = Xs ∞ − X s τ Xs .
(G. 16) In Table G13 we compare the parameters we used to generate the synthetic data and the 1006 parameters we identified as producing the best fit to these synthetic data (which included added 1007 noise). We see that the maximum likelihood parameters are very similar to those from which 1008 the synthetic data were produced, demonstrating the theoretical capability of this approach to 1009 parameterise a model of I Ks . The practicality of this approach with real data is still to be explored.
1010 Table G13: Table comparing parameter values used to simulate the synthetic data for the I Ks model and the maximum likelihood fit to these data.