Mechanisms Underlying Interactions Between Low-Frequency Oscillations and Beat-to-Beat Variability of Celullar Ventricular Repolarization in Response to Sympathetic Stimulation: Implications for Arrhythmogenesis

Background and Objectives: Enhanced beat-to-beat variability of ventricular repolarization (BVR) has been linked to arrhythmias and sudden cardiac death. Recent experimental studies on human left ventricular epicardial electrograms have shown that BVR closely interacts with low-frequency (LF) oscillations of activation recovery interval during sympathetic provocation. In this work human ventricular computational cell models are developed to reproduce the experimentally observed interactions between BVR and its LF oscillations, to assess underlying mechanisms and to establish a relationship with arrhythmic risk.

Materials and Methods: A set of human ventricular action potential (AP) models covering a range of experimental electrophysiological characteristics was constructed. These models incorporated stochasticity in major ionic currents as well as descriptions of β-adrenergic stimulation and mechanical effects to investigate the AP response to enhanced sympathetic activity. Statistical methods based on Automatic Relevance Determination and Canonical Correlation Analysis were developed to unravel individual and common factors contributing to BVR and LF patterning of APD in response to sympathetic provocation.

Results: Simulated results reproduced experimental evidences on the interactions between BVR and LF oscillations of AP duration (APD), with replication of the high inter-individual variability observed in both phenomena. I_CaL, I_Kr and I_K1 currents were identified as common ionic modulators of the inter-individual differences in BVR and LF oscillatory behavior and were shown to be crucial in determining susceptibility to arrhythmogenic events.

Conclusions: The calibrated family of human ventricular cell models proposed in this study allows reproducing experimentally reported interactions between BVR and LF oscillations of APD. Ionic factors involving I_CaL, I_Kr and I_K1 currents are found to underlie correlated increments in both phenomena in response to sympathetic provocation. A link to arrhythmogenesis is established for concomitantly elevated levels of BVR and its LF oscillations.

1. Background and Objectives

Beat-to-beat variability of repolarization (BVR) is an inherent property of ventricular electrical function (Thomsen et al., 2006; Baumert et al., 2016). When enhanced, this temporal variability has been associated with arrhythmia vulnerability in patients with structural heart disease (Tereshchenko et al., 2009), drug-induced long QT syndrome (Hinterseer et al., 2008), heart failure (Hinterseer et al., 2010), and catecholaminergic polymorphic ventricular tachycardia (Paavola et al., 2015). A link between increased BVR and arrhythmogenesis has been established in a range of animal models as well (Thomsen et al., 2006; Gallacher et al., 2007; Wijers et al., 2018). Various approaches have been proposed in the literature to quantify BVR at the level of the body surface electrocadiogram (ECG), including measurements of QT interval variability (Baumert et al., 2016), T-wave alternans (Verrier et al., 2011), or T-wave morphology variations (Ramirez et al., 2017).

Recent studies have shown that BVR presents a clear low-frequency (LF) oscillatory pattern that can be quantified from the ECG by measuring LF oscillations of the T-wave vector, so-called Periodic Repolarization Dynamics (PRD) (Rizas et al., 2014, 2016). PRD has been shown to be unrelated to heart rate variability or respiratory activity and has been postulated to most likely reflect the effect of phasic sympathetic activity on the ventricular myocardium. Increases in PRD have been associated with destabilization of repolarization leading to ventricular arrhythmias and sudden cardiac death (Rizas et al., 2014, 2017). The described T-wave oscillations have been suggested to reflect oscillations of the ventricular action potential (AP) duration (APD) (Hanson et al., 2014; Rizas et al., 2016; Porter et al., 2018). In in vivo studies on heart failure patients, APD has been shown to indeed oscillate at the same LF range (Hanson et al., 2014). Additional studies have demonstrated that both LF oscillations of APD and BVR are significantly augmented in response to physiologically-induced increased sympathetic activity, with a close interaction between both observed increments (Porter et al., 2017, 2018).

The mechanisms underlying the interactions between BVR and LF patterning of APD in response to sympathetic provocation (SP) and its potential link to arrhythmogenesis remain to be investigated. Regarding BVR, a growing number of studies, both experimental and computational, have provided evidence on the role of ion channel stochasticity and Ca²⁺ cycling variations as underlying mechanisms of temporal variability at different scales, covering from isolated cells (Lemay et al., 2011; Pueyo et al., 2011, 2016a; Antoons et al., 2015; Kistamas et al., 2015; Nánási et al., 2017) to coupled cells / tissue (Zaniboni et al., 2000; Pueyo et al., 2011; Lemay et al., 2011; Magyar et al., 2015; Nánási et al., 2017) to whole heart (Yamabe et al., 2007; Baumert et al., 2016). Furthermore, the action of adrenergic stimulation in modulating those BVR mechanisms and facilitating arrhythmia initiation by the formation of afterdepolarizations and triggered activity has been reported in single cells (Johnson et al., 2010, 2013; Heijman et al., 2013; Szentandrássy et al., 2015; Hegyi et al., 2018) and in the whole heart (Gallacher et al., 2007). In respect of LF oscillations of APD, computational investigations in single cells have suggested that sympathetic nerve activity promotes their generation by both a direct β-adrenergic (βA) action and through the intermediary of mechano-electric feedback (Pueyo et al., 2016b). In the presence of disease-related conditions, like Ca²⁺ overload and reduced repolarization reserve (RRR), these oscillations have been shown to contribute to pro-arrhythmia (Pueyo et al., 2016b).

In the present study, which builds on the work published in Pueyo et al. (2016b), a set of stochastic human ventricular AP models are developed to reproduce the sympathetically-mediated interactions between BVR and LF patterning of APD observed experimentally, to investigate their underlying mechanisms and to establish a link to arrhythmic risk. The developed models are representative of a whole range of AP characteristics and include biophysically detailed descriptions of the electrophysiology, Ca²⁺ dynamics, βA signaling and mechanics of human ventricular cells in health and disease. Stochastic gating of ion channels are incorporated into major currents active during AP repolarization. An approach based on the Automatic Relevance Determination technique (MacKay, 1996) is adopted to unravel the major ionic contributors to augmented BVR and LF oscillations of APD in response to SP, with subsequent analysis of the involved mechanisms. The relationship between the unraveled mechanisms and arrhythmogenesis is established by a methodology grounded on Canonical Correlation Analysis (Hotelling, 1936).

2. Materials and Methods

2.1. Human Data

Previously acquired human data has been described in detail elsewhere (Porter et al., 2018). Briefly, eleven heart failure patients with cardiac resynchronization therapy defibrillator devices had activation recovery intervals (ARIs) recorded from left ventricular epicardial electrodes alongside simultaneous non-invasive blood pressure and respiratory recordings. Heart rate was clamped by right ventricular pacing. Recordings took place during resting conditions and following an autonomic stimulus (Valsalva maneuver). The study was approved by the West London Ethics Committee and conformed to the standards set by the Declaration of Helsinki (latest revision: 64th WMA General Assembly). Informed consent was obtained in writing from all subjects.

2.2. Stochastic Human Ventricular Models

2.2.1. Models of Electrophysiology

The ORd human ventricular epicardial cell model (O'Hara et al., 2011) served as a basis to construct a set of AP models covering a range of experimentally observed electrophysiological characteristics. Each AP model in the dataset, which represents a different virtual cell, was obtained by varying the ionic conductances of the following currents: rapid delayed rectifier K⁺ current, I_Kr; slow delayed rectifier K⁺ current, I_Ks; transient outward K⁺ current, I_to; L-type Ca²⁺ current, I_CaL; inward rectifier K⁺ current, I_K1; sodium current, I_Na; sodium-K⁺ pump current, I_NaK; and sodium-Ca²⁺ exchanger current, I_NaCa. A total of 500 models were initially generated by sampling the nominal conductance values of the ORd model in the range ±100% using the Latin Hypercube Sampling method (McKay et al., 1979; Pueyo et al., 2016b).

Out of all the generated models, only those satisfying the calibration criteria shown in Table 1 were retained. Such criteria were based on experimentally available human ventricular measures of steady-state AP characteristics taken from O'Hara et al. (2011), Guo et al. (2011), Britton et al. (2017), Jost et al. (2008), and Grandi et al. (2010). These characteristics included: APD_90|50, denoting 1 Hz steady-state APD at 90%|50% repolarization (expressed in ms); RMP, standing for resting membrane potential (in mV); V_peak, measuring peak membrane potential following stimulation (in mV); and ΔAPD₉₀, calculated as the percentage of change in APD₉₀ with respect to baseline when selectively blocking I_Ks, I_Kr or I_K1 currents (measured in ms). After applying the described calibration criteria, the initial set of 500 models was reduced to a set of 161 selected models. In addition, models leading to pro-arrhythmic events at baseline conditions were excluded because they did not allow quantification of BVR or LF oscillations of APD, thus resulting in a final population of 123 models. For each of those models, the parameters θ_Ks, θ_Kr, θ_to, θ_CaL, θ_K1, θ_Na, θ_NaCa, and θ_NaK were defined to take the values of the factors multiplying the nominal conductances of I_Ks, I_Kr, I_to, I_CaL, I_K1, I_Na, I_NaK, and I_NaCa, respectively, with respect to the original ORd model, i.e., I_j = θ_jI_{j, ORd}, where I_{j, ORd} represents current j in the ORd model, with j being one of the elements in the set {Ks, Kr, to, CaL, K1, Na, NaCa, NaK}.

TABLE 1

Table 1. Calibration criteria applied onto human ventricular cell models.

Stochasticity was incorporated into the equations describing the ionic gating of four major currents active during AP repolarization, namely I_Ks, I_Kr, I_to, and I_CaL, following the approach described in Pueyo et al. (2011) For a gating variable x, the temporal evolution of the probability of this gate being open was calculated as in Equation (1), where the variance of the stochastic term introduced to formulate the Stochastic Differential Equation (SDE) describing ionic fluctuations was inversely proportional to the number of channels of each species. In Equation (1), x_∞ and τ_x represent the steady-state value of x and the time constant to reach that steady-state value, with x, x_∞ and τ_x being functions of voltage, while w is a Wiener process. The number of channels N associated with each species j were obtained for each virtual cell by multiplying the ionic factor θ_j of that cell by the corresponding number of channels in the ORd model, i.e., N_j = θ_jN_{j, ORd}. Further details on estimation of channel numbers for the ORd model are presented in the Supplementary Material (section 1.1 and Table S1).

$\begin{array}{ll}dx=\frac{x_{\infty }-x}{{\tau }_x}dt+\frac{\sqrt{x_{\infty }+(1-2x_{\infty })x}}{\sqrt{{\tau }_{x}N}}dw.& (1)\end{array}$

2.2.2. Models of PKA Phosphorylation

βA stimulation (βAS) effects were modeled as in Pueyo et al. (2016b) by using a modified version of the Xie et al. (2013) model, with definition of graded and dynamic phosphorylation levels of cellular protein kinase A (PKA) substrates. This model was updated from the original βA signaling formulation proposed in Soltis and Saucerman (2010) to slow down the I_Ks phosphorylation and dephosphorylation rate constants to fit experimental observations. PKA-mediated phosphorylation of phospholemman (PLM) was accounted for by increasing the Na⁺-K⁺-ATPase (NKA) affinity for the intracellular Na⁺ concentration, as in Xie et al. (2013). RyR phosphorylation was described in this study following the formulation proposed in Heijman et al. (2011).

2.2.3. Models of Electromechanical Coupling

An extended version of the Niederer model (Niederer et al., 2006), adjusted to human cell characteristics, as in Weise and Panfilov (2013) and Pueyo et al. (2016b), was used for the electromechanical coupling model. The current through stretch-activated channels (SACs), I_SAC, was introduced as in Pueyo et al. (2016b), with the total current obtained as the sum of the current through K⁺-selective and non-specific cationic SACs. Further details can be found in the Supplementary Material.

2.2.4. Simulation of Baseline and Sympathetic Provocation

A 0.1 Hz periodic stepwise dose of the βA agonist isoproterenol (ISO) was simulated, in accordance with the pattern of muscle sympathetic nerve activity in humans (Pagani et al., 1997). For the first half of the simulated ISO period, the ISO dose was set to either 0.01 μM, for simulated baseline conditions, or 1 μM, for simulated SP, while it was 0 μM for the second half in both cases. Additionally, phasic changes in hemodynamic loading accompanying enhanced sympathetic activity were simulated at the same 0.1 Hz frequency by varying the stretch ratio following a sinusoidal waveform with a maximum change of 1% for baseline conditions and 10% for SP. Sympathetically induced changes in βAS and hemodynamic loading were considered to be in-phase with each other. A total of 640 beats (320 for baseline and 320 for SP) were simulated while pacing the cells at 1 Hz frequency. Figure S1 illustrates simulation of βAS and stretch effects at baseline and in response to sympathetic provocation, while Figure S2 illustrates the APD time series of a cell in the generated population in response to the simulated protocol. For comparison purposes, additional simulations were run under constant βAS and/or hemodynamic loading.

2.2.5. Simulation of Disease-Related Conditions

On top of simulating physiological conditions, models describing disease conditions were built by including representations of: Reduced Repolarization Reserve (RRR), defined by simultaneous blockades of I_Kr and I_Ks currents; and Ca²⁺ overload, defined by increases in the extracellular Ca²⁺ levels. In both cases, an approach like the one described in Pueyo et al. (2016b) was used. Mild disease conditions were simulated by a 1.5-fold increment in the extracellular Ca²⁺ concentration and 7.5% and 20% inhibitions of I_Kr and I_Ks currents, respectively. Moderate disease conditions were simulated by a 2.5-fold increment in the extracellular Ca²⁺ concentration and 22.5% and 60% inhibitions of I_Kr and I_Ks currents, respectively. Severe disease conditions were simulated by a 4-fold increment in the extracellular Ca²⁺ concentration, 30% and 80% inhibitions of I_Kr and I_Ks currents, respectively, and by additionally increasing the conductance of non-specific cationic SACs as described in Isenberg et al. (2003) (G_{SAC, ns} changed from 0.006 nS/pF for physiological, mild and moderate disease conditions to 0.01 nS/pF for severe disease conditions). Table S2 summarizes how physiological as well as mild, moderate and severe disease conditions were simulated in this study.

2.3. Measurements of Repolarization Variability

For each of the developed AP models, APD at 90% repolarization, denoted as APD in the following, was calculated for every beat of the stochastic realizations. A triangulation measure (T1) was calculated as the difference between APD at 90% and 50% repolarization. The last L = 120 beats of each condition (baseline and SP) were used for evaluation of measures describing BVR and LF oscillatory behavior. Averages of those measures over stochastic realizations were computed.

2.3.1. Beat-to-Beat Variability of Repolarization

The following BVR measures were evaluated:

• Standard deviation of APD over the last L beats:

$\begin{array}{ll}{m}_{\text{SD}}=\sqrt{\frac{1}{L-1}{\displaystyle \sum _{l=1}^L}{(\text{APD}(l)-\overline{\text{APD}})}^2}& (2)\end{array}$

where $\overline{\mathrm{\text{APD}}}$ is the average APD over those L beats.

• Normalized variance of APD over the last L beats:

$\begin{array}{ll}{m}_{\text{NSD}}=\frac{m_{\text{SD}}^2}{{\overline{\text{APD}}}^2}.& (3)\end{array}$

• Short-Term Variability (STV) of APD, defined as the average distance perpendicular to the identity line in the Poincaré plot, computed as the average over windows of L_win = 30 beats sliding every one beat along the last L = 120 simulated beats:

$\begin{array}{ll}{m}_{\text{STV}}=\frac{1}{L-L_{win}+1}\left({\displaystyle \sum _{l=1}^{L-L_{win}+1}}{\displaystyle \sum _{i=l}^{l+L_{win}-1}}\frac{|\text{APD}(i+1)-\text{APD}(i)|}{(L_{win}-1)\sqrt{2}}\right).& (4)\end{array}$

• Normalized STV:

$\begin{array}{ll}{m}_{\text{NSTV}}=\frac{m_{\text{STV}}^2}{{\overline{\text{APD}}}^2}.& (5)\end{array}$

2.3.2. Low-Frequency Repolarization Variability

Spectral analysis was performed to compute LF variability measures following the methodology described in Porter et al. (2018). The APD time series of the last L = 120 beats, for either baseline or SP, was linearly detrended. Power Spectral Density (PSD) was estimated after fitting an autoregressive model to the detrended APD time series using the Yule-Walker method. The optimal order of the autoregressive model was chosen in the range between L/3 and L/2 to minimize Akike's Information Criterion, with a requisite on the residuals to pass a whiteness test. Two measures were extracted from the estimated PSD:

• LF power (m_PLF), calculated as the integral of the PSD over the [0.04, 0.15] Hz band.

• Normalized LF power (m_NPLF): LF power normalized by the total power in the [0.04, 0.5] Hz frequency band.

2.4. Contributors to BVR and LF Oscillations

Automatic Relevance Determination (ARD) was used to unravel individual and common factors, in the form of ionic conductance levels, contributing to BVR and LF oscillations of APD in response to SP. ARD is a Bayesian sparsity method, first proposed in the context of neural network models (MacKay, 1996), which has been successfully used to determine the relevance of various input features to given measures (see e.g., Rasmussen and Williams, 2006).

In a regression problem where an output variable (in this case, a BVR or LF oscillatory measure) is aimed to be predicted by several input variables (in this case, the conductances of ionic currents), it commonly happens that some of the variables are irrelevant to the prediction. However, when a finite dataset is analyzed, random correlations between the irrelevant inputs and the output are always obtained, diminishing the capability of the techniques employed for the prediction. A method like ARD, able to infer which input variables are relevant and prune all the irrelevant ones, is advantageous. ARD works by adjusting multiple weight constants, one associated with each input, which are inferred from the data and automatically set to be large for the relevant features and small for the irrelevant ones. The fact that ARD renders a sparse set of explanatory variables makes its results more interpretable than for other correlation-based methods (see e.g., Gunn and Kandola, 2002 for the relation between sparsity and interpretability).

Each virtual cell n out of the N simulated models was considered as a data point determined by its D = 8 parameters (factors multiplying ionic conductances). Those factors were stacked in a row vector x⁽ⁿ⁾ = [${\theta }_{\mathrm{\text{Ks}}}^{(n)},{\theta }_{\mathrm{\text{Kr}}}^{(n)},{\theta }_{\mathrm{\text{to}}}^{(n)},{\theta }_{\mathrm{\text{CaL}}}^{(n)},{\theta }_{\mathrm{\text{K1}}}^{(n)},{\theta }_{\mathrm{\text{Na}}}^{(n)},{\theta }_{\mathrm{\text{NaCa}}}^{(n)},{\theta }_{\mathrm{\text{NaK}}}^{(n)}$], representing the feature vector of each data point. All data were stacked in the feature matrix X, i.e., X = [x⁽¹⁾; … ;x^(N)]. Hence an element of X, denoted as x_{n, i}, was the value of the i-th conductance parameter of virtual cell n. In addition, we used y as a wildcard to denote the column vector with the values of the analyzed variability measure for the data points. Hence, the values in y can either correspond to a temporal BVR measure or a measure of the magnitude of APD LF oscillations: m_SD, m_NSD, m_STV, m_NSTV, m_PLF and m_NPLF. To simplify the training process of the algorithm, the values of y were standardized to zero mean and unit variance. Using this input-output definition we posed the following regression model.

$\begin{array}{ll}{y}^{(n)}=f({\text{x}}^{(n)})+r^{(n)}& (6)\end{array}$

where r⁽ⁿ⁾ is additive random Gaussian noise with variance ${\sigma }_r^2$ and f is a function linking the inputs and the outputs. Typical choices for f include linear, polynomial or neural network functions, with the ones most extensively used by the Bayesian learning community being Gaussian Processes (Rasmussen and Williams, 2006), which represent a powerful and flexible non-parametric option:

$\begin{array}{ll}f({\text{x}}^{(n)})~{GP}\left(m({\text{x}}^{(n)}),c({\text{x}}^{(n)},{\text{x}}^{(n^{\prime })})\right)& (7)\end{array}$

where m(x⁽ⁿ⁾) is the mean function and c(x⁽ⁿ⁾, x^(n′)) is the covariance function between data points n and n′. In its simplest form, m(x⁽ⁿ⁾) = 0 and all the complexity of the model is captured by the covariance function. The covariance is commonly described by linear, polynomial or radial basis functions, or other more complicated functions (see e.g., Rasmussen and Williams, 2006). In this work, a linear function was used for the covariance:

$\begin{array}{ll}c({\text{x}}^{(n)},{\text{x}}^{(n^{\prime })})={\displaystyle \sum _{i=1}^D}{\sigma }_{d,i}^2{x}_{n,i}{x}_{n^{\prime },i}.& (8)\end{array}$

Considering this choice, f(x⁽ⁿ⁾) can be shown to define a set of linear functions with respect to x⁽ⁿ⁾, where directions (i.e., the different factors contained in each x⁽ⁿ⁾) are weighted according to ${\sigma }_{d,i}^2$.

ARD was applied to optimize type II Maximum Likelihood (ML-II) with respect to ${\sigma }_{d,i}^2$ and ${\sigma }_r^2$. Specifically, a quasi-Newton method (in the case of our implementation, L-BFGS, see e.g., Boyd and Vandenberghe, 2004) was used to find the values of the hyperparameters leading to maximization of the following function:

$\begin{array}{l}{L}({\sigma }_{d,1}^2,\cdots ,{\sigma }_{d,8}^2,{\sigma }_r^2)\hfill \\ =\frac{1}{2}logdet{\text{C}}_{\text{ext}}({\displaystyle {\sigma }_{d,1}^2},\cdots ,{\sigma }_{d,8}^2,{\sigma }_r^2)\hfill \\ +\frac{1}{2}{\text{y}}^T{\text{C}}_{\text{ext}}{({\sigma }_{d,1}^2,\cdots ,{\sigma }_{d,8}^2,{\sigma }_r^2)}^{-1}\text{y}+\frac{N}{2}log(2\pi )\hfill & (9)\end{array}$

where ${\text{C}}_{ext}({\sigma }_{d,1}^2,\cdots ,{\sigma }_{d,8}^2,{\sigma }_r^2)=\text{C}({\sigma }_{d,1}^2,\cdots ,{\sigma }_{d,8}^2)+{\sigma }_r^2\text{I}$, with I being the identity matrix and $\text{C}({\sigma }_{d,1}^2,\cdots ,{\sigma }_{d,8}^2)$ being the matrix obtained by evaluating the covariance function c(x⁽ⁿ⁾, x^(n′)) for every pair of data points in X. To avoid overfitting, ten-fold cross validation was applied. Results are presented after averaging the ten corresponding values for each ${\sigma }_{d,i}^2$. The higher the value of ${\sigma }_{d,i}^2$, the more relevant the i-th factor (input parameter) is for the prediction.

This methodology allows establishing which factors are more relevant to predict a given output measure (i.e., a BVR or LF oscillatory measure). In the following, these relevance values are presented as normalized values so that they add up to one to facilitate assessment of the relative relevance of each factor. Since relevance factors do no account for the sign of the contribution, that is, whether an increase in the BVR or LF oscillation measure corresponds to upregulation or downregulation of an ionic current, the Gaussian Process regression was interpreted as a linear regression where the covariance matrix is C_ext and the sign of each contribution was calculated as

$\begin{array}{ll}{\text{s}}_{{\theta }_i}=\text{sign}({({\text{C}}_{\text{ext}}^{-1}\mathbf{\text{X}})}^T\mathbf{\text{y}})& (10)\end{array}$

where θ_i is each of the conductance parameters and T denotes matrix transposition.

Finally, to address the fact that a factor may only be relevant in association with another one, the same methodology was applied after removing one factor (ionic conductance) at a time. If after removing a specific factor, the relevance associated with another factor was found to be significantly changed, a tight relationship between the effects of the two factors was postulated and common mechanisms underlying such a relationship were explored.

This method is implemented in Python 3 using the GPy, Gaussian Process Toolbox (see Sheffield ML group, 2012) and is available in Data Sheet 2 of the Supplementary Material (section 1.6).

2.5. Contributors to Arrhythmogenesis

Canonical Correlation Analysis (CCA) (Hotelling, 1936; Hardoon et al., 2004) was used to identify the ionic conductances with the largest contribution to the occurrence of arrhythmogenic events under simulated diseased conditions. This method has been widely used in several different applications (see e.g., Torres et al., 2007; Kaya et al., 2014; Zhu et al., 2014 for some representative examples).

Similarly to the description of ARD above, the data were stacked in the feature matrix X, with x_{n, i}, being the value of the i-th factor for virtual cell n. A binary vector z of length N was generated, which contained a value of 1 in the positions corresponding to virtual cells for which pro-arrhythmic events were observed following SP and 0 otherwise.

Given X and z, CCA was applied to compute the values of the canonical variables w_x and w_z such that:

$\begin{array}{ll}({\text{w}}_x^{*},w_z^{*})=arg{\displaystyle \underset{{\text{w}}_x,w_z}{max}}\text{corr}(\text{X}{\text{w}}_x,\text{z}{w}_z)& (11)\end{array}$

with corr being the linear correlation between the projected versions of X and z, i.e., Xw_x, zw_z. The elements of vector ${\text{w}}_x^{*}$ represent the projection of ionic factors into a subspace common with $\text{z}{w}_z^{*}$ and can be interpreted as the correlations of each of these factors with the presence of pro-arrhythmic events. Hence, the higher the value of an element in ${\text{w}}_x^{*}$, the higher the relevance of such factor to the events in z.

3. Results

3.1. Sympathetic Provocation Increases BVR and LF Oscillations of APD

Figure 1 shows representative examples of zero-mean time series of experimental ARI (ARI - $\overline{\mathrm{\text{ARI}}}$, with $\overline{\mathrm{\text{ARI}}}$ denoting temporal mean of ARI, left panel) and simulated APD (APD - $\overline{\mathrm{\text{APD}}}$, with $\overline{\mathrm{\text{APD}}}$ denoting temporal mean of APD, right panel) and corresponding PSDs at baseline and following SP. In both experiments and simulations, a remarkable increase in BVR in response to SP can be clearly appreciated from the APD series. Also, the experimental and simulated spectra corresponding to SP show notably more marked peaks in the LF band as compared to baseline.

FIGURE 1

Figure 1. (Left panel) Experimental zero-mean ARI series (ARI - $\overline{\mathrm{\text{ARI}}}$) and corresponding spectra at rest (left) and following Valsalva maneuver (right). (Right panel) Simulated zero-mean APD series (APD - $\overline{\mathrm{\text{APD}}}$) and corresponding spectra at baseline (left) and following sympathetic provocation (right). The LF region of the spectra is shadowed in red and the high frequency region in green.

Of note, the peaks in the high frequency band present in the experimentally recorded data were not analyzed in this study, as vagal or respiratory effects were not included in our simulations for being out of the scope of the present study. The simulated results presented in this and the next sections correspond to simulation of mild disease conditions, since these are compared with experimental results obtained from heart failure patients (see section 2.1). Results for physiological conditions remained qualitatively unchanged with respect to those shown for mild disease conditions.

Figure 2 shows relative measures of BVR and LF oscillations at baseline and following SP for each individual of the experimental and simulated datasets (the cases shown in Figure 1 are highlighted in blue). For the vast majority of individuals, m_NSD and m_NPLF increased in response to augmented sympathetic activity. Importantly, both the level of BVR and LF oscillations as well as the magnitude of change in response to SP presented a high degree of variation between individuals, as shown in Figure 2. As expected, the m_NSD values in the simulations were higher than in the experiments, as simulations correspond to single epicardial cells while experimental data is from left ventricular epicardial electrograms and, thus, includes the effects of intercellular coupling acting to mitigate cell-to-cell variability.

FIGURE 2

Figure 2. Left: Normalized variance m_NSD (top) and Normalized LF power m_NPLF (bottom) at rest and following Valsalva maneuver calculated from experimental ARI series. Right: m_NSD (top) and m_NPLF (bottom) at baseline and following SP calculated from simulated APD series. The cases presented in Figure 1 are highlighted in blue.

In both experiments and simulations, the sympathetically-mediated increases in BVR and LF oscillations were confirmed either when quantified in absolute terms by m_SD, m_STV and m_PLF or in relative terms by m_NSD, m_NSTV and m_NPLF.

3.2. There Is a Close Interaction Between BVR and LF Oscillations of APD, Particularly in Response to Sympathetic Provocation

Table 2 shows correlation values between measures of BVR and LF oscillations of APD, both calculated using absolute and normalized indices. As can be seen in Table 2, the LF power of APD, m_PLF, was highly correlated with BVR measured by the short-term variability of APD, m_STV, and, even to a larger extent, by the standard deviation of APD, m_SD. This observation held true when the correlation was evaluated both at baseline and in response to SP. The strong association found between BVR and LF oscillations of APD in our SP simulations was in line with the one measured experimentally, where the Spearman correlation coefficient between m_PLF and m_SD was 0.679.

TABLE 2

Table 2. Spearman correlation coefficients between simulated BVR and LF oscillation measures.

When normalized measures were considered, Table 2 shows that the correlation between the normalized LF power of APD, m_NPLF, and the normalized BVR measures, m_NSTV and m_NSD, was notably reduced. This highlights the relevance of absolute APD values in modulating the interactions between BVR and LF oscillations of APD. The reduction in correlation after considering normalized measures was particularly so for baseline conditions, while following SP there was still a high interaction between normalized BVR and LF oscillations of APD.

Figure S3 illustrates the simulated relationships between the absolute measures m_PLF and m_SD and between the relative measures m_NPLF and m_NSD at baseline and in response to SP.

Based on the fact that the two ways of evaluating BVR, i.e., by standard deviation and by short-term variability of APD, led to very similar outcomes in terms of the relationship with LF oscillations of APD, the results in the next sections will be shown for m_STV and its normalized counterpart m_NSTV. For APD oscillatory behavior, m_PLF and m_NPLF will be used.

3.3. K⁺ and Ca²⁺ Current Densities Are Common Modulators of BVR and LF Oscillations of APD

Figure 3 illustrates the major contributors to the values of m_STV, m_NSTV, m_PLF, and m_NPLF found in our simulated population in response to SP. The sign of the relationship between the contributing ionic current conductances and the evaluated BVR or LF oscillation measurements was negative in all relevant cases, meaning that downregulation of the ionic current density led to an increment in the analyzed measurement. Note that each bar in the graphs of Figure 3 represents relative relevance with respect to the other evaluated factors, all adding up to one. According to the results in Figure 3, m_STV and m_PLF shared the same major contributors to their observed values following SP. Specifically, the three ionic conductances with the most relevant role in determining the values of m_STV and m_PLF were those of I_Kr, I_K1, and I_CaL currents. For the normalized measurements m_NSTV and m_NPLF, a substantial reduction in the relevance of I_Kr conductance was observed with respect to that quantified for the non-normalized measurements. I_K1 and I_CaL current conductances remained as the two most relevant contributors to the values of m_NSTV and m_NPLF following SP.

FIGURE 3

Figure 3. Relevance of ionic current conductances to m_PLF (A), m_NPLF (B), m_STV (C) and m_NSTV (D), calculated from simulated APD series under SP.

To assess potential associations between ionic conductances in their contributions to the evaluated BVR and LF oscillations measures, the same ARD technique was applied after removing one ionic conductance at a time. For the majority of cases, the computed relevance levels were highly similar after such removals, meaning that there is no co-dependency in the contribution of the different ionic conductances. However, when I_K1 conductance was removed from the analysis, the relevance of other repolarization currents, like I_Kr and I_Ks, in their contribution to m_NSTV was notably increased. This increment reveals common mechanisms in the contributions of all these repolarization currents to the m_NSTV values following SP.

Since the same ionic conductances were found to modulate BVR and LF oscillations of APD following SP, simulations were ran in which βAS and stretch were modeled as constant, with assigned values corresponding to the maximal effects in the above simulations. As can be seen in Figure S4, in those cases I_Kr and I_K1 were still the major modulators of BVR whereas the contribution of I_CaL was drastically decreased. Thus, I_CaL modulation of BVR was mediated by the increment in the LF oscillations of APD, while the role of I_Kr and I_K1 as modulators of BVR did not present such a strong dependence.

For healthy conditions, results were essentially the same as those shown in Figure 3 for mild disease conditions, with only a slight decrease in the relevance of I_NaCa contribution to m_PLF. This is illustrated in Figure S5.

3.4. Modulation of BVR and LF Oscillations of APD by K⁺ and Ca²⁺ Current Densities Is Explained by Their Effects on Ionic Gating Stochasticity, βAS, and Hemodynamic Loading

Before describing the mechanisms by which I_Kr, I_K1, and I_CaL current densities modulate BVR and LF oscillatory measures following SP, the differential effects of the two components associated with enhanced sympathetic activity, namely βAS and mechanical stretch, to such measures were analyzed. Figure 4 illustrates the variations in BVR and LF oscillation measurements in the simulated population for different scenarios, including combined phasic βAS and mechanical stretch, only phasic βAS, only phasic mechanical stretch and only phasic mechanical stretch without SACs. Results showed that the largest contribution to LF oscillations, measured either by m_PLF or m_NPLF, was caused by phasic mechanical stretch, particularly when SACs were included in the models. Regarding BVR, both effects contributed to m_STV and m_NSTV, even if not in an additive manner and with the contribution of βAS being larger than that of mechanical stretch. Additional effects associated with stochastic ionic gating of currents active during AP repolarization added to the BVR values presented in Figure 4.

FIGURE 4

Figure 4. Distributions of BVR and LF oscillation measurements for simulated scenarios including individual and combined βAS and mechanical stretch effects, with and without the contribution of SACs.

3.4.1. Mechanisms Underlying the Role of I_K1 as a Modulator of BVR and LF Oscillations of APD

The role of I_K1 current density as a modulator of APD oscillatory behavior following SP was only relevant when phasic mechanical stretch was simulated and particularly so when SACs were included in the models. The mechanism of action was as follows. Downregulation of I_K1 increased resting membrane potential (Figure 5A) and this increment was associated with an enhancement of the total I_SAC current in the zenith of the oscillation, where phasic stretch reached maximal values (Figures 5B,C). These effects altered the AP shape at the end of the repolarization phase (Figure 5D) and this, in turn, had an impact on the calculated APD. In particular, the magnitude of the APD oscillations was amplified (Figure 5E), which led to increases in both m_PLF and m_NPLF (Figure 5F).

FIGURE 5

Figure 5. (A) Resting membrane potential vs. I_K1 current conductance in the population of virtual cells; (B) I_SAC current for two examples corresponding to upregulated and downregulated I_K1 while keeping all the other currents at their default values in the ORd model; (C) Minimum I_SAC current value vs. I_K1 current conductance in the population; (D) AP traces and (E) zero-mean APD series (APD - $\overline{\mathrm{\text{APD}}}$) for the examples in (B); (F) m_NPLF values vs. I_K1 current conductance in the population.

Furthermore, I_K1 current density had an impact on modulating BVR following SP, especially when including the effects of SACs. Specifically, the above described alterations in AP morphology induced by I_K1 downregulation, manifested as a slowing down of the final part of AP repolarization, rendered the AP more sensitive to the effects of stochastic ionic gating. This led to increased variability in APD values of consecutive beats, thus enlarging m_STV and m_NSTV.

3.4.2. Mechanisms Underlying the Role of I_Kr as a Modulator of BVR and LF Oscillations of APD

The impact of I_Kr current density on the magnitude of BVR and LF oscillations of APD was related to modulation of AP repolarization duration. This is evidenced by the fact that the contribution of I_Kr conductance was very relevant in the modulation of m_PLF and m_STV but was notably reduced for their normalized counterparts m_NPLF and m_NSTV.

In the case of m_PLF, the mechanism of action was as follows. I_Kr downregulation led to AP prolongation, which in our simulations including phasic βAS and stretch could be seen as an increase in both the minimum and the average APD within each oscillation period (Figure 6A). The observed AP lengthening correlated with an increment in the magnitude of the APD oscillations, quantified by the APD range (Figure 6B). This was the result of amplified effects of βAS and stretch on the prolonged AP. In relation to the amplified oscillation amplitude, m_PLF was increased. Representative examples are shown in Figure 6C, where the case with longer APD induced by downregulated I_Kr was associated with larger LF oscillations.

FIGURE 6

Figure 6. (A) Minimum APD vs. I_Kr current conductance in the population of virtual cells; (B) APD range vs. minimum APD; (C) zero-mean APD series (APD - $\overline{\mathrm{\text{APD}}}$) for two examples corresponding to downregulated and upregulated I_Kr while keeping all the other currents at their default values in the ORd model.

In the case of m_STV, the lengthening of AP repolarization induced by I_Kr downregulation led to more accentuated temporal voltage variations. This occurred under phasic βAS, stretch and the combination of both effects associated with enhanced sympathetic activity.

3.4.3. Mechanisms Underlying the Role of I_CaL as a Modulator of BVR and LF Oscillations of APD

The contribution of I_CaL to BVR and LF oscillations was relevant under both simulated βAS and mechanical stretch, with an important role of SACs in explaining I_CaL modulation of APD oscillations.

I_CaL downregulation shortened the AP plateau, leading to more triangular APs (Figure 7A). This, in turn, magnified the effects of phasic βAS and accentuated the APD differences within each simulated oscillation period. This change produced an increase in the magnitude of LF oscillations of APD, associated with increments in both m_PLF and m_NPLF (Figure 7B). Representative examples of low and high BVR and LF oscillations of APD related to up- and downregulation of I_CaL current are presented in Figure 7C. In close correspondence with the above described mechanisms, the more triangular AP induced by I_CaL downregulation facilitated larger voltage fluctuations. This was seen as increased m_STV and m_NSTV.

FIGURE 7

Figure 7. (A) Average triangulation vs. I_CaL current conductance and (B) m_PLF values vs. average triangulation for the population of virtual cells; (C) zero-mean APD series (APD - $\overline{\mathrm{\text{APD}}}$) for two examples corresponding to downregulated and updownregulated I_CaL while keeping all the other currents at their default values in the ORd model.

Under simulated mechanical stretch on top of βAS, there was an additional change in the amplitude and duration of intracellular and subspace Ca²⁺ concentrations as well as in the I_SAC current. All these effects modified the AP repolarization morphology, enhancing the differences within each simulated oscillation period. As a consequence, m_PLF and m_NPLF were further increased and, correspondingly, m_STV and m_NSTV too.

3.5. Severe Disease Conditions Accentuate Both BVR and LF Oscillations of APD, Leading to Electrical Instabilities

Disease conditions simulated by Ca²⁺ overload and RRR had an impact on sympathetically-mediated BVR and LF oscillations of APD. Specifically, when severe disease conditions were simulated, including also an associated increase in the conductance of non-specific cationic SACs, pro-arrhythmic events could be observed. These occurred in 35% of the cases in our population and took the form of early afterdepolarizations (EADs), EAD bursts and spontaneous beats. Examples are presented in Figure 8.

FIGURE 8

Figure 8. Pro-arrhythmic events observed following SP in cells under simulated severe disease conditions.

For those cases where arrhythmogenic events were observed under severe disease conditions (denoted as subpopulation A), BVR and LF oscillations of APD were increasingly accentuated for higher levels of disease conditions, as illustrated in Figure 9. As can be noted from the figure, m_NSTV and m_NPLF took larger values for progressively higher levels of Ca²⁺ overload and RRR. Similarly occurred for the non-normalized indices m_STV and m_PLF. Those cases not presenting arrhythmogenic events under severe disease conditions (denoted as subpopulation NA) showed lower values of BVR and LF oscillation measures for both mild and moderated disease conditions. This can be appreciated in Figure 9 as well.

FIGURE 9

Figure 9. Violin plots representing the distributions of log(m_PLF), m_NPLF, m_STV, and m_NSTV for mild and moderate disease conditions. The whole population of models is divided into two subpopulations: the set of cells presenting (denoted by A) and not presenting (denoted by NA) pro-arrhythmic events following SP under severe disease conditions.

The results of Canonical Correlation Analysis (CCA) performed to assess major contributors to pro-arrhythmic events under severe disease conditions are presented in Figure 10. According to these results, the ionic currents with a major involvement in pro-arrhythmicity were I_Kr, I_CaL, I_K1, and I_NaK, the first three being major modulators of BVR and LF oscillations of APD. The sign of the relationship between ionic conductances and pro-arrhythmicity was negative (i.e., current downregulation facilitating pro-arrhythmic events) in all cases except for I_CaL.

FIGURE 10

Figure 10. Relevance of ionic current conductances to pro-arrhythmic events.

The role of I_Kr, I_CaL, and I_K1 in contributing to pro-arrhythmicity is further illustrated in Figure 11, which shows the distribution of virtual cells as a function of their I_Kr, I_CaL, and I_K1 conductances (θ_Kr, θ_CaL, and θ_K1, respectively). As can be appreciated, pro-arrhythmic cells were most commonly located in regions with low θ_Kr and θ_K1, thereby exemplifying how I_Kr and I_K1 downregulation contribute to pro-arrhythmicity. The effect of I_CaL was only significant in the region where θ_Kr < 1, implying that the role of I_CaL was dependent on I_Kr expression. The information needed to reproduce Figure 11 is available in Data Sheet 1 of the Supplementary Material (section 1.5).

FIGURE 11

Figure 11. Location of cells presenting (red) and not presenting (gray) pro-arrhythmic events under simulated severe disease conditions as a function of relevant ionic current conductance values.

4. Discussion

A population of human ventricular stochastic AP models was built and shown to reproduce a range of responses in terms of BVR and LF oscillations of APD following enhanced sympathetic activity, as reported experimentally (Porter et al., 2018). The models included descriptions of electrophysiology, βA signaling, mechanics and ionic gating stochasticity and served to investigate the interactions between the two investigated phenomena, namely temporal variability and LF oscillatory behavior of APD, following sympathetic provocation. Ionic mechanisms underlying inter-individual differences in those phenomena were ascertained and individual characteristics associated with concomitantly large beat-to-beat variability and LF oscillations of repolarization were established. These were linked to higher susceptibility to electrical instabilities in the presence of disease conditions like Ca²⁺ overload and RRR.

4.1. Relationship Between Sympathetically-Mediated BVR and LF Oscillations of APD in a Human Ventricular Population

Increases in LF oscillations of repolarization in response to enhanced sympathetic activity have been described at the level of the electrocardiographic T-wave and QT interval in humans and animals (Negoescu et al., 1997; Rizas et al., 2014, 2016) and at the level of the ventricular APD in ambulatory patients (Hanson et al., 2014; Porter et al., 2018). A direct effect related to enhanced activity of the sympathetic nerves innervating ventricular myocardium, rather than just an effect attributable to heart rate variability, has been proved (Negoescu et al., 1997; Rizas et al., 2014; Porter et al., 2018). In this study, phasic βAS and mechanical stretch were simulated in association with muscle sympathetic nerve activity patterns during enhanced sympathetic activity (Pagani et al., 1997). Pacing at a constant rate was applied to the models. In accordance with experimental observations, increments in absolute and normalized LF power of APD have been overall measured in our population. Nevertheless, there is a high degree of inter-individual variability, with some individual cases showing no change or even a decrease in LF oscillations of APD in response to SP, which is in line with experimental reports as well.

Additionally, clinical and experimental studies have reported that enhanced sympathetic activity leads to increased BVR in patients with the long QT syndrome type 1 (Satomi et al., 2005) and animal models of this disease (Gallacher et al., 2007) as well as in heart failure patients (Porter et al., 2017). Our human ventricular AP models, by including stochastic expressions of ionic current gating, allowed investigation of BVR at baseline and in response to SP. In agreement with experimental evidences, most of the models in our diseased population have shown sympathetically-mediated increments in BVR. The increase in BVR in the referred experimental/clinical studies as well as in our simulations of disease could be explained by βAS effects under conditions of reduced IKs, which is indeed the case in our simulations and in long QT syndrome type 1 investigations and could also be the case in heart failure following previous reports suggesting downregulation of this current in failing hearts (Long et al., 2015). Also, mechanical effects associated with increased sympathetic activity could synergistically enhance BVR. Furthermore, in our simulations, a wide range of individual behaviors in terms of BVR patterns could be characterized following SP, in line with experimental data.

The interactions between BVR and LF oscillations of APD have been recently investigated in ambulatory patients with heart failure following a standard sympathetic provocation maneuver (Porter et al., 2018). In the present study, a strong correlation between BVR and LF oscillation measures has been measured as well by simulation of SP through phasic βAS and mechanical stretch in human ventricular myocytes. This holds true for physiological conditions and for disease conditions, simulated by Ca²⁺ overload and RRR, which are characteristic of diseased hearts like those of heart failure patients. In both simulations and experiments the variability measurements m_SD, m_NSD, m_PLF and m_PLF were quantified. In addition, the BVR measurement m_STV, which accounts for information on the APD variation between consecutive beats and has been extensively used for arrhythmic risk prediction (Thomsen et al., 2004; Hinterseer et al., 2010), was included in this study together with its APD-normalized version m_NSTV.

The strong correlation between m_STV and m_PLF found in simulations and experiments can be explained in light of our simulation outcomes. On the one hand, an increment in temporal APD variability associated with random ionic gating directly augments the LF power of APD, as it induces a rise in the power of APD at all frequencies. Although the measurement m_NPLF normalizes m_PLF by the total power, this marker turns out to be more insensitive to the amplitude of the LF oscillations of APD than m_PLF, while still indicative of the presence or absence of such oscillatory behavior. In the case of BVR, the normalized measurement m_NSTV has been quantified on top of m_STV to correct for the dependence on the APD. Even if the applied APD correction is able to reduce the correlation between APD and m_NSTV to a good extent, it does not abolish it completely. The very strong correlation between m_PLF and m_STV, both at baseline and following SP, dropped to very low correlation when m_NPLF and m_NSTV were evaluated at baseline. Following SP, the correlation between m_NPLF and m_NSTV was still remarkable, which can be explained by the fact that the presence of a marked LF oscillatory pattern directly impacts the temporal APD variability by increasing beat-to-beat APD differences.

4.2. Main Contributors to Increased BVR and LF Oscillations of APD Following Enhanced Sympathetic Activity

The tight relationship between BVR and LF oscillations of APD following enhanced sympathetic activity suggests there could be common modulators of both phenomena. By building a population of virtual cells representing a range of experimentally reported characteristics, in this study it was possible to elucidate the ionic current conductances with a major contribution to inter-individual differences in absolute (m_STV and m_PLF) and normalized (m_NSTV and m_NPLF) BVR and LF oscillation markers. For such elucidation, an approach based on the Automatic Relevance Determination (ARD) technique was developed. Similar approaches have been proposed in the context of magnetoencephalography (Nummenmaa et al., 2007) and wireless communications (Jacobs, 2012), among others, but to the best of our knowledge this is the first time an ARD-based technique is used to identify ionic modulators of cardiac electrophysiological phenomena.

In Pueyo et al. (2016b) the mechanisms underlying LF oscillations of ventricular APD were investigated by simulating phasic βAS and mechanical stretch in association with enhanced sympathetic activity. Differential I_Ks and I_CaL phosphorylation and dephosphorylation kinetics in response to βAS together with variations in Ca²⁺ cycling and SACs in response to stretch were found to synergistically underlie LF oscillatory behavior under SP. While that study provided meaningful insights into the bases for LF oscillations of ventricular repolarization, only an average cell was modeled, which did not allow investigation of inter-individual differences in LF oscillations of APD as in the present study. Also, the models of the population built here are stochastic, as opposed to the deterministic models employed in Pueyo et al. (2016b), thus allowing to quantify BVR at baseline and its change in response to SP. This is of major relevance for investigation of the interactions between BVR and LF oscillations of APD and of their modulators in a whole population.

Our results highlighted the relevance of I_Kr, I_CaL, and I_K1 conductances in modulating inter-individual differences in both BVR and LF oscillatory pattern of APD under SP. Regarding I_Kr, its downregulation was shown to be a key factor for augmentation of m_STV and m_PLF but less important when considering their normalized counterparts m_NSTV and m_NPLF. Concerning LF oscillations of APD, there is little investigation in the literature into factors acting to modulate their magnitude. In Pueyo et al. (2016b), a reduction in the repolarization current was shown to amplify APD oscillatory behavior. Our results are in line with such observations. Considering the fact that m_NPLF does not reflect the magnitude of the oscillations but mostly its presence or absence, this normalized marker was found not to be modulated by I_Kr. Regarding BVR, a variety of experimental, clinical and computational studies have addressed the role of ionic current conductances in modulating beat-to-beat temporal variability quantified by markers such as m_STV or m_SD. In accordance with the results presented in Pueyo et al. (2011) for baseline conditions and Heijman et al. (2013) for βAS, our study has shown I_Kr downregulation to act as a contributor of BVR magnification. Since such a contribution is to a large extent mediated by APD lengthening, it becomes importantly reduced when measured by markers that include APD normalization, such as m_NSTV or m_NSD.

Another very relevant current in the modulation of BVR and LF oscillatory behavior of APD was I_CaL. Although no previous studies in the literature have investigated the role of I_CaL as a modulator of LF oscillation amplitude, there have been a number of studies addressing its role as a modulator of BVR. In Lemay et al. (2011), I_CaL downregulation was shown to increase the random channel fluctuation effects in guinea pig models, which is in good agreement with our presented results. On top of the contribution of I_CaL, a role for I_Ks and persistent I_Na currents in enhancing BVR was also demonstrated in Lemay et al. (2011). We could not find such a role for those two currents, which could be due to differences between species [guinea pig in Lemay et al. (2011) and human in this study] and to the fact that this study investigated conditions of enhanced sympathetic activity rather than baseline conditions.

Regarding I_K1 regulation, this is, to the best of our knowledge, the first study identifying its relevance to BVR and LF oscillations of APD. In our results, I_K1 downregulation appears as a relevant contributor when SACs are incorporated into the models to simulate mechanical stretch changes associated with SP. Under donwregulated I_K1, SACs contribute to alter the AP shape during the last part of repolarization in a phasic manner, leading to increments in both BVR and LF oscillations.

As chronotropic effects of sympathetic provocation have been well documented in in vivo studies, computational simulations were additionally carried out while pacing the virtual cells at higher frequencies. The main ionic contributors I_Kr, I_K1, and I_CaL are confirmed to remain very relevant to explain inter-individual differences in BVR and LF oscillatory behavior in response to SP. Of note, the relevance of I_NaK in determining LF oscillations of APD increases when the analysis is performed for pacing frequencies above 1 Hz.

4.3. Pro-arrhythmic Events Associated With Increased BVR and LF Oscillations of APD Under Severe Disease Conditions

Ca²⁺ overload and RRR are properties commonly present in diseased hearts, like those of patients with heart failure, ischemic heart disease or post-myocardial infarction (Dhalla and Temsah, 2001; Sridhar et al., 2008; Varró and Baczkó, 2011; Guo et al., 2012; Nissen et al., 2012; Gorski et al., 2015). In this study, BVR and LF oscillations of APD have been found to become increasingly accentuated in response to disease progression. These results are in line with those reported in previous clinical, experimental and theoretical studies of the literature. In isolated myocytes and animal models of diseases like diabetes, heart failure or post-myocardial infarction, exaggerated temporal APD variability has been observed in association with Ca²⁺ overload and RRR (Maltsev et al., 2007; Sridhar et al., 2008; Wu et al., 2008; Meo et al., 2016). In the long QT syndrome type 1, involving loss of I_Ks function, elevated ventricular repolarization variability in response to βAS has been documented and mechanisms have been proposed based on animal models, isolated myocytes and computer simulation research (Gallacher et al., 2007; Johnson et al., 2010, 2013; Heijman et al., 2013). In chronic atrioventricular block dogs, where ventricular remodeling importantly compromises repolarization reserve, beat-to-beat APD variability has been found to be augmented with respect to healthy dogs (Stams et al., 2016); an observation also confirmed at the level of ventricular myocytes (Antoons et al., 2015). A mechanical challenge in the form of preload variability has been reported to be essential in that augmentation, with mechano-electrical feedback through stretch-activated channels (SACs) postulated as a major mechanism (Stams et al., 2016). In Pueyo et al. (2016b), the presence of disease conditions has been reported to lead to notably augmented LF oscillations of APD.

Under severe disease conditions, arrhythmogenic manifestations have been found to arise in individual cases of our population presenting large temporal repolarization variability, either quantified at the LF band (LF oscillations) or at all frequencies (BVR). These observations are in agreement with studies relating disproportionate APD fluctuations, particularly in response to enhanced sympathetic activity, and the generation of afterdepolarizations and arrhythmias. In Gallacher et al. (2007) the authors used an in vivo canine model of the long QT syndrome type 1 to demonstrate that βAS enhanced temporal and spatial variability of ventricular repolarization, which precipitated Torsades de Pointes (TdP) arrhythmias. The association between increased BVR and the onset of TdP arrhythmias has also been demonstrated in dogs with chronic atrioventricular block (Thomsen et al., 2006; Wijers et al., 2018). In ventricular myocytes and wedge preparations from human end-stage failing hearts, βAS has been shown to generate electrical abnormalities that result in EADs and delayed afterdepolarizations (DADs) (Veldkamp et al., 2001; Lang et al., 2016). Using a rabbit model mimicking electrophysiological and contractile alterations in human HF, βAS has been reported to be a key factor in inducing DADs and increasing the propensity for triggered arrhythmias (Pogwizd et al., 2001). At the level of the surface ECG, increased BVR and LF oscillations of repolarization have been shown to be risk predictors of ventricular arrhythmias and sudden cardiac death (Wu et al., 2008; Rizas et al., 2014, 2017; Baumert et al., 2016).

Provided the tight relationship between magnification of BVR and LF oscillations of APD and pro-arrhythmic risk, the existence of common modulators has been explored in the present study. Canonical Correlation Analysis has been proposed to identify ionic factors contributing to pro-arrhythmic risk following enhanced sympathetic activity. CCA revealed the important role of I_K1, I_Kr, and I_CaL in the development of pro-arrhythmic events. These same factors are those primarily involved in modulation of sympathetically-mediated BVR and LF oscillations of APD. The role of I_K1 in contributing to arrhythmogenesis has been reported in a rabbit model of heart failure, where the combination of upregulated I_NaCa, downregulated I_K1 and residual βA responsiveness has been shown to increase the propensity for triggered arrhythmias (Pogwizd et al., 2001). In our study, the contribution of I_K1 downregulation to pro-arrhythmicity in association with elevated temporal variability might have been more prominent if our population of stochastic AP models had been built based on an electrophysiological model more likely producing delayed afterdepolarizations under downregulated I_K1, and possibly upregulated I_NaCa, as compared to the ORd model. The role of I_Kr in arrhythmogenesis has been well established in a variety of previously published investigations. In Sridhar et al. (2008), the loss of repolarizing currents, including I_Kr, has been described to lead to increased BVR, repolarization instability and afterdepolarizations in myocytes from dogs susceptible to sudden cardiac death. In Pueyo et al. (2016b) reduced I_Kr and I_Ks have been reported to cause AP irregularities associated with enhanced LF oscillations of APD induced by sympathetic provocations. The implications of I_Kr inhibition in promoting ventricular arrhythmias associated with increased temporal APD dispersion has been further demonstrated in animal models of disease (Stams et al., 2016). On top of K⁺ currents, the present work has identified I_CaL current as another relevant contributor to pro-arrhythmia associated with elevated BVR and LF oscillations of APD, even if conditioned to the presence of reduced I_Kr. In line with these results, increased I_CaL has been demonstrated to facilitate electrical abnormalities in the form of EADs in ventricular myocytes from human failing hearts (Veldkamp et al., 2001). The contribution of increased I_CaL to arrhythmogenesis during βAS has been also shown in Johnson et al. (2013) under reduced I_Ks.

In this study, other currents, like I_Ks and I_NaL, were found to have minor relevance as contributors to arrhythmogenesis in association with temporal dispersion of repolarization. This is contrast to previous studies showing major roles of I_Ks downregulation and I_NaL upregulation (Undrovinas et al., 2006; Gallacher et al., 2007; Maltsev et al., 2007; Wu et al., 2008; Johnson et al., 2010, 2013; Heijman et al., 2013). This discrepancy may be explained by differences between species, modeling characteristics and, importantly, investigated conditions, since this study has focused on the investigation of arrhythmic events occurring following enhanced sympathetic activity.

4.4. Limitations

The stochastic models built in this study included random gating descriptions for major ionic currents active during AP repolarization like I_Ks, I_Kr, I_to, and I_CaL, as in previous studies of the literature (Pueyo et al., 2016a; Tixier et al., 2017). Future studies could include stochasticity in other currents like I_NaL, whose contribution to BVR has been reported in canine ventricular models (Heijman et al., 2013).

In the present work the ORd ventricular AP model has been used, which was developed based on extensive undiseased human data. In this model the effect of varying the I_Ks current on AP is significantly smaller than in other human ventricular cell models, like the ten Tusscher-Panfilov model (ten Tusscher and Panfilov, 2006). The low relevance of I_Ks as an ionic modulator of BVR and LF oscillations of APD found in this work may have to do with it. In Pueyo et al. (2016b), which served as a starting point for the present work, several electrophysiological, mechanical and adrenergic signaling models were tested and only some quantitative differences could be found, while the conclusions remained qualitatively the same for all models. Nevertheless, the role of certain ionic currents in modulating inter-individual differences in BVR and LF oscillatory behavior, as investigated in this study, might still be different if another AP model were used as a basis. This should be addressed in future works.

Also in relation to the use of the ORd model as a basis for the development of the population of models in our study, it should be noted that other ventricular AP models with updated mechanisms of Ca²⁺ induced Ca²⁺ release could provide additional insight into the occurrence of spontaneous Ca²⁺ release and delayed afterdepolarizations in association with elevated BVR and LF oscillations of APD. Indeed, previous studies have been shown that Ca²⁺ handling abnormalities are a major driver of BVR during βAS (Johnson et al., 2013) and a link between Ca²⁺ handling and arrhythmia liability during increased sympathetic activity has been demonstrated, particularly in the setting of heart failure (Johnson and Antoons, 2018).

The population of human ventricular cells used in this work was generated by varying the conductances of eight ionic currents. Ionic parameters other than maximal current conductances might also represent relevant mechanisms underlying the interactions between BVR and LF oscillations of APD. In particular for the I_CaL current, previous studies have proved that modulation of other biophysical properties, like a reduction in the amplitude of the non-inactivating pedestal component of I_CaL, allows to effectively suppress EADs without blocking peak I_CaL, thus preserving excitation-contraction coupling (Madhvani et al., 2015). Future work could address the investigations of the present study by generating a population of virtual cells where biophysical ionic parameters other than maximal conductances were varied, which could eventually lead to findings that developed into more clinically useful therapeutic approaches.

The present study has focused on single cells, while the available experimental data on the interactions between BVR and LF oscillations of human APD are from in vivo measurements in ambulatory heart failure patients. Simulated results qualitatively reproduced the behavior observed in the experiments. Future work could include assessment of those interactions in tissue and whole-heart models. Nonetheless, cell-to-cell coupling has been shown to be remarkably reduced in heart failure and other disease conditions, which would render cell and tissue results close to each other.

Statistical approaches based on ARD and CCA have been used in this study. Future works could investigate generalization of these techniques to consider nonlinear relationships by using kernel functions, even if a larger number of simulations would be required to avoid overfitting.

5. Conclusions

Human ventricular models including descriptions of cell electrophysiology, ion channel stochasticity, β-adrenergic signaling and mechanical stretch were developed. These models reproduced experimentally reported interactions between beat-to-beat variability and low-frequency oscillations of repolarization in response to enhanced sympathetic activity. Ionic factors underlying correlated increments in both phenomena were investigated, which included downregulation of the inward and rapidly activating delayed rectifier K⁺ currents and the L-type Ca²⁺ current. Concomitantly elevated levels of beat-to-beat repolarization variability and its low-frequency oscillations in diseased ventricles led to electrical instabilities and arrhythmogenic events. This investigation serves as a basis for future studies aiming at improving arrhythmic risk stratification and guiding the search for more efficient anti-arrhythmic therapies.

Author Contributions

EP and PT devised the project, the main conceptual ideas and proof outline, and were responsible for overseeing the research and providing critical insight and recommendations regarding the focus, structure, and content of the paper. DS-P and JF-B performed computational simulations and analyzed the data results. BP and SvD contributed with technical details and analysis support. All authors participated in writing and proofreading throughout the publication process.

Funding

This work was supported by the European Research Council through grant ERC-2014-StG 638284, by MINECO (Spain) through project DPI2016-75458-R, by MULTITOOLS2HEART-ISCIII, by Gobierno de Aragón (Reference Group BSICoS T39-17R) cofunded by FEDER 2014-2020, by European Social Fund (EU) and Gobierno de Aragón through project LMP124_18 and a personal predoctoral grant to DS-P. Computations were performed by ICTS NANBIOSIS (HPC Unit at University of Zaragoza).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2019.00916/full#supplementary-material

Data Sheet 1. Population of virtual cells indicating the pro-arrhythmic models distribution (Figure-11-master.zip). This Data Sheet contains the ionic factors associated with each virtual cell, both for the group presenting and the group not presenting pro-arrhythmic events (as shown in Figure 11).

Data Sheet 2. Automatic Relevance Determination (automatic-relevance-master.zip). This Data Sheet presents the Python code used to unravel individual and common factors, in form of ionic conductance levels, contributing to Beat-to-beat Variability of Repolarization and Low-Frequency Oscillations.

Data Sheet 3. This document contains additional information about the stochastic ORd model and the stretch-activated channels and βAS formulation. In addition, details on the simulation of healthy and disease conditions are provided. Additional figures are included to facilitate understanding.

References

Antoons, G., Johnson, D. M., Dries, E., Santiago, D. J., Ozdemir, S., Lenaerts, I., et al. (2015). Calcium release near L-type calcium channels promotes beat-to-beat variability in ventricular myocytes from the chronic AV block dog. J. Mol. Cell. Cardiol. 89, 326–334. doi: 10.1016/j.yjmcc.2015.10.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Baumert, M., Porta, A., Vos, M. A., Malik, M., Couderc, J.-P., Laguna, P., et al. (2016). QT interval variability in body surface ECG: measurement, physiological basis, and clinical value: position statement and consensus guidance endorsed by the European Heart Rhythm Association jointly with the ESC Working Group on Cardiac Cellular Electrophysiology. Europace 18, 925–944. doi: 10.1093/europace/euv405

PubMed Abstract | CrossRef Full Text | Google Scholar

Boyd, S., and Vandenberghe, L. (2004). Convex Optimization. Cambridge: Cambridge University Press.

Google Scholar

Britton, O. J., Bueno-Orovio, A., Virág, L., Varró, A., and Rodriguez, B. (2017). The electrogenic Na+/K+ pump is a key determinant of repolarization abnormality susceptibility in human ventricular cardiomyocytes: a population-based simulation study. Front. Physiol. 8:278. doi: 10.3389/fphys.2017.00278

CrossRef Full Text | Google Scholar

Dhalla, N. S., and Temsah, R. M. (2001). Sarcoplasmic reticulum and cardiac oxidative stress: an emerging target for heart disease. Emerging Ther. Targets 5, 205–217. doi: 10.1517/14728222.5.2.205

PubMed Abstract | CrossRef Full Text | Google Scholar

Gallacher, D. J., Van de Water, A., van der Linde, H., Hermans, A. N., Lu, H. R., Towart, R., et al. (2007). In vivo mechanisms precipitating Torsades de Pointes in a canine model of drug-induced long-QT1 syndrome. Cardiovas. Res. 76, 247–256. doi: 10.1016/j.cardiores.2007.06.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Gorski, P. A., Ceholski, D. K., and Hajjar, R. J. (2015). Altered myocardial calcium cycling and energetics in heart failure a rational approach for disease treatment. Cell Metab. 21, 183–194. doi: 10.1016/j.cmet.2015.01.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Grandi, E., Pasqualini, F. S., and Bers, D. M. (2010). A novel computational model of the human ventricular action potential and Ca transient. J. Mol. Cell. Cardiol. 48, 112–121. doi: 10.1016/j.yjmcc.2009.09.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Gunn, S. R., and Kandola, J. S. (2002). Structural modelling with sparse kernels. Mach. Learn. 48, 137–163. doi: 10.1023/A:1013903804720

CrossRef Full Text | Google Scholar

Guo, D., Liu, Q., Liu, T., Elliott, G., Gingras, M., Kowey, P. R., et al. (2011). Electrophysiological properties of HBI-3000: a new antiarrhythmic agent with multiple-channel blocking properties in human ventricular myocytes. J. Cardiovasc. Pharmacol. 57, 79–85. doi: 10.1097/FJC.0b013e3181ffe8b3

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, X., Gao, X., Wang, Y., Peng, L., Zhu, Y., and Wang, S. (2012). Iks protects from ventricular arrhythmia during cardiac ischemia and reperfusion in rabbits by preserving the repolarization reserve. PLoS ONE 7:e31545. doi: 10.1371/journal.pone.0031545

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanson, B., Child, N., Van Duijvenboden, S., Orini, M., Chen, Z., Coronel, R., et al. (2014). Oscillatory behavior of ventricular action potential duration in heart failure patients at respiratory rate and low frequency. Front. Physiol. 5:414. doi: 10.3389/fphys.2014.00414

PubMed Abstract | CrossRef Full Text | Google Scholar

Hardoon, D. R., Szedmak, S., and Shawe-Taylor, J. (2004). Canonical correlation analysis: an overview with application to learning methods. Neural Comput. 16, 2639–2664. doi: 10.1162/0899766042321814

PubMed Abstract | CrossRef Full Text | Google Scholar

Hegyi, B., Bányász, T., Izu, L. T., Belardinelli, L., Bers, D. M., and Chen-Izu, Y. (2018). β-adrenergic regulation of late Na⁺ current during cardiac action potential is mediated by both PKA and CaMKII. J. Mol. Cell. Cardiol. 123, 168–179. doi: 10.1016/j.yjmcc.2018.09.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Heijman, J., Volders, P. G., Westra, R. L., and Rudy, Y. (2011). Local control of β-adrenergic stimulation: effects on ventricular myocyte electrophysiology and Ca²+-transient. J. Mol. Cell. Cardiol. 50, 863–871. doi: 10.1016/j.yjmcc.2011.02.007

CrossRef Full Text | Google Scholar

Heijman, J., Zaza, A., Johnson, D. M., Rudy, Y., Peeters, R. L., Volders, P. G., et al. (2013). Determinants of beat-to-beat variability of repolarization duration in the canine ventricular myocyte: a computational analysis. PLoS Comput. Biol. 9:e1003202. doi: 10.1371/journal.pcbi.1003202

PubMed Abstract | CrossRef Full Text | Google Scholar

Hinterseer, M., Beckmann, B.-M., Thomsen, M. B., Pfeufer, A., Ulbrich, M., Sinner, M. F., et al. (2010). Usefulness of short-term variability of qt intervals as a predictor for electrical remodeling and proarrhythmia in patients with nonischemic heart failure. Am. J. Cardiol. 106, 216–220. doi: 10.1016/j.amjcard.2010.02.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Hinterseer, M., Thomsen, M. B., Beckmann, B.-M., Pfeufer, A., Schimpf, R., Wichmann, H.-E., et al. (2008). Beat-to-beat variability of QT intervals is increased in patients with drug-induced long-QT syndrome: a case control pilot study. Eur. Heart J. 29, 185–190. doi: 10.1093/eurheartj/ehm586

PubMed Abstract | CrossRef Full Text | Google Scholar

Hotelling, H. (1936). Relations between two sets of variates. Biometrika 28, 321–377. doi: 10.1093/biomet/28.3-4.321

CrossRef Full Text | Google Scholar

Isenberg, G., Kazanski, V., Kondratev, D., Gallitelli, M. F., Kiseleva, I., and Kamkin, A. (2003). Differential effects of stretch and compression on membrane currents and [Na⁺]_c in ventricular myocytes. Prog. Biophys. Mol. Biol. 82, 43–56. doi: 10.1016/S0079-6107(03)00004-X

CrossRef Full Text | Google Scholar

Jacobs, J. P. (2012). Bayesian support vector regression with automatic relevance determination kernel for modeling of antenna input characteristics. IEEE Trans. Antennas Propag. 60, 2114–2118. doi: 10.1109/TAP.2012.2186252

CrossRef Full Text | Google Scholar

Johnson, D. M., and Antoons, G. (2018). Arrhythmogenic mechanisms in heart failure: linking β-adrenergic stimulation, stretch, and calcium. Front. Physiol. 9:1453. doi: 10.3389/fphys.2018.01453

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnson, D. M., Heijman, J., Bode, E. F., Greensmith, D. J., van der Linde, H., Abi-Gerges, N., et al. (2013). Diastolic spontaneous calcium release from the sarcoplasmic reticulum increases beat-to-beat variability of repolarization in canine ventricular myocytes after β-adrenergic stimulation. Circ. Res. 112, 246–256. doi: 10.1161/CIRCRESAHA.112.275735

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnson, D. M., Heijman, J., Pollard, C. E., Valentin, J.-P., Crijns, H. J., Abi-Gerges, N., et al. (2010). Iks restricts excessive beat-to-beat variability of repolarization during beta-adrenergic receptor stimulation. J. Mol. Cel. Cardiol. 48, 122–130. doi: 10.1016/j.yjmcc.2009.08.033

CrossRef Full Text | Google Scholar

Jost, N., Varro, A., Szuts, V., Kovacs, P. P., Seprényi, G., Biliczki, P., et al. (2008). Molecular basis of repolarization reserve differences between dogs and man. Circulation 118:S342.

Google Scholar

Kaya, H., Eyben, F., Salah, A. A., and Schuller, B. (2014). “Cca based feature selection with application to continuous depression recognition from acoustic speech features,” in 2014 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP) (Florence: IEEE), 3729–3733. doi: 10.1109/ICASSP.2014.6854298

CrossRef Full Text | Google Scholar

Kistamas, K., Szentandrassy, N., Hegyi, B., Vaczi, K., Ruzsnavszky, F., Horvath, B., et al. (2015). Changes in intracellular calcium concentration influence beat-to-beat variability of action potential duration in canine ventricular myocytes. J. Physiol. Pharmacol. 66, 73–81.

PubMed Abstract | Google Scholar

Lang, C. N., Menza, M., Jochem, S., Franke, G., Perez Feliz, S., Brunner, M., et al. (2016). Electro-mechanical dysfunction in long QT syndrome: role for arrhythmogenic risk prediction and modulation by sex and sex hormones. Prog. Biophys. Mol. Biol. 120, 255–269. doi: 10.1016/j.pbiomolbio.2015.12.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Lemay, M., de Lange, E., and Kucera, J. P. (2011). Effects of stochastic channel gating and distribution on the cardiac action potential. J. Theor. Biol. 281, 84–96. doi: 10.1016/j.jtbi.2011.04.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Long, V. P., Bonilla, I. M., Vargas-Pinto, P., Nishijima, Y., Sridhar, A., Li, C., et al. (2015). Heart failure duration progressively modulates the arrhythmia substrate through structural and electrical remodeling. Life Sci. 123, 61–71. doi: 10.1016/j.lfs.2014.12.024

PubMed Abstract | CrossRef Full Text | Google Scholar

MacKay, D. J. (1996). “Bayesian methods for Backpropagation networks,” in Models of Neural Networks III (New York, NY: Springer), 211–254.

Google Scholar

Madhvani, R. V., Angelini, M., Xie, Y., Pantazis, A., Suriany, S., Borgstrom, N. P., et al. (2015). Targeting the late component of the cardiac L-type Ca²⁺ current to suppress early afterdepolarizations. J. Gen. Physiol. 145, 395–404. doi: 10.1085/jgp.201411288

PubMed Abstract | CrossRef Full Text | Google Scholar

Magyar, J., Banyasz, T., Szentandrassy, N., potential, K., P Nanasi, P., and Satin, J. (2015). Role of gap junction channel in the development of beat-to-beat action potential repolarization variability and arrhythmias. Curr. Pharm. Des. 21, 1042–1052. doi: 10.2174/1381612820666141029102443

PubMed Abstract | CrossRef Full Text | Google Scholar

Maltsev, V. A., Silverman, N., Sabbah, H. N., and Undrovinas, A. I. (2007). Chronic heart failure slows late sodium current in human and canine ventricular myocytes: implications for repolarization variability. Eur. J. Heart Fail. 9, 219–227. doi: 10.1016/j.ejheart.2006.08.007

PubMed Abstract | CrossRef Full Text | Google Scholar

McKay, M. D., Beckman, R. J., and Conover, W. J. (1979). A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 21:239. doi: 10.2307/1268522

CrossRef Full Text | Google Scholar

Meo, M., Meste, O., Signore, S., Sorrentino, A., Cannata, A., Zhou, Y., et al. (2016). Reduction in Kv current enhances the temporal dispersion of the action potential in diabetic myocytes: insights from a novel repolarization algorithm. J. Am. Heart Assoc. 5:e003078. doi: 10.1161/JAHA.115.003078

PubMed Abstract | CrossRef Full Text | Google Scholar

Nánási, P. P., Magyar, J., Varró, A., and Ördög, B. (2017). Beat-to-beat variability of cardiac action potential duration: underlying mechanism and clinical implications. Can. J. Physiol. Pharmacol. 95, 1230–1235. doi: 10.1139/cjpp-2016-0597

PubMed Abstract | CrossRef Full Text | Google Scholar

Negoescu, R., Dinca-Panaitescu, S., Filcescu, V., Ionescu, D., and Wolf, S. (1997). Mental stress enhances the sympathetic fraction of QT variability in an RR-independent way. Integr. Physiol. Behav. Sci. 32, 220–227. doi: 10.1007/BF02688620

PubMed Abstract | CrossRef Full Text | Google Scholar

Niederer, S. A., Hunter, P. J., and Smith, N. P. (2006). A quantitative analysis of cardiac myocyte relaxation: a simulation study. Biophys. J. 90, 1697–1722. doi: 10.1529/biophysj.105.069534

PubMed Abstract | CrossRef Full Text | Google Scholar

Nissen, J. D., Thomsen, M. B., Bentzen, B. H., Diness, J. G., Diness, T. G., Jespersen, T., et al. (2012). Attenuated ventricular β-adrenergic response and reduced repolarization reserve in a rabbit model of chronic heart failure. J. Cardiovasc. Pharmacol. 59, 142–150. doi: 10.1097/FJC.0b013e318238727a

PubMed Abstract | CrossRef Full Text | Google Scholar

Nummenmaa, A., Auranen, T., Hämäläinen, M. S., Jääskeläinen, I. P., Sams, M., Vehtari, A., et al. (2007). Automatic relevance determination based hierarchical Bayesian MEG inversion in practice. Neuroimage 37, 876–889. doi: 10.1016/j.neuroimage.2007.04.021

PubMed Abstract | CrossRef Full Text | Google Scholar

O'Hara, T., Virág, L., Varró, A., and Rudy, Y. (2011). Simulation of the undiseased human cardiac ventricular action potential: model formulation and experimental validation. PLoS Comput. Biol. 7:e1002061. doi: 10.1371/journal.pcbi.1002061

PubMed Abstract | CrossRef Full Text | Google Scholar

Paavola, J., Väänänen, H., Larsson, K., Penttinen, K., Toivonen, L., Kontula, K., et al. (2015). Slowed depolarization and irregular repolarization in catecholaminergic polymorphic ventricular tachycardia: a study from cellular Ca²⁺ transients and action potentials to clinical monophasic action potentials and electrocardiography. Eurospace 18, 1599–1607. doi: 10.1093/europace/euv380

PubMed Abstract | CrossRef Full Text

Pagani, M., Montano, N., Porta, A., Malliani, A., Abboud, F. M., Birkett, C., et al. (1997). Relationship between spectral components of cardiovascular variabilities and direct measures of muscle sympathetic nerve activity in humans. Circulation 95, 1441–1448. doi: 10.1161/01.CIR.95.6.1441

PubMed Abstract | CrossRef Full Text | Google Scholar

Pogwizd, S. M., Schlotthauer, K., Li, L., Yuan, W., and Bers, D. M. (2001). Arrhythmogenesis and contractile dysfunction in heart failure: roles of sodium-calcium exchange, inward rectifier potassium current, and residual β-adrenergic responsiveness. Circ. Res. 88, 1159–1167. doi: 10.1161/hh1101.091193

PubMed Abstract | CrossRef Full Text | Google Scholar

Porter, B., Bishop, M. J., Claridge, S., Behar, J., Sieniewicz, B. J., Webb, J., et al. (2017). Autonomic modulation in patients with heart failure increases beat-to-beat variability of ventricular action potential duration. Front. Physiol. 8:328. doi: 10.3389/fphys.2017.00328

PubMed Abstract | CrossRef Full Text | Google Scholar

Porter, B., Van Duijvenboden, S., Bishop, M. J., Orini, M., Claridge, S., Gould, J., et al. (2018). Beat-to-beat variability of ventricular action potential duration oscillates at low frequency during sympathetic provocation in humans. Front. Physiol. 9:147. doi: 10.3389/fphys.2018.00147

PubMed Abstract | CrossRef Full Text | Google Scholar

Pueyo, E., Corrias, A., Virág, L., Jost, N., Szél, T., Varró, A., et al. (2011). A multiscale investigation of repolarization variability and its role in cardiac arrhythmogenesis. Biophys. J. 101, 2892–2902. doi: 10.1016/j.bpj.2011.09.060

PubMed Abstract | CrossRef Full Text | Google Scholar

Pueyo, E., Dangerfield, C., Britton, O., Virág, L., Kistamás, K., Szentandrássy, N., et al. (2016a). Experimentally-based computational investigation into beat-to-beat variability in ventricular repolarization and its response to ionic current inhibition. PLoS ONE 11:e0151461. doi: 10.1371/journal.pone.0151461

PubMed Abstract | CrossRef Full Text | Google Scholar

Pueyo, E., Orini, M., Rodríguez, J. F., and Taggart, P. (2016b). Interactive effect of beta-adrenergic stimulation and mechanical stretch on low-frequency oscillations of ventricular action potential duration in humans. J. Mol. Cell. Cardiol. 97, 93–105. doi: 10.1016/j.yjmcc.2016.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Ramirez, J., Orini, M., Tucker, J. D., Pueyo, E., and Laguna, P. (2017). Variability of ventricular repolarization dispersion quantified by time-warping the morphology of the T-waves. IEEE Trans. Biomed. Eng. 64, 1619–1630. doi: 10.1109/TBME.2016.2614899

PubMed Abstract | CrossRef Full Text | Google Scholar

Rasmussen, C. E., and Williams, C. K. (2006). Gaussian Processes for Machine Learning (Cambridge, MA: MIT Press).

Google Scholar

Rizas, K. D., Hamm, W., Kääb, S., Schmidt, G., and Bauer, A. (2016). Periodic repolarisation dynamics: a natural probe of the ventricular response to sympathetic activation. Arrhythm. Electrophysiol. Rev. 5:31. doi: 10.15420/aer.2015:30:2

PubMed Abstract | CrossRef Full Text | Google Scholar

Rizas, K. D., McNitt, S., Hamm, W., Massberg, S., Kääb, S., Zareba, W., et al. (2017). Prediction of sudden and non-sudden cardiac death in post-infarction patients with reduced left ventricular ejection fraction by periodic repolarization dynamics: MADIT-II substudy. Eur. Heart J. 38, 2110–2118. doi: 10.1093/eurheartj/ehx161

PubMed Abstract | CrossRef Full Text | Google Scholar

Rizas, K. D., Nieminen, T., Barthel, P., Zürn, C. S., Kähönen, M., Viik, J., et al. (2014). Sympathetic activity–associated periodic repolarization dynamics predict mortality following myocardial infarction. J. Clin. Invest. 124, 1770–1780. doi: 10.1172/JCI70085

PubMed Abstract | CrossRef Full Text | Google Scholar

Satomi, K., Shimizu, W., Takaki, H., Suyama, K., Kurita, T., Aihara, N., et al. (2005). Response of beat-by-beat QT variability to sympathetic stimulation in the LQT1 form of congenital long QT syndrome. Heart Rhythm 2, 149–154. doi: 10.1016/j.hrthm.2004.11.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Sheffield ML group (2012). GPy: A Gaussian Process Framework in Python. Available online at: http://github.com/SheffieldML/GPy

Soltis, A. R., and Saucerman, J. J. (2010). Synergy between CaMKII substrates and β-adrenergic signaling in regulation of cardiac myocyte Ca2+ handling. Biophys. J. 99, 2038–2047. doi: 10.1016/j.bpj.2010.08.016

CrossRef Full Text | Google Scholar

Sridhar, A., Nishijima, Y., Terentyev, D., Terentyeva, R., Uelmen, R., Kukielka, M., et al. (2008). Repolarization abnormalities and afterdepolarizations in a canine model of sudden cardiac death. Am. J. Physiol. Regul. Integr. Comp. Physiol. 295, R1463–R1472. doi: 10.1152/ajpregu.90583.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Stams, T. R., Oosterhoff, P., Heijdel, A., Dunnink, A., Beekman, J. D., van der Nagel, R., et al. (2016). Beat-to-beat variability in preload unmasks latent risk of Torsade de Pointes in anesthetized chronic atrioventricular block dogs. Circulation 80, 1336–1345. doi: 10.1253/circj.CJ-15-1335

PubMed Abstract | CrossRef Full Text | Google Scholar

Szentandrássy, N., Kistamás, K., Hegyi, B., Horváth, B., Ruzsnavszky, F., Váczi, K., et al. (2015). Contribution of ion currents to beat-to-beat variability of action potential duration in canine ventricular myocytes. Pfluegers Arch. Eur. J. Physiol. 467, 1431–1443. doi: 10.1007/s00424-014-1581-4

PubMed Abstract | CrossRef Full Text | Google Scholar

ten Tusscher, K. H., and Panfilov, A. V. (2006). Alternans and spiral breakup in a human ventricular tissue model. Am. J. Physiol. Heart Circ. Physiol. 291, H1088–H1100. doi: 10.1152/ajpheart.00109.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Tereshchenko, L. G., Fetics, B. J., Domitrovich, P. P., Lindsay, B. D., and Berger, R. D. (2009). Prediction of ventricular tachyarrhythmias by intracardiac repolarization variability analysis. Circ. Arrhyth. Electrophysiol. 2:276. doi: 10.1161/CIRCEP.108.829440

PubMed Abstract | CrossRef Full Text | Google Scholar

Thomsen, M. B., Verduyn, S. C., Stengl, M., Beekman, J. D., de Pater, G., van Opstal, J., et al. (2004). Increased short-term variability of repolarization predicts d-sotalol–induced torsades de pointes in dogs. Circulation 110, 2453–2459. doi: 10.1161/01.CIR.0000145162.64183.C8

PubMed Abstract | CrossRef Full Text | Google Scholar

Thomsen, M. B., Volders, P. G., Beekman, J. D., Matz, J., and Vos, M. A. (2006). Beat-to-beat variability of repolarization determines proarrhythmic outcome in dogs susceptible to drug-induced Torsades de Pointes. J. Am. Coll. Cardiol. 48, 1268–1276. doi: 10.1016/j.jacc.2006.05.048

PubMed Abstract | CrossRef Full Text | Google Scholar

Tixier, E., Lombardi, D., Rodriguez, B., and Gerbeau, J.-F. (2017). Modelling variability in cardiac electrophysiology: a moment-matching approach. J. R. Soc. Interface 14:20170238. doi: 10.1098/rsif.2017.0238

PubMed Abstract | CrossRef Full Text | Google Scholar

Torres, D. A., Turnbull, D., Sriperumbudur, B. K., Barrington, L., and Lanckriet, G. R. (2007). “Finding musically meaningful words by sparse cca,” in Neural Information Processing Systems (NIPS) Workshop on Music, the Brain and Cognition (Vancouver, BC).

Google Scholar

Undrovinas, A. I., Belardinelli, L., Undrovinas, N. A., and Sabbah, H. N. (2006). Ranolazine improves abnormal repolarization and contraction in left ventricular myocytes of dogs with heart failure by inhibiting late sodium current. J. Cardiovas. Electrophysiol. 17, S169–S177. doi: 10.1111/j.1540-8167.2006.00401.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Varró, A., and Baczkó, I. (2011). Cardiac ventricular repolarization reserve: a principle for understanding drug-related proarrhythmic risk. Br. J. Pharmacol. 164, 14–36. doi: 10.1111/j.1476-5381.2011.01367.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Veldkamp, M. W., Verkerk, A. O., Van Ginneken, A. C., Baartscheer, A., Schumacher, C., de Jonge, N., et al. (2001). Norepinephrine induces action potential prolongation and early afterdepolarizations in ventricular myocytes isolated from human end-stage failing hearts. Eur. Heart J. 22, 955–963. doi: 10.1053/euhj.2000.2499

PubMed Abstract | CrossRef Full Text | Google Scholar

Verrier, R. L., Klingenheben, T., Malik, M., El-Sherif, N., Exner, D. V., Hohnloser, S. H., et al. (2011). Microvolt T-wave alternans: physiological basis, methods of measurement, and clinical utility consensus guideline by International Society for Holter and Noninvasive Electrocardiology. J. Am. Coll. Cardiol. 58, 1309–1324. doi: 10.1016/j.jacc.2011.06.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Weise, L. D., and Panfilov, A. V. (2013). A discrete electromechanical model for human cardiac tissue: effects of stretch-activated currents and stretch conditions on restitution properties and spiral wave dynamics. PLoS ONE 8:e59317. doi: 10.1371/journal.pone.0059317

PubMed Abstract | CrossRef Full Text | Google Scholar

Wijers, S. C., Sprenkeler, D. J., Bossu, A., Dunnink, A., Beekman, J. D. M., Varkevisser, R., et al. (2018). Beat-to-beat variations in activation-recovery interval derived from the right ventricular electrogram can monitor arrhythmic risk under anesthetic and awake conditions in the canine chronic atrioventricular block model. Heart Rhythm 15, 442–448. doi: 10.1016/j.hrthm.2017.11.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, L., Guo, D., Li, H., Hackett, J., Yan, G.-X., Jiao, Z., et al. (2008). Role of late sodium current in modulating the proarrhythmic and antiarrhythmic effects of quinidine. Heart Rhythm 5, 1726–1734. doi: 10.1016/j.hrthm.2008.09.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, Y., Grandi, E., Puglisi, J. L., Sato, D., and Bers, D. M. (2013). β-adrenergic stimulation activates early afterdepolarizations transiently via kinetic mismatch of PKA targets. J. Mol. Cell. Cardiol. 58, 153–161. doi: 10.1016/j.yjmcc.2013.02.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamabe, M., Sanyal, S. N., Miyamoto, S., Hadama, T., Isomoto, S., and Ono, K. (2007). Three different bradycardic agents, zatebradine, diltiazem and propranolol, distinctly modify heart rate variability and QT-interval variability. Pharmacology 80, 293–303. doi: 10.1159/000107103

PubMed Abstract | CrossRef Full Text | Google Scholar

Zaniboni, M., Pollard, A. E., Yang, L., and Spitzer, K. W. (2000). Beat-to-beat repolarization variability in ventricular myocytes and its suppression by electrical coupling. Am. J. Physiol. Heart Circ. Physiol. 278, H677–H687. doi: 10.1152/ajpheart.2000.278.3.H677

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, X., Suk, H.-I., and Shen, D. (2014). “Multi-modality canonical feature selection for alzheimer's disease diagnosis,” in International Conference on Medical Image Computing and Computer-Assisted Intervention (Boston, MA: Springer), 162–169.

Google Scholar