Model-driven optimal experimental design for calibrating cardiac electrophysiology models

Models of the cardiomyocyte action potential (AP) have contributed immensely to the understanding of heart function, pathophysiology, and the origin of heart rhythm disturbances. However, AP models are nonlinear, complex, and can contain more than a hundred differential equations, making them difficult to parameterise. Therefore, cellular cardiac models have been limited to describing ‘average cell’ dynamics, when cell-specific models would be ideal to uncover inter-cell variability but are too experimentally challenging to be achieved. Here, we focus on automatically designing experimental protocols that allow us to better identify cell-specific maximum conductance values for each major current type—optimal experimental designs—for both voltage-clamp and current-clamp experiments. We show that optimal designs are able to perform better than many of the existing experiment designs in the literature in terms of identifying model parameters and hence model predictive power. For cardiac cellular electrophysiology, this approach will allow researchers to define their hypothesis of the dynamics of the system and automatically design experimental protocols that will result in theoretically optimal designs.

Two types of experimental protocols (u) are optimised: voltage clamp and current clamp (see § 4.2). 102 We defined voltage clamp protocols as sequences of steps with a magnitude and duration that could 103 be varied by OED as the design variables d. Our current clamp protocols consisted of short bursts 104 of current injection with a fixed magnitude and duration, here the interval between burst was used as 105 design variable (see Figure 1). OED can also work within a variety of constraints, such as restrictions 106 in duration of the experiment or applied voltage ranges, to ensure the optimal designs are practically   Each of the protocols appears to elicit rich and varied dynamics. Since these optimal protocols 120 were designed to maximise the information for all model parameters for inference, we evaluate the 121 performance of these optimised protocols and compare them against some protocols from the literature.  In terms of reducing the uncertainty and error in the inferred parameters, the two most successful We may also expect the best protocol from those OEDs would be a robust good 'all-round' protocol.

147
Each of the OED is by definition optimal under a given design measure, whilst a robust good all-round

161
Examples of cross-measure assessments for the OEDs are shown in Figure 6A for voltage clamp 162 and Figure 7A for current clamp. Figure 6A shows our best voltage clamp protocol design, the OED  Figure 6C for voltage clamp and Figure 7C for current clamp.
169 Figure 6C and Figure 7C show that the protocols based on multiple models (bottom row) are Normalised score (%)  for the 1 Hz pacing protocol: not only is the cross evaluation matrix consistently worse than the OED 207 protocols ( Figure 7B), but the posterior obtained in the practical assessment is also wider than all of 208 the OED protocols ( Figure 5). Moreover, the results from using the biomarkers ( Figure 5A) were not 209 only more uncertain than the other approaches but also biased in some of the parameter estimates, e.g.

210
I N a and I to . The results suggest that the OED protocols were much better in estimating the model 211 parameters than those in the literature. compared to some of the conventional protocols. 223 We explored various designs adapted from the OED literature. Most of them were able to improve 224 the uncertainty in the inferred parameters compared to the benchmark protocols. We tried to seek the 225 best protocol amongst the OED generated protocols. However, we were not able to identify the best 226 one, as each of them perform slightly differently under different situations. For example Protocol K 227 was the best in the cross-measures comparison, whilst Protocol D was better in reducing parameter 228 uncertainty and error. We found that the protocols based on averaged model were consistently (slightly) 229 better than the others, which may be a better way for designing patch-clamp protocols, which implicitly 230 makes some allowance for differences in ion current kinetics.

231
Fitting conductances to whole-cell recordings while using ion channel kinetics models taken (or

242
Amongst all the methods that we compared, the biomarker approach performed the worst, even 243 though it has been one of the most popular approaches. We further note that some of the biases in the 244 biomarker posteriors were due to the biases introduced in estimating the biomarkers with noisy data.

245
Such an issue is illustrated in Figure 8A, where biomarkers were estimated for 1000 noise realisations

262
The inset of Figure 8B shows a magnification around the apex of the AP, highlighting the biases of the guarantee that the designs will work better, or even work well, as discussed in (Lei et al., 2020b).

275
The results, however, have generated experimental protocols that can readily be tested experimentally,  In the design of experiments, the results of OED are experimental designs that are optimal with respect to some statistical criteria-optimal design measures. Formally, we consider a (controlled) nonlinear differential equation model of the form where x is the vector of model states, u is the vector of system inputs, θ is the vector of model stimulus, d could be the durations between each stimulus, as described in details in the next sections.

298
The optimisation of the experimental design procedure is then defined as where the function Φ represents the optimal design measure to be minimised.

300
All optimal design measures discussed in this paper assume the model calibration process to be where g j is the maximum conductance or permeability and x j (V ) is some nonlinear function of the 317 voltage V m for the current of type j, and we consider the major currents only (Eq. (6)). These nonlinear is the output defined in Eq.
(2). The vector of system inputs u describes the applied voltage-clamp 322 protocol V (t) which we assume to be the same as V m .

323
Alternatively, under the current-clamp configuration, the observed membrane voltage V m follows where I obs (t) is the same as Eq. (4). In this case, the membrane voltage V m is the output defined in 325 Eq.
(2). The vector of system inputs u describes the applied current-clamp protocol, i.e. the stimulus 326 current I stim (t). 327 We assume that these kinetics terms are known (and correct). We are interested in finding all g j ,

328
which determine the magnitude of the currents when the channels are fully open, therefore we have 329 θ = {g Na , g CaL , g Kr , g Ks , g to , g NaCa , g K1 , g NaK } (6) to be inferred, which are the maximum conductance or permeability of the fast sodium current, the L-330 type calcium current, the rapid-delayed rectifier potassium current, the slow-delayed rectifier potassium 331 current, the transient potassium current, the sodium-calcium exchanger current, the inward rectifier 332 potassium current, and the sodium-potassium pump current, respectively. We note that this ignores 333 some smaller background currents in the model. Each OED is optimal for the model used to perform the design. However, we may want to obtain 338 a design that is optimal for a range of biological assumptions instead of focusing on one particular 339 model, so that the same data can be used to study multiple similar models. Therefore, we also consider 340 using multiple models for calculating the cost function; instead of using only one model to optimise Here we use five adult ventricular (epicardial) AP models to represent an averaged behaviour: ten   The voltage-clamp protocol V is the system input u(d), and we choose it to be a piecewise function 353 defined as for a N -step protocol; equivalently we can define ∆T i = T i − T i−1 with T 0 = 0, then we obtain In order to ensure the resulting voltage-clamp protocol V is feasible to be run experimentally, we In the current-clamp mode, the stimulus current I stim is the system input u(d).
We further constrain ∆T i ∈ [100, 2000] ms, ∀i, to make the designed experiments practically usable. 374 We applied only the LSA designs for this experimental setting; to apply GSA designs for current clamp 375 experiments, we would need to first explore the boundaries (not necessarily being a hypercube) of the 376 parameter space that generate APs. The sensitivity matrix based on LSA, S L (d | θ =θ), is defined as Here the subscript i (row) in y i runs through all model outputs in y and all sampled time points t k * , 381 and the subscript j (column) θ j goes over all model parameters in θ. The Fisher information matrix 382 (FIM) is given by where Σ is the covariance matrix of the measurement data (noise), and we assume that data are a 384 single output (a single current trace) with i.i.d. Gaussian noise through time such that Σ = σ 2 I with 385 σ constant and I the identity matrix.

386
Many local design criteria are based on the parameter (co)variance matrix defined as optimal, D-optimal, and E * -optimal designs, Atkinson and Donev, 1992) are given in Table 1, where 399 λ max (X) and λ min (X), respectively, denote the maximum and minimum eigenvalues of a matrix X.

400
Local design criteria Cost functions Table 1: Local design criteria for OED. Here we have assumed the parameter covariance matrix C θ for some given model parameters θ can be approximated as σ −2 [S T L S L ] −1 , where σ −2 is dropped from the design criteria (cost functions) as it is assumed to be constant. E * is the modified E criterion.
The design criteria in Table 1  improving these properties of the confidence ellipsoid, the uncertainty of the inferred parameters using 405 the data measured under these optimal design protocol would (in theory) be reduced.

406
One obvious limitation of the local designs is that they are local, i.e. they depend on a particular 407 choice of the model parameters (θ =θ in the equations above). Hence the design is only guaranteed 408 to be optimal for that particular choice of parameters. However, we do not know the true model where p(θ j ) is the probability density function of the j th parameter, defined over a domain Θ (a  Global sensitivity analysis-based designs 417 One of the most straight forward ways of doing GSA-based design is by replacing the S L in Table 1  Sobol sensitivity index is given by where Var (y i ) is the total variance of the model output y i , Var θ j (X) = E θ j (X 2 ) − (E θ j (X)) 2 denotes 430 the conditional variance only over parameter θ j , and E θ −j [ · ] indicates the expected value taken over 431 all parameters θ except parameter θ j . Criteria using S G are given in  Cross-measure evaluations 454 We defined a cross-measure matrix X for each optimised protocol, such that each entry X i,j was the 455 criterion j score for this protocol with currents generated by model i, with 1 < j < N measure and The posterior distribution of the parameters was