We employed mathematical models to understand variability in the response to HERG-blocking drugs using a technique of parameter randomization followed by multivariable regression described in the Methods section and elsewhere.^{17,18} respectively illustrate the procedure and show an especially relevant result obtained with the TNNP model.^{19} Two sets of randomly-chosen parameters (Trials 97 and 270) result in baseline APDs extremely similar to that produced by the control set of parameters (Trial 1). However, these changes in parameters greatly influence the response of the model myocytes to HERG blockade. Compared with the AP prolongation produced by HERG block in the control myocyte (23.7 ms), Trial #97 showed a smaller response (ΔAPD = 15.9 ms) whereas Trial #270 exhibited a more dramatic response (ΔAPD = 45.3 ms). Thus, variations in parameters can influence both the behavior of the unperturbed model^{17,18,24} and the response of the model myocyte to drugs.

The response of the population of 300 model cells is illustrated in . shows that ΔAPD in the population of models varies over a wide range (8.4 to 53.9 ms). illustrates that in the TNNP model, ΔAPD produced by HERG blockade is uncorrelated with the baseline APD before HERG blockade, similar to results observed clinically.^{6} Thus, model myocytes that display an especially deleterious response to HERG block cannot be identified simply from their baseline APs.

The multivariable regression procedure^{17,18} computes a matrix **B** that relates the input matrix of parameters **X** to the matrix of outputs **Y**, i.e. **Y**≈**X*B**. Thus, parameter sensitivities in **B** indicate how much changing a parameter affects a model output. displays the values in **B** corresponding to ΔAPD; i.e. the length of each bar shows how much that parameter influences the cell's response to HERG block, and parameters are ranked from the most to the least important. Green and red bars indicate parameters for which an increase in the parameter increases or decreases ΔAPD, respectively. therefore provides a comprehensive, systems-level understanding of how the model cell's electrophysiological properties influence its repolarization reserve. Some results are expected based on prior knowledge. For instance, important parameters include: the density of I_{Kr} itself (G_{Kr}; rank 1), the density of I_{Ks} (G_{Ks}, rank 2), and the kinetics of I_{Ks} activation (p_{xs}, rank 5). These are consistent with the obvious notion that I_{Kr} density will influence the response to I_{Kr} block, and with the conventional wisdom that I_{Ks} provides the main reserve current that allows HERG block to be tolerated.^{25,26} Conversely, parameters related to fast Na^{+} current (I_{Na}) and transient outward K^{+} current (I_{to}), which activate at the beginning of the AP and are considered unimportant for repolarization, rank relatively low. One advantage of the present strategy, however, is that the unbiased analysis can identify surprising behaviors. For instance, the voltage dependence of HERG inactivation (V_{xr2}) is the third most important contributor whereas the voltage dependence of HERG activation (V_{xr1}) is the least important parameter. shows additional simulations that confirm these predictions. A +20 mV shift in V_{xr1} has virtually no effect on the response to HERG block whereas a +20 mV shift in V_{xr2} greatly increases ΔAPD produced by HERG block (44.8 ms).

To identify additional counterintuitive findings, we plotted the regression coefficients for APD versus the corresponding regression coefficients for ΔAPD (). With this representation, parameters located in the 2^{nd} and 4^{th} quadrants have opposing effects on APD and ΔAPD and are therefore somewhat surprising. For instance, an increase in G_{Ca}, located in the 4^{th} quadrant, will prolong the baseline APD but improve the response of the cell to HERG block by reducing ΔAPD. Additional simulations were performed to determine the mechanism underlying this paradoxical result. As shown in , Increasing G_{Ca} does indeed lead to greater APD but reduced ΔAPD, as predicted by the regression model. This occurs because greater L-type Ca^{2+} current elevates the AP plateau, which leads to increased activation of I_{Ks}. Larger I_{Ks} during the plateau makes the cell less susceptible to block of I_{Kr}, and ΔAPD is consequently smaller. To express the idea of counterintuitive results in more quantitative terms, we calculated the absolute difference between the parameter sensitivities corresponding to APD and ΔAPD. The 5 parameters for which this quantity was greatest, indicating that a parameter influences the baseline AP and the drug response differentially, are listed in . Thus, the regression method can identify surprising model behaviors that can be understood in greater detail through more mechanistic simulations.

| **Table 1**Ranked list of TNNP model parameters that showed the greatest difference in the magnitude between the parameter sensitivities corresponding to APD and ΔAPD |

We extended our regression method to four other ventricular myocyte models.

^{20–23} After performing, with each model, an analysis similar to that shown in (

Figs. S1–S4), we ranked the parameters in each model from most to least important. The ranks of the common parameters were then averaged to derive a consensus list of the 10 that were consistently important in determining repolarization reserve (). Besides the density of HERG itself (G

_{Kr}), important parameters included G

_{K1}, G

_{Ca}, G

_{Ks}, and K

_{NaK}, the density of the Na

^{+}-K

^{+} pump. The appearance on the list of G

_{Kr}, G

_{Ks} and G

_{K1}, the three main K

^{+} currents that contribute to repolarization, makes sense physiologically. For instance, it has been demonstrated that I

_{Ks} can limit AP prolongation when this occurs by other means such as I

_{Kr} or I

_{K1} block.

^{25,26} Similarly, block of I

_{K1} has been shown to influence the response to HERG-blocking drugs,

^{25} and has been hypothesized in simulations to play an important role in the development of early afterdepolarizations.

^{27} In contrast, the role of the Na

^{+}-K

^{+} pump in determining repolarization reserve has been less well-explored and therefore suggests a potential avenue for future research.

| **Table 2**Consensus list of model parameters that were consistently important in the determining repolarization reserve across several mathematical models. |

The matrix framework of our approach allows one to understand complicated phenotypes, in which several electrophysiological properties are altered, as superimposed combinations of the individual effects. Since **Y** ≈ **X**·**B**, the approximate change in outputs in a given cell (**Ŷ**) can be computed as the vector describing that cell's parameter changes (**x**) multiplied by the matrix of parameter sensitivities **B**. shows examples of how this framework can provide insight. To compute an element of **ŷ**, corresponding for instance to ΔAPD in a particular cell, one multiplies each change in a parameter (element of **x**) by the corresponding element of **B**, then sums the individual products to calculate the overall response. For instance, G_{Ca} and G_{Ks} parameter sensitivities have opposite signs with regards to how they influence APD but the same sign with respect to ΔAPD. This means that if G_{Ca} or G_{Ks} are increased or decreased in parallel, the effects on APD will offset one another but the effects on ΔAPD will be additive. Thus, in the TNNP model, a cell with decreased G_{Ca} and G_{Ks} exhibits an especially dramatic increase in APD after HERG block, even if the baseline APD is essentially normal (). We can similarly understand clinically relevant responses to HERG blockade through matrix multiplication. For instance, Yang et al identified a silent mutation in KCNQ1, the gene encoding the α-subunit of I_{Ks}, that causes a single amino acid substitution (R583C) and predisposes individuals to drug-induced arrhythmia.^{28} It was not known, however, if the increased susceptibility was due primarily to the decrease in G_{Ks} or to the positive shift in the voltage--dependence of activation (parameter V_{xs} in the model). The regression model predicts that these two changes contribute roughly equally to the increased ΔAPD (). Similarly, we can simulate heart failure (HF) by changing four model parameters, as done previously by Winslow et al.^{29} illustrates how each change contributes to reduced repolarization reserve in HF, as simulated with the TNNP model. The regression model predicts that even though the decrease in G_{to} is large, this does not alter repolarization reserve because this parameter's regression coefficient is near zero. In contrast, the small decrease in G_{K1} markedly decreases repolarization reserve due to a larger value in **B**. Changes in SERCA and NCX contribute mildly to the increased ΔAPD.