Quantifying the efficacy of voltage protocols in characterising ion channel kinetics: A novel information‐theoretic approach

Voltage‐clamp experiments are commonly utilised to characterise cellular ion channel kinetics. In these experiments, cells are stimulated using a known time‐varying voltage, referred to as the voltage protocol, and the resulting cellular response, typically in the form of current, is measured. Parameters of models that describe ion channel kinetics are then estimated by solving an inverse problem which aims to minimise the discrepancy between the predicted response of the model and the actual measured cell response. In this paper, a novel framework to evaluate the information content of voltage‐clamp protocols in relation to ion channel model parameters is presented. Additional quantitative information metrics that allow for comparisons among various voltage protocols are proposed. These metrics offer a foundation for future optimal design frameworks to devise novel, information‐rich protocols. The efficacy of the proposed framework is evidenced through the analysis of seven voltage protocols from the literature. By comparing known numerical results for inverse problems using these protocols with the information‐theoretic metrics, the proposed approach is validated. The essential steps of the framework are: (i) generate random samples of the parameters from chosen prior distributions; (ii) run the model to generate model output (current) for all samples; (iii) construct reduced‐dimensional representations of the time‐varying current output using proper orthogonal decomposition (POD); (iv) estimate information‐theoretic metrics such as mutual information, entropy equivalent variance, and conditional mutual information using non‐parametric methods; (v) interpret the metrics; for example, a higher mutual information between a parameter and the current output suggests the protocol yields greater information about that parameter, resulting in improved identifiability; and (vi) integrate the information‐theoretic metrics into a single quantitative criterion, encapsulating the protocol's efficacy in estimating model parameters.

(vi) integrate the information-theoretic metrics into a single quantitative criterion, encapsulating the protocol's efficacy in estimating model parameters.

K E Y W O R D S
experimental design, identifiability, information theory, ion channel kinetics, parameter estimation

| INTRODUCTION
2][3][4] By developing robust models for physiological processes, the effects of pharmacological agents can be predicted and evaluated for safety before proceeding with lengthier and more expensive trials.However, these models usually have a large number of unknown parameters, which are typically estimated through experimental measurements.
One physiological process prone to interference from pharmaceuticals is the hERG potassium ion channel. 5This channel carries the rapid delayed rectifier potassium current, I Kr , which, among other things, helps co-ordinate the beating of the heart.There is great interest in the modelling of I Kr channel dynamics within cardiac cells, as this channel is prone to blockade by certain pharmacological agents, which can lead to potentially life-endangering arrhythmic side-effects. 6It is worth mentioning that this is not a phenomenon exclusive to hERG; Na þ and Ca 2þ ion channels also play fundamental roles in generating and propagating the cardiac action potential, 7 and are also susceptible to blockade by pharmaceutical agents, either as a side-effect 8 or therapeutically to treat certain cardiac disorders. 9 Kr channel dynamics are voltage-dependent and time-dependent, and characterised by a number of parameters that determine the relationship between the channel current and the time-varying voltage signal applied to the channel.The experimental data with which to estimate these parameters are often gathered through voltage-clamp experiments, 10,11 wherein a cell is stimulated by a known voltage signal, and the resultant current is measured.This voltage signal is referred to as the voltage protocol.
The protocols typically last several minutes, and results in the gradual damage and destruction of the cell.This results in a decrease in data quality as the experiment progresses, and often times means that the experiment will have to be carried out across multiple cells.In Reference 12, six classically used voltage protocols are analysed.Additionally, a novel short voltage protocol is proposed, which primarily consists of a sum of three sinusoidal waves.This new protocol has the advantage of being far shorter than previously used voltage protocols, which improves data quality, allows models to be characterised for individual cells, and makes investigations into cell-to-cell variability practical.
The results of the investigation in Reference 12 concluded that when model parameters were estimated with the novel protocol, the model was able to reproduce the experimental measurements for not only the novel voltage protocol, but also the other six classical protocols.This demonstrated the superiority of the novel protocol.
However, the focus of the investigation was on estimating the parameter values based on a Bayesian inference scheme, rather than the design of the voltage protocol itself.The novel protocol was largely constructed through experience and broad understanding of how variations in voltage affect the measured current through the model parameters.This is to say that a systematic method for construction of novel protocols and comparing their efficacy quantitatively, while being computationally practical and highly desirable, is lacking a conclusive approach in literature, although a systematic method for designing and optimising voltage protocols using Bayesian inference was described in Reference 13.When faced with multiple options voltage protocols, such as during the optimal design of a novel protocol, the extension of the methodology in Reference 12 would be to carry out the entire process of solving the inverse problem for each new candidate protocol-which would be costly and time-consuming due to the experimental procedure and subsequent Markov Chain Monte Carlo (MCMC) analysis for Bayesian inference.
In this article, a new computational approach based on information theory is used to assess the ability of different voltage-clamp protocols in estimating the parameters of a model describing ion channel kinetics. 12Emphasis is also given on computational efficiency-this pipeline is easily parallelisable and is less computationally intensive than several other methods for electrophysiological simulations. 14,15This would allow candidate protocols to be assessed in terms of individual parameter identifiability without experimental data, only by using simulated data and estimates of key metrics central to information theory.Essentially, this framework would allow the experiment designer to compare the information contained by the candidate protocols for each parameter, and their aggregates.As a result, the candidate protocol which yields the highest information for the parameters can be given a higher preference for selection.Optimising over a space of candidate protocols would, in theory, yield a protocol with the maximum possible information gain for each parameter, resulting in an optimally designed experiment. 16he information-theoretic metrics are demonstrated for the seven protocols presented. 12For each voltage protocol, the core steps of the framework are: 1. Generate random samples of the parameters from chosen prior distributions; 2. Run the model to generate model output (current) for all samples; 3. Construct reduced-dimensional representations of the time-varying current output using Proper Orthogonal Decomposition; 4. Estimate information-theoretic metrics (between the current output and the parameters) such as mutual information, entropy equivalent variance, and conditional mutual information using non-parametric methods; 5. Interpret the metrics; for example, a higher mutual information between a parameter and the current output suggests the protocol yields greater information about that parameter, resulting in improved identifiability; and 6.Integrate the information-theoretic metrics into a single quantitative criterion, encapsulating the protocol's efficacy in estimating model parameters.
The remaining paper is organised as follows.In the next section, the methodology used in the present work is explained in detail, and organised into sections dedicated to each of the aforementioned steps of the framework.In Section 3, the results of an investigation using the framework to compare voltage protocols found in literature are presented and discussed, as well as discussions surrounding the characteristics and general applicability of the proposed framework.Finally, additional details and supplementary investigations that do not fit within the main body of the article are presented and discussed in the Supporting Information.

| METHODOLOGY
This section presents methodological details of the adopted model for ion channel kinetics, description of voltage protocols, parameter priors and associated sampling, data preprocessing, dimensionality reduction, and finally the information-theoretic metrics and their estimation.The overall framework is shown in Figure 1.

| Model for ion channel kinetics
The mathematical model used to characterise the I Kr channel is a Markov model 17 based on the Hodgkin-Huxley set of equations 18,19 that approximate the cellular ionic kinetics.In this implementation, the model used is a deterministic mean-field approximation for the stochastic ion channel model.The channel is described as four states: closed ( ), and closed and inactivated ([IC]), and as such, the values of these states represent the probability of an individual I Kr channel being in that corresponding state, or equivalently, the proportion of the channels found within a cell membrane being in that state.The system of ordinary differential equations that describes these channel states is as follows: where k 1 -k 4 are the transition rates associated with the states, and the eight unknown parameters P 1 -P 8 characterise the exponential voltage-dependence of these rates.These relationships are shown below: The fourth channel state, [IC], is not represented by an ODE and is instead constrained by the sum of probabilities being equal to one: Solving the system of ODEs for a time-dependent voltage input, V t ð Þ, yields the values for each of the channel states over time.An example of such a solution (using protocol Pr7 from Reference 12 and using reference parameter values; see Section 2.3) is shown in Figure 2. Since the I Kr current can only flow when the channel is open and activated ([O]), the current can be modelled as a standard Ohmic expression, dependent on the proportion of channels that are open: where G Kr is the conductance of the channel, V is the voltage applied across the cell membrane, and E K is the resting potential.G Kr is unique as a parameter as it does not affect the system of ODEs being solved.Instead, it simply scales the magnitude of the current curve according to the maximum conductance of the channel.It is important to note however, that this does not mean that an estimate for G Kr can be obtained independently from the other parameters, as the values of [O] are unknown a priori.As such, it may be possible to fit a I Kr curve identically to different sets of parameters; those with a low conductance but high [O] values and those with a high conductance and low [O] values.Due to this unique identifiability challenge, specific measures have to be taken to estimate G Kr , such as those described in Reference 12 and Appendix 4.2 in Data S1, which are independent from the main protocol design.For example, the large voltage steps present at the start of Pr6 and Pr7 used to estimate the lower bound of G Kr can be implemented in any voltage protocol, and it would therefore stand to reason that any candidate protocol could be adapted in this manner to improve the identifiability of G Kr .For this reason, while G Kr was sampled and treated in the same fashion as the other unknown parameters, it will be excluded from the information theoretic analysis of the project, as the ability of each protocol to identify the ODE model parameters is of principal interest.The resting potential for potassium ions, E K , was determined in Reference 12 by measuring the concentration gradient of K þ and evaluating the Nernst equation for each experiment.Since E K is not influenced by the voltage protocol, and only a minor degree of variability was observed, an average of these results was used for evaluating Equation (9) in this work.
F I G U R E 2 Voltage trace of protocol Pr7, along with the ordinary differential equation solutions for the reference parameters P !ref (described in Section 2.3) at each time step.The characteristic sum-of-sinusoids section occurs between t = 3000 and t = 6500 ms.The [IC] channel state is omitted, as it is simply defined as in Equation (8).

| Voltage protocols
A total of seven candidate voltage protocols, referred to as Pr1-Pr7, are analysed.A full description of each protocol, and the rationale behind their designs can be found in the original work. 12The first five protocols are square step signals, based on previous experiments in literature, 20,21 and designed to induce activated (Pr1-Pr3), inactivated (Pr4), and deactivated (Pr5) ion channel states.It is worth mentioning that these protocols are typically used to characterise specific processes, such as the activation or inactivation of the ion channel, in relative isolation.For a full explanation of the design of Pr2-Pr5, the reader is referred to. 22The sixth voltage protocol (Pr6) is a combination of simulated cardiac action potentials, including a range of different cells, species, and some signals exhibiting arrhythmic or otherwise pathological behaviour, and has hitherto been used exclusively as a validation protocol.The seventh protocol (Pr7) is the novel short-duration sinusoidal protocol proposed in Reference 12, and unlike Pr1-Pr5, is designed to characterise the entire model using the results obtained from one experiment.Voltage protocol Pr7 is presented in Figure 2, and the six remaining protocols are shown in Appendix 4.1 in Data S1.An additional protocol, consisting of Pr2-Pr5 appended to each other, is introduced in the brief investigation outlined in Appendix 4.5 in Data S1.

| Parameter sampling and reference values
The ODE system is characterised by eight parameters P For the purposes of this investigation, the results of the Bayesian inference scheme in Reference 12 will serve as reference parameters, P !ref .These parameters were found to be the closest-fitting parameters to the experimental current traces, 12 and as such yield physiologically relevant results.
Suitable prior distributions are chosen for random sampling in the following manner.Initially, the prior distributions were identical to the one used for the Bayesian inference scheme in Reference 12, but this was found to result in many samples yielding physiologically unlikely I Kr curves.By examining the samples where this occurred, it was found that this behaviour occurred for samples with parameter values significantly higher than P !ref , implying that the upper bounds are too high.The prior distribution for P !was, thus, redefined with a decreased upper bound.The lower and upper bounds for all the parameters for the adopted uniform prior distributions are shown in Table 1.
The prior distributions are randomly sampled to generate n ¼ 10,000 parameter vectors P !i , where i ¼ 1, 2,ÁÁÁ, n f g .For each protocol individually, the model is simulated for every parameter sample P !i to generate the corresponding current output I Kr,i .The choice of n ¼ 10,000 is justified by analysis of the I Kr curve modes following principal component analysis (PCA), which is described in Section 2.4.
As an additional step, for each parameter vector sample, the transition rates k 1k 4 are computed, and if any of these rates are outside the physiologically feasible range of 1:67 Â 10 À5 to 1000 ms À1 , the parameter vector is resampled.This step was taken to ensure the same sampling procedure as in Reference 12 was followed; transition rates below 1:67 Â 10 À5 ms À1 would result in a practically redundant level of voltage dependence, and the upper bound of 1000 ms À1 extends beyond what is physiologically expected from hERG.

| Data processing and dimensionality reduction
For every P !i , the current output I Kr,i obtained by numerical solution of the model is function of time and obtained at discrete time-intervals t 1 , t 2 , ÁÁÁ, t m f g .The dimension of this current output vector is quite large.For example, Pr7 is 8 s long, and if a discrete time interval of 0.1 ms is utilised, then m ¼ 80,000 (the total number of time steps).
The non-parametric methods available for the estimation of information-theoretic quantities (Section 2.5) provide poor estimates when vector sizes are in such high dimensions. 23Thus, dimensionality reduction of the I Kr,i vector is proposed.This is performed through the use of PCA. 24,25The data matrix M ℝ nÂm , whose ith row contains the current output I Kr,i corresponding to the parameter P !i , is factorised through singular value decomposition (SVD) as: where Σ ℝ nÂm is a diagonal (rectangular) matrix containing the singular values of M, and U ℝ nÂn and V ℝ mÂm are orthogonal matrices.This is equivalent to the columns of UΣ being the principal component scores, and the columns of V being the eigenvectors, or 'modes'.Dimensionality reduction from m-dimensions to k-dimensions (with k ( m) is achieved by only considering the first k dominant modes and discarding other modes.For all the protocols, k ¼ 3 was found to be sufficient, which captured an average of 98% of the variance across all protocols, with a minimum explained variance of 95.4% for Pr7, thus preserving most of the data content in the reduced dimensional space.Following this, the m-dimensional I Kr,i vectors are represented as three-dimensional vectors, which is a suitable size for reliable estimation of the information-theoretic quantities.SVD was implemented by a randomised truncated method described in Reference 26.The distribution of the first k principal components (the reduced dimensionality representation) is shown in Figure 3. Analysis of the eigenvectors/modes was used to justify the number of samples used; increasing the sample size beyond n ¼ 10,000 did not significantly change the shape of the modes, as can be seen in Figure 4.
Finally, the distributions of the parameters and PCs were transformed to standard normal distributions (zero mean and unit variance), in order to improve the performance of the non-parametric estimators for the information-theoretic F I G U R E 3 Probability density functions for the transformed data sets.While the distributions for each principal component are significantly more Gaussian than they would otherwise be, it is not possible to make them perfectly Gaussian using homeomorphic transformations.Note that the width of the transformed distribution has been scaled by a factor of 100 to allow for presentation on the same axes.
quantities.This is because, as a result of the local Gaussian approximations used in the entropy estimators, 27,28 applying data that resembles a Gaussian distribution would result in a more accurate approximation between the estimator and the true distribution, thereby lowering the error.
For the transformed principal component scores, this was implemented using a simple power transform of y ¼ x λ , where the optimal value of λ is determined by a minimisation algorithm 29 employed in tandem with a test for normalcy. 30This is implemented by finding the λ that yields the lowest squared kurtosis and skew, since normal distributions are characterised by possessing close to zero skew and kurtosis.For the parameter samples P !i , a log transform was found to be sufficient to achieve normalcy.Note that the application of homeomorphic transformations such as these to entire datasets does not affect the information-theoretic quantities of interest. 31

| Information-theoretic quantities of interest
The information-theoretic quantities of interest are mutual information and conditional mutual information, which are both derived from the central concept of entropy.For a continuous random variable X with probability density function where S is the support set of p x ð Þ. Entropy quantifies the uncertainty in a random variable. 32,33By observing a related random variable Y , the probability density of X would update from the prior distribution p x ð Þ to the posterior distribution p xjy ð Þ.The posterior uncertainty can be quantified again by the entropy H XjY ð Þ.Thus, the amount of information gained about X by observing Y is the difference between the prior and the posterior entropies.In this sense, gain in Comparison of the modes/eigenvectors from principal component analysis with 10,000 and 20,000 samples.The parameter sample sets were generated separately, that is, do not contain overlapping data.The lack of significant change in these modes indicates that a sample size of 10,000 is sufficient.
information is equivalent to reduction in uncertainty.This difference of entropies is also known as mutual information MI X;Y ð Þbetween X and Y , and can be written as In the context of parameter estimation and identifiability, the mutual information between the parameters and the measurements is relevant, as it quantifies the information contained in the measurements about the parameters.
To quantify correlations between parameters, the information-theoretic quantity of conditional mutual information (CMI) is relevant.It is defined as the expected value of the mutual information (MI) between two variables, X and Y , when a third variable Z is known 34,35 : If X and Y are two parameters and Z represents the measurements, then CMI X, Y ; Z ð Þquantifies the amount of information that would be gained about X, in addition to what has already been gained by observing Z, if Y was also known.In this sense, CMI X, Y ; Z ð Þ, quantifies the dependence/correlation between the parameters X and Y .For a more detailed discussion on the information-theoretic quantities presented above, the reader is referred to References  31,33,36,37.In an optimally designed experiment, it is desirable for there to be a strong correlation between the results of the experiment, and the underlying parameters the experimenter is trying to identify.This is equivalent to maximising the MI between the output of the model (which would be the measurement in an experiment) and the input parameters.Additionally, to improve individual parameter identifiability it is usually desirable for the CMI between the combinations of input parameters (given model output) to be minimised.This is because a high correlation between parameters implies they may only be estimated jointly, not separately, despite the objective of the experiment being identifying each individual parameter.It is for these reasons that MI and CMI are both used as the basis for evaluating protocol performance.
In our setting, only samples (Section 3) from the joint density function of the parameters and the measurement (current output from the model) are available.Since the underlying distributions are unknown, MI and CMI must be estimated directly from the samples.Such estimation algorithms are either parametric, where the form of the density function is known (or assumed) and its parameters estimated from the samples, or non-parametric, where such assumptions are not made and the density function is approximated through samples (kernel density estimation or nearest neighbours); an overview of these techniques can be found in Reference 38.In this work, assuming the form of the posterior pdf would be risky due to the stiffness of the ODE system, which results in a somewhat non-uniform distribution of 'acceptable' parameters, that is, parameter vectors who yield physiologically viable I Kr curves.While it would be possible to carry out an additional study to approximate the pdf, and use this to parametrise a parametric entropy estimator, non-parametric entropy estimators have been shown to be accurate and computationally efficient for similar applications. 27,28,37For this reason, the non-parametric k nearest-neighbours based estimators described in References 28,39 are used to estimate MI and CMI.
Once the estimates for MI and CMI are obtained, it is desirable to combine these into a single metric that can be used to assess the individual parameter identifiability, as well as overall information content, of each protocol.Multiple candidate criteria are considered.The first statistical criterion S 1 is based on the criterion in Reference 37, that is maximised by increasing MI and decreasing average CMI.It is equivalent to the average MI across all parameters, minus the average of all pairwise CMI values: where P i represents the random variable corresponding to the ith parameter, ζ represents the random variable corresponding to the reduced dimensional current output, and the total number of parameters is N ¼ 8.The intent behind this proposal is to create a measurement analogous 37 to D-optimality, 16,40 which is the process of maximising the determinant of the Fisher information matrix of the experiment.In this analogy, the elements on the main diagonal of the information matrix could be considered the MI, and the remaining elements each of the CMI combinations.This can be applied to other classical optimisation metrics; an analogue to A-optimality can be found by simply summing the values of MI P i ; ζ ð Þ, and calculating the eigenvalues of the information matrices for E-optimality.Caution must be exercised here, as while Fisher and Mutual information can be thought of as quantifying the same relationship, and have been found to produce congruent results, 37 they are still different quantities and are not interchangeable.A more detailed explanation on the similarities and differences can be found in Appendix 4.7 in Data S1.
Relying on the mean MIs and CMIs has a drawback that the measure S 1 averages information gains across all the parameters, and thus may be overly influenced by few parameters or fail to identify parameters with very low information gain.While this is unavoidable to a certain extent, due to the fact that eight MIs and 8 2 ¼ 28 pairwise CMI values need to be combined into a single metric, it is possible to adapt S 1 to better describe the consistency of the protocol's ability to identify individual parameters.One such metric could be to use the standard deviation of the MIs as a criterion: where μ is the average MI for that protocol.To meet the objective of having a single criterion by which to compare protocols, Equations ( 14) and ( 15) can be combined to yield a third possible criterion: where η is a weighting factor for the standard deviation.Maximising S 3 implies maximising MIs (thus maximising individual parameter information gains), minimising CMIs (thus minimising pairwise correlations in the parameter posteriors), and minimising the spread of MIs (thus penalising large variations in MIs).
It is worth considering that, despite being a measure of how evenly distributed the MIs are, S 2 (and consequently S 3 ) still do not provide information about the minimum level of performance.Thus, another suitable metric can be the minimum information content across all the parameters: Depending on what the desired outcome of the protocol designer is, other metrics combining MIs and CMIs may also be developed.Furthermore, a designer may even consider multiple metrics in a multi-objective design optimisation framework to balance multiple metrics.
A drawback of MIs is that their magnitudes (measured in nats or bits) are difficult to interpret as they have little physical interpretation aside from when the magnitudes are extremely high or low, implying identifiability or unidentifiability.This is important for comparing voltage protocols or other experiment designs, as one of the objectives of the analysis is finding the degree to which protocols out-perform one another on the basis of information content, and balancing this with other considerations.To assist in the interpretation of MI as the degree of parameter identifiability, a concept known as the entropy equivalent variance (EEV) 31 is also calculated and employed in the comparison of voltage protocols.
The EEV is derived from the estimated MI and the variance (σ 2 ) of the prior distribution for each variable: and represents the variance (error) of a hypothetical device that can directly measure each parameter in lieu of estimating the parameter from the measurements through an inverse problem, while providing identical MI.A protocol whose EEV for a particular parameter is double that of another protocol's indicates that the distribution of the noise that corrupts the hypothetical direct measurement of the parameter has twice the variance, which is a far more meaningful interpretation than simply concluding that the first protocol contains x more nats or bits of information.Therefore, it would follow reason to redefine S4 as S 4 in terms of the EEVs as: Note that unlike for S 1 -S 3 and S4 , where higher magnitudes imply better protocols, for S 4 lower magnitudes imply better protocols.

| Summary of computer codes
The following pseudocodes provide a summary of the above sections, describing the structure and process of obtaining information estimates for each protocol from the initial prior parameter distributions.In this work, the CVODE solver as part of SUNDIALS 41,42 was used due to the stiff nature of the ODE systems for certain parameter combinations, and all codes were implemented in Python and C.
By far, the most computationally expensive process of the framework is the invocation of the ODE solver (Algorithm 1, line 27) for each of the 10,000 samples.The rest of the algorithms, such as dimensionality reduction and entropy 0:0,0:0,0:0 ½ , t 0:1, j 1 9: h max 0:1 ► The output timestep is set to 0.1 ms, equal to the resolution of the protocols and original experimental data.In reality, much smaller timesteps will be used by the ODE solver to achieve convergent results.I Kr,i ODESolve(P !i ,PR) 28: end for 29: Let M ℝ nÂm be a data matrix for the current protocol, where m is equal to the number of time steps for that protocol, and each row i corresponds to the calculated current I Kr,i 30: end for estimation, incur fractions of the computational cost due to primarily comprising of simple, vectorized operations.The principle advantage in cost offered by this methodology is the relatively low sample-size required when compared with MCMC-based methods.Using this method, a full Bayesian analysis using MCMC simulation would be carried out to generate posterior distributions of the parameters for each protocol.Protocols with a higher efficacy/information content would yield posterior distributions with a narrower variance; there would be less uncertainty in the parameter values.This is in essence the same conclusion as can be drawn from a higher MI content estimated in this method, but with MCMC requiring a far greater number of samples.For example, the MCMC simulation in Reference 12 ran for 250,000 samples, which is a 25-fold increase in computation.Additionally, other steps involved in Bayesian inference schemes, such as global optimisation and adaptive sampling, can also incur high computational costs.A final advantage is that simply simulating n samples is easily parallelisable and can be shared across several machines if required, compared with a MCMC analysis which is characteristically state-dependent.
It is also worth considering that, depending on the number of samples and timesteps, a large amount of RAM is used by the dimensionality reduction process (Algorithm 2, line 1) in loading the matrix UΣV T , which would restrict this process to workstations and HPC clusters with sufficient memory capacity.
An overview of the computational cost of this method, and a comparison with a MCMC-sized workload, derived from the single-core computational time of a single sample (using reference parameter values) is presented in Table 2.The actual values may vary depending on parameter samples used.The computation times were obtained on an Intel Core i9-13900K CPU with 64 GB DDR5 5200 MHz RAM, running SUNDIALS 2.4.0 compiled with GCC 12.3.0.
Algorithm 2 Dimensionality reduction and normalisation of M ℝ nÂm , and estimation of information-theoretic metrics and derivative criteria 1: UΣV T factorised form of M, where the columns of UΣ consist of principal component scores, and V contains the eigenvectors 2: M is truncated after the first k ¼ 3 dominant modes, resulting in a reduced dimensionality data matrix M DR ℝ nÂk 3: 4: Let P !x be a parameter vector that contains all n samples for that specific parameter 5: for x 1 to 8 do 6: μ arithmetic mean of P The transformed parameters are used henceforth in place of the original parameters, but will also be referred to as P for notational simplicity.8: end for 9: 10: for PR 1 to 7 do 11: Let M !DR,x be a vector containing all n scores of the principal component x for the current protocol 12: for x 1 to k do 13: The estimated information content of each protocol, in the form of the MI for each parameter, is presented in Table 3.
As explained in Section 2.4, the information content related to conductance is not considered a factor in the comparison of voltage protocols, and a more detailed justification of this decision can be found in Appendix 4.2 in Data S1.The results in Table 3 show that, at a glance, the MI content of each protocol is comparable, with only a few entries that deviate greatly.For example, the MI of Pr1 for P 6 stands out the most, at just 0.006 nats, implying poor identifiability for this parameter.Pr1 also suffers with identifying P 8 , but to a lesser degree than P 6 .Other areas of poor identifiability are Pr2 for P 3 , and Pr4 for P 4 (0.03 and 0.06 nats, respectively).Such low values suggest that the design of these protocols lacks features that are affected by these parameters, and variations in these parameters do not significantly affect the simulated I Kr curves.The remaining protocols are, overall, more consistent in providing information about all the parameters.However, as previously mentioned, analysis of the information quantities expressed in nats (or bits) carry a limited amount of physical interpretability, thereby making direct comparisons between them difficult and potentially misleading.By applying Equation 18, the MIs can be transformed to the more interpretable EEVs, which are presented in Table 4.
Note that, while EEV is far more interpretable than information values expressed as nats or bits, the overall ranking of the protocols is equivalent for the two metrics.In other words, the hypothetical conclusion 'Protocol A contains more information than Protocol B' would hold true regardless of whether the information is expressed as MI or EEV; what would be difficult to ascertain without the EEV is the extent to which A is better than B. This means that the three criteria S 1 -S 3 , despite being expressed in nats, would yield the same ranking of protocols if they were converted to EEV-based metrics in the fashion of S 4 .
As can be seen from the MI values in Table 3, the ranking of each entry in terms of EEV and MI is equivalent, as expected from the definition in Equation (18).In Table 5, the results of the proposed criteria, S 1 -S 4 , are presented.At a glance, it is clear that these metrics assign different rankings across the candidate protocols-sometimes dramatically different, in the case of S 4 (whose value is sought to be minimised whereas S 1 -S 3 are to be maximised in an optimal protocol).While it is possible that each of these criteria are valid in their own right, due to them appraising different aspects of an optimally designed protocol, to fulfil the requirement of this framework producing a single framework, further analysis of these results is required to select the most appropriate criterion S, and verify that it is indeed an accurate representation of the information content of the candidate voltage protocols.In order to facilitate interpretation of these results, data from Tables 3 and 5 were collated into information heat maps for each of the protocols.Examples of this for Pr1 and Pr 7 are shown in Figure 5.The information content of each protocol is given in a single figure, consisting of the MI and CMI estimates, as well as the S criterion for reference.Heat maps for the remaining five protocols are presented in Appendix 4.4 in Data S1.
Figure 5 clearly shows that Pr1 possesses extremely high and low values across both the MI and CMI, with the maximum and minimum values across all the protocols for MI and CMI being found in its heat map.Recall that low MIs and high CMIs are both detrimental to parameter identifiability.Thus, P 1 and P 7 appear to most identifiable with this protocol.However, the CMI plot shows that P 7 , unlike P 1 , has high correlations with respect to other parameters.This implies that most information about P 7 is contained in the joint distributions of its combinations with other parameters, that is, a high correlation exists in the posterior of P 7 and its combinations with other parameters.Alternatively viewed, this implies that if P 7 was known or fixed, then other parameters would individually become much better identifiable.On the other hand, Pr7 does not show large variations in MI and CMI content, implying moderately good information content (MI) for all the parameters and moderate correlations between all pairwise combinations.In this sense, Pr7 does seem superior to Pr1.
The heat maps also show that there does not seem to be a strong trend in identifiability across parameters between protocols-that is to say, certain parameters do not seem to be intrinsically more or less identifiable than others.This implies that it is possible to create 'hybrid' protocols that may yield results with consistently high information content across parameters.A large number of different protocols could be simulated, and protocols with high information content for specific parameters (such as Pr1 for P 1 and P 7 ) could be selected to derive a new set of information-rich protocols from.This could be achieved, in its simplest form, by appending truncated versions of the protocols to one another, in a similar fashion to the method outlined in Appendix 4.5 in Data S1, and optimising the degree to which each protocol is shortened in order to minimise protocol length.Another possibility would be identifying the most information-rich features of the eigenvectors/modes of the current curves, in order to determine what type of current variations yield high information about individual parameters.
The proposed single quantitative criteria of S 1 -S 4 , can also be compared for each protocol, and the most appropriate one selected.Such objective scalar metrics can also be used in an optimisation/optimal-design framework with parameterised voltage protocols to devise new and better protocols.The S 1 -S 4 values for the seven protocols are provided in Table 5, and show a general lack of concordance in how they rank each protocol.This is to be expected, as most of the criteria are proposed from separate observations or rationalisations made about the data.It is also worth highlighting that S 1 and S 3 , based on information aggregates, are sought to be maximised, while S 2 and S 4 , the average variance and the worst-performing EEV, are sought to be minimised.According to S 1 , Pr1 is the best performing protocol, with Pr7 trailing by a significant margin (À0.06 for Pr1 and À0.14 for Pr7).From the MI results in Table 3 and Figure 5, this is due to some of the high MI values present for Pr1, although the shortcomings associated with relying on average information discussed in Section 2.5 immediately become apparent as Pr1 is not sufficiently penalised for its extremely low identifiability for P 6 and P 8 .While S 2 does identify those protocols with the most obvious gaps in parameter information (Pr1 and Pr4), it does not provide any information beyond this, and therefore makes it hard to compare the remaining similar-performing protocols.Furthermore, while in this case the outlier information values are very low, if a protocol were to have outliers in the form of extremely high information, they would be penalised regardless.This makes S 2 unsuitable as a standalone metric, although it could be used for other analyses alongside more complete data.Its incorporation into S 3 seems to have improved the performance of S 1 , but still suffers from the aforementioned problems.
S 4 , while arguably simpler than the aggregate criteria due to being the worst-performing EEV for each protocol, does meet the objective of the desired criterion, which is to provide a metric by which to compare protocols for their ability to characterise the I Kr channel model.Since the model cannot simulate reliable results without a good estimate of all parameters, there is a need to ensure a minimum level of identifiability for each parameter.This requirement is met by S 4 , as it becomes immediately apparent which protocols are the most consistently information-rich.Therefore, as a singular metric to rate the performance of each protocol, S 4 is the most appropriate for this task out of those proposed, and will be used moving forwards in a design optimisation context.
Upon examining the results in Table 5, several conclusions can be drawn from the suitability of each S by contextualising these results, and Pr1-Pr7 themselves, with the original investigation. 12The rationale behind basing the design of Pr7 on a sum of three sinusoidal waves is stated as being to 'sweep' the time and voltage-dependent characteristics of the I Kr channel, which came from the experience and intuition of the experiment designer.The square-wave based protocols, Pr1-Pr5, would therefore lack this element, and it would therefore not be surprising if these protocols would fail to explore one or multiple characteristics of the ion channel.This can be seen in the apparently very low, or 'missing', information these protocols contain about certain parameters.The observation that these square-wave protocols carry deficits in different parameters to one another is supported by the differences inherent in each protocol.For example, Pr2 varies the length of each peak voltage, whereas Pr3 varies the magnitude of the peak voltage, and are otherwise identical.
Pr7, on the other hand, does not appear to suffer from this flaw, which is reflected in its values for S 4 , but not necessarily for S 1 -S 3 .The fact that 12 concludes that the model parameterised with Pr7 outperforms literature models is further proof that S 4 is an appropriate metric for judging the efficacy of a voltage protocol.It is worth keeping in mind that the full Bayesian investigation in Reference 12 was only carried out for Pr7; replicating this with another protocol and comparing the results with the predicted rankings for each S would be the surest way to confirm both the accuracy of the information estimates, and the ranking scheme derived from them.
To verify that these observations are in fact due to the interactions between the voltage protocols and the parameters of the cellular model, and not due to any conceptual errors in the framework and/or numerical errors in MI/CMI/EEV, the veracity of the results can be tested by observing the change in the model output as the parameters of interest are changed.When a protocol that contains a low amount of information for particular parameters is simulated or otherwise carried out, one would expect the resulting I Kr curve not to be significantly affected by perturbations of those parameters.In other words, a simplified form of sensitivity analysis is carried out, which has previously been successfully implemented in the context of electrophysiological parameterisation. 15,43The expected outcome can be observed in Figure 6, which shows (simulated) I Kr curves being virtually indistinguishable in Pr1, even when P 6 is increased by as much as 25%, but the same change in P 6 yielding a far more apparent change in Pr7.
It is also of note that Pr7 and Pr3 have highly similar results in terms of S 2 , S 3 , S 4 and average (C)MI values, indicating that these two voltage protocols contain very similar amounts of information regarding the model parameters, and is also consistent across parameters.However, this has to be considered along with the fact that Pr7 was designed to be as short as possible in order to allow full sets of results to be obtained from individual cells.Indeed, Pr3 has a far greater length, running for 57.8 s compared with just 8 s for Pr7.This can be interpreted as Pr7 having a higher information density, or average rate of gain in information, than Pr3.
It might therefore follow that the resulting S scores are made less interpretable by the different lengths of the simulated protocols.In other words, if a protocol is assigned a lower S than a longer protocol, would this difference in scores be compensated if the shorter protocol were extended to the length of the higher scoring protocol?Although the MI estimator calculates these results based on the protocols reduced to the same number of dimensions, the principal components are by definition affected by the shape of their original protocol, including its length.There does not appear to be an obvious solution to this, as lengthening shorter protocols would defeat the advantages that they provide for experimentalists, and truncating longer protocols would eliminate features which are required to identify certain parameters.
To illustrate this, Figure 7 shows how the EEV of P 1 -P 8 varies across the length of the protocol for Pr3, Pr5, and Pr7.These three protocols were selected for the broad range of protocol lengths they represent, and were divided into F I G U R E 6 Model output for Pr1 and Pr7.A section of each plot, shown as the shaded area, is enlarged to highlight the extent to which model output has changed for each protocol.Input parameters were equal to the reference parameters, with P 6 increased by up to 25%.P 6 was selected due to its low identifiability for Pr1 and high identifiability for Pr7, interpreted from estimated entropy equivalent variance values.
five 'slices', with the first slice comprising the first 20% of the protocol, the second slice comprising the first 40%, and so on.While the overall trend for all parameters shows identifiability increasing as the protocol progresses, there are instances where the EEV for certain parameters increases, implying that there are features that increase uncertainty for these values.It is important to note that this is not 'negative MI', which cannot exist, but rather that the correlation between output and parameter is weakened upon observing additional data, as each slice of a protocol is effectively a different protocol derived from the original.
It would be prudent for experimentalists to consider this limitation of measuring overall information content, rather than information density, if they wished to employ this technique to compare experiment designs.For example, if an experiment designer wished to compare the information content of two different sinusoidal protocols, it would be recommended to either make the protocols of the same duration, or to accept a possible reduction in parameter identifiability in exchange for the advantages of a shorter protocol, such as the ability to complete experiments on a single cell, and reducing the effect of data quality degradation.
Another limitation to this framework is the requirement of an accurate numerical model, which could pose problems if the framework were to be applied to a system which has not been modelled in existing literature, and would be especially problematic if there is no way to propose a reasonably accurate model from theoretical knowledge.While this may be a significant drawback in newly emerging fields where the literature is sparse, this problem would also be F I G U R E Graphs showing the approximate variation in entropy equivalent variance (EEV) across the duration of Pr3, Pr5, and Pr7.The trends for each protocol show that EEV decreases as the protocol progresses, which is expected.As each 'slice' of the protocol was processed in the framework as a separate protocol, some of these results show EEV increasing for specific parameters, as the rest of the parameters decrease.This suggests that some features present in the protocols decrease the correlation between those parameters and the output.Speculating further, this may be a result of information about some parameters coming at a cost of information to others, implying a high conditional mutual information and low degree of independent identifiability.If the time-dependence of MI estimates is considered important to an experimentalist, for example if they wished to know how identifiability is distributed in time so that they can design their protocol to be information rich in the beginning, when the cell is at its healthiest, it might be useful to confirm this by comparing more protocols and taking additional slices.
shared by other techniques requiring a numerical model, such as machine learning.To confirm the framework's reliance on accurately simulated synthetic data, a cursory investigation was carried out (for full details please refer to Appendix 4.6 in Data S1).In this, the entire methodology was carried out for two protocols using an altered version of the model described in Section 2.1, with major errors introduced into the governing equations.As expected, this had a significant effect on the resulting MI and S criterion estimates, with some parameters having a near-zero estimated MI level, which can be interpreted as the erroneous model effectively ignoring these parameters.It is worth mentioning that, while the results are clearly invalid, they were very similar for the two protocols, and obtaining results such as these could provide an indication to the experimentalist that the underlying numerical model possesses a serious flaw.Interestingly, this means that a secondary purpose of the framework could be to confirm the veracity of a proposed numerical model, whereby the framework yielding nonsensical MI/EEV estimates signifies the presence of high structural model error.
In addition to the proposed use of S 1 -S 4 to rank and compare protocols in general terms, the results obtained through this framework could be used to draw conclusions about the information content for individual parameters, which could be useful for an experiment designer who wishes to design a voltage protocol to measure a specific parameter or combination of parameters.Indeed, Reference 22 explains how Pr2-Pr5 are designed in such a way to parametrise-specific transition rates for the ion channel model, as opposed to the entire system, as is the case for Pr7.For example, Pr2 is designed to characterise the time constant of activation for the channel, which corresponds to the transition between [C] and [O] and [I] and [IC] in this paper.As such, one would expect there to be higher information gain for the parameters associated with this state transition (P 1 -P 4 ).However, the data in Table 4 appears to contradict this, as there is no clear correlation between information content and whether the parameter in question belongs to the P 1 -P 4 group.Similar observations can be made for the rest of the protocols mentioned in Reference 22.
It seems likely that this discrepancy might be explained by differences in methodology between References 12 and 22; three of the four model fitting methods described in the latter paper focus on identifying the steady-state and timedependent transition rates, as opposed to the estimation of the values for P 1 -P 8 in the former.Since the transition rates are expressed as functions of combinations of parameters, it stands to reason that while a protocol may be able to yield results that can be used to accurately fit the transition rates, the identifiability for individual parameters may still be low, as this is what the entropy estimators are attempting to quantify and nothing else.Another possible source for this unexpected observation could be model error; perhaps the eight parameters used in this model could be expanded to a greater number of parameters, hinting at relationships between parameters and channel states that are lost in the model.It would be worth investigating this further, focusing on comparing the information content about model parameters within a protocol of interest, rather than comparing protocols.For example, if the model were to be characterised by the steady-state and time-dependent activation and inactivation transition rates, would the mutual information of those parameters tally with what the protocol was designed to capture?With solely the results obtained from this work, it is impossible to definitively state whether or not this is the case, but a brief investigation was carried out in Appendix 4.6 in Data S1 in order to explore the effects of appending the results of multiple protocols, to reflect the intended use for these protocols being to identify certain parameters in relative isolation.
Finally, we anticipate that this framework for evaluating voltage protocols would be applicable to a broad range of Markovian models describing voltage-dependent ion channel dynamics, such as those describing Na v 1:5 and L-type Ca 2þ channels.This is because there are no specific conditions or features unique to the dynamics presented in hERG that would render the analysis invalid for other channels, and the overall framework does not depend on the model selected for Section 2.1.However, it is possible that other models might exhibit more complex or inconclusive patterns of results that would require the addition of a fourth (or fifth) POD element in Section 2.4, which in turn might reduce the accuracy of the entropy estimators, who struggle in higher dimensions.While it would be useful to explore this, only one ion channel was explored in this work for brevity, and I Kr was selected both for its relation to the original work, 12 and its propensity for blockade as a side-effect for a very wide range of drugs. 44hat simply employing the average information content across all parameters neglected the crucial requirement of a minimum level of identifiability for each parameter.Therefore, the criterion describing the minimum level of identifiability, S 4 , is proposed as being most appropriate.It is worth noting that other combinations of MIs and CMIs (and perhaps even higher than pairwise conditional mutual information) may also be considered based on the desires of the experiment designer.For example, if the requirements were to be changed so that information regarding one parameter is no longer important (e.g., it is already known), this can easily be taken into account, due to the available breakdown of information content per parameter.Nevertheless, the ideas of MIs and CMIs provide suitable metrics, which carry intuitive meaning of 'information' contained in a protocol.Furthermore, the proposed framework allows for efficient computation of the information-theoretic metrics and can be adapted.
Plotting the MIs and CMIs as heat maps has proven to be a valuable tool for analysing the information content of each protocol, as it allows for a quick and intuitive comparison between MI values, combinations of CMI, and the S score of that protocol.When used in conjunction with S, these plots of information metrics fulfil the requirement of providing a means through which to rank several candidate voltage protocols/experiment designs, while the criterion S is appropriate for use when a single numerical value is required.The latter is useful to be employed in optimisation frameworks, where S could be optimised over a space of parameterised candidate protocols.
Lastly, the results of the information-theoretic analysis are verified with sensitivity analysis and align with the expected outcome (as argued in Reference 12) of Pr7 being the most optimally-designed voltage protocol for characterising the ion channel model.This is based on the comparisons of average information gain, and minimum level of identifiability across all the unknown parameters, along with the fact that it is the shortest protocol. 12

Algorithm 1 7 :
Obtaining samples of I Kr 1: for i 1 to n do ► Where n is the number of samples i.e. 10,000 2: P !i a random sample from the (uniform) prior distribution of P !3: end for ► For each of these samples, ion channel current I Kr is calculated for each protocol PR 4Defines voltage and time data V ! and T max depending on the value of PR 8: Y !0

F I G U R E 5
Heat maps displaying the information content (mutual information [MI], conditional mutual information [CMI]) for the two chief voltage protocols of interest.Pr1 was selected as the main point of comparison for Pr7 due to it out-performing Pr7 in terms of S criterion, which is based on average information content rather than more comprehensive techniques.Being the novel protocol of interest, Pr7 scored second highest in terms of S, but has visibly far more consistent information content across all parameters when compared to Pr1.The remaining five voltage protocols, shown in Appendix 4.3 in Data S1, show similar behaviour to Pr7, albeit with poorer S x scores and slightly higher variance in MI and CMI.
T A B L E 2 Computation times associated with each voltage protocol.
T A B L E 4 Entropy equivalent variance of each parameter for each voltage protocol.Summary of proposed information criteria for each voltage protocol.