Accounting for variability in ion current recordings using a mathematical model of artefacts in voltage-clamp experiments

Mathematical models of ion channels, which constitute indispensable components of action potential models, are commonly constructed by fitting to whole-cell patch-clamp data. In a previous study, we fitted cell-specific models to hERG1a (Kv11.1) recordings simultaneously measured using an automated high-throughput system, and studied cell-cell variability by inspecting the resulting model parameters. However, the origin of the observed variability was not identified. Here, we study the source of variability by constructing a model that describes not just ion current dynamics, but the entire voltage-clamp experiment. The experimental artefact components of the model include: series resistance, membrane and pipette capacitance, voltage offsets, imperfect compensations made by the amplifier for these phenomena, and leak current. In this model, variability in the observations can be explained by either cell properties, measurement artefacts, or both. Remarkably, by assuming that variability arises exclusively from measurement artefacts, it is possible to explain a larger amount of the observed variability than when assuming cell-specific ion current kinetics. This assumption also leads to a smaller number of model parameters. This result suggests that most of the observed variability in patch-clamp data measured under the same conditions is caused by experimental artefacts, and hence can be compensated for in post-processing by using our model for the patch-clamp experiment. This study has implications for the question of the extent to which cell-cell variability in ion channel kinetics exists, and opens up routes for better correction of artefacts in patch-clamp data. This article is part of the theme issue ‘Uncertainty quantification in cardiac and cardiovascular modelling and simulation’.


Introduction
Mathematical modelling and computational simulations have been remarkably successful in providing mechanistic insight into many electrophysiological phenomena. Quantitative models of the action potential have demonstrated their usefulness in basic research and are beginning to be used in safety-critical applications [1][2][3]. Mathematical models of ion channels constitute indispensable components of these action potential models. Even when models are fitted to the best available data, uncertainty in their parameter values remains, which can be due to measurement uncertainty and/or physiological variability [4]. Thus, identifying and quantifying sources of uncertainty enables informed decision making when using models in safety-critical applications [5].
Whole-cell patch-clamp experiments (in voltage-clamp configuration) are a common source of data for calibrating ion channel models. To study the dynamics of ion channels, currents through the cell membrane are often measured with a patch-clamp amplifier. In voltage-clamp mode, a patch-clamp amplifier is a sensitive feedback amplifier that rapidly calculates, applies and reports the small currents necessary to maintain a given voltage across a cell's membrane (and vice versa for current-clamp mode) [6]. Typical peak current magnitudes are of the order of pA to μA, depending on the size and type of cell; the voltage across the cell membrane (potential inside minus outside) typically is within the range −140 to +60 mV.
We use the term kinetics to describe the voltage-dependent opening and closing of ion channels in response to changes in membrane voltage. Maximal conductance determines the magnitude of the current that would flow if all the channels were open, a value proportional to the number of channels in the membrane. High levels of variability in kinetics have been observed in many studies [7][8][9][10][11][12][13]. A previous study [14] analysed hERG current kinetics recorded simultaneously in 124 cells using an automated high-throughput patch-clamp machine. The experiments used Chinese hamster ovary (CHO) cells stably expressing hERG1a. Since the cells expressed the same genes and were measured at the same time under highly similar conditions, one might expect the resulting current kinetics to be very similar across different cells. However, a high level of variability was seen, similar to that in previous studies using manual patch-clamp experiments conducted one cell at a time over several days [15]. This raises a question-What is the origin of the observed variability? Figure 1a shows an idealized voltage-clamp experiment, where the cell is connected directly to an ammeter which records the current of interest, I ion , while clamping the membrane to the command voltage, V cmd (its equivalent circuit is shown in the electronic supplementary material, figure S1). In other words, a perfect patch-clamp experiment can be represented with  (1.1) This is the first of three ways of modelling the patch-clamp system that will be discussed in this manuscript (each enclosed in boxes as above).
In realistic voltage-clamp experiments, as illustrated in figure 1b, the cell membrane acts as a capacitor in parallel to the ion currents. Between the pipette electrode and the cell, there is a series resistance, which induces a mismatch between the membrane voltage V m and the command voltage V cmd , causing a shift in measured current-voltage relationships [16,17]. Furthermore, there is a voltage offset introduced between the pipette electrode and the cell [18]; and the wall of the pipette (or the well plate in an automatic system) behaves as a capacitor. Finally, a finite seal resistance can cause a substantial leak current that contaminates the recording of the current of interest. All of these can be partially compensated by the patch-clamp amplifier using real-time hardware adjustments, or addressed in post-processing, and the remainder are what we term experimental artefacts. Although real-world experiments suffer from these remaining artefacts, many studies assume the experimental apparatus perfectly compensates for any discrepancies, so the idealized assumptions (equations (1.1)-(1.2)) are frequently used when analysing experimental data. Often the existence of a leak current is acknowledged and an estimate for it is subtracted in a post-processing step to modify the I out in equation (1.2), as we will discuss later.
A previous microelectrode study [19] estimated and subtracted a transient cell membrane capacitive current post hoc using a mathematical model that included many of these components. In this paper, we relax further the typical, ideal voltage-clamp assumptions by introducing a mathematical model that allows and accounts for artefacts, imperfect amplifier compensations for artefacts, together with any residual uncompensated leak current. We then validate the mathematical model experimentally using electrical model cells, for which we designed a new type of electrical model cell which exhibits simple dynamic behaviour. Using this new mathematical model, variability in observations can be attributed to measurement artefacts as well as current properties. We develop a method to optimize the ion current maximal conductances, a set of current kinetic parameters, and measurement artefact parameters for more than one cell at the same time. Finally, we compare fits and assess the predictions of models calibrated with either ideal voltage-clamp assumptions or realistic voltage-clamp assumptions. While different cells have varying maximal current conductance, the study has implications for the significance or even existence of cell-cell variability in ion channel kinetics, and opens up routes for better correction of artefacts in patch-clamp data. Figure 2 presents an expanded circuit for the more realistic voltage-clamp experiment shown in figure 1b, including the amplifier components that compensate for the additional artefacts in the realistic experiment. Our goal is to observe the ion current across the cell membrane, I ion . This current is present in the 'Cell Model' in figure 2 (shown with black, dashed box). Between the Cell Model and the Headstage (green, dashed box) is where the pipette (or the well plate in an automatic system) sits, separating the cell membrane and the pipette electrode, which includes many of the undesired artefact components shown in figure 1b. There are five main undesired effects in this voltage-clamp set-up: (i) parasitic/pipette capacitance (a capacitance effect induced by the pipette wall), (ii) (cell) membrane capacitance, (iii) series resistance (a lumped term for all resistances between the pipette electrode and the cell), (iv) voltage offset (due to amplifier offsets, electrode offsets, junction potentials, etc.), and (v) leak current (a current leak through the pipette-cell seal and any leak through the membrane). Table 1 contains a glossary of symbols and parameters used for these quantities throughout this paper.

A detailed mathematical model of a voltage-clamp experiment
To model the dynamics of the current of interest, I ion , the assumptions behind equations (1.1)-(1.2) can be removed to instead model the entire voltage-clamp experiment and amplifier compensations with the equations below. Their derivation is given in the electronic supplementary material, §S2, where each of the undesired effects and how they are typically compensated is modelled based on published circuitry [17,18,[20][21][22][26][27][28][29][30][31].
ion channel model , l e a k c u r r e n t dV m dt In these equations, α is the requested proportion of series resistance compensation (a machine setting, typically 70-85%), and g leak and E leak are the conductance and the reversal potential of the leak current. The meaning of the remaining symbols is given in  Table 1. Glossary of symbols and parameters. We also denote a (machine or post-processing) estimate of a parameter X as X * , and the error in the estimate of the same parameter as X † . The range of typical values are taken from [18,20- [14]; 100-400 pF has been measured for individual cardiomyocytes in various species [23] (∼ 190 pF for human [24]); and can be as high as 100 nF for large cells such as Xenopus oocytes [25] Figure 2. A more realistic voltage-clamp experiment equivalent circuit. This includes undesired factors such as voltage offset (V off ), series resistance (R s ) between the pipette electrode and the cell, cell capacitance (C m ), pipette capacitance (C p ) and leakage current (I leak ), which can introduce artefacts to the recordings. The circuit also includes the components within a typical amplifier that are designed to compensate the artefacts. The blue (a) and orange (b) components are two idealized multiplying digital-to-analogue converters that control the amount of compensation. We assume that these, and the transimpedance amplifier and differential amplifier (green), are ideal electrical components. Please refer to table 2 for a description of the symbols. (Online version in colour.) In the real cell data analysed later, a leak compensation was applied to the data itself, prior to model fitting. In whole-cell recordings, typically estimates (denoted by * ) of the leak parameters in equation (2.2) are made in post-processing by assuming it is the only current active during a 'leak step' between voltages where I ion ≈ 0. Finally, I out is adjusted by subtracting estimated leak and to give I post as a better approximation of I ion .

Validating the mathematical model with electrical model cell experiments
Before applying the voltage-clamp experiment model to currents from a real biological cell, we test the performance of our mathematical model with electrical model cells-circuits made of electrical hardware components that mimic real cells. Some of these model cells are commercially available, e.g. HEKA MC 10, HEKA TESC, Axon or Molecular Devices PATCH-1U, and used to test and calibrate patch-clamp amplifiers. We are interested in both the current readout of a voltage-clamp experiment and the membrane voltage that the cell experiences. Therefore, we designed a circuit that connects the model cell to two amplifiers, one in voltage-clamp mode and one in current-clamp mode, as shown in figure 3. This set-up can simultaneously perform the conventional voltage-clamp procedure on the model cell with one amplifier, while using the other in current-clamp mode (clamped to zero) to measure the clamped voltage at the terminal corresponding to the membrane via the currentclamp. Effectively, this set-up allows the membrane voltage V m of the electrical model cell to be recorded while performing the conventional voltage-clamp current measurement.  (a) Electrical model cell design In figure 3a,b is a circuit which is equivalent to commercially available model cells (the design was based on the HEKA TESC), when under 'whole-cell' mode. In this study, we call this a Type I Model Cell. It consists of a capacitor and a resistor in parallel to mimic the membrane capacitance, C m , and membrane resistance, R m . Unlike real ion channels, this simple electrical model cell lacks any current dynamics in the R m resistor representing ion currents. We therefore developed a new type of model cell, termed a Type II Model Cell, which exhibits simple current dynamics when stepping to different voltages. Figure 3 c shows the equivalent circuit of our Type II Model Cell. In addition to the usual C m and R m connected in parallel, this model cell has an extra component (R k in series with C k ) connected in parallel, to mimic the addition of another ion current with some kinetic properties. The time constant (τ k = R k C k ) for this extra component was chosen to be O(100) ms, which is of the same order of magnitude as I Kr dynamics. The dynamics of the Type II Model Cell allow us to test and understand the effects of series resistance, etc., and enable us to validate our mathematical voltage-clamp experiment model experimentally to check we have correctly modelled amplifier compensations.

(b) Validation of the mathematical model
The experimental recordings using the simultaneous voltage-clamp and current-clamp set-up are shown in figure 4 (solid lines), with Type I (figure 4a,c) and Type II (figure 4b,d) Model Cells. The measurements were performed with a holding potential of 0 mV.
We performed two sets of experiments, firstly with no amplifier compensation, and secondly with automatic amplifier compensation using a computer controlled amplifier (HEKA EPC 10 Double Plus) where amplifier settings could be set with high precision. Here, automatic adjustment of the compensation settings, including V off , C p , C m and R s , was performed using the HEKA Patchmaster software, which closely follows the compensation procedure operators perform by hand on many manual patch-clamp amplifiers (pp. 80-84 of the HEKA manual [32]).
For the simulations, parameters were set in equations (2.1)-(2.7) to correspond to each set of experiments: for no amplifier compensation {V † off , C * p , C * m , R * s , g leak } were set to zeros; while for the automatic amplifier compensation those parameters were set to the amplifier's estimates. Results are shown in figure 4 for part of the voltage clamp protocol to be used for I Kr measurements later in this study (current under the full protocol is shown in the electronic supplementary material, figure S4). Our model (dashed lines) is able to capture both the current and the membrane voltage very closely for all these experiments. Note the differences between the membrane voltage V m (grey) and the command voltage V cmd (orange/red). For example, the uncompensated case, due to the voltage drop across R s in the Type I Model Cell, V m shows a simple offset-like deviation that depends on the current magnitude. But in the Type II Model Cell V m exhibits nonlinear dynamics while V cmd does not. When the amplifier was actively compensating, the differences between V m and V cmd were successfully reduced. All of these details are captured by the mathematical model so we are confident that equations (2.1)-(2.7) are a good representation of the model cells, amplifier compensations and artefacts that occur when compensations are disabled or imperfect.

(c) Parameter inference without compensations
Next, we attempt to use only the uncompensated, raw voltage-clamp measurements (i.e. only I out in figure 5a, and V cmd ) to infer the underlying membrane voltage, V m , and the parameters of the model cells. We then compare the model V m predictions with the current-clamp measurements. Here our focus is on the Type II Model Cell, since this should be the more challenging of the two and more similar to a real ionic current (similar results for the Type I Model Cell are in the electronic supplementary material, §S5). To optimize the model parameters root-meansquared error (RMSE) between the simulated and recorded I out was minimized using a global  optimization algorithm [33]. All optimization was done with an open source Python package, PINTS [34], and simulations were performed in Myokit [35]. All codes and data are freely available at https://github.com/CardiacModelling/VoltageClampModel. Figure 5a shows the fitted model I out (bottom, orange dashed line) and its corresponding prediction of the membrane voltage, V m (top, red dashed line), compared against the experimental recordings (solid lines). Figure 5b further shows that the fitted model is able to predict measurements under an independent, unseen voltage-clamp protocol-a series of action potential-like waveforms (grey lines in the first panel). Note the excellent predictions for V m and its deviation from V cmd . The series of action potential-like waveforms were concatenated from those in a previous study [14], they are composed of linear ramps and steps to permit use on highthroughput machines that cannot clamp to arbitrary waveforms. Magnifications of the protocols are shown later in figure 7c-f ; they include a delayed afterdepolarization (DAD)-like protocol, an early afterdepolarization (EAD)-like protocol, and action-potential-like protocols with different beating frequencies.
Table 2 also shows a comparison of the values of the hardware components (typical manufacturing tolerances for these are ±1 to 2%) in figure 3, the amplifier's estimation using a standard test pulse, and the fitted values using the mathematical model. Our model-inferred values are much closer to the component labels than the amplifier estimates because the Type II Model Cell exhibits nonlinear dynamics (as do real cells); whereas the amplifier uses a simple square-wave test pulse and assumes a simple resistor-capacitor model cell (i.e. a Type-I Model Cell) when estimating the parameters. That is, there is a difference between the electrical model cell we attached and the circuit the amplifier is designed to compensate, thus leading to inaccurate estimation of some components. For example, even though we did not apply any voltage offset, V off , in the experiment, the amplifier incorrectly estimated an offset of −1.2 mV. Proceeding with this amplifier-estimated value would lead to a voltage offset artefact of V † off = −1.2 mV in all recordings. The idea is that using the full model could capture these effects and take them into account when calibrating current kinetic parameters.
Thus far, we considered a realistic voltage-clamp experiment and developed a detailed mathematical model for such a setting, where imperfect compensations made by the amplifier and imperfect leak current subtraction are included. We then validated this mathematical model via electrical model cell experiments, demonstrating that our model captures the main effects of the voltage-clamp artefacts and amplifier compensations. In the next part of the study, we apply our mathematical model to experimental CHO-hERG1a data recorded previously.

Application to variability in CHO-hERG1a patch-clamp data
After experimentally validating the mathematical model of the full voltage-clamp experiment with two electrical model cells, it is now applied to experimental data from real cells. Here, a high-throughput dataset from a previous publication is used [14]. The dataset contains 124 voltage-clamp recordings of the potassium current that flows through hERG (Kv11  measurements were performed on CHO cells stably transfected with hERG1a at 25 • C, using the Nanion SyncroPatch 384PE, a 384-well automated patch-clamp platform with temperature control. All the data were leak-corrected (see equation (2.9)), E-4031 subtracted, and passed a semi-automated quality control. The recordings, and parameter values derived from them, showed good agreement with earlier work using manual patch; for details on the methods please see Lei et al. [14].
In the previous study [14], all of the observed variability in the post-processed data was assumed to be due to biological variability in I Kr conductance and kinetics, an assumption which we term 'Hypothesis 1', and so 124 cell-specific variants/parameterizations of the I Kr model were created. Kinetic parameters fitted to cell-specific data using the staircase calibration protocol enabled very good predictions for eight independent validation protocols. However, covariance in the inferred parameters across cells led the authors to speculate about a voltage offset (V † off ) being responsible for much of the variability [14].
An alternative hypothesis is that all cells have the same I Kr kinetics (functional properties of the channel proteins), but different maximal conductances (expression levels of the proteins), and that the observed variability in currents is due to differences in the patch-clamp artefacts and compensations for each cell. We term this set of assumptions 'Hypothesis 2'. Figure 6 shows a schematic overview of the two hypotheses.
We will now introduce the I Kr model and discuss the extent to which patch-clamp artefacts and imperfect compensations can explain the variability in the biological recordings.   Figure 6. A schematic overview of the two hypotheses we explore. Hypothesis 1 assumes the collected data represent a perfect, idealized voltage-clamp experiment where any experimental artefacts have been perfectly compensated by the amplifier and leak subtraction, so all the observed variability is rooted in biological variability (varying conductance and model kinetics in every cell), we term these 'independent kinetics models' . Hypothesis 2 assumes that the observed variability is due to differences in the voltage-clamp experimental artefacts, and that all of the cells share identical ion channel kinetics (although the maximum conductance is allowed to vary across cells), we term these 'identical kinetics models' . (Online version in colour.) (a) A mathematical model of I Kr I Kr is represented with a model used in previous studies [14,15,36,37]. The current is described with two Hodgkin and Huxley-style gating variables (a for 'activation' and r for 'recovery' from inactivation) and a standard Ohmic expression, where g Kr is the maximal conductance. E K is the reversal potential (Nernst potential) for potassium ions, which can be calculated directly from concentrations either side of the membrane using the Nernst equation (equation 2 in Lei et al. [14]). The gates a and r are governed by the ordinary differential equations and

(b) Inference with the full voltage-clamp experiment model
Earlier we saw how the full voltage-clamp experiment model (equations (2.1-2.7)) had parameters that could be successfully inferred from data for the electrical model cells. We tested how this inference scheme performs when combined with the I Kr model of equation (4.1) in a synthetic data study. All of the parameters could be identified (model parameters g Kr and θ together with artefact parameters C m , C p , R s , V † off and g leak ) from a single cell recording (full results are shown in the electronic supplementary material, §S6). We then tested this approach with real experimental data. However, values of the series resistance, R s , were consistently estimated at the lower bound we had imposed, which we consider unrealistic; results and code for this analysis can be found in the online repository. This behaviour could be due to imperfections in the representation of I Kr by the model (model discrepancy [38]), as all parameters in the full voltage-clamp experiment model were successfully inferred in the cases of electrical model cells or synthetic I Kr data. We therefore propose a simplification of the voltage-clamp experiment model while capturing the principal causes of variability.

(c) A simplified voltage-clamp experiment model
Modelling the entire voltage-clamp machine may not be necessary, because the timescales of the components in the system span multiple orders of magnitude; ranging from the order of 0.1 μs (e.g. τ clamp ) to tens of ms (e.g. C-slow) or many seconds for ion channel gating (e.g. activation of I Kr at room temperature [39]).
For these I Kr investigations the two fastest processes, τ clamp and τ z , can be approximated as instantaneous responses. That is, equations (2.4) and (2.7) can be approximated as V p ≈ V clamp and I out ≈ I in , respectively. After analysing the local sensitivity of the voltage-clamp experiment model (electronic supplementary material, §S3), we found that the effects of V † off and imperfect R s compensation are most apparent in the observed current (on timescales relevant to I Kr ). As a result, we further assume that (i) τ sum , part of the fast amplifier processes, is instantaneous; (ii) the effects of C p and C m are negligible; and (iii) C * m , R * s ≈ C m , R s . Finally, the data were leak subtracted, where the leak current parameters (g * leak and E * leak ) were estimated by fitting equation (2.2) to current during the −120 to −80 mV leak-ramp at the beginning of the measurements, to yield zero current at holding potential [14]. We allow this leak subtraction to be imperfect by retaining a small residual leak current with parameters g † leak and E † leak . With these assumptions, equations (2.1)-(2.7) become ion channel model For all symbols refer to table 1, and I * leak is given by equation (2.8). The effective reversal potential of the residual leak current, E † leak , is chosen to be the holding potential (−80 mV) because the primary leak subtraction (fit of g * leak and E * leak ) ensured approximately zero current at holding potential. Therefore, we have only two voltage-clamp model parameters (V † off and g † leak ) to infer along with the ion current parameters (g Kr and θ). Note that the other parameters (α, C * m , and R * s ) can be obtained from the amplifier settings without performing any inference.
This simplified voltage-clamp experiment model was successfully applied to the model cell experiments (results in electronic supplementary material, §S7) and, like the full model, was able to correct imperfect R m estimates made by the amplifier.

(d) Optimization of model parameters
We now present the parameter optimization schemes used to test each hypothesis. (i) Hypothesis 1: Cell-specific kinetics with no artefacts Under Hypothesis 1, any experimental artefacts are perfectly compensated by the amplifier and leak subtraction in post-processing. Any variability in the kinetic parameters obtained in current models fitted to the resulting data would reflect underlying variability in the biological function of the channels. To test this hypothesis, the parameter inference scheme detailed in Lei et al. [14] was employed. In summary, a Bayesian inference scheme was used which resulted in very narrow posterior distributions. This scheme uses some parameter transforms so that the optimization algorithm [33] searches in log-transformed space for certain parameters [37]. Here, we look for the most likely parameter set under that scheme, which is identical to that given by a least squareerror fit. So the log-likelihood of a given parameter set for cell i is proportional to Under this hypothesis, I model i is a function of just conductance g i and kinetic parameters θ i , and is given by equation (4.1) while assuming V m = V cmd . So we performed an optimization to maximize L i by finding g i and θ i for each cell i independently. The resulting models under this hypothesis are termed independent kinetics models.
(ii) Hypothesis 2: Identical kinetics for all cells, with cell-specific artefacts Hypothesis 2 assumes that the observed variability is due to the imperfect voltage-clamp experiments. Under this assumption, models fitted to the data should have the same kinetic parameters, θ * , across all cells, that is, θ i = θ * for any cell i. But there will be a cell-specific I Kr conductance, g i , and different patch-clamp experiment parameters for each cell too, φ i = {V † off,i , g † leak,i }. These models are termed identical kinetics models. To impose the assumption that all N cells have the same kinetics, and that the observed variability arises only from the experimental artefacts, the log-likelihood becomes where L i , still defined by equation (4.13), is the log-likelihood for the ith cell. But under this hypothesis, I model i is the post-processed current I post in equation (4.12) and hence L i is a function of the artefact parameters φ i too.
Optimizing L is a high-dimensional optimization problem, which is computationally expensive. This burden is reduced with a Gibbs-sampling style scheme; breaking the optimization problem into two: (i) optimizing the common kinetics parameters, θ; and (ii) optimizing the cellspecific parameters, {g i , φ i } i=1,...,N . To evaluate the maximum likelihood of θ, we nest optimization schemes. That is, for any single estimate of θ , we optimize {g i , φ i } for each cell i independently to compute an approximate likelihood for θ. The estimate of θ is then updated by running a single iteration of the outer optimization loop, followed again by optimization to convergence for {g i , φ i } in each cell. This overall scheme converges to the full set of optimal parameters: θ * , ..,N . The algorithm is detailed in the electronic supplementary material, §S8.
(e) Variability in ion channel kinetics or variability in patch-clamp artefacts?
The performance of the models arising from the two hypotheses is compared next. For Hypothesis 1 the results from Lei et al. [14] are used, where all the variability was assumed to be the result of kinetic variability (giving 9 parameters × 124 cells = 1116 parameters in total). For Hypothesis 2, we use the results from optimizing L(θ) in equation (4.14). A table of optimized parameter values, steady state-voltage relations, and time constants-voltage relations for the new identical kinetics models can be found in the electronic supplementary material, §S10.
As in Lei et al. [14], fits and predictions are quantified using relative root mean square error (RRMSE), defined as the root mean square error between the model simulation and the  Using this RRMSE quantification, the difference in the absolute size of the current across cells (due to varying conductance) is eliminated such that the scores are comparable between cells. Figure 7 shows for the independent kinetics models (red), and the corresponding raw current traces are shown in the three panels above. The model predictions for the same three cells with the identical kinetics model are shown in green. These protocols were used in the previous study [14] where they were intended to explore physiologically relevant predictions for I Kr behaviour. The same analysis applied to the remaining three protocols is shown in the electronic supplementary material, figure S7. The median, 10th percentile and 90th percentile RRMSE values for each protocol are shown in the electronic supplementary material, table S3.
For the calibration protocol shown in figure 7a, the RRMSE histogram of independent kinetics models (Hypothesis 1, shown in red) clearly shows lower errors (and hence better fits) than the identical kinetics models (Hypothesis 2, in green). This is expected, as under the independent kinetics Hypothesis, nine parameters can be varied per cell (leading to 9 × 124 = 1116 degrees of freedom overall), while under the identical kinetics hypothesis there are only three (for 8 + 3 × 124 = 380 degrees of freedom). This increased freedom allows Hypothesis 1 to achieve better fits, but could also lead to overfitting [40].
In overfitting, an increased number of degrees-of-freedom leads to improved fits to calibration data, but at the expense of a loss of predictive power on new datasets. Examples in figure 7b hint that this may be happening. Here, the 'median quality' Hypothesis 1 prediction (with nine cell-specific parameters) over-predicts the peak current under the Validation #3 protocol. By contrast, the Hypothesis 2 prediction (with three cell-specific parameters) is much closer to the data, suggesting that the cell-specific Hypothesis 1 parameter set may have varied kinetic parameters to fit details that are better explained by patch artefacts. To summarize such fits across all 124 recordings we can examine the histograms of RRMSE beneath each protocol in figure 7. Further evidence that overfitting might be occurring under Hypothesis 1 is shown by figure 7cf, where the red RRMSE histogram is shifted right, indicating that validation predictions show higher error for Hypothesis 1 than Hypothesis 2. In summary, Hypothesis 2 provides predictions in new situations that overall are at least as good as, if not better than, Hypothesis 1. Coupled with the greatly reduced number of degrees-of-freedom in Hypothesis 2 (736 fewer parameters), this is a strong indication that identical kinetics is the preferred hypothesis for these data.
In the identical kinetics models, all the variability must be explained by the voltageclamp artefact parameters φ i and the cell conductance g Kr,i in equation (4.14). Figure 8 shows the histograms and the pairwise scatter plots of the obtained {g * Kr,i , φ * i } i=1,...,124 . The artefact parameter values are within reasonable ranges (∼ ±5 mV for V † off and g † leak g Kr ) and are not strongly correlated, which lends further credibility to the identical kinetics models (Hypothesis 2).
Finally, one might expect that the mean of the kinetic parameters identified under Hypothesis 1 would be similar to the single set of kinetic parameters identified under Hypothesis 2. Kinetic parameters and model predictions from the two hypotheses are very similar, but not exactly the same (full results in electronic supplementary material, §S10). The difference is due to the experimental artefacts that are unaccounted for in Hypothesis 1 being included in the (model) kinetic parameters. If Hypothesis 2 is correct, the kinetic parameters derived using it are the more physiologically relevant.
We also examined a 'Hypothesis 3' where conductance, kinetics and artefacts were all cellspecific. Results are shown in electronic supplementary figure S10. Fits improved slightly (as expected with more parameters) but predictions were of similar quality to Hypothesis 1, and not as good as Hypothesis 2. Figures and code for this analysis can be found in the online repository.

Discussion
In this paper, we have introduced a new mathematical model for voltage-clamp experiments that allows and accounts for experimental artefacts and imperfect compensations of these artefacts, as well as imperfections in leak current subtraction. The model gives insight into the process of how experimental artefacts are introduced and the performance of the compensation procedure. For example, often the largest currents have the highest signal-to-noise ratio and could be interpreted  as the 'cleanest' recordings. But the model (equation (2.5)) shows that the effect of R s in distorting the applied voltage is proportional to the current size. So, for a given series resistance and a fixed percentage compensation (α), the bigger the current then the larger this artefact.
We validated the mathematical model through experiments using two types of electrical model cells, where we showed that our mathematical model is able to rectify imperfect amplifier estimations; then applied it to account for artefacts whilst inferring parameters of an I Kr model from hERG1a CHO cell data. This is, to our knowledge, the first time a detailed voltage-clamp experiment model has been used for parameter optimization in ion channel modelling.
Whole-cell patch-clamp data show a high degree of variability [14,15,36,39,41], with differences in observed current kinetics between experiments. As recordings measure wholecell 'macroscopic' current one might expect variation due to stochastic 'microscopic' ion channel opening. But here we have observed currents to be reproducible within a given cell [14,36]that is, the apparent kinetic variability upon protocol repeats in the same cell at different times is much smaller than cell-to-cell variability. We expect the maximum conductance of the current to differ between cells, due to varying cell sizes and gene expression levels. Since each CHO cell over-expresses the same channel gene, there is no obvious reason for the current kinetics to vary cell-to-cell, especially if the cells are from the same culture and recorded simultaneously under very similar conditions, as they were in our high-throughput data. Each voltage-clamp experiment has a different cell membrane capacitance, pipette (or wellplate in automated clamp) capacitance, voltage offset, series resistance and leak current. This led us to propose two competing hypotheses (as shown in figure 6): Hypothesis 1 is that there are cellspecific kinetics with no artefacts; and Hypothesis 2 is that there are identical kinetics for all cells with cell-specific artefacts. A parameter optimization technique was developed for Hypothesis 2, so that we were able to optimize all the model parameters at once (cell-specific conductances and measurement artefact parameters, with a single set of kinetic parameters shared by all the cells).
Models based on cell-specific conductances and kinetics (or cell-specific conductances, kinetics and artefacts) included many parameters/degrees-of-freedom. The simpler model of just cellspecific conductances and artefacts made better predictions, on average, over all 124 cells that were analysed. Considering Occam's razor, this makes the 'identical kinetics' hypothesis, in which fewer assumptions are made, the favourable explanation. So it is more plausible that the observed variability in the current traces arose from patch-clamp experimental artefacts, that is imperfect amplifier compensations and imperfect leak current subtraction, than from varying current kinetics. One might expect this result, since the current is carried by thousands of identical ion channel proteins, especially in our over-expression cell line.
It should be noted that these results and this interpretation are specific to our preparation, and cell-cell biological variability in kinetics in general cannot be ruled out. For instance, in native myocytes one might speculate that differing subunit expression and other signalling-related changes in proteins' states in the membrane (e.g. by PI3K signalling [42]) could also confer real biological variability in kinetics from cell to cell. But this study suggests that the major factor causing apparent variability in expression system current kinetics will be artefacts introduced by the patch-clamp experimental procedure, and offers an approach to address it.
Determining the origin of variability is particularly important for forward propagation of uncertainty in cardiac modelling. If each ion current biologically exhibited a high degree of variability in kinetics from cell to cell, then a description of this variability for all currents would be necessary to make predictions at the whole-cell level that accounted for it. However, if the majority of the variability is due to experimental artefacts, then a single set of kinetic parameters could accurately describe the physiology of each current type. If our findings also apply to myocytes, the standard approach to building action potential models with a single set of kinetic parameters for each current is appropriate [43][44][45][46], and the observed cell-cell variability in kinetics data is not required to be propagated forward in action potential simulations (as some studies have examined [47]). Nevertheless, we may still need an approach like the one demonstrated here to determine unbiased ion channel kinetic parameters to build the most physiologically relevant action potential models; as opposed to taking the mean of biased recordings which could accidentally include some experimental artefacts within the kinetics models.
Identifying ion current kinetics and separating out experimental artefacts is crucial for many cardiac electrophysiology studies. For example, ion channel mutation studies often conclude with a statement such as 'there is a 5-10 mV shift in the half-activation potential' [48,49]. However, given the variability that we observe in patch-clamp data, often the cell to cell variability in halfactivation potential can be in the range of 10-15 mV. Therefore, it is important to separate out experimental artefacts from real biological effects. The same principle applies to other cardiac electrophysiology studies, such as drug studies and the basic characterization of ion channel kinetics.
The findings raise the question of whether a mathematical model of the voltage-clamp experiment could remove the need for amplifier compensations altogether. This approach works well when we have an almost perfect model of the current, as in the electrical model cell cases (figure 5). But we do not recommend it in general, because unlike the reduced model, the full voltage-clamp experiment model (equations (2.1)-(2.7)) with an I Kr model appeared to produce spurious fits to real cell data. This behaviour is possibly due to I Kr model discrepancy-extra artefact components in the full model may attempt to 'mop up' real I Kr current which the approximations in this simple I Kr model do not explain. We envision that a more detailed I Kr model, or techniques for accounting for discrepancy [38], could allow us to infer parameters for the full patch-clamp model. Yet our simplified voltage-clamp model (equations (4.7)-(4.12)) appears to account for imperfect amplifier compensation even with an inevitably imperfect mathematical model of I Kr kinetics.
In this study, we have demonstrated that a reduced voltage-clamp experiment model can unify the kinetics of I Kr measured across many cells. Our model represents the voltage-clamp experimental set-up and could be applied to any voltage-clamp data gathered with a standard patch-clamp amplifier. Our full model of the voltage-clamp experiment could be particularly useful for currents like the fast sodium current (I Na ) where the time constants are similar to those in the capacitance artefacts.
Finally, it should be possible to generalize our method to current-clamp experiments by introducing an extra feedback circuit to the full voltage-clamp model for standard patch-clamp amplifiers [21,22]. In current-clamp mode, the leak current correction is usually performed during the measurement, therefore equation (4.12) needs to be incorporated into equation (2.6). However, note that conventional microelectrode amplifiers, such as the Axoclamp and the Axoprobe, have a different headstage design and a re-analysis of their circuits to build a similar model would need to be undertaken [22].

Conclusion
A mathematical model was derived to describe the entire voltage-clamp experiment including artefacts, imperfect amplifier compensations, and imperfect leak current subtraction. Using this model, variability in a set of experimental observations could be explained either by varying current properties or varying measurement artefacts. After comparing model predictions in each case, the results suggest that most of the observed variability in a set of expression system patchclamp data measured in high-throughput under the same conditions is caused by experimental artefacts. These varying experimental artefacts can be compensated in post-processing by fitting the mathematical model for the patch-clamp experiment at the same time as fitting ion channel kinetics. This study raises questions for the biological significance of any cell-cell variability in macroscopic ion channel kinetics, and provides for better correction of artefacts in patch-clamp data.