Reducing complexity and unidentifiability when modelling human atrial cells

Mathematical models of a cellular action potential (AP) in cardiac modelling have become increasingly complex, particularly in gating kinetics, which control the opening and closing of individual ion channel currents. As cardiac models advance towards use in personalized medicine to inform clinical decision-making, it is critical to understand the uncertainty hidden in parameter estimates from their calibration to experimental data. This study applies approximate Bayesian computation to re-calibrate the gating kinetics of four ion channels in two existing human atrial cell models to their original datasets, providing a measure of uncertainty and indication of potential issues with selecting a single unique value given the available experimental data. Two approaches are investigated to reduce the uncertainty present: re-calibrating the models to a more complete dataset and using a less complex formulation with fewer parameters to constrain. The re-calibrated models are inserted back into the full cell model to study the overall effect on the AP. The use of more complete datasets does not eliminate uncertainty present in parameter estimates. The less complex model, particularly for the fast sodium current, gave a better fit to experimental data alongside lower parameter uncertainty and improved computational speed. This article is part of the theme issue ‘Uncertainty quantification in cardiac and cardiovascular modelling and simulation’.

S1 Datasets and simulations S1.1 Data sources Pg. H233 [6] inactivation † , f f ∞ Fig. 2B [6] Pg. H230 [6] I to activation, r r ∞ Fig Table 1: Summary of patch clamp experimental datasets used in modelling papers in human atrial myocytes for each channel studied. Ticks and crosses are used to indicate which datasets are included in the original model calibration and which compose the unified dataset. * There are some differences in the terminology used in each model. The N model refers to the inactivation gates of I Na as h 1 and h 2 , and the I Kur channel as I sus . † The L-type calcium current has a calcium-dependent inactivation process which is not calibrated in this experiment (discussed further in results). ‡ In some cases it was not clear from the modelling paper where the comparison data plotted were obtained from. In these cases, the data points from the modelling paper itself are used and a protocol assumed based on the experimental paper cited. § In some cases it was not explicit which figure among a number of possibilities within a cited data source was used; the choice was inferred from the modelling paper.
2 S1.2 Temperature adjustment The N model was created to simulate a human atrial action potential at 306.15K (33 • C), whereas the C model was created to simulate at 310K (∼37 • C). Time constant measurements of rise and decay rates are temperature-dependent, and thus it was important to account for this during the calibration. During ABC, time constant measurements from the experimental sources were adjusted to the temperature of the model being calibrated using a Q10 factor from an experimental source. The Q10 factors used are: I Na : 2.79 [11]; I CaL : 1.7 (activation), 1.3 (inactivation), both calculated from values in [6]; I {to,Kur} : 2.2 [9]. In the original publication for the C model, I to and I Kur were fitted at room temperature (295.15K) and then adjusted by the authors to the model temperature (310K) by dividing time constants by a factor of 3 [2]. To maintain consistency, we also calibrate the C model at room temperature rather than make any adjustment to the experimental data. The time constants for this channel are adjusted before being used to simulate a full action potential. The S model was always calibrated by adjusting experimental data to 310K. Figure 1 shows the temperature-adjusted datasets for all experiments across channel models. When comparing channel models on the same figure throughout this work, their time constants are adjusted for the same temperature of 310K. For full action potential simulations, the N and C model time constants are kept at their model temperature, and only the S model time constant adjusted from 310K to 306.15K when it is added to the N model.

S1.3 Voltage-clamp protocols and summary statistic functions
Throughout this section, the lettering in the headers refers to the figure describing voltage protocols for the channel. All curve fitting for finding time constants was carried out using the scipy.optimize Python library.

S1.3.1 I Na
All experiments by Sakakibara et al. [3] used equal extracellular and pipette solution sodium concentration of 5mM at a temperature of 290K (17 • C). The time constant of activation experiment from Schneider et al. [4] used sodium concentration of 120mM in the extracellular solution and 70mM in the pipette solution, and the experiment was conducted at 297K (24 • C). A: Steady-state activation (p.538 [3]). The standard protocol to measure steady-state activation of the channel holds the membrane at a sub-threshold potential of −140mV and then steps the channel to a series of voltage steps between −100mV and 20mV, with intervals of 10mV. The steps last for 1s and there are 10s between each step.
The degree of steady-state activation is measured by recording the peak current during the voltage step (normalised to cell capacitance). The conductance is calculated by dividing the peak current by the forcing term, usually assumed to be the potential difference to the Nernst potential of the primary ion carrier, g = I V −EX . In the computational protocols, we directly measure the conductance of the channel to bypass this calculation. To plot the activation curve, the conductance is plotted against the voltage step normalised to its maximum value in any voltage step.
B: Time constant of activation (p.85 [4]). The time constant constant of inactivation is measured by fitting an equation to the current trace from a standard steady-state activation protocol as described in the previous protocol. In this case, the holding potential is −135mV and the steps are from −65mV and 15mV in steps of 10mV. The time at holding potential between each step was not given and so assumed to be 10s (more than enough for the I Na channel to return to steady-state) and each test pulse lasted for 12ms.
The activation time constant was measured by fitting the entire current trace at a pulse to the equation I Na = I Na,max 1 − e −t/τm 3 e −t/τ h + constant (p.87 [4]). In this equation, t is the time in ms, I Na,max is the peak of the current trace, τ m is the activation time constant and τ h is the inactivation time constant.
C: Steady-state inactivation (p.541 [3]). The protocol used to measure steady-state inactivation is often also referred to as an availability protocol. The membrane is held at a holding potential of −140mV for 10s, then stepped to a conditioning potential for 1s to activate and inactivate the channel. The membrane potential is returned to the holding potential for 2ms before stepping to a test pulse at −20mV for 30ms. A series of different conditioning pulses between −140mV and −40mV in steps of 10mV is used to test the amount of inactivation of the channel at different voltages.
The steady-state inactivation is measured by recording the peak current during the test pulse, normalised to the current in a test pulse when no conditioning step is applied (usually the maximum current amplitude). In the virtual voltage clamp experiment, we measure conductance directly in this step (which is equivalent as the forcing is the same during each test pulse and eliminated during the normalisation).
D: Fast/slow inactivation time constant (p.539 [3]). The protocol to determine inactivation time constants is a simple steptrain of test pulses from a holding potential of −140mV for 10s to a series of 100ms test pulses to voltages from −50mV to −20mV in steps of 10mV.
The fast and slow inactivation constants are determined by fitting the decay part of the current trace (after the peak current) to the equation I Na = A 1 e −t/τ f +A 2 e −t/τs +A 0 (p.538 [3]) where A are amplitude variables, t is the time and τ f and τ s are the fast and slow time constants of inactivation respectively. E: Fast/slow recovery time constant (p.538 [3]). The time constants of recovery from inactivation were determined using a double pulse protocol. The first pulse is a conditioning pulse for 1000ms to −20mV followed by a recovery period to a holding potential between −140mV and −90mV in steps of 10mV, of varying length between 2 − 1000ms (this was not specified and we assumed a series of t r = 2 i where i = 1, 2...10). The recovery period is followed by a test pulse identical to the conditioning pulse. We assumed 10s at holding potential between each pair of pulses.
The recovery time constant is measured through a series of processing steps. Firstly, the peak current in the test pulse is normalised to the peak current in the preceding conditioning pulse to give a measurement of proportion of the channel recovery for each recovery time period. This recovery measure is plotted against the recovery time period and the resulting curve is fit to a double exponential equation where r is the proportion of recovery, A are amplitude parameters, t 5 is the recovery time period and τ r(f) and τ r(s) are the fast and slow recovery time constants respectively. These values are calculated for each holding potential.

S1.3.2 I CaL
Experiments by Mewes and Ravens [5] were carried out at room temperature (assumed by authors to be 295K) and external solution with calcium concentration of 1.8mM. Those by Li and Nattel [6] were conducted at 309K with 2.0mM external calcium concentration. Experiments by Sun et al. [7] were completed at room temperature (assumed by authors to be 296K) and external calcium concentration of 1mM. It is difficult to estimate the level of intracellular calcium concentration during these experiments due to the mechanisms of the intracellular calcium stores and buffering. In our virtual voltage clamp experiments, the level of the intracellular calcium was kept constant at the resting value from the published N and C models (72.5nM and 101.3nM respectively) and, for the S model, set to the same value as the C model (101.3nM). A, B: Steady-state activation (p.1309 [5], p.H228 [6]). In Mewes and Ravens [5], steady-state activation was assessed using a conventional steptrain protocol from a holding potential of −40mV to 450ms steps between −35mV to 15mV in intervals of 5mV, with 10s between each pulse. In Li and Nattel [6], the activation curve was generated from the IV curve dataset which used a similar step train protocol. This time the holding potential was −80mV and the 300ms test pulses ranged from −80mV to 20mV in steps of 10mV. Both activation curves were generated as described in I Na steady-state activation.
B: Activation time constant (p.H233 [6]). The single activation time constant value was determined from the current trace evoked during the 10mV test pulse in the activation protocol from [6]. The activation time constant was estimated by fitting the upstroke of the normalised current trace to I CaL = 1 − Ae −t/τa where A is an amplitude parameter, t is the time and τ a is the activation time constant.
C: Steady-state voltage-dependent inactivation (p.H229 [6]). The voltage-dependent steadystate inactivation was assessed using a standard availability protocol. Conditioning pulses to a series of voltages between −80mV and 50mV in steps of 10mV were followed immediately by a test pulse to 10mV for 300ms. In [6], there are three datasets using different length of conditioning pulses of either 150ms, 300ms or 1000ms. We use the 1000ms prepulse dataset as this was used for calibration in both modelling papers. The inactivation curve was calculated as in I Na steady-state inactivation.
D, E: Fast/slow voltage-dependent inactivation time constant (p.H229-H230 [6], p.H1628 [7]). Fast and slow inactivation time constants were calculated by fitting a biexponential equation to the decay portion of the current trace during test pulses from a holding potential. In [6], three different holding potentials were used and we use the same as the modelling papers: −80mV. Test pulses lasted 300ms to steps between −10mV and 30mV in intervals of 10mV (we assume 10s between pulses). In [7], the holding potential was −80mV, test pulses had duration 1000ms with the same levels as above, and were preceded by a 500ms pulse to −40mV. In both cases, time constants of inactivation were calculated by fitting the decay portion of the current trace during the test pulses to I CaL = A 0 + A f e −t/τ f + A s e −t/τs were A are amplitude parameters, t is the time and τ f and τ s are the fast and slow time constants respectively.
F: Fast/slow voltage-dependent recovery time constant (p.H229-H230 [6]). Recovery time constants were assessed with a two-pulse protocol as in I Na recovery experiments. In this case, the holding potentials were −80mV, −60mV and −40mV with conditioning and test pulses both to 10mV for 300ms. Recovery periods were generated from t r = 2 i where i = 1, 2, ..., 11 based on the range of data points in Fig 4B [6]. The data was processed as in I Na recovery experiments with the exception that a single exponential function was used to fit the −80mV recovery curve to give a single (slow) recovery time constant. A: steady-state activation [5], B: steady-state and time constant of activation [6], C:steady-state inactivation [6], D: fast and slow time constants of inactivation [6], E: fast and slow time constants of inactivation [7], F: fast and slow time constants of recovery from inactivation [6]. The recovery protocol is repeated at multiple holding potentials (only one shown). See text for details of how the current traces are processed into summary statistics.

S1.3.3 I to
Experiments by Shibata et al. [8] and Wang et al. [9] were conducted at room temperature (assumed to be 295K). Firek  Some data for I to was extracted directly from the N and C modelling papers [1,2] as the source of the comparison experimental data plotted was not clear from the text. In these cases, conditions were assumed to be the same as the lab which produced the model. Thus, conditions for the assumed experimental data from Nygren et al. [1] was set to the same conditions as Firek and Giles [10] above, and data from Courtemanche et al. [2] was assumed to have been collected at the conditions in Wang et al. [9].
A, B: Steady-state activation (p. H1776 [8], p.1065 [9]). In Shibata et al. [8], steady-state activation was determined by holding at a potential of −60mV for 20s, then depolarising for 15ms to a step between −30mV to 80mV, and finally stepping to −40mV for 100ms. The activation curve was generated by recording the peak current amplitude in the final step (the 'tail' current). Wang et al. [9] used a standard steady-state activation protocol with a holding potential of −80mV to a series of test potentials between −40mV and 50mV and measured the peak current during the 1000ms depolarising step (before processing as above for I Na ).
C: Activation time constant (p.H305 [2]). Activation time constant data corresponding to that plotted in the modelling paper could not be found in the cited experimental source [9]. A simple steptrain protocol as described in [9] was assumed with a holding potential of −50mV and 100ms test pulses. The activaton time course of I to was fitted to a single exponential equation: I to = A 0 − Ae −t/τa where A parameters are amplitudes, t is the time course and τ a is the activation time constant. D: Deactivation time constant (p.H305 [2]). The deactivation time constant data source is also uncertain and it is assumed to use a similar protocol as in Fig. 9A in Wang et al. [9]. This protocol holds the membrane potential at −50mV for 20s before applying a 10ms conditioning pulse to 50mV. This is followed by a test pulse to elicit a tail current, which is fit to a single exponential (same as the activation time constant above) to determine the deactivation time constant.
E, F: Steady-state inactivation (p.34 [10], p.1065 [9]). In Firek and Giles, a standard steady-state availability protocol was applied as described in more detail in I Na steady-state inactivation. The holding potential is −80mV, followed by a 400ms conditioning pulse to levels between −80mV and 16mV, and finally a 400ms test pulse to 0mV to activate the outward current. Wang et al. [9] similarly uses a standard availability protocol with the same holding potential followed by 1000ms conditioning pulse to a range of voltages between −90mV and 30mV followed by a 1000ms test pulse to 60mV. Output was processed into summary statistics as described in I Na steady-state inactivation.
G, H: Inactivation time constant (p.66 [1], p.H305 [2]). It was not clear in either N or C model where experimental comparison data of time constants of inactivation for I to were obtained from. Consequently, simple protocols were assumed based on single-pulse protocols in experimental papers originating from the same labs. For [1], a single pulse protocol from a holding potential of −80mV to 400ms test pulses between 0mV and 40mV was applied and the decay phase of the current trace fit to a single exponential equation as above for activation time constants. For [2], a similar protocol was assumed based on [9] with a holding potential of −50mV and 100ms test pulses to between −40mV and 50mV.
I, J: Recovery time constant (p.66 [1], p.H305 [2]). Similarly to the inactivation time constant data, it was unclear where the recovery data was obtained from for both N and C models. For the N model, we assume the recovery protocol in [8] was used. This is a standard two-pulse recovery protocol with holding potentials of −100mV, −80mV and −60mv and 100ms test pulses to 20mV. For the C model, we assume the protocol in [9] was used with holding potential between −60mV and 40mV and 200ms test pulses to 50mV.

S1.3.4 I Kur
Experiments by Wang et al. [9] and Firek and Giles [10] use the same conditions stated above for I to .
As before, in some cases it was unclear where the comparison data in the modelling papers was obtained from. Details are given in the following sections.
A: Steady-state activation, activation time constant (p.1069 [9]). I Kur was measured from a holding potential of −50mV followed by a 1s prepulse to 50mV (experimentally used to inactivate the I to current which would otherwise interfere with isolating I Kur ). The potential is returned to −50mV for 20ms before being stepped to a range of 100ms test pulses between −40mV to 50mV each followed by a repolarising pulse to −10mV. I Kur was measured as the peak current in the final repolarising pulse and the activation curve determined as previously described for I Na steady-state activation. The activation time constants were determined by fitting the time course of I Kur trace during the test pulses to a single exponential function as described above.
B: Deactivation time constant (p.305 [2]). As it was not clear where the data points were obtained from, the protocol from Fig. 9 in [9] was assumed to have been used. This is the same protocol as described in deactivation time constant for I to .
C, D: Steady-state inactivation (p.34 [10], p.1068 [9]). For Firek and Giles [10], the protocol was the same as the steady-state inactivation protocol described for steady-state inactivation of I to with an increase of the length of the conditioning pulse from 400ms to 2500ms. I Kur was measured as the steady-state current at the end of the test pulse. In Wang et al. [9], steady-state inactivation is measured from the steady-state current at the end of a 2000ms test pulse to 40mV after a 1000ms conditioning pulse from a range of voltages. The holding potential was −60mV.
E, F: Inactivation time constant (p.66 [1], p.H305 [2]). It was unclear where the experimental data points from the N model paper came from. We assume they are generated from the protocol used in [10] from the same lab. This protocol is a simple 400ms step to a range of potentials from a holding potential of −80mV. In [10], the current decay is fit to a double exponential equation in order to separate I to and I Kur decay rates. As in the virtual voltage clamp there is only I Kur , we fit the current decay to a single exponential.
Similary for the C model, we assume a similar protocol as in [9] was used. This is a simple 2000ms step to a range of test potentials from a holding potential of −50mV. The decay portion of the current trace is fit to a single exponential function as above.
G: Recovery time constant (p.66 [1]). As previously for I to , we assume that the recovery protocol from [8] was used to determine a recovery time constant for I Kur . The protocol is as described for I to recovery time constant.

S2 Approximate Bayesian computation
Formally, ABC approximates the true parameter posterior distribution P (λ|D) by P (λ|ρ(D, D) ≤ ) where D are the experimental data, ρ is the chosen distance function,D is the model summary statistics and is the threshold value. We use the Toni ABC sampler based on sequential Monte Carlo to infer our parameter posterior distributions [12]. In this sampler, the ABC process above is repeated through a number of iterations with reducing . A population of parameter samples, referred to as 'particles', are propagated through each iteration and represent a discrete surrogate to the continuous posterior distribution.
At each iteration of the algorithm, the previous population of particles are perturbed slightly, by a multivariate Gaussian kernel, and used as the prior distribution for the current iteration. is reduced over iterations and chosen as the median distance of samples from the previous iteration. The particle population number is set according to the number of parameters being constrained in the experiment by considering the size of the parameter sampling hyperspace and assuming at least two particles in each dimension. A limit of 10000 particles is enforced due to computational demands. The algorithm terminates when less than 1% of parameter samples are accepted in a given iteration indicating the algorithm is struggling to improve on the current optimum. This criterion is chosen over termination at an absolute value of the distance metric termination used in other studies due to differences in number of model parameters and availability of data between experiments.
To compare our model summary statistics to the experimental data (which include error measurements at each point), we use a weighted Euclidean distance function where w i are the weights for each data point, σ i is the standard deviation associated with the error at each data point, δ is a regularisation factor applied when weights are close to zero, and n exp|i is the number of data points within the experiment. In this work δ was set to 0.05. Once generated, the weights are also mean-normalised to improve convergence of the ABC algorithm. The choice of distance function reduces the weighting of data points based on the magnitude of experimental uncertainty. This allows us to use information of which experimental data points we are relatively more confident about during the model calibration, propagating some of the experimental uncertainties through to this stage in the model development [13]. In some cases for I to and I Kur , it was not clear in the modelling paper where the data were obtained. In these cases, the points displayed in the figure from the modelling paper were digitised and 10% standard deviation error assumed (based on the error of similar measurements given in [9]).
Additionally, we account for the fact that we are simultaneously calibrating to multiple datasets by weighting according to the number of data points in a specific experiment. This provides balance between the different types of channel behaviour, rather than favouring an experiment with a greater number of data points. Each dataset is normalised to the maximum value in that experiment to avoid preference towards datasets with measurements at a larger scale in the ABC loss function.
A uniform prior distribution is used for each model parameter in the first iteration. The width of this prior is set to be sufficiently wide to cover the range of physiological possibilities. After convergence, the parameter posterior distributions are inspected and, if observed to be restricted by the lower or upper limit of the prior, the calibration is re-started with a wider prior. For the S model, prior ranges are set as in [14].      Table 3: Gating kinetics in Nygren model of I Na channel current(see Table 6 in [1] Table 7: Summary of results for parameters of standardised I Na model using unified dataset. * Parameters were searched in log 10 space and are presented in linear space. The prior for these parameters is still in the original log 10 space.    Table 13: Gating kinetics in Nygren model of I to channel current (see Table 8 in [1]