For visualization of combined effects of the channels on EAD induction, we performed EAD analysis for varying combinations of block of the three most sensitive non-hERG channels (ICaV, INaL, and IKs) regulating EAD generation, each ranging from 0 to 100% at Cdrug = IC60,hERG. **Figure 2A** represents the EAD test for the models. The red region represents the set of combinations of blocks for these currents that resulted in the EADs in the OHR model (EAD+ region). The region in green covers the parameter subspace where no EADs were observed in the OHR model (EAD− region). Separation of the EAD+ region and EAD− region is outlined by the blue and yellow surfaces for the OHR and OHRmv model, respectively (**Figure 2A**). Among the non-hERG channels, block of ICaV had the highest modulatory effects on EAD generation under normal

simulation of six drugs in the OHR model. Red and green traces show that the drug is defined as TdP+ and TdP−, respectively.

conditions. Under critical hERG block (60%), block of ICaV by more than 30% resulted in suppression of the EADs in both the OHR and OHRmv models (**Figure 2A**). Block of INaL in the absence of the block of ICaV did not result in suppression of the EADs in the OHR model. Even in the OHRmv model with improved INaL formulation, block of more than 60% of the late sodium current was required to suppress pause-induced EADs. Block of IKs currents increased the amount of block of INaL and ICaV currents required to prevent EADs. A similar parametric space for EAD generation was also analyzed under enhanced late sodium currents (**Figure 2B**), associated with LQT3. In the OHRmv model the conductance of the late sodium current was doubled. Block of the hERG channel by 30% was enough to induce EADs in the OHRmv model under enhanced late sodium currents. EAD generation was examined at this 30% hERG block and is shown in **Figure 2B**. For the simulations with enhanced late sodium currents, the difference in the effects of INaL and ICaV was significantly reduced. At 30% hERG block in the OHRmv model, block of greater than 20 and 30% was required for EAD suppression by ICaV and INaL, respectively (**Figure 2B**).

**Table 3** lists the accuracy of TdP risk prediction using EADs as the classification criteria. As an illustration, the data points for the drugs in one of the datasets (Dataset 5) are overlaid on the the parametric space in **Figure 2A** at hERG blocks of 60% to visualize the agreement between drugs EAD and TdP risks. The drugs with positive TdP risk are shown in red, while the drugs with negative TdP risk are shown in green. Actual AP profiles were also simulated for drugs in Datasets 5 and 7 at 60% hERG block concentrations taking into account the block of all seven channels. **Figure 2C** gives representative examples of the AP profiles for the multi-channel block of six of the drugs from Dataset 7 simulated in the OHR model. Our results show a good concordance between drugs torsadogenic risk and its propensity to induce EADs (i.e., most of the torsadogenic drugs resulted in pause-induced EADs in the models while no EADs are observed for majority of the non-torsadogenic drugs at 60% hERG block drug concentrations) across all the datasets. Majority of the datasets (Datasets 1, 2, 4, 6) contains values of blocks for two non-hERG channels (ICaV, INa,fast) of which only ICaV had an impact on EAD occurrence in the OHR and OHRmv models. Hence, for these datasets, OHR

and OHRmv models give identical accuracies as the drugs with greater than 30% ICaV block would result in EAD suppression in both of these models. For the Datasets 3, 5, and 7, druginduced blocks of multiple ion channels are reported. However, very few drugs among these datasets are located in the region between the EAD risk decision surfaces for the OHR (blue surface) and OHRmv (yellow surface) models (i.e., the region in the parameter space that would result in a different prediction between the OHR and OHRmv models, **Figure 2A**). Hence, similar accuracies are observed for both of the models across all datasets. For the Dataset 7, Ranolazine (TdP−) was the only drug that was located on the negative side of the EAD risk decision surface for the OHRmv and to the positive side of the decision surface for the OHR, and hence predicted correctly by the OHRmv model but not by the OHR model. On the contrary for Dataset 5, Quinine (TdP+) ended up on the negative side of the EAD risk separating surface for the OHRmv and to the positive side of the EAD risk decision surface of the OHR model resulting in its incorrect prediction using the OHRmv model (**Figure 2A**).

#### 3.2. Binary TdP Risk Discrimination from Direct Features

Although the biophysical models can provide mechanistic insights underlying TdP genesis, the benefits of using biophysical models in terms of classification is unclear. Here, we wanted to examine the performance of the classifiers built on direct features using the proposed method. The predictive power of the TdP risk classifiers built on the direct features using the proposed method (MCB@EAD) is shown in **Table 4** (two-step classifier column). Classification scores on the direct features at EFTPC are reported for comparison (one-step classifier column). The predictive ability of the two-step classifier was comparable or better than those for the classifiers built on the various features in the original datasets and also better than those for the one-step classifier. Most datasets comprise in-vitro assay data for druginduced block of IKr, INa,fast, and ICaV (Datasets 1, 2, 4, and 6). Drug-induced blocks of additional channels were reported in three datasets (Datasets 3, 5, and 7). The classifiers built using the block of IKr and ICaV as inputs provided high TdP risk prediction scores (**Table 4**). Utilizing the block of INaL as an additional input feature improved the prediction for the two datasets (Dataset 3 and Dataset 7) by classifying correctly one more drug, Ranolazine, as TdP−. On the contrary, in Dataset 5, addition of the block of INaL to the features reduced the number of correctly classified drugs by one (classifying Ritonavir incorrectly as TdP−). Taking into account the block of additional ion channels (INa,fast for Datasets 1, 2, 4, 6 and INa,fast, IKs, Ito, IK1 for Datasets 3, 5 ,7) did not provide any further improvement in the performance of the classifier for all datasets.

For visualization of our two-step approach we presented two-dimensional risk maps. Drug-induced block of the ICaV or the sum of ICaV and INaL blocks is plotted against the hERG ratio ( IC60,hERG EFPTCdrug ) for each drug in a two-dimensional risk map (**Figure 3**). The values of the hERG ratio (x-axis) and the block of ICaV (y-axis) that provide the best discrimination vary widely

TABLE 4 | Binary TdP classifier scores of one-step and two-step classification on the direct features for the eight datasets.

*The total number of drugs in the particular dataset are reported in the NDrugs column. Table also lists the feature used for classification in the original methods and corresponding scores if applicable. The last column shows a list of channel currents used to construct the classifiers. For comparision with the classification accuracy obtained with the direct features, we reported the highest accuracy obtained for the classifiers from the 13 derived features in isolation at different pacing rate/cell type in the bracket in One step and Two step classifier columns. \*The derived features were extracted from simulations of the drug-induced block of all the reported ion channel currents in the in-vitro assay datasets.*

across the datasets (**Figure 3**). Using the threshold of 22 and 100 for the block of ICaV and hERG ratio, respectively, resulted in the perfect classification for Dataset 1 (Mirams et al., 2011). However, for Dataset 2 the hERG ratio threshold of 200 and threshold for ICaV of 57 provided the best classification. Looking at the risk maps in **Figure 3** we can see that the thresholds that provide best classification accuracy vary across datasets resulting in variations in the obtained high risk zones across datasets. For datasets 3, 5, 7 where we also considered block of INaL as one of the input feature in addition to the block of ICaV the value on y-axis of the risk map were reported using Equation (5). A two-dimensional risk map for Dataset 5 for two alternate TdP definitions is shown in **Figure 4**. On the two-dimensional risk map with only ICaV block on the y-axis, we also highlight, using a blue rectangular outline, the separation between the region with EAD presence and absence observed in the in-silico simulation at varying ICaV blocks. For example, at critical hERG block simulations in the model at block of ICaV less than 30% would result in EADs. Block of ICaV by more than 30% results in suppression of the observed EADs. There is approximate correspondence between the EAD observance region in the in-silico model and the high torsadogenic region obtained from classifying the direct-features using machine learning.

# 3.3. Derived vs. Direct Features as Predictors for TdP Risk

We examine the performance of classifiers from various derived features extracted from the simulated action potential and calcium transient in predicting TdP risk. Drug-induced multichannel block was simulated in the OHR model for all the compounds in the merged dataset (Dataset 8) considering the block IKr and ICaV channels. Simulations were also carried out taking into account drug-induced block of all the channels with available IC50 values. Taking into account the block of additional channels resulted in similar accuracies and are reported in the Supplemental Material. **Figure 5** shows the performance of the logistic-regression classifier to discriminate TdP+ and TdP− drugs (**Figures 5A,B**) and to predict drug-induced EADs (**Figures 5C,D**). The classifiers were built on 13 features extracted from the steady-state APs and Ca2+ transients (for 3 cell types at 3 pacing rates). The derived features were extracted at either EFTPC of the drugs (**Figures 5A,C**) or at concentrations at which each drug would produce 60% hERG block (Cdrug = IC60,hERG) (**Figures 5B,D**). Among the various derived features obtained at Cdrugs = EFTPC, diastolic Ca2+ levels provided the best discrimination score between the TdP+ and TdP− drugs (∼79% accuracy at 1 Hz in epi cell type, **Figure 5A**) in agreement with a

et al., 2017), (G) Dataset 7 (*IKr*, *ICaV* ) (Li et al., 2017), (H) Dataset 7 (*IKr*, *ICaV* , *INaL*) (Li et al., 2017), and (I) Dataset 8 (merged dataset). The regions in green are low risk regions and the regions in red are high risk areas. Red dots (•) indicate TdP+ drugs and green dots (•) indicate TdP− drugs. For comparsion purposes, we superimpose a blue rectangular outline that shows the separation between EAD+ and EAD− regions of parameter space. For binary classification the high and intermediate risk drugs in Dataset 7 were assigned TdP+ and the low risk drugs were assigned to TdP−.

previous report (Lancaster and Sobie, 2016). Classifier was also built using APD50 and Diastolic Ca2+ together as inputs (the combination that provided the best prediction in Lancaster and Sobie, 2016) but did not give an improved classification for the merged dataset. Several derived features performed well (>90% accuracy) for EAD risk prediction (**Figure 5**). However, at drug concentrations equal to IC60,hERG (Cdrug = IC60,hERG), each of the 13 features from the Ca2+ transient and AP provided high classification scores (∼85% maximum) for TdP risk assessment (**Figure 5C**). The highest accuracy for each of the features was obtained at different pacing rates (0.5, 1, and 2 Hz) and for different cell types (endo, mid and epi). The maximum classification scores to discriminate the drugs that induced EAD in the model from the drugs that did not induce EADs was 100% (**Figure 5D**) for each of the features. Our results suggest that at fixed hERG block concentrations where trigger events such as EADs arise, several derived features obtained from the model including features from the Ca2+ transient can provide good TdP risk prediction. These derived features also highly correlate with EADs (**Figure 5**). However, the derived features extracted from simulations of the drug-induced effects in ventricular myocyte model did not result in much improvement in TdP risk assessment over the classifiers built on the direct-features using the proposed method (**Table 4**). Moreover, combining the direct and derived features to build the TdP risk classifiers also did not improve the classification performance (results are not shown).

#### 3.4. Direct Features Perform Well to Allow Tertiary Risk Classification

Our results show that classifiers built on the direct features serve as excellent predictors of TdP risk of the drugs categorized into binary risk groups. However, a working committee under the CiPA initiative led by FDA has recently categorized 28 drugs into

tertiary risk categories (low, medium, and high risk compounds) (Colatsky et al., 2016; Fermini et al., 2016). Hence, we test the predictive capabilities of the classifier based on direct features to classify the drugs categorized into tertiary risk categories. **Figures 6A,B** show a two-dimensional risk map for the 12 drugs in Li et al. (2017). These drugs have been analyzed previously using modified OHR model that incorporate dynamic hERG channel interactions. The 12 drugs are a subset of the 28 drugs categorized into three risk categories (CP1, CP2, and CP3) under the CiPA initiative. As a first step, we developed a TdP risk map using only the block of ICaV channel. The hERG ratio threshold of 150 and the threshold of ICaV block of 45%, the values that provided best classification for the merged dataset, were utilized for classifying the low high and intermediate risk drugs from low risk drugs **Figure 6A**. Arbitrary value of hERG ratio and ICaV block of 25 and 15, respectively (a value greater than the maximum hERG ratio and maximum block of ICaV among the high risk drugs in Datasets 7 and 9) was utilized to separate high risk drugs from the low and intermediate risk drugs. Three of the four drugs in low risk (CP3) category were classified correctly. Ranolazine was the only misclassified drug. The boundaries of red zone were defined to include all high risk drugs and hence all the drugs from CP1 category were correctly classified. However, several drugs in intermediate risk were incorrectly classified. Next, we built regression classifier using as input metric the sum of block of ICaV, INaL channels and the degree of drug-trapping parameter which was shown to be essential to improve risk prediction of intermediate risk drugs in the original dataset. **Figure 6B** shows the risk map built using this metric. The threshold along the y-axis was obtained from the regression coefficients. The hERG ratio threshold of 25 and 150 were utilized as before. Including the degree of drug trapping characterized by open-bound/closed-bound ratio for the drugs at steady-state (Li et al., 2017) as one of the features in addition to the blocks of CaV and NaL channels, resulted in the perfect separation of the 12 drugs in 3 categories (**Figure 6B**).

We employed the same approach to test all of the 28 drugs (Fermini et al., 2016) categorized in CP1, CP2, and CP3 categories under the CiPA initiative. In-vitro assays for 12 of these 28 drugs were reported in Crumb et al. (2016) and analyzed in Li et al. (2017). To augment this dataset, we extract the IC50 values for hERG and CaV blocks from Datasets 1, 2, 3, 5, 6, 7 resulting in characterization of 26 of the 28 drugs categorized under the CiPA initiative. Two drugs (Azimilide and Loratidine) were absent in all of the datasets analyzed here and hence were not taken into consideration. We used the mean value of the block if the drug was present in more than one

dataset. The final dataset of the 26 drugs is reported in the Supplemental Material. **Figure 6C** show the two-dimensional risk maps for the 26 drugs with hERG ratio on the x-axis and the block of CaV on the y-axis. The 45% ICaV block threshold and hERG ratio threshold of 150 yielded almost perfect binary classification (high and intermediate vs. low risk drugs) with only 2 (Ranolazine and Tamoxifen) drugs of the 8 from the CP3 category and one drug (Clozapine) of the 16 drugs from CP1 and CP2 cateogry (high and intermediate risk drugs) classifying incorrectly. However, no clear separation was observed amoung the drugs in CP1 and CP2 categories, with several drugs from CP2 category ending up in the high risk region. Considering additional features such as INaL and the degree of drug trapping (if either of the features were not available their value was set to zero) and utilizing the threshold obtained from training the Dataset 7 (**Figure 6B**) resulted in only 3 intermediate risk drugs (Clarithromycin, Domperidone, Droperidol) and 1 low risk drug (Tamoxifen) of the 26 drugs to be misclassified (**Figure 6D**). The classifier already performs well in spite of testing data from heterogeneous sources with some missing values. Further refinement of the method may be possible when a dataset is available with all 28 drugs characterized with a consistent methodology.

#### 3.5. Diverse Definition of Drugs Torsadogenicity Lead to Different Prediction Accuracies

Different binary definitions for the drug's torsadogenic risk have been used across the literature (Lancaster and Sobie, 2016; Ando et al., 2017; Wi´sniowska and Polak, 2017). **Tables 5**, **6** list the different classification accuracy scores obtained for four different binary TdP definitions (Datasets 5 and 8, respectively). Classifiers were constructed on the block of one or multiple-ion channels as inputs at critical hERG block concentrations (Cdrug = IC60,hERG). The various definitions not only resulted in variability of classification scores (**Tables 5**, **6**), but also changed the role of different ion channels in accurate TdP risk classification

dots (•) indicate intermediate risk drugs. Metric = % block CaV + % block NaL + degree of drug trapping.

TABLE 5 | Accuracy scores of TdP classifiers on the direct features for Dataset 5 under four different TdP definitions.

*Target1: TdP*+ = *CM1, CH1, R1, R2, R3, and FDA label QT prolongation and torsade warnings. Target 2: TdP*+ = *CM1 and CM2. Target 3: TdP*+ = *CM1. Target 4: TdP*+ = *CM1 and CM3. Row reporting the maximum accuracies obtained using the derived features is reported for comparision with the direct features.*

(**Table 5**). Using the block of ICaV currents provided the best accuracy scores under two of the four definitions (Target 1 and Target 2) for Dataset 5. Including the effects on additional TABLE 6 | Accuracy scores of TdP classifiers on the direct features for Dataset 8 under four different TdP definitions.

*Target1: TdP*+ = *CM1, CH1, R1, R2, R3, and FDA label QT prolongation and torsade warnings. Target 2: TdP*+ = *CM1 and CM2. Target 3: TdP*+ = *CM1. Target 4: TdP+ = CM1 and CM3. Row reporting the maximum accuracies obtained using the derived features is reported for comparision with the direct features.*

ion-channels did not improve classification scores for these two definitions. Taking into account the block of late sodium currents in addition to CaV channels, provided the best classification accuracy for remaining two of the four TdP definitions (Target 3 and 4) (**Table 5**).

Step2: Derived features

86 76 74 73

# 4. DISCUSSION

Evaluation of drug-induced alterations in multiple cardiac ionchannel currents to determine the drug's torsadogenic potential is currently under investigation through initiatives like CiPA (Comprehensive In-vitro Proarrhythmia Assay) (Sager et al., 2014; Fermini et al., 2016). We have developed a novel twostep method (MCB@EAD) for classification of drugs according to their torsadogenic risk. Using the proposed method, we examined the drug effects at fixed hERG block (i.e., 60% block) concentrations for all tested compounds. This approach allows to isolate the effects of hERG and non-hERG channels in the classification problem. The proximity of the drug's EFTPC to the concentration that results in the critical hERG block provides one of the metrics for determining the drug's TdP risk. For the drugs that induce this critical hERG block at concentrations below a set threshold, the drug-induced effects on non-hERG channels provide an additional metric to determine pro-arrhythmic risk independently of the drug's EFTPC. Our classifier shows improved or equivalent prediction to existing methods. However, one of the advantages compared to previous studies is that the direct and derived features based MCB@EAD classifiers were tested on several in-vitro assay datasets reported previously, as well as on a large composite dataset obtained by merging the different datasets together. One of the important findings of the study was that MCB@EAD TdP classifiers from the direct features provides excellent TdP risk prediction and performs identical to the TdP classifiers from the derived features, which are extracted from complex biophysical models. Although the derived features provided by the biophysical models did not improve the predictive capability for TdP risk assessment, the biophysical models helped determine the amount of block that generates EADs (i.e., the concentration at which the direct features are analyzed using the MCB@EAD classifier). The proposed method not only performs comparably or better than the previous classifiers (**Table 4**) across various in-vitro assay datasets published previously, but also highlights the link between direct and derived feature based classifiers. The results also show strong correlation between the drugs that generate EADs and the drugs with positive TdP risk.

#### 4.1. Ion-Channels Critical for TdP Risk Prediction

Although the role of multiple ion-channels have been suggested for improved TdP risk prediction, classifiers have been primarily built on the blocks of IKr, ICaV, and INa,peak currents (Mirams et al., 2011; Christophe, 2013, 2015; Kramer et al., 2013; Lancaster and Sobie, 2016). A recent assay reports the drug-induced effects on seven ion channels (Crumb et al., 2016) providing an opportunity to identify the ion channels that are important for pro-arrhythmic risk assessment. The results of our parametric simulations of EAD indicate a potential role of block of ICaV, INaL, and IKs currents, in addition to IKr, for determination of torsadogenic risk of the drugs (**Figure 2**). It should be noted that the EAD simulations results are highly dependent on the ventricular myocyte model. For example, block of late sodium current in the OHRmv plays a more prominent role in regulation of AP sensitivity to EADs as compared to the OHR model (**Figure 2A**). The IKs plays a much bigger role in Ten Tusscher and Panfilov model (Ten Tusscher and Panfilov, 2006) than OHR model in regulation of APD as shown in Mirams et al. (2014). The present datasets have limited examples of block of ICaV, INaL, and IKs currents in the same compounds. Moreover, among the non-hERG channels, the regulatory effect of ICaV block was the highest with the block of ICaV by only 30% resulting in EAD suppression at critical hERG current block (**Figure 2**). The classifiers constructed on the block of IKr and ICaL provided the best discrimination between torsadogenic and non-torsadogenic drugs for the majority of the datasets tested here, including the dataset where drug-induced effects on seven ion channels were reported (**Figure 2**, **Tables 4**, **5**). Our results suggest that among different channels, examination of block of ICaV and IKr might be the most critical for TdP risk prediction. Relative block of ICaV and IKr (among the three currents measured in the in-vitro assay) was shown to provide the best risk prediction, with no role of peak/fast sodium currents in improving the classification in Kramer et al. (2013).

Examination of late sodium block can be important for the drugs with low to moderate ICaV block as these drugs would be predicted TdP+ if only IKr and ICaV block are considered for risk prediction. Moderate to high block of INaL by these drugs can result in suppression of EADs (**Figure 2**) indicating lower TdP risk. Earlier datasets did not report values for INaL block. Dataset 5 reports the value of drug-induced block of seven ion-channels, including the block of INaL. Among the drugs with low to moderate ICaV in Dataset 5, only three drugs [Ranolazine(CM3), Toremifene(CM2) and Quinine(CM3)] have greater than 30% INaL block at critical hERG block concentrations. For the limited data with inconsistent risk categorization, taking into account the INaL block did not improve predictive power of the classifiers (**Tables 4**, **5**). A small improvement in TdP prediction was observed for Datasets 3 and 7 (Okada et al., 2015; Li et al., 2017) when considering drug-induced block of INaL as one of the input features by correctly classifying Ranolazine (the only drug with high late sodium block in absence of ICaV block) in both the datasets (**Table 4** and **Figure 3**). The limited data and inconclusive/minor improvement in torsadogenic risk classification make it difficult to ascertain the role of ion channels such as INaL and IKs in predicting TdP risk.

## 4.2. Predictive Power of Direct vs. Derived Features

In-silico biophysical models can be thought of as a complex nonlinear transfer function, which translates the drug-induced multichannel block effects at channel level (input) to alterations in APs and calcium transients at cellular/tissue levels (output). Several in-silico electrical biophysical models of human ventricular cell models have been published over the last decade (e.g.,Ten Tusscher and Panfilov, 2006; Grandi et al., 2010; O'Hara et al., 2011; Himeno et al., 2015). TdP risk classification on features extracted from the drug-induced responses in isolated cell (Mirams et al., 2011; Christophe, 2013, 2015; Lancaster and Sobie, 2016), tissue (Trenor et al., 2013; Kubo et al., 2017) or organ level (Okada et al., 2015) computational models can provide physiological/mechanistic insights. Moreover, in-silico models serve as an excellent tool for evaluation of drug-safety in diseased conditions (Trenor et al., 2013; Kubo et al., 2017). Our simulations in the OHRmv model under pathological conditions (enhanced late sodium currents) reveal that EAD can appear at significantly lower drug concentrations as only 30% hERG block was required to induce pause-induced EADs under pathological conditions compared to 60% hERG block under normal conditions. Moreover, the modulatory effects of nonhERG channels on EAD induction was also significantly different for simulations under pathological conditions (**Figure 2**). On the other hand, biophysical models show considerable differences in their formulations and can lead to different predictions based on the chosen model. Simulations of EAD generation show a significant difference between the surfaces separating EAD+ region from EAD− region, that are obtained from the OHR and OHRmv models (**Figure 2**). Recently, several efforts have been carried out for optimization of in-silico cardiac cell models for pro-arrhythmia risk assessment (Dutta et al., 2016; Mann et al., 2016; Li et al., 2017).

Statistical/machine learning classifiers that use measured invitro block of multiple cardiac channels (direct features) as their input (Kramer et al., 2013; Mistry et al., 2015) demonstrated comparable accuracy as compared to TdP risk classifiers built on derived features, questioning the need of additional complexity provided by the in-silico models. On the contrary, the study by Okada et al. suggested the need of highly detailed threedimensional cardiac models for pro-arrhythmic risk assessment and showed relatively low predictive ability using the direct features and also certain derived features from the in-silico lumped parameter (zero-dimensional) cellular models (Okada et al., 2015). Derived features from in-silico simulations that incorporate dynamic drug-hERG channel interactions were shown to improve prediction of TdP risk (Li et al., 2017). For all the datasets tested here, including the datasets in Okada et al. (2015) and Li et al. (2017), we showed that the classifiers built on the direct features performed equally or better than the previously developed classifiers on the derived features (**Table 4**). Our results show that for currently available in-vitro assay datasets simple models based on the direct features can provide similar accuarcy to more complex models based on derived features. It should be noted that our two-dimensional risk classifiers on the direct features also utilized insights gained from the computational models (the direct features are examined at critical hERG block concentrations where EAD can arise in the in-silico models). Our parametric simulation for EAD induction highlights one of the possible reasons for the insignificant improvement in predictive power of classifiers built on the derived features from the in-silico models. Although a nonlinear surface is obtained from the in-silico models separating the EAD+ and EAD− regions (**Figure 2**), a hyperplane

$$a \times b \\ \text{lock}\_{I\_{\text{CaV}}} + b \times b \\ \text{lock}\_{I\_{\text{NaL}}} + c \times b \\ \text{lock}\_{I\_{\text{Ks}}} + d = 0 \quad (6)$$

constructed using direct features can result in nearly identical separation, where a, b, c, and d are the parameters of the hyperplane and blockICaV , blockINaL , and blockIKs are the values of block of ICaV INaL, and IKs, respectively. Moreover, with most of the datasets comprising values for block of few channels (Mirams et al., 2011; Kramer et al., 2013) and the much higher incidence of drug-induced block of particular ion channels (INaL, ICaV, and IKr) even when drug-induced modulation of several channels are examined (Crumb et al., 2016), the result is a congregation of majority of the data in a small region of the plausible highdimensional risk space (e.g., see **Figure 2**). For example, the data in Mirams et al. (2011) and Kramer et al. (2013) would fall on a single edge of the 3D EAD space in **Figure 2** in the absence of values for drug-induced block on INaL and IKs in these datasets. This allows risk classification to be performed by a hyperplane with a single parameter, such as blockICaV . Here, we utilized an additional metric, i.e., the hERG ratio, to further improve the classification performance of the direct-feature based classifiers (**Figures 3**, **4**, **6**). For the limited data currently available, risk classification using simple statistical models built on the direct features as the one presented here may suffice.

#### 4.3. Diversity in the Proposed Derived Features

The classifiers built on derived features obtained from the insilico models are based on certain underlying physiological phenomenon (APD, increase in calcium levels, etc.). Hence, derived features are thought to allow better extrapolation to examine drug targets other than those in the training set. However, diverse derived features from the in-silico models have been suggested as possible candidate metrics. Several features from the biophysical models, such as APD50, APD90, calcium level peak, and CaD90 provided the best classification depending on the selected in-silico model (Mirams et al., 2011). Other derived features (EADs, TDR, change in ICaV & INaL) extracted from the AP and calcium transient (Christophe, 2013, 2015; Li et al., 2017) have also been suggested as possible candidate metrics for TdP risk prediction. Rather than examining the individual features separately, a recent study performed a comprehensive feature selection among 331 metrics and determined that two metrics, APD50 and diastolic Ca2+ in the OHR model at 1 Hz pacing, provided the best discrimination between torsadogenic and non-torsadogenic drugs (Lancaster and Sobie, 2016). The overall diversity in reported plausible candidate metrics for TdP risk classification can be attributed to different simulation protocols, drug concentrations and biophysical models. We showed that several derived features obtained from the in-silico models may track together and provide equal predictive power for risk classification when examined independent of drug EFTPC (**Figure 5C**). Equal predictive ability of several features makes it difficult to determine the underlying causal mechanism. In addition, the identical performance of several derived features limits the extensibility of the classifier to untrained targets, as classification results depend on the specific set of features chosen to perform the classification. For example, examination of untrained ion-channel targets using a classifier with diastolic Ca2+ level as the primary risk discriminating feature predicts a decrease in torsadogenic risk for increased Na+ − Ca 2+ currents (Lancaster and Sobie, 2016). On the contrary, TdP risk prediction under Na+ − Ca2+ modulation using a classifier with APD50 or APD90 as the primary discriminating feature would predict opposite effects, with decreased Na+ − Ca2+ exchanger current being associated with decreased TdP risk. Moreover, the derived features obtained from the highly complex biophysical models did not result in improved prediction over the classfiers built on the direct features using the proposed method.

#### 4.4. Limitations

One of the primary limitations is the quality of the datasets itself. The variability in the IC50 values among the several datasets can be one of the reasons for the observation of different thresholds for the hERG ratio and ICaV block that resulted in the best discrimination between TdP+ and TdP− drugs (**Figure 3**). Quantification of the uncertainties in the in-vitro channel screening data and their effects on risk prediction are presented in Johnstone et al. (2016). Inconsistencies in risk definition presents another important challenge for torsadogenic risk assessment. Wi´sniowska and Polak (2017) reports a comprehensive list of compounds that have been inconsistently defined as TdP+ or TdP− in different studies to develop torsadogenic risk classifiers. The different categorizations can lead to different interpretations and accuracy scores for TdP risk determination (**Table 5**). Standardization of torsadogenicity definition, which would allow comparison of the performance of different classifiers/features, is required. Certain steps in this direction have been started. Based on a general consensus, a working group formed under CiPA initiative picked 28 compounds and categorized each into three groups (Colatsky et al., 2016; Fermini et al., 2016)

#### REFERENCES

for testing/training of the classifiers. In-silico simulations of dynamic drug-channel interactions might be essential to further improve the TdP risk assessment (Li et al., 2017). Inclusion of the drug-binding parameter, in addition to the amount of block of ion-channels, resulted in 100% prediction using our approach. Sufficient IKr block was assumed to be necessary for TdP generation in our method. The effects of non-hERG channels are thought to enhance or mitigate the torsadogenic effects of IKr block. The method resulted in excellent predictive performance across several datasets that report drug-induced block of various ion channels only. However, drug-induced enhancement of ion-channel currents such as INaL can result in increased TdP risk in the absence of hERG block (Lacerda et al., 2008; Yang et al., 2014). The method could be further extended to examine such effects when more data are available. The present work not only provides a new method for invitro ion-channel screening based TdP risk classification but also highlights several important issues in regards to the use of drug-induced multi-channel blockage for torsadogenic risk prediction.

#### AUTHOR CONTRIBUTIONS

JP, VG, and JR wrote the manuscript. JP, VG, and JR designed the research. JP and VG performed the simulations. JP, VG, and JR analyzed the data.

#### SUPPLEMENTARY MATERIAL

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

Pharmacol. Therapeut. 129, 109–119. doi: 10.1016/j.pharmthera.2010. 08.008

combining a patch clamp and heart simulator. Sci. Adv. 1:e1400142. doi: 10.1126/sciadv.1400142

**Conflict of Interest Statement:** All authors are employees of IBM Research.

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.

The reviewer EP and handling Editor declared their shared affiliation.

Copyright © 2017 Parikh, Gurev and Rice. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Synergistic Anti-arrhythmic Effects in Human Atria with Combined Use of Sodium Blockers and Acacetin

Haibo Ni 1, 2, 3, Dominic G. Whittaker 1 , Wei Wang1 , Wayne R. Giles 4 , Sanjiv M. Narayan5 and Henggui Zhang1, 2, 3, 6 \*

*1 Biological Physics Group, University of Manchester, Manchester, United Kingdom, 2 Space Institute of Southern China, Shenzhen, China, 3 Key Laboratory of Medical Electrophysiology, Ministry of Education, Collaborative Innovation Center for Prevention and Treatment of Cardiovascular Disease/Institute of Cardiovascular Research, Southwest Medical University, Luzhou, China, 4 Faculties of Kinesiology and Medicine, University of Calgary, Calgary, AB, Canada, 5 Department of Medicine, Stanford University School of Medicine, Stanford, CA, United States, 6 School of Computer Science and Technology, Harbin Institute of Technology, Harbin, China*

#### Edited by:

*Eleonora Grandi, University of California, Davis, United States*

#### Reviewed by:

*Christopher Huang, University of Cambridge, United Kingdom Andrew G. Edwards, Simula Research Laboratory, Norway*

\*Correspondence: *Henggui Zhang<EMAIL_ADDRESS>
#### Specialty section:

*This article was submitted to Integrative Physiology, a section of the journal Frontiers in Physiology*

Received: *31 July 2017* Accepted: *08 November 2017* Published: *23 November 2017*

#### Citation:

*Ni H, Whittaker DG, Wang W, Giles WR, Narayan SM and Zhang H (2017) Synergistic Anti-arrhythmic Effects in Human Atria with Combined Use of Sodium Blockers and Acacetin. Front. Physiol. 8:946. doi: 10.3389/fphys.2017.00946* Atrial fibrillation (AF) is the most common cardiac arrhythmia. Developing effective and safe anti-AF drugs remains an unmet challenge. Simultaneous block of both atrial-specific ultra-rapid delayed rectifier potassium (K+) current (IKur) and the Na+ current (INa) has been hypothesized to be anti-AF, without inducing significant QT prolongation and ventricular side effects. However, the antiarrhythmic advantage of simultaneously blocking these two channels vs. individual block in the setting of AF-induced electrical remodeling remains to be documented. Furthermore, many IKur blockers such as acacetin and AVE0118, partially inhibit other K+ currents in the atria. Whether this multi-K+-block produces greater anti-AF effects compared with selective IKur-block has not been fully understood. The aim of this study was to use computer models to (i) assess the impact of multi-K+-block as exhibited by many IKur blokers, and (ii) evaluate the antiarrhythmic effect of blocking IKur and INa, either alone or in combination, on atrial and ventricular electrical excitation and recovery in the setting of AF-induced electrical-remodeling. Contemporary mathematical models of human atrial and ventricular cells were modified to incorporate dose-dependent actions of acacetin (a multichannel blocker primarily inhibiting IKur while less potently blocking Ito, IKr, and IKs). Rate- and atrial-selective inhibition of INa was also incorporated into the models. These single myocyte models were then incorporated into multicellular two-dimensional (2D) and three-dimensional (3D) anatomical models of the human atria. As expected, application of IKur blocker produced pronounced action potential duration (APD) prolongation in atrial myocytes. Furthermore, combined multiple K+-channel block that mimicked the effects of acacetin exhibited synergistic APD prolongations. Synergistically anti-AF effects following inhibition of INa and combined IKur/K+-channels were also observed. The attainable maximal AF-selectivity of INa inhibition was greatly augmented by blocking IKur or multiple K+-currents in the atrial myocytes. This enhanced anti-arrhythmic effects of combined block of Na+- and K+-channels were also seen in 2D and 3D simulations; specially, there was an enhanced efficacy in terminating re-entrant excitation waves, exerting improved antiarrhythmic effects in the human atria as compared to a single-channel block. However, in the human ventricular myocytes and tissue, cellular repolarization and computed QT intervals were modestly affected in the presence of actions of acacetin and INa blockers (either alone or in combination). In conclusion, this study demonstrates synergistic antiarrhythmic benefits of combined block of IKur and INa, as well as those of INa and combined multi K+-current block of acacetin, without significant alterations of ventricular repolarization and QT intervals. This approach may be a valuable strategy for the treatment of AF.

Keywords: atrial-selective block, atrial fibrillation, sodium and potassium current block, multiscale simulation, synergistic antiarrhythmic effect

#### INTRODUCTION

Despite recent advances in the management of Atrial fibrillation (AF), the world's most common cardiac arrhythmia (Dobrev et al., 2012; Nattel and Dobrev, 2017), developing effective and safe antiarrhythmic drugs for treatment of AF remains challenging (Aguilar-Shardonofsky et al., 2012; Aguilar et al., 2015). Frequently these antiarrhythmic agents promote ventricular arrhythmias (Dobrev et al., 2012; Woods and Olgin, 2014; Voigt and Dobrev, 2016) by prolonging cellular action potential durations (APDs). The associated QT-interval prolongation can lead to life-threatening consequences. Developing atrial-selective drugs is acknowledged to be a current strategy for the treatment of AF (Burashnikov et al., 2007).

Atrial and ventricular tissues show intrinsic regional differences in their cellular ion channel properties, thus suggesting a basis for developing atrial-selective drugs. For example, the atrial and ventricular fast sodium (Na+) channel currents (INa) exhibit different voltage-dependent inactivation properties, opening the opportunity for atrial-selective Na+ channel blockade (Burashnikov et al., 2007; Antzelevitch and Burashnikov, 2009; Zygmunt et al., 2011). Previous simulation studies have demonstrated that by optimizing state-dependent Na+-channel blocking dynamics (i.e., drug-channel interaction parameters), atrial-selective block of INa could be achieved and that could maximize pharmaceutical effects on the atria while minimizing their proarrhythmic actions in the ventricles (Aguilar-Shardonofsky et al., 2012; Aguilar et al., 2015).

Another tissue-specific difference between the atria and ventricles is that the ultra-rapid delayed rectifier potassium current (IKur, carried by the KV1.5 channel) contributes to repolarization in the atria but plays little role in the ventricles (Tamargo et al., 2009; Ravens and Wettwer, 2011). Recent studies suggest that atrial-selective blockade of IKur may be an effective pharmacological treatment of AF (Li et al., 2008; Pavri et al., 2012; Loose et al., 2014; Ford et al., 2016). Although the efficacy of IKur block in the treatment of AF remains controversial (Burashnikov and Antzelevitch, 2008), multiple IKur blockers have been developed (Tamargo et al., 2009; Loose et al., 2014; Wettwer and Terlau, 2014; Ford et al., 2016). Interestingly, these IKur blockers actually target multiple channels, and are known to inhibit other K+ currents including Ito and IK,ACh in the atria (Burashnikov and Antzelevitch, 2008). Examples of such blockers include AVE0118 (Gögelein et al., 2004), AVE1231 (Wirth et al., 2007), AZD7009 (Persson et al., 2005), and acacetin (Li et al., 2008). Among these channel blockers, acacetin, a natural flavone initially isolated from a traditional Chinese medicine Xuelianhua, potently blocks IKur, Ito, and IK,ACh, and has a smaller potency in inhibiting IKr and IKs (Li et al., 2008), similar to AVE0118 (Gögelein et al., 2004; Haan et al., 2006). Acacetin is regarded as a promising atrial-selective agent for the treatment of AF (Li et al., 2008). However, the actions of acacetin on atrial electrophysiology, especially its effects following AF-induced electrical remodeling of atrial electrophysiological properties (Dobrev et al., 2012), remain to be elucidated. Furthermore, since most IKur blockers inhibit other K+ channels, the question whether the "additional" inhibitive actions produce favorable antiarrhythmic effects has not been addressed thoroughly. A better understanding of these effects of modulating multiple ion channels on atrial excitation and recovery/repolarization may provide insights into evaluating and developing antiarrhythmic drugs.

Interestingly, simultaneous multiple-channel blocking of both depolarization and repolarization currents is attracting more attention since empirical observations suggest that such multichannel blockers generally mediate more effective antiarrhythmic effects (Kirchhoff et al., 2015; Reiffel et al., 2015; Hartmann et al., 2016). A recent numerical and experimental study on the canine heart (Aguilar et al., 2015) suggested that blocking K+ currents enhanced the anti-arrhythmic effects and AF-selectivity of INa blockade. In their study, IKur block was modeled using a simple pore block scheme by reducing the conductance of the channel. As the kinetics of drug action plays an important role in the effects of IKur blockers (Scholz et al., 2013; Ellinwood et al., 2017), in simulating IKur block a state-dependent block model reproducing a realistic blocker is more favorable. Once again, the effects of combined INa and IKur block on the human atria, especially in the setting of AF-induced electrical remodeling which reduced IKur, remain to be elucidated. It is also unclear how multiple-channel blockade may affect QT interval.

In the present study, it was hypothesized that combined block of INa and K+-currents (predominantly IKur) could produce antiarrhythmic benefits compared with the application of either blocker alone in the setting of AF-induced electrical remodeling. We have tested the hypothesis with the following three aims: (i) to identify and illustrate the effects of the realistic IKur

blocker, acacetin, on atrial electrophysiology following AFrelated remodeling; (ii) to assess whether combined INa and IKur block produce synergistic antiarrhythmic effects; and (iii) to investigate the action of such drug combinations on ventricular electrophysiology.

#### METHODS

### Modeling Electrophysiology of the Human Heart

To simulate human atrial electrophysiology, an updated Colman et al. model for atrial electrophysiology (Colman et al., 2013, 2017) was used. For in silico study of effects of chronic AF- (cAF) induced electrophysiological remodeling on the atria, we incorporated the cAF model parameters from our previous study (Colman et al., 2013) into the updated atrial single cell model (for details please see Online Supplement Material 1.1).

To assess the effects of the anti-AF drugs on the human ventricles, simulations were performed to investigate the actions of the anti-AF drugs on the ventricular AP, INa and QT intervals in the electrograms. In these simulations, the mathematical model developed by O'Hara et al. (2011) was used to represent the ventricular electrophysiology. Additionally, the INa formulation in the model was replaced by the one in the Luo–Rudy model (Luo and Rudy, 1994), which enabled electrical excitation to propagate in the tissue model.

More detailed descriptions of the electrophysiological models of human atrial and ventricular cells are given in Online Supplementary Material 1.1.

# Modeling State-Dependent INa Block

As in previous studies (Aguilar-Shardonofsky et al., 2012; Aguilar et al., 2015), INa block was simulated using a guarded receptor model with dynamical drug-channel interactions. This approach allows for investigations of the role of the specified parameters for selected INa blockers, and effects of combined IKur block on the atrial selectivity of Na+-channel block. The guarded receptor model considers the binding and unbinding kinetics of the drug to INa channels in a drug concentration-dependent manner. They can be described by first-order transition equations (Aguilar-Shardonofsky et al., 2012; Aguilar et al., 2015). It was also assumed that the drug predominantly binds to the activated and/or inactivated states of INa. The blockade of INa is given by Aguilar-Shardonofsky et al. (2012) and Aguilar et al. (2015):

$$I\_{\rm Na} = \mathcal{g}\_{\rm Na} \left( 1 - B\_A - B\_I \right) m^3 hj \left( V\_m - E\_{\rm Na} \right) \tag{1}$$

$$\frac{dB\_A}{dt} = K\_A \left[ D\_{Na^+} \right] m^3 h j \left( 1 - B\_A - B\_I \right) - L\_A B\_A \tag{2}$$

$$\frac{dB\_I}{dt} = K\_I \left[ D\_{Na^+} \right] \left( 1 - h \right) \left( 1 - B\_A - B\_I \right) - L\_I B\_I \tag{3}$$

where gNa is the maximum conductance of INa; BA and BI are the fractional blockade of activation and inactivation channels; m is the activation gate state variable, h and j are the inactivation gate state variables; Vm the transmembrane potential; ENa the reversal potential of Na+; KA, KI the binding constants and LA, LI the unbinding constants; [DNa+ ] is the concentration of a Na+-blocker. As in previous studies (Aguilar-Shardonofsky et al., 2012; Aguilar et al., 2015), a concentration of 60µM was utilized unless otherwise stated; this concentration was chosen based on previous experimental and modeling studies (Zhu et al., 2006; Moreno et al., 2011; Aguilar-Shardonofsky et al., 2012; Aguilar et al., 2015); a parameter set (KA = 100 ms−1 · M−1 , KI = 100 ms−1 · M−1 , LA = 1 ms−1 , LI = 0.01 ms−1 ) was first used to represent the kinetics of an INa-selective blocker.

In our investigations of the AF-selectivity of INa block following AF-remodeling, the binding and unbinding constants of the INa blockers were varied to evaluate the dependence of INa block on these parameters, and whether an atrial-selective anti-AF action could be achieved in cAF-remodeled myoctes. The AFselectivity of Na+-channel blockade was defined as the product of atrial-selectivity, rate-selectivity and block efficacy. With fractional block (Bf ) by Na+-channel blockers being measured as the relative reduction in the peak of INa, the rate-selectivity was defined as the ratio of Bf measured in an atrial myocyte paced at 6 Hz to that paced at 1 Hz (Aguilar-Shardonofsky et al., 2012; Aguilar et al., 2015). Atrial-selectivity was used to determine the extent of atrial-ventricular difference in response to each drug. This was represented by the ratio of Bf observed from an atrial myocyte to that of a ventricular cell both paced at 1 Hz. In this study, we defined block efficacy (E) as:

$$E = \frac{1.0}{1.0 + \left(\frac{0.5}{B\_{\text{f},6\text{Hz}}}\right)^4} \tag{4}$$

where Bf,6Hz is the fractional block of INa measured in an atrial cell paced at 6 Hz. Different from Aguilar et al. (2015), we introduced block efficacy to constrain the measure of AFselectivity when the fractional block observed in a ventricular cell paced was minimal (and could result in a great atrialselectivity), which otherwise could give a great value in AFselectivity regardless of a small Bf,6Hz.

To assess the dependence of the AF-selectivity of INa block on the drug action kinetics, the unbinding constants LA and LI were first varied over a parameter space from 10−5 to 100 ms−1 , while KA and KI were fixed (see Figure S6 of online Supplementary Materials for more details). The resultant unbinding constants were used in subsequent optimizations varying KA and KI . The parameter space was {1, 10, 100, 500, 2,500, 10,000} for KA and {1, 10, 100, 200, 500, 2,500} for KI . The parameter space fell into a likely range of INa blockers as summarized in Aguilar-Shardonofsky et al. (2012).

#### Modeling Effects of Acacetin on Atrial and Ventricular Electrophysiology

Acacetin was the chosen IKur blocker in the present study. To reveal the functional effects of (i) pure IKur block vs. (ii) the effects of combined K+ currents block by acacetin on human atrial electrophysiology, the actions of acacetin were modeled by considering its effects on (a) IKur only, and (b) all the respective K+ currents as detailed in **Table 1**. This approach allows for modeling the effects of the selective IKur block as well as uncovering the role of "additional" inhibitory effects of acacetin on other K+ currents.

TABLE 1 | Concentration-dependent block of K+-currents by acacetin (Li et al., 2008; Wu et al., 2011).

#### Modeling Effect of Acacetin on IKur

Previous modeling studies have demonstrated the important role of the kinetic properties of drug actions in IKur block (Tsujimae et al., 2008; Almquist et al., 2010; Scholz et al., 2013). In addition, the pharmaceutical effects of acacetin on IKur are characterized by use- and rate-dependencies (Wu et al., 2011), which have also been observed in other IKur blockers (Pavri et al., 2012; Ford et al., 2016). Therefore, it was necessary to adopt a statedependent block model (Brennan et al., 2009) for simulating the blockade of IKur by acacetin. Similar to our approach for modeling INa block, the binding and unbinding kinetics of a drug was described by a first-order transition equation, in contrast to simulating IKur block by reducing its conductance in Aguilar et al. (2015). Experimental studies revealed that acacetin binds to both the open and closed gates of KV1.5 (Wu et al., 2011). Therefore, following the guarded receptor formulas given in Equations (1–3), the formulation of inactivation-state binding and unbinding kinetics in INa block was modified to simulate the closed-state block of IKur by acacetin. The guarded receptor model of IKur block by acacetin is given by:

$$I\_{\rm Kur} = \mathcal{g}\_{\rm Kur} \left(1 - B\_O - B\_C\right) ai \left(V\_m - E\_K\right) \tag{5}$$

$$\frac{dB\_O}{dt} = K\_O \exp\left(Z\_{KO}\frac{V\_m F}{RT}\right) \left[D\_{K^+}\right] ai \left(1 - B\_O - B\_c\right)$$

$$-L\_O \exp\left(-Z\_{LO}\frac{V\_m F}{RT}\right) B\_O \tag{6}$$

$$-L\_O \exp\left(-Z\_{LO}\frac{V\_m F}{RT}\right) B\_O \tag{6}$$

$$\frac{dB\_c}{dt} = K\_c \exp\left(Z\_{Kc}\frac{V\_m F}{RT}\right) [D\_{K^+}] \left(1 - a\right) i \left(1 - B\_O - B\_c\right)$$

$$-L\_c \exp\left(-Z\_{Lc}\frac{V\_m F}{RT}\right) B\_c \tag{7}$$

where gKur isthe conductance of IKur; Bo and BC are the fractional block on open and closed state variables, respectively; a and i are the activation and inactivation gate variables; EK is the reversal potential of potassium; F, R, and T are the Faraday's constant, universal gas constant and temperature respectively. KO and Kc are the binding constants; LO and Lc are the unbinding constants; ZKO, ZLO, ZKc and ZLc are the drug charge parameters for the corresponding binding or unbinding processes; [DK+ ] is the concentration of acacetin applied. The binding and unbinding parameters were obtained by fitting the model to the experimental data on the rate-dependent blockade of IKur by acacetin (Wu et al., 2011), as detailed in Online Supplementary Material 1.1.

**Figure 1** shows a simulated frequency-dependent block of IKur by acacetin, and this is compared to the experimental data (**Figures 1A,B**). As shown, repeating the voltage command (**Figure 1B**, insert) at 0.5 Hz resulted in an approximately 50% blockade in this current after application of 3µM acacetin. Increasing the voltage command rate to 4 Hz significantly increased the relative fractional block to approximately 63% (**Figure 1C**).

#### Modeling Effect of Acacetin on Ito, IKr, and IKs

In addition to inhibiting IKur in the atria, acacetin potently blocks both Ito and IK,ACh, and also modulates IKr and IKs, exhibiting multiple K+-current block. The parameters of Hill equations describing use-dependent inhibitions of these channels by acacetin are shown in **Table 1**. In the simulations, the effects of acacetin on these channels were modeled using a simple pore block model (Yuan et al., 2015). In the present study, we did not simulate the effects of acacetin on IK,ACh inhibitions as the role of autonomic regulation on AF is beyond the scope of the study.

#### Simulations of the Effects of Acacetin on Human Ventricle

The effects of acacetin on human ventricular APs are unknown, although experimental data demonstrated that acacetin at 30µM did not affect the heart rate and QT interval in isolated rabbit hearts (Li et al., 2008). In the present study, it was assumed that similar effects on the K+ currents (Ito, IKr, IKs) in atrial myocytes could be extrapolated to the ventricular myocytes. We acknowledge that IKur is negligible in ventricles (Ravens and Wettwer, 2011), therefore in simulations of blocking IKur alone, the ventricular electrophysiology was not affected.

#### Tissue Models

The effects of acacetin and INa blockers on atrial and ventricular electrophysiology were further evaluated using tissue models. The monodomain equation (Clayton et al., 2011) was employed to simulate the excitation wave propagation in the myocardium. 1D models of human atrial strands were used to quantify the effects of channel blockers on atrial conduction velocity and APD restitution properties. Changes in ventricular depolarization and repolarization in response to these drugs were evaluated using a 1D model representing a transmural strand of ventricular tissue. In order to evaluate the antiarrhythmic effects of the channel blockers on re-entrant excitations in atria in the setting of cAF-induced remodeling, both idealized 2D models representing an isotropic slab of atrial tissue and an anatomically accurate 3D model of the human atria (Aslanidi et al., 2011; Colman et al., 2013, 2017; Whittaker et al., 2017) were employed to simulate the behavior of re-entrant excitations in atrial tissue. Pseudo-ECGs (pECGs) (Gima and Rudy, 2002; Baher et al., 2006) were computed as a measure of the excitation rates of intissue with sprial excitation waves. Detailed descriptions of these tissue models and pECGs are given in Online Supplementary Material 1.2, 1.3.

#### RESULTS

The updated Colman et al. human atrial model was first used to simulate effects cAF-induced remodeling on the

frequencies plotted against the pulse number of the voltage step. The simulated data (lines) were compared with experimental values (squares). The relative fraction was obtained by normalizing the end-step current measured from each pulse following application of acacetin to that of control. Experimental data were digitalized from Wu et al. (2011).

action potential (AP) and calcium transient (CaT). Details are presented in Online Supplementary Material 2.1. The resultant changes in APD, APD restitution and CaT following cAF-induced remodeling as compared to those under the normal condition showed good agreement with previous experimental (Bosch et al., 1999; Osaka et al., 2000; Workman et al., 2001; Dobrev and Ravens, 2003; Voigt et al., 2012) and simulation studies (Zhang et al., 2005; Grandi et al., 2011; Colman et al., 2013, 2017; Wilhelms et al., 2013).

# Effects of Application of Acacetin on Human Atrial Cells

To reveal the roles of inhibition of individual channels by acacetin in modulating cellular AP by acacetin, both the individual and combined block of Ito, IKr, IKur, and IKs by acacetin (3.2µM) in simulated SR at cycle length 1,000 ms without (normal) or with cAF-related electrical remodeling were simulated. **Figures 2A,B** illustrates the effects of individual and combined K+ channel block by acacetin on AP waveform. The alterations to APD relative to the control are summarized in **Figures 2C,D**.

In the absence of electrical remodeling, in normal myocytes at a cycle length of 1,000 ms, simulated IKs or IKr block by acacetin (3.2µM) presented no significant alterations to the atrial AP: although the atrial repolarization was delayed by 1.3 and 5.9 ms, respectively, the plateau phase was not affected, which is consistent with the minimal potency of acacetin on these channels (**Table 1**). Similar effects were also obtained from our simulated IKr and IKs block by the compound in the cAFremodeled atrial cells.

Selective block of Ito (alone) by acacetin elevated the atrial plateau potential in both normal and cAF-remodeled myocytes paced at 1 Hz, and this led to modest prolongations in APD30 (by 6.4 and 3.6 ms for normal and AF-remodeling myocytes, respectively). The changes in APD90 due to Ito block varied between the two conditions: under normal conditions the atrial APD90 was shortened by 2.4 ms, whereas it was prolonged by 3.7 ms following cAF-remodeling at a stimulus rate of 1 Hz.

In contrast, blocking IKur alone by acacetin resulted in a pronounced alteration to the shape and duration of the AP in both normal and cAF-remodeled myocytes. The inhibition in IKur significantly elevated the plateau potential of atrial AP (by 7.1 and 5.7 mV in normal and cAF-remodeled myocytes, respectively), and this was accompanied by marked prolongations in APD30 (by 105.9 and 23.6 ms in normal and cAF-remodeled myocytes, respectively). The prolongation in APD90 induced by the IKur block was 9.8 ms for normal atrial cells, and was more pronounced (23.6 ms) in cAFremodeled myocytes, despite that IKur was down-regulated by cAF-remodeling.

We note that combined effects of acacetin (3.2µM) on multiple K+-currents produced greater alterations to the AP than those of any individual blocking effect. We have quantified effects produced by the combined block and compared it with the sum of the changes seen in each individual block. Synergistic effects were observed in the changes in APD90, represented by a further prolongation of 9.3 and 1.1 ms in APD90 in normal and cAFremodeled myocytes, respectively. Additionally, the effects of IKur block dominated the AP-modulation by the compound, which is consistent with the high potency of acacetin on the channel (**Table 1**).

# Effects of Sodium Blocker and Acacetin on cAF-Remodeled Atrial Myocytes and Ventricular Cells

#### Effects on Single Myocyte AP and INa

Individual and combined effects of Na+-block (indicated by Bl·INa) and K+-block by acacetin (3.2µM) on human atrial electrophysiology after cAF-remodeling were simulated to assess any anti-AF benefits. Effects of acacetin (representing K+-block) were simulated in different settings: (i) IKur block alone (denoted by Bl·IKur) and (ii) combined block of all K+-currents in **Table 1** (denoted by Comb·Bl·IX). In addition, the effects of Na+- and K +- block on human ventricular myocytes were also studied to assess the atrial-selectivity of the block. The results are shown in **Figures 3A–C** and quantitative measurements are shown in **Figures 4A–C**.

For atrial myocytes paced at 1 Hz, Bl·IKur, and Comb·Bl·IX prolonged atrial APD (**Figure 4Ai** and as presented in Effects of Application of Acacetin on Human Atrial Cells) whilst their effects on peak INa were minimal (reducing peak INa by less than 0.3%). Bl·INa alone slightly reduced the peak INa by 1.63% without affecting the APD. The fractional inhibition in INa by INa-block was slightly increased by the addition of Bl·IKur or Comb·Bl·IX (**Figure 4Aii**). In the ventricles, the simulated application of acacetin induced a prolongation of 19.5 ms in APD90 compared with that in control (drug-free) condition (**Figures 3B**, **4Ci**). Bl·INa alone showed a negligible inhibitory effect on the ventricular INa (by 0.42%), which was also not affected by combining Bl·INa and Comb·Bl·IX (**Figure 4Cii**).

In atrial myocytes paced at 6 Hz, AP alternans were observed under the drug-free condition (**Figure 3Ci**): the APD varied between 100.1 and 88.1 ms. In the presence of AP alternans, the changes in APD by the Na+- and K+- blockers were quantified by comparing the corresponding big APs at baseline and after drug actions. These values were selected based on the characteristics of AP (APD prolongations seen in the long AP, and reduced INa for the short AP) that may be anti-arrhythmic. The fractional reductions in peak INa were calculated from the INa associated with the shorter APs. The results showed that applying Bl·IKur or Comb·Bl·IX alone both abolished the AP alternans while prolonging the APD to 116.3 and 126.3 ms and reducing the peak INa by 5.9 and 20.1%, respectively. The application of Bl·INa alone produced a minor APD prolongation (3.5 ms) and a reduction of 16.2% in peak INa. Combining block of INa with Bl·IKur or Comb·Bl·IX promoted the genesis of AP alternans, resulting in substantial prolongations in the APD of the big APs (by 35.4 ms for Bl·IKur + Bl·INa, and 55.6 ms for Comb·Bl·IX + Bl·INa) and dramatic decreases in the peak INa (by 57.5% for Bl·IKur + Bl·INa and 88.2% for Comb·Bl·IX + Bl·INa) in the corresponding small APs. These results suggest that the combined block of Bl·INa and Comb·Bl·IX/Bl·IKur exhibited synergistic antiarrhythmic effects manifested by prolongation in APD and reduction in peak INa. However, an increased susceptibility to AP alternans was observed at a fast pacing rate of 6 Hz, which may be potentially proarrhythmic at fast heart rates.

#### Effects on Steady-State Restitutions of APD and Conduction Velocity

Steady-state APD restitutions of cAF-remodeled human atria were simulated at both the cellular and tissue levels. In single myocyte simulations (Figure S5 in Online Supplementary Material 2.2), APD was prolonged over the entire range of simulated basic cycle lengths (BCL) for Bl·IKur and Comb·Bl·IX as compared to the control (drug-free) conditions. The reduction in peak INa in Bl·INa was rate-dependent and significantly greater at fast pacing rates. K+-block alone (Bl·IKur or Comb·Bl·IX) slightly shifted the rate-dependence of peak INa to larger BCLs. In comparison to the effects of individual current block scenarios, synergistic reductions in peak INa were observed following combined blocks of Bl·INa with Bl·IKur or Comb·Bl·IX over a wide range of BCLs. As compared to the drug-free conditions, AP alternans were observed at greater BCLs after K+-block, and this was further increased by combined Na+- and K+-block (Figure S5 in Online Supplementary Material 2.2).

.

FIGURE 3 | Simulated AP and INa traces of cAF-remodeled atrial myocytes and ventricular cells in response to Na+- and K+- block regimens. (Ai) APs from

#### Effects of Combined Na+- and K +- Block on the AF-Selectivity of INa block

AF-selectivity of Na+ blockers in cAF-remodeled hearts was examined by varying the drug binding and unbinding constants over wide parameters spaces to provide information concerning drug-Na+-channel interactions for various drug candidates. This was done by independently changing LA and LI for fixed {KA, KI} (Online Supplementary Material 2.3); and then varying KA and KI for fixed {LA, LI}. In this way we obtained the maximum AF-selectivity over the parameter space of drug binding and unbinding kinetics. Simulations with varied KA and KI were repeated for Bl·INa + Bl·IKur and Bl·INa + Comb·Bl·IX.

**Figures 6A–D** illustrates the block efficacy (defined in Equation 4), rate-selectivity, atrial-selectivity and the resultant AF-selectivity for Bl·INa alone and the combined block as a function of KA and KI . For Bl·INa alone, the block efficacy increased with increase of KI , whereas the rate-selectivity was reduced by increasing KA or KI . The AF-selectivity reached a maximum value of 9 at KA = 1 and KI = 500 ms−1 · M−1 . The combined blocks achieved significantly greater AF-selectivity than Bl·INa alone: the maximum attainable AF-selectivity was increased by nearly 6-fold for Bl·IKur + Bl·INa and more than 14-fold for Comb·Bl·IX + Bl·INa as compared to Bl·INa alone (**Figure 6C**). These dramatic increases were attributed to the significantly greater values in all metrics contributing to the AF-selectivity. The maximal block efficacy achieved by Bl·INa alone was 0.77, and was increased to 0.81 and 0.94 for Bl·IKur + Bl·INa and Comb·Bl·IX + Bl·INa, respectively. A more appreciable increase in the maximal rate-selectivity was observed by the combined blocks as compared to Bl·INa alone (8-fold for Bl·IKur + Bl·INa and nearly 10-fold for Comb·Bl·IX + Bl·INa). Additionally, the atrial-selectivity was also increased by the combined block, although to a lesser extent. Bl·IKur + Bl·INa exhibited a greater atrial-selectivity than that of Comb·Bl·IX + Bl·INa since Bl·IKur was assumed to have no effect on the ventricles.

Furthermore, these simulations revealed that the block efficacy, rate-selectivity and atrial-selectivity were strongly

dependent on the inactivation state binding rate KI . These measures were also dependent on the open-state binding kinetics KA, but to a much lesser extent. The block efficacy was mainly determined by KI : an increase in KI led to a significant increase in the block efficacy. In combined block, the maximal rate- and atrial-selectivity were observed for KA = KI = 1 ms−1 · M−1 and increases in KI resulted in substantial reductions in the rate- and atrial-selectivity. In Bl·INa alone, the parameter set KA = 1 ms−1 · M−1 , KI = 200 ms−1 · M−1 produced a maximal value in atrial-selectivity. Collectively, the optimal KI that maximized AF-selectivity was 200 ms−1 · M−1 for Bl·INa and Bl·IKur + Bl·INa, and smaller (100 ms−1 · M−1 ) for Comb·Bl·IX + Bl·INa. The optimal KA = 1 ms−1 · M−1 was seen for all conditions. Increasing KA consistently resulted in a smaller rate- and atrial-selectivity and therefore reduced AF-selectivity. These results suggest that the inactivation-state binding rate might be a more favorable targeting parameter than the open-state binding kinetics in optimizing AF-selectivity of Na+-blockers.

#### Effects of INa and IKur Block on Spiral Excitation Events in cAF-Remodeled Atria Two-Dimensional Simulations

Using the cross-shock protocol, spiral waves were initiated in a 2D model representing a tissue slab of cAF-remodeled human atria. For each condition, a 10-s episode of electrical activity was simulated. Representative snapshots of the re-entrant waves in control (drug-free) and following application of drugs are presented in **Figure 7A**. The trajectories of the tips of re-entrant rotors under these conditions were traced and are shown in **Figure 7B**. The number of rotors during the time course of wave evolution was also measured (**Figure 7C**). The simulated pseudo-ECGs, membrane potential traces extracted from a representative

myocyte and the corresponding fractional block of INa and IKur are detailed in Online Supplementary Material 2.4.

Under the control (drug-free) condition, a single rotor was formed at approximately t = 630 ms; this broke into two spiral waves at t = 830 ms. These two rotors were stably anchored with star-shaped tip trajectories at the bottom half of the slab and persisted throughout the rest of the simulated 10-s episode (**Figures 7Ai,Bi,C**). In simulating drug actions, each drug was applied at t = 2,500 ms. For Bl·INa the dual rotors progressively became unstable, and the tips of the spiral waves meandered out of the tissue at approximately t = 6,000 ms, leading to selftermination of the re-entrant waves (**Figures 7Aii,Bii,C**). The dual rotors persisted throughout the period of the simulation after applying Bl·IKur alone (**Figures 7Aiii,Biii,C**). For Bl·INa + Bl·IKur, the rotor at the bottom left corner of the slab became unstable and meandered out of the tissue at approximately t = 3,800 ms, whereas the trajectory of the second rotor was confined to a small tissue area until t = 7,000 ms and then gradually became chaotic, forming up to 3 transient rotors that self-terminated at t = 8,334 ms (**Figures 7Aiv,Biv,C**). A

tissue slab in drug-free condition or after applying the drugs. (A) Snapshots of simulated re-entrant excitation events; The time sequence (ms) is indicated at the top left corner of each screenshot. (B) Tip trajectories of re-entrant waves. (C) Temporal evolution of total number of rotors represented by the number of spiral wave tips in tissue. In both (A,B), (i) Drug-free (CTL) condition, (ii) Bl·INa alone, (iii) Bl·IKur alone, (iv) Combined Bl·INa and Bl·IKur, (v) Applying Comb·Bl·IX alone, and (vi) Combined Bl·INa and Comb·Bl·IX. Rate constants for INa block: *KA* = 1 ms−1 · M−1 , *KI* = 100 ms−1 · M−1 , *LA* = 1 ms−1 , *LI* = 0.01 ms−1 .

similar but more marked effect was seen in the simulations that addressed the aggregate effects of acacetin: the bottom left rotor quickly meandered out of the tissue at t = 3,510 ms whilst the tip trajectory of the other rotor became chaotic and terminated at t = 5,439 ms (**Figures 7Av,Bv,C**). We note that the combined Bl·INa and Comb·Bl·IX exerted the strongest potency in terminating re-entrant excitations in these simulations: the two rotors transiently turned unstable and chaotic and selfterminated at t = 3,555 ms, with up to 5 rotors during the excitation in the slab (**Figures 7Avi,Bvi,C**).

The anti-arrhythmic benefits of combined Na+- and K+ block were clearly revealed by additional simulations assuming the use of a different Na+-blocker (KA = 1 ms−1 · M−1 , KI = 200 ms−1 · M−1 , LA = 1 ms−1 , LI = 0.01 ms−1 ) and at a reduced dose (75%) of both Na+-blocker and acacetin. The life span of reentrant excitations was measured and shown in **Figure 8A**. Also, the pECG was computed and the segment from t = 3,000 ms to 500 ms before the termination of re-entries (or the end of the simulation if the rotor sustained) was analyzed using the Fast Fourier Transform to obtain the dominant frequency (DF) of the re-entrant excitations, which is illustrated in **Figure 8B**.

At the simulated doses of acacetin, applying Bl·IKur alone did not lead to termination of re-entrant waves within the duration of the simulation (7,500 ms after TDrug), whereas the rotors were terminated in the simulations for Comb·Bl·IX at both doses (**Figure 8Ai–iv**), thus demonstrating enhanced anti-AF benefits of combined K+-channel block. For the simulated Bl·INa alone, the Na+-blocker with KI = 100 ms−1 · M−1 led to termination of AF at the control dose (lifespan of 4,102 ms) but not at the reduced dose; increasing KI of the Na+-blocker to 200 ms−1 ·M−1 resulted in a reduced lifespan (1,305 and 1,459 ms for [DNa+ ] = 60 and 45 µM, respectively). The lifespan for the combined Bl·INa + Comb·Bl·IX was consistently shorter than that of any individual applications of Bl·INa or Comb·Bl·IX alone in all cases. A similar augmented anti-arrhythmic effect (shown as shortened lifespan of re-entry) was also observed for the combined Bl·INa + Bl·IKur for [DNa+ ] = 45 µM and KI = 200 ms−1 · M−1 (**Figure 8Aiv**) but not for the rest of the cases.

A consistent decrease in the DF was observed in the drugmodulated re-entrant excitations as compared to those in the drug-free condition (**Figure 8Bi–iv**). In the drug-free condition, the DF extracted from the pECG was 8.63 Hz, which is within the range of similar clinical data (Jarman et al., 2012). Applying Na+- or K+- block individually resulted in slowing of the rate of the rotors, and this was dependent on the concentrations and parameters of the blockers. For Bl·IKur the DF was 8.25 Hz

comparison. In the left column (i,iii) of (A,B), Bl·INa was simulated with *KI* = 100 ms−1 · M−1 ; in the right column (ii,iv), *KI* = 200 ms−1 · M−1 . In the top panels (i,ii) of (A,B), [*DK*+ ] = 3.2µM, [*DNa*+ ] = 60µM; in the bottom panels (iii,iv), [*DK*+ ] = 2.4µM, [*DNa*+ ] = 45 µM. In all panels, *KA* = 1 ms−1 · M−1 , *LA* = 1 ms−1 , *LI* = 0.01 ms−1 .

with the control dose and 8.30 Hz for the reduced dose. In the simulations with Comb·Bl·IX, the DF was 6.63 Hz and was not affected by the 25% reduction in the dose of the compound. For Bl·INa alone the DF was 7.81 Hz and substantially smaller (6.24 Hz) for Na+-blockers of KI = 100 ms−1 · M−1 and KI = 200 ms−1 · M−1 , respectively. An enhanced deceleration of the rotors was observed for Bl·INa + Bl·IKur in all cases. The DF for simulations with Comb·Bl·IX + Bl·INa was not computed due to the short lifespan in these events.

#### 3D Simulations

The antiarrhythmic effects of acacetin and Na+-block on the reentrant waves in the cAF-remodeled atria were also evaluated using our 3D anatomical model of the human atria (Colman et al., 2013, 2017). A 10-s episode of sustained re-entrant excitation was first initiated in the cAF-remodeled atria in the drugfree condition; this produced the initial conditions of the 3D model for additional 10-s episode simulations. Next, the behavior of electrical waves of another 10-s episode simulation in the drug-free condition and after applying the selected blockers was analyzed and compared. **Figure 9A** shows snapshots of excitation wave evolution following the 10-s episode of re-entrant excitation events. pECGs computed from the excitation waves are shown in **Figure 9Bi–iv**. **Figures 9C,D** illustrates power spectrum analyses of the pECGs and a comparison of the lifespan of the electrical waves in the drug-free condition and after applying the drugs/blockers. In the drug-free condition, stable reentrant waves around the pulmonary veins and left atrium were observed; the power spectrum density (PSD) manifests a singlefocused peak around 8.16 Hz. Following applying Bl·INa, the reentrant waves became less organized and also decelerated. This was characterized by a smaller dominant frequency (6.58 Hz) and less focused PSD distribution; the excitation was not terminated. Following the application of Comb·Bl·IX, the re-entrant wave soon became unstable and eventually disappeared after t = 3,279 ms. Note that the peak PSD amplitude was much smaller and its distribution was much broader as compared to the drug-free condition. The combined drugs further destabilized the re-entrant waves and reduced the lifespan to approximately 1,120 ms.

# Simulated Effects of Acacetin and Na+-Current Blocker on the QT Interval

Further simulations were performed using a 1D ventricular transmural strand model to evaluate the changes in the waveform of electrograms in consequence of applying the drugs. A comparison of the computed electrogram waveforms are illustrated in **Figures 10A,B**, and the QT intervals are quantified in **Figure 10C**. Blocking small fractions of IKr, IKs and Ito in the human ventricles by acacetin slightly increased the QT interval by 21 ms. The electrograms were not noticeably affected by applications of the Na+-blocker with the simulated parameters. Results from both cases did not show dramatic QT prolongation, indicating no dramatic effects affecting ventricular repolarization process which might promote ventricular arrhythmias.

# DISCUSSION

Even decades after goal-directed work, successful development of effective and safe antiarrhythmic drugs for treating AF has not been accomplished and remains a major unmet clinical need. In a recent study on the canine heart, enhanced anti-arrhythmic effects and AF-selectivity of INa blockade by additional IKur block was demonstrated (Aguilar et al., 2015). Whether similar effects could be obtained in the human atria, especially following cAF-induced electrical remodeling which reduces IKur, remained unclear. How the combined Na+ and K+-block modulates the QT interval also remained incompletely understood. In this study, the effects of IKur (combined with a modest block in Ito, IKr, and IKs as presented by acacetin, a compound shown to be effective in anti-AF treatment) and INa block (two potentially effective atrial-selective block on human atrial electrophysiology) were investigated in silico using multiscale models of the human atria and state-dependent block scheme. The simulation results demonstrate that both Na+-block and K+-block exhibited anti-arrhythmic effects in the atria following cAF-remodeling, despite reduced IKur by the remodeling. The present study highlighted that in addition to combined Na+- and K+-block, combined multi-K+-channels also exerted beneficial synergistic antiarrhythmic effects when compared with single channel block whilst having modest impact on ventricular repolarization (QT interval). This study suggests that multi-channel block (either combined Na+-K+-block, or combined multi-K+-block) may be a favorable strategy for the development of novel pharmaceutical therapies for AF.

# Effects of INa block

An atrial-ventricular difference in the properties of INa, especially in the voltage dependence of steady-state inactivation, has been reported (Li et al., 2002; Burashnikov et al., 2007; Chen et al., 2016; Fan et al., 2016; Caves et al., 2017). In these studies, the voltage dependent steady-state inactivation curves for INa were found to be negatively shifted (by 5–16 mV) in atrium as compared to the ventricular parameters. This difference gave rise to an on-going interest in developing an atrial-selective blocker of INa as a strategy in terminating AF (Burashnikov et al., 2007; Antzelevitch and Burashnikov, 2009; Zygmunt et al., 2011; Morotti et al., 2016; Caves et al., 2017). As was done in previous studies (Aguilar-Shardonofsky et al., 2012; Aguilar et al., 2015), in this study, the kinetic parameters in drug actions of Na+ blockers were varied over wide parameter spaces to reveal AFselectivity of Na+-blockers in the ventricles and fibrillating atria. Our results demonstrated that in the presence of AF-remodeling, an atrial-selective block of INa could produce different effects between atrial and ventricular cells (**Figures 3**, **4**) and that the AF-selectivity could be maximized by optimizing the binding and unbinding rates of the Na+-blocker (**Figure 6**). Note also that, the fractional inhibition of INa by the Na+-blocker exhibited a substantial dependence on the rate of pacing (**Figures 3**, **4** and Online Supplementary Material 2.2), which was quantified using the rate-selectivity (**Figure 6**).

At the cellular level, Na+-block resulted in a significant inhibition in INa at fast pacing rates, but minor effects within

the range of normal heart rates in the atria (**Figures 5**, **6**). The antiarrhythmic effects of these changes were demonstrated in simulations of multicellular atrial tissue. In a 1D atrial strand model, applying Na+-block progressively enhanced the rateadaptations of Vmax and CV over a larger range of BCLs, whereas the atrial APD was not affected. At fast pacing rates, Vmax and CV were decreased significantly, suggesting reduced excitabilities of atrial myocytes (**Figure 5**). These results are in concordance with the recent study (Aguilar et al., 2015) where similar effects of Na+-block on the canine atria were demonstrated in silico and experimentally in coronary-perfused hearts.

In 2D tissue simulations, this Na+-block shortened the lifespan and caused slowing in the excitation rate of the spiral waves (**Figures 7**, **8**). In the 3D anatomical model, applying Na+-block alone produced antiarrhythmic effects by slowing the re-entrant excitations (**Figure 9**). Furthermore,

in our 1D model of transmural ventricular strand, the simulations suggested that the Na+-block had minimal impact on ventricular repolarization, as judged by modest QT interval prolongation. These results demonstrated that Na+-block could be beneficial in suppressing re-entrant activities in the cAF-remodeled atria, with modest impact on ventricular repolarization.

# Effects of K+-Current Block

K +-current blockers delay the repolarization phase of the AP and thus prolong atrial APD and refractory period. This can cause disruptions and eventually termination of the re-entrant circuits (Hancox et al., 2016). However, K+-channel blockers such as dofetilide and sotalol (which potently inhibit IKr) have a substantial risk of prolonging QT interval and promoting Torsades de pointes arrhythmias (Hondeghem and Snyders, 1990; Yap and Camm, 2003). In principle, blocking atrial-specific K+ channels may exert antiarrhythmic effects in the atria while minimizing potential risks of adverse effects in the ventricles.  Interestingly, many existing IKur blockers potently block other K +-channels including Ito and IK,ACh (Gögelein et al., 2004; Wirth et al., 2007; Burashnikov and Antzelevitch, 2008; Li et al., 2008). The additional blockades of these channels may contribute to the antiarrhythmic effects of those drugs, which warrant further investigations.

In this study, acacetin, a compound initially isolated from the traditional Chinese medicine Xuelianhua, was selected as a representative IKur blocker. The effects of acacetin on atrial electrophysiology were evaluated in two ways: (a) the effects of acacetin blocking IKur only; and (b) the full actions of acacetin on the targeting channels (Ito, IKur, IKr, and IKs) (Li et al., 2008). This approach allowed for investigations into the effects of IKur block alone as well as the potential benefits of additional-but-modest inhibition of other K+-currents in the human atria.

#### Selective IKur Block

Blocking IKur with 3.2µM acacetin exerted APD prolongation (9.8 ms) under the baseline/normal conditions (**Figure 2**). Experimental data show that dependent on the baseline AP waveform the effect of IKur block on human atrial APD70−90 under normal (SR) conditions can manifest as prolongation or shortening in the APD, (Workman et al., 2001; Wettwer et al., 2004; Schotten et al., 2007; Burashnikov and Antzelevitch, 2008; Loose et al., 2014). Additionally, the prolongation in APD by IKur block observed in the present study is similar to our previous paper (Colman et al., 2017) concerning the effects of genetically down-regulated IKur. Moreover, inhibiting IKur under normal conditions elevated the AP plateau potential and prolonged APD30 (**Figure 2**). Both effects matched well with experimental studies (Workman et al., 2001; Wettwer et al., 2004; Schotten et al., 2007; Burashnikov and Antzelevitch, 2008; Loose et al., 2014) and our simulation study (Colman et al., 2017).

We note that in the cAF-remodeling cells, a more pronounced prolongation in APD (by 23.6 ms for 1 Hz and 16.2 ms at 6 Hz, **Figure 4**) was observed in the presence of 3.2µM acacetin, despite that this current was down-regulated by cAF-remodeling (Wagoner et al., 1997; Brandt et al., 2000; Christ et al., 2008). These results are in accordance with previous experimental results of blocking IKur with MK-0448 (Pavri et al., 2012; Loose et al., 2014). In addition, IKur block exhibited enhanced ratedependent adaptations in APD both at the cellular (Online Supplementary Material 2.2) and 1D strand models (**Figure 5**). Importantly, the CV restitution curve shifted toward higher BCLs, indicating that this tissue is less capable of conduction of atrial excitation waves at high rate while maintaining conduction of slow waves (**Figure 5**). In 2D tissue simulations, applying IKur block alone (3.2µM acacetin) destabilized the cores of rotors (i.e., potential organizing centers for AF), and slightly slowed their excitation rates, but failed to terminate them (**Figure 7**), suggesting a limited efficacy of terminating AF by IKur block alone. Similarly, a recent modeling study by Aguilar et al. (2017) suggested that the antiarrhythmic efficacy of IKur block was substantially decreased in the presence of AF-induced electrical remodeling. Also, the experimental study (Burashnikov and Antzelevitch, 2008) showed that block of IKur by 4-AP of small doses had limited efficacy in suppressing AF in canine atria. This may represent the fact that IKur is reduced at high frequencies (as discussed/suggested in Feng et al., 1998a; Burashnikov and Antzelevitch, 2008; Wu et al., 2011) and shown in **Figure 1**) as well as by cAF-induced remodeling (Wagoner et al., 1997; Brandt et al., 2000; Christ et al., 2008). In addition, IKur is primarily active during phase 2 of AP, and hence pure IKur block exerted a relatively greater prolongation in APD30 than APD90 (**Figure 3**), in contrast to other K+-block including dofetilide which mediates anti-AF effects by prolonging the terminal phase of the AP (Roukoz and Saliba, 2007).

In this study, IKur block was simulated using a state-dependent block model, which successfully reproduced the use- and ratedependent inhibition of acacetin (**Figure 1**). The rate-dependent block of IKur exerted a higher fractional inhibition in the current at faster pacing rates, which likely produces greater anti-AF effects in the presence of high-frequency excitations as seen during AF. Along with the previous modeling studies on investigating effects of IKur block (Almquist et al., 2010; Scholz et al., 2013; Ellinwood et al., 2017), this study demonstrated the importance of explicitly considering the kinetic properties of the block in computational efforts of understanding the consequences and underlying mechanisms of IKur block.

#### Effects of Combined K+-Current Block

The combined K+-block (as exhibited by acacetin and many other IKur blockers) resulted in synergistic APD prolongation as well as an increased efficacy in terminating re-entry in tissue as compared to the pure IKur block.

Note that at the single myocyte level, the combined actions of acacetin produced greater prolongation in atrial APD than the sum of changes due to drug-induced block of individual channel in normal and cAF-remodeled myocytes (**Figure 2**). Additionally, the combined K+-block increased the rateadditivity of APD as compared to the pure IKur block (Online Supplementary Material 2.2). This was also consistently observed in the 1D simulation (**Figure 5**). In the setting of pure IKur block, the elevated and prolonged plateau phase of the AP could promote the activation of IKr/IKs, which in return may accelerate the repolarization of AP-phase 3 (Colman et al., 2017). Therefore, additional inhibition in IKr by an identical fraction is expected to result in a greater APD prolongation than a pure IKur or IKr/IKs block.

In 2D simulations, the combined K+-block produced an enhanced efficacy in suppressing AF compared with the pure IKur block: promoting meandering of rotor tips (**Figure 7B**), shortening the lifespan of re-entries (**Figure 8A**) and slowing of spiral wave excitations (**Figure 8B**). Rotor meandering is one mechanism by which spiral waves may meet non-conducting boundaries to extinguish re-entry (Narayan et al., 2013; Pandit and Jalife, 2013; Rappel et al., 2015).

The effects of acacetin (3.2µM) on the ventricular AP and QT interval was assessed in a single cell model and 1D transmural strand model by assuming similar blockade effects of the compound on the human ventricles and atria. It was shown that following applying acacetin, the ventricular repolarization and QT interval was both preserved with slight prolongations around 21 ms (**Figure 10**). Our results are close to the previous experimental study (Li et al., 2008) showing that QT intervals were not prolonged by acacetin in isolated rabbit hearts and anesthetised dogs.

The synergistic effects demonstrated by the combined K+ blocks have implications on developing novel pharmaceutical anti-AF therapies. Given that Ito, IKr and IKs contribute to the repolarizations of ventricular APs, inhibitions in these channels may promote risks of side effects in the ventricles. In this regard, combined block of atrial-specific K+ channels may be favorable. Recently, another two families of K+-channels that are dominantly expressed in the atria have been acknowledged: the small-conductance Ca2+-activated K+ (SK) channels (ISK) (Qi et al., 2014), and the two-pore K+ (K2P3.1) channel (ITASK−1) (Schmidt et al., 2015), further to the well-known constitutively active acetylcholine-activated K+ current (IK,ACh). Combined block of these atrial-specific channels may exert greater and safer antiarrhythmic effects in the atria, warranting future investigations.

#### Synergistic Effects of Combined Na+- and K +- Block

The present study reveals novel and significant synergistic effects of combined block of Na+- and K+-currents (INa and pure-IKur/multi-K+-block) and demonstrates the additional synergistic anti-arrhythmic effects derived from the multi-K+ channel block in cAF-remodeled atria.

In cAF-remodeled atria, combined Na+- and K+-block significantly increased the fractional INa inhibition and APD prolongation (**Figures 3**, **4**) and promoted pronounced AP alternans at 6 Hz, with complex effects in human AF (Narayan et al., 2011). In the simulations varying the blockade kinetics of INa block, the combined block dramatically augmented the attainable maximal AF-selectivity in consequence of enhanced atrial-selectivity and rate-selectivity as compared to the pure Na+-block (**Figure 6**).

In the 1D model of an atrial strand, combined Na+ and K +-block produced synergistic reductions in Vmax and CV; the threshold of BCL allowing a 1:1 conduction was increased as compared to the control conditions (**Figure 5**). In simulated re-entrant waves in 2D and 3D atria, the combined Bl·INa + Comb·Bl·IX exhibited a greater efficacy in suppressing AF, with a decreased lifespan of rotors as compared to that by either individual block (**Figures 7**–**9**). Although the combined Bl·INa + Bl·IKur did not further reduce the lifespan of spiral waves as compared to the Bl·INa alone, the combination did lead to the extinction of one of the two rotors (**Figures 7Aiv,Biv**) and deceleration of re-entrant activations (**Figure 8B**). Follow-up simulations showed that consistent synergistic antiarrhythmic effects could be obtained with reduced doses of Na+- and K+ blockers.

The non-specific multi-channel blockade is increasingly recognized as a strategy for pharmaceutical therapy of AF both experimentally (Sicouri et al., 2010; Aguilar et al., 2015; Kirchhoff et al., 2015; Hartmann et al., 2016) and clinically (Koskinas et al., 2014; Reiffel et al., 2015). In a previous study (Aguilar et al., 2015), synergistic anti-arrhythmic effects were demonstrated both in silico and experimentally in healthy canine hearts. Additionally, the favorable synergistic antiarrhythmic effects have also been reported in combined block of ISK and INa in an experimental atrial-fibrillated guinea pig model (Kirchhoff et al., 2015). Also, the recent HARMONY trial (Reiffel et al., 2015) revealed synergistic AF-suppressing effects for combined use of ranolazine and dronedarone. While revealing the synergistic effects of combined Na+- and K+- block in cAF-remodeled human atria, this study supports and adds insights into the ongoing efforts in developing multi-channel block as a strategy for the treatment of AF.

#### Limitations and Future Work

In the absence of the required detailed experimental data, when simulating effects of acacetin on Ito, IKs and IKr, the dosedependence block of acacetin was assumed to be identical in both human atrial and ventricular cells. This assumption may warrant further investigations. In addition, the parameters of atrial Ito have been reported to be different from those of ventricular Ito in human (Amos et al., 1996). Previous studies reported that the IC50 of 4-AP block of atrial Ito was one-third of that of ventricular Ito (Amos et al., 1996; Nattel et al., 2000). If a similar atrial-vs.-ventricular difference in the IC50 of Ito and/or IKr/IKs could exist for acacetin, the effects of acacetin on the ventricular electrophysiology would be less significant than our simulations, which might result in to a smaller change in ventricular INa and APD for the combined block of Bl·INa and Comb·Bl·IX, and thus enhance the computed atrial-selectivity and AF-selectivity of the combined block. Given that applying acacetin in vivo did not prolong QT intervals in isolated rabbit hearts and anesthetised dogs (Li et al., 2008), any significant prolongation of the ventricular APD and QT interval is unlikely (**Figures 3B**, **4C**, **10**). Therefore, our assumption of no atrial-ventricular difference in the potency of acacetin on K+-currents may not affect our conclusions concerning the atrial-selectivity of combined Na+ and K+-block.

Additionally, in the absence of detailed experimental data for state-dependent block of Ito, IKr, and IKs by acacetin, the block of these channels was modeled using a single pore block model. The IC50 values (**Table 1**) were determined by fitting the concentration-response relation of the step current at 40 mV in previous experimental studies (Li et al., 2008; Wu et al., 2011). A recent study suggests that the IC50 values may be dependent on the voltage protocols applied, and this cannot be reflected by single pore models. In future studies, the pore block model for Ito, IKr, and IKs can be replaced by a state-dependent block model when such experimental data become available. Also, the present work did not attempt to model the effects of acacetin on IK,ACh, although the study shows the current is potently blocked by the compound. The 2D and 3D simulations of atrial tissue, while validated, may not fully capture the complexity of fibrosis-tissue interfaces which are seen in structurally remodeled atria and were not simulated in these monodomain experiments.

Thirdly, our simulation results showed a moderate QT prolongation of around 20 ms following applying both INa blocker and acacetin. While a QT prolongation of less than 5 ms does not raise a regulatory concern (Committee for Medicinal Products for Human Use, 2012), implications of QT prolongations between 5 and 20 ms remain inconclusive (Committee for Medicinal Products for Human Use, 2005). In the present study, though the extent of QT prolongation of 20 ms is far less than the threshold of discontinuation criteria of 60 ms as indicated in Committee for Medicinal Products for Human Use (2005), it would indeed raise a positive flag in thorough QT tests and necessitate extended safety assessment and intensive patient monitoring during late stages of trials (Committee for Medicinal Products for Human Use, 2012). On the other hand, the approach we used in accounting for the effects of acacetin on ventricular myocytes may result in upper bound of QT prolongation, since the potency of acacetin was assumed to be identical in atria and ventricles.

Fourthly, the threshold in BCL inducing AP alternans was increased by K+-block (**Figure 5**, Figure S5). However, AP alternans seen at slower pacing rates has been linked with occurrence of AF (Narayan et al., 2002, 2011). Therefore, the increased threshold in BCL developing AF by K+-block can be potentially proarrhythmic. The safety of K+-block and its proarrhythmic potential in the atria should be addressed in future studies.

Fifthly, there are limitations in the approaches used in simulating the INa blockers in single myocytes and tissue. Similar to previous studies (Aguilar-Shardonofsky et al., 2012; Aguilar et al., 2015), in our simulations, the drug action on INa was modeled through a state-dependent block assuming drugs binding to both activated and inactivated states of INa, and the gating variables of INa were modeled using an Hodgkin-Huxley scheme. The limitations in this approach outlined in Aguilar-Shardonofsky et al. (2012) therefore apply in the present study. The results of the use-dependent block may be affected by the models used (Aguilar-Shardonofsky et al., 2012). However, the previous study (Aguilar-Shardonofsky et al., 2012) compared this modeling scheme with simulations using a Markov model, showing qualitative agreement in major findings. Therefore, the major conclusions drawn from this study may not be affected by the selected modeling approach for INa and drug interactions. Furthermore, in optimizing the AFselectivity of the INa and K+-current blockers, the concentrations of Na+ and K+ blockers were fixed at 60 uM. This may potentially impose limitations in discomposing the role of the binding parameters in the modulatory effects of the blockers because of the very slow kinetics of the drug binding to its targeted channel at this high concentration. It warrants further studies by varying the concentration of blockers to simulate the optimized effects of the AF-selectivity of INa blocker. In tissue simulations, effects of drugs were modeled by increasing their doses homogeneously, simultaneously and instantaneously. The realistic actions of INa blockers in tissue, however, may be different. Also, in tissue simulations a homogenous cell model was used. As previous study (Feng et al., 1998b) showed atria are electrically heterogeneous, future work is needed to

assess how tissue heterogeneities affect the efficacy of atrialselective pharmaceutical interventions. Furthermore, the current simulations did not take the cardiac autonomic regulation into account in order to take into considerations of acacetin on IKACh. Future studies on interactions of atrial-selective anti-arrhythmic drug actions and autonomic systems may also render valuable findings.

It is important to acknowledge that administration of class Ic agents for Na+-block can cause cardiac arrhythmia and increased mortality (Echt et al., 1991). Further investigations are therefore warranted to assess the safety of the simulated Na+-block in the heart, especially in the ventricles.

# CONCLUSIONS

By using state-dependent drug block models and our mathematical models of the human atria, the antiarrhythmic effects of atrial selective Na+- and K+-blockers on the cAFremodeled atria were evaluated. The combined block of multiple K +-currents as well as simultaneous block of Na+- and K+ currents produced synergistic antiarrhythmic effects. Our results suggest that developing multi-channel (multiple K+ currents and/or combined Na+- and K+-current) block is a potentially valuable strategy for the treatment of AF.

## REFERENCES

# AUTHOR CONTRIBUTIONS

HZ and HN conceived the study. HN designed experiments, developed and validated computational models, and performed numerical experiments. HN, DW, and WW analyzed data. HN, DW, WW, WG, SN, and HZ interpreted data and wrote the manuscript.

# FUNDING

This work was supported by grants from the British Heart Foundation FS/14/5/30533, EPSRC (UK) (EP/J00958X/1; EP/I029826/1), MC-IRSES CORDIS3D (317766), NSFC (61179009), Shenzhen Science and Technology Innovation Committee (JCYJ20151029173639477; JSGG20160229125049615). SN is supported by grants from the National Institutes of Health (R01 Hl83359).

## SUPPLEMENTARY MATERIAL

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

London: EMEA. Available online at: http://www.ema.europa.eu/docs/en\_GB/ document\_library/Scientific\_guideline/2009/09/WC500002879.pdf

isolated heart model of atrial fibrillation. Heart Rhythm 12, 409–418. doi: 10.1016/j.hrthm.2014.12.010

and clinical studies. Circ. Arrhythm. Electrophysiol. 8, 1325–1333. doi: 10.1161/CIRCEP.115.002956

**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.

Copyright © 2017 Ni, Whittaker, Wang, Giles, Narayan and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Trigger vs. Substrate: Multi-Dimensional Modulation of QT-Prolongation Associated Arrhythmic Dynamics by a hERG Channel Activator

Michael A. Colman1 \*, Erick A. Perez Alday 2 , Arun V. Holden1 and Alan P. Benson1

*1 School of Biomedical Sciences and Multidisciplinary Cardiovascular Research Centre, University of Leeds, Leeds, United Kingdom, 2 Division of Cardiovascular Medicine, Oregon Health and Science University, Portland, OR, United States*

#### Edited by:

*Eleonora Grandi, University of California, Davis, United States*

#### Reviewed by:

*Sara Dutta, United States Food and Drug Administration, United States Kazuharu Furutani, Osaka University, Japan*

\*Correspondence:

*Michael A. Colman<EMAIL_ADDRESS>
#### Specialty section:

*This article was submitted to Computational Physiology and Medicine, a section of the journal Frontiers in Physiology*

Received: *02 August 2017* Accepted: *19 September 2017* Published: *04 October 2017*

#### Citation:

*Colman MA, Perez Alday EA, Holden AV and Benson AP (2017) Trigger vs. Substrate: Multi-Dimensional Modulation of QT-Prolongation Associated Arrhythmic Dynamics by a hERG Channel Activator. Front. Physiol. 8:757. doi: 10.3389/fphys.2017.00757* Background: Prolongation of the QT interval of the electrocardiogram (ECG), underlain by prolongation of the action potential duration (APD) at the cellular level, is linked to increased vulnerability to cardiac arrhythmia. Pharmacological management of arrhythmia associated with QT prolongation is typically achieved through attempting to restore APD to control ranges, reversing the enhanced vulnerability to Ca2+-dependent afterdepolarisations (arrhythmia triggers) and increased transmural dispersion of repolarisation (arrhythmia substrate) associated with APD prolongation. However, such pharmacological modulation has been demonstrated to have limited effectiveness. Understanding the integrative functional impact of pharmacological modulation requires simultaneous investigation of both the trigger and substrate.

Methods: We implemented a multi-scale (cell and tissue) *in silico* approach using a model of the human ventricular action potential, integrated with a model of stochastic 3D spatiotemporal Ca2+ dynamics, and parameter modification to mimic prolonged QT conditions. We used these models to examine the efficacy of the hERG activator MC-II-157c in restoring APD to control ranges, examined its effects on arrhythmia triggers and substrates, and the interaction of these arrhythmia triggers and substrates.

Results: QT prolongation conditions promoted the development of spontaneous release events underlying afterdepolarisations during rapid pacing. MC-II-157c applied to prolonged QT conditions shortened the APD, inhibited the development of afterdepolarisations and reduced the probability of afterdepolarisations manifesting as triggered activity in single cells. In tissue, QT prolongation resulted in an increased transmural dispersion of repolarisation, which manifested as an increased vulnerable window for uni-directional conduction block. In some cases, MC-II-157c further increased the vulnerable window through its effects on *I*Na. The combination of stochastic release event modulation and transmural dispersion of repolarisation modulation by MC-II-157c resulted in an integrative behavior wherein the arrhythmia trigger is reduced but the arrhythmia substrate is increased, leading to variable and non-linear overall vulnerability to arrhythmia. Conclusion: The relative balance of reduced trigger and increased substrate underlies a multi-dimensional role of MC-II-157c in modulation of cardiac arrhythmia vulnerability associated with prolonged QT interval.

Keywords: QT interval, action potential duration, arrhythmia trigger, arrhythmia substrate, hERG activators, computational modeling

#### INTRODUCTION

Conditions in which the QT interval of the electrocardiogram (ECG) is prolonged, such as heart failure (Hart, 1994) and inherited or acquired long QT syndromes (LQTS) (Schwartz et al., 2012), are associated with an increased risk of ventricular arrhythmias (Tomaselli and Zipes, 2004; Moss and Kass, 2005). The prolonged QT interval reflects prolongation of the ventricular cellular action potential duration (APD), which can result in arrhythmias through an increase in cell-level arrhythmia triggers and/or modification of tissue-level arrhythmia substrates (Kalin et al., 2010; Benson et al., 2011a; Tse, 2016).

Cell-level triggers arise because delayed ventricular repolarisation modifies membrane and subcellular Ca2+ handling (Clusin, 2003; Neme ˇ c et al., 2016), leading to either re-activation of the L-type Ca2+ current and early afterdepolarisations (Lankipalli et al., 2005), or Ca2+ overload of the sarcoplasmic reticulum (SR), causing spontaneous SR Ca2+ release events and delayed-afterdepolarisations through activation of the forward-mode Na+-Ca2+ exchange current (INaCa), which in turn can result in triggered activity if the delayed afterdepolarisation is of sufficient magnitude (Janse, 2004). Cardiomyocytes typically exhibit a threshold dependence of the occurrence of spontaneous release events as a function of SR Ca2+ load, wherein the probability of a spontaneous release event rapidly rises from 0 to 1 within a critical region of SR Ca2+ (Venetucci et al., 2007; Campos et al., 2015). Although originating at the cell level, triggers need to be coordinated at the tissue level in order to develop into arrhythmias: a critical compact region of tissue simultaneously exhibiting triggered behavior is required to initiate propagation of the trigger (Noble, 1972; Clayton et al., 2011; Bezekci et al., 2015; Campos et al., 2015).

An increase in the tissue-level substrates for arrhythmias (that is, the necessary conditions for triggered activity to propagate and develop into arrhythmias) arise because APD prolongation is rarely homogenous in and between the different regions of the ventricles (e.g., transmurally, or from base to apex) (Antzelevitch, 2005; Glukhov et al., 2010). This heterogeneous APD prolongation increases the spatial dispersion of repolarisation, potentially leading to regions of recovered (i.e., excitable) tissue partially bordered by still refractory (i.e., unexcitable) tissue. A propagating trigger event occurring in such a location can be partially blocked by the refractory tissue, leading to re-entrant arrhythmias (Pandit and Jalife, 2013). The spatiotemporal region where such partial conduction block could occur is termed the "vulnerable window" (VW) (Starmer et al., 1993; Shaw and Rudy, 1995; Benson et al., 2008, 2011a). It follows that heterogeneous APD prolongation increases not only the spatial dispersion of repolarisation, but also the VW, i.e., the arrhythmia substrate.

The interaction of triggers and substrates determines the initiation of arrhythmia: an arrhythmia cannot be initiated without both a suitably-sized and -timed trigger and the necessary substrate to allow that trigger to propagate in a reentrant manner (Kalin et al., 2010). Arrhythmias, therefore, are not cellular events, but tissue-level events.

Management of arrhythmias associated with QT prolongation can be achieved by attempting to restore the ventricular APD to control ranges (Nachimuthu et al., 2012), thus reversing the increases in arrhythmogenic triggers and substrates associated with APD prolongation. One such strategy is the use of human ether-a-go-go-related-gene channel (hERG) activators that enhance the repolarising rapid delayed rectifier K+ current (IKr), thus reducing the APD (Grunnet et al., 2008; Wu and Sanguinetti, 2016). However, many anti-arrhythmic drugs have pro-arrhythmic effects (Kumar and Zimetbaum, 2013); such drugs can shorten the action potential and the QT interval (reducing arrhythmias associated with a prolonged QT interval), but they may have additional and unintended effects that increase (rather than reduce) the propensity for arrhythmias under certain conditions. For example, we have shown in a previous experimental and computational study that the hERG activator NS1643, one of the most effective and best characterized hERG activators (Hansen et al., 2006), successfully restores APD toward healthy durations and reduces arrhythmia triggers, but is associated with an increase in the VW, i.e., the substrate for arrhythmias, due to effects on the post-repolarisation refractory period (Peitersen et al., 2008). Furthermore, while low concentrations of NS1643 activates IKr and shortens APD, it has been shown that higher concentrations of NS1643 blocks (rather than activates) IKr (Bilet and Bauer, 2012).

The potentially pro-arrhythmic effects of hERG activators, such as NS1643 has prompted the search for novel hERG activators that do not display these effects. One recentlyidentified compound is MC-II-157c, an NS1643 analog. MC-II-157c activates IKr at low concentrations and, unlike NS1643, it continues to activate IKr at high concentrations; it may also block the sodium current (INa) (Guo et al., 2014). However, it remains unknown how the IKr activation and INa block seen with this new compound affect arrhythmia triggers and substrates, and importantly, how these (increased or decreased) triggers and substrates interact to induce arrhythmias (if at all).

Computational models provide powerful research tools to understand the intricacies of arrhythmia trigger and substrate interaction, as they allow us to modify parameters under precisely controlled conditions and quantify the resultant tissue-level arrhythmic behavior, and can predict how these arrhythmias will manifest in a clinical setting (e.g., changes to the ECG). We therefore used a computational approach to quantify the interaction of pharmacologically-modified arrhythmia trigger and substrate, using the novel hERG activator MC-II-157c as an example.

We wanted to quantify the modified triggers and substrate that result from MC-II-157c ion channel actions. To this end, we use detailed single cell and tissue level models to study the effect of QT prolongation and its modulation by MC-II-157c on: (i) APD heterogeneity in isolated cells; (ii) SR Ca2+ loading and subsequent spontaneous SR Ca2+ release events; (iii) the probability of spontaneous SR Ca2+ release manifesting as triggered action potentials in single cell and ectopic activity in tissue; and (iv) the vulnerability to the initiation of re-entrant like conduction patterns.

#### METHODS

We implemented a multi-scale in silico approach to study the interactions between arrhythmia trigger and substrate in conditions associated with prolonged QT intervals, and their modulation by the hERG activator MC-II-157c.

# Isolated Cell Models – Intracellular Ca2+ Handling

In order to simulate triggered and ectopic activity underlain by spontaneous Ca2+ release events, a spatial model of intracellular Ca2+ handling which explicitly accounts for stochastic state transitions and spatial coupling is required. We therefore implemented an efficient, idealized reduction of a previously developed and validated general model of spatio-temporal Ca2+ handling with realistic structure (Colman et al., 2017), using a similar approach implemented by other groups (e.g., Restrepo et al., 2008). Briefly, 15 × 20 × 65 spatially-discrete individual calcium release units (CRUs) were modeled throughout the geometry of the cell (**Figure 1**). Each CRU comprises of five compartments with associated Ca2+ concentrations: the intracellular spaces of the dyadic cleft space ([Ca2+]ds), subspace ([Ca2+]SS) and bulk-cytosolic space ([Ca2+]i), and the network and junctional SR ([Ca2+]nSR, ([Ca2+]jSR). The bulk cytosol, subspace and network SR are diffusively coupled to neighboring CRUs; the dyadic cleft space and junctional SR are not spatially coupled to neighbors. The fundamental model equations describing this system are:

$$\begin{aligned} \frac{\mathbf{d}[\text{Ca}^{2+}]\_{\text{i}}}{\text{dt}} &= \beta \left(\mathbf{D}\nabla^{2}[\text{Ca}^{2+}]\_{\text{i}} + \varphi\_{\text{i}} + (\upsilon\_{\text{SS}}/\upsilon\_{\text{i}})f\_{\text{SS}}\right) \\ \frac{\mathbf{d}[\text{Ca}^{2+}]\_{\text{SS}}}{\text{dt}} &= \mathbf{D}\nabla^{2}[\text{Ca}^{2+}]\_{\text{SS}} + \varphi\_{\text{SS}} - f\_{\text{SS}} + (\upsilon\_{\text{ds}}/\upsilon\_{\text{SS}})f\_{\text{ds}} \\ \frac{\mathbf{d}[\text{Ca}^{2+}]\_{\text{nSR}}}{\text{dt}} &= \beta\_{\text{nSR}} \left(\mathbf{D}\nabla^{2}[\text{Ca}^{2+}]\_{\text{nSR}} + \varphi\_{\text{nSR}} - \left(\upsilon\_{\text{j}\text{SR}}/\upsilon\_{\text{nSR}}\right)f\_{\text{j}\text{SR}}\right) \\ \frac{\mathbf{d}[\text{Ca}^{2+}]\_{\text{ds}}}{\text{dt}} &= \varphi\_{\text{ds}} - f\_{\text{ds}} \\ \frac{\mathbf{d}[\text{Ca}^{2+}]\_{\text{j}\_{\text{SR}}}}{\text{dt}} &= \varphi\_{\text{j}\text{SR}} - f\_{\text{j}\text{SR}} \end{aligned}$$

where transfer between compartments is given by:

$$\begin{array}{rcl} J\_{\rm SS} & = \left( \left[ \text{Ca}^{2+} \right]\_{\rm SS} - \left[ \text{Ca}^{2+} \right]\_{\rm i} \right) \text{\tiny\tau}\_{\rm SS}^{-1} \\ J\_{\rm ds} & = \left( \left[ \text{Ca}^{2+} \right]\_{\rm ds} - \left[ \text{Ca}^{2+} \right]\_{\rm SS} \right) \text{\tiny\tau}\_{\rm ds}^{-1} \\ J\_{\rm jSR} & = \left( \left[ \text{Ca}^{2+} \right]\_{\rm jSR} - \left[ \text{Ca}^{2+} \right]\_{\rm nSR} \right) \text{\tiny\tau}\_{\rm jSR}^{-1} \end{array}$$

where vx is the volume of compartment x, βx is an instantaneous buffering term, and ϕx is the reaction term. Stochastic dynamics are modeled for the RyRs and LTCCs, which are part of the dyadic cleft space reaction term. All parameters and reaction terms are given in the Supplementary Material.

#### Isolated Cell Models - Action Potential Model

Ion currents were described by a simplified version of the O'Hara-Rudy dynamic (ORd) human ventricular cell model (O'Hara et al., 2011), wherein the major currents only (INa, IKr, IKs, IK1, Ito, INaK) were included, without the further details of phosphorylation included in the original study; Ca2+ currents (ICaL, INaCa, ICap, ICab) are described by the Ca2+ handling model. These simplifications were implemented to improve computational efficiency and for integration with the general stochastic Ca2+ handling model described above. The model provides specific formulations to describe the heterogeneity in ionic currents of the transmural cell types (endocardial, midmyocardial, and epicardial) found in the left ventricular free wall of the human heart (see O'Hara et al., 2011) for details.

This integrated stochastic framework captures spontaneous Ca2+ release events that could lead to triggered activity (see details in Simulating Spontaneous Release Events in Tissue Models below). Default model action potentials and cytosolic Ca2+ transients, for the three transmural cell types of the simplified ORd model with stochastic Ca2+ handling, are shown in **Figure 1B**. The updated cell model exhibits action potential and Ca2+ transient properties similar to the original cell model and within the range of experimental data presented in the original study (O'Hara et al., 2011): APD = 272–360 ms and intracellular Ca2+ transient magnitude of ∼0.6µM for the three cell-types during control pacing; note that the present study does not consider heterogeneity in the intracellular Ca2+ handling system and thus the Ca2+ transient is more homogeneous between the cell types than in the original study (see Limitations). The model is therefore considered suitable for the mechanistic study undertaken.

#### Prolonged QT and Pharmacological Modulation

We were interested in general cases of QT prolongation rather than modeling the kinetics of very specific conditions, while still making our results broadly applicable to clinical conditions, such as LQTS. We therefore simulated QT interval prolongation (i.e., prolongation of the ventricular APD) in one of three ways: (i) downregulation of the slow delayed rectifier K+ current (IKs) maximal conductance by 50%, similar to LQTS1, which we term "prolonged QT variant a" (PQTa) in the remainder of the manuscript; (ii) downregulation of IKr maximal conductance by

50% (PQTb), similar to LQTS2; and (iii) upregulation of ICa,L maximal conductance by 50% (PQTc), similar to LQT8 (Bohnen et al., 2017).

Effects of 10µM of the hERG activator MC-II-157c on IKr were modeled by modifying the IKr formulation according to experimental data (Guo et al., 2014): maximal conductance was decreased by 12%, activation was shifted by −14 mV and inactivation by +14 mV, and deactivation kinetics were slowed 3.3-fold (note that, although the maximal conductance of IKr is reduced by MC-II-157c, its kinetic effects enhance the activity of the current; see The hERG Activator MC-II-157c Partially or Fully Reverses APD Prolongation Heterogeneously). Effects of MC-II-157c on blocking INa (which are not as well characterized as its effects on IKr; Guo et al., 2014) were simulated by reducing the maximal conductance of INa by 0, 40, and 80%.

# Tissue Models

We used a 20 mm 1D virtual tissue strand (Kléber and Rudy, 2004) for quantifying transmural propagation and vulnerability, with equal spatial distributions of endocardial, midmyocardial and epicardial cells. A 20 × 40 mm 2D tissue sheet (Clayton et al., 2011) was used for simulations examining propagation of triggered activity, with equal distributions of endocardial, midmyocardial and epicardial tissue in the x direction. For examining intramural propagation and ectopic activity in 3D, we used an anatomically detailed 3D ventricular wedge model, obtained by diffusion tensor MRI, and used equal proportions of endocardial, midmyocardial and epicardial tissue in the transmural direction (see Benson et al., 2007, 2011b; Walton et al., 2013 for details).

All tissues were isotropic, i.e., conduction velocity was set to be equal in all directions: We used an electrical diffusion coefficient of D = 0.048 mm2ms−1 , to give a conduction time along the 1D strand (i.e., a transmural activation time) of 40 ms (cf. Glukhov et al., 2010), and a plane wave conduction velocity of 0.5 m.s−1 in all tissues. The body surface potential was computed by placing the ventricular wedge model in a human torso mesh; the forward problem was solved by a boundary element method, as has been described in previous studies (Perez Alday et al., 2015, 2016, 2017). ECGs were derived from the body surface potential by selecting elements of the torso mesh which correspond to the ECG electrodes.

# Simulating Spontaneous Release Events in Tissue Models

Performing tissue level simulations using the fully detailed spatial Ca2+ handling model to describe individual cells is computationally intractable due to the large number of equations that need to be solved in such situations. Furthermore, the focus of this study was not to dissect the mechanisms of spontaneous Ca2+ release in single cell, but rather to understand the considerations determining the manifestation of sub-cellular Ca2+ release as triggered action potentials and propagating electrical excitation in the presence of prolonged APD and its pharmacological modulation.