The O'Hara-Rudy dynamics (ORd) model (O’Hara et al., 2011) was utilized to simulate the electrophysiology of human ventricular myocytes in this study. The ORd model is a comprehensive human cell model that was created using human experimental data. To overcome its unphysiologically slow conduction velocity (Elshrif and Cherry, 2014), the original I _Na in the ORd model was substituted with that in the Tusscher et al. biophysically detailed model (TNNP06 model) (Ten Tusscher and Panfilov, 2006).

A conventional Hodgkin-Huxley model of a cardiac cell was implemented at the cellular level, with the model equation being: ∂ V m ∂ t = − ( I i o n + I s t i m ) C m (1) where V _m is the membrane potential, I _ion is the sum of all transmembrane ionic currents, and I _stim is the externally applied stimulus current. C _m is the membrane capacitance.

The cell model of heart failure (HF) used in this study was based on Elshrif et al.’s research (Elshrif et al., 2015), where a collection of HF-induced ion channel remodeling effects were incorporated into the ORd model. Similarly, the effects of CORM-2 (i.e., a CO-releasing molecule) were modeled based on previous research by Al-Owais et al. (Al-Owais et al., 2021) and were incorporated into the healthy and HF cell models. The reason we chose CORM-2 rather than CO is that CORM-2 is one of the most common CO-releasing molecules in biological research, and is safer and more controllable than CO. More details can be found in Sections SII and SIII in the Supplementary Material.

Available experimental data regarding the effects of ranolazine and HW-0168 from previous studies have been summarized in Table 1 (Antzelevitch et al., 2004; Rajamani et al., 2008; Beyder et al., 2012; Moreno et al., 2013; Dong et al., 2019). Specifically, ranolazine has been shown to exert dose-dependent blocking effects on I _Na (Beyder et al., 2012), I _NaL (Antzelevitch et al., 2004), I _NaCa (Antzelevitch et al., 2004), I _CaL (Antzelevitch et al., 2004), I _Kr (Rajamani et al., 2008). Dose-response curves for ranolazine-affected ion channels were fitted using the following Hill functions:

I _Na f N a RAN = 1.0 1.0 + ( [ RAN ] / 53.6 ) 2.4 (2)

I _NaL f N a L RAN = 1.0 1.0 + ( [ RAN ] / 6.23 ) 1.0 (3)

I _NaCa f N a C a RAN = 1.0 1.0 + ( [ RAN ] / 91.0 ) 1.48 (4)

I _CaL f C a L RAN = 1.0 1.0 + ( [ RAN ] / 296.0 ) 1.0 (5)

I _Kr f N a RAN = 1.0 1.0 + ( [ RAN ] / 12.0 ) 1.0 (6)

where [RAN] is the dose of ranolazine used in experiments.

The fitting results are illustrated in Figure 1A. For the I _Kr activator, HW-0168, only I _Kr was reported to be affected by the drug (Dong et al., 2019); therefore, the data were fitted using Eq. 7:

I _Kr f K r H W = 1.8 1.0 + ( 0.41 / [ H W ] ) 0.73 + 1 (7) where [HW] is the dose of HW-0168 used in the experiment.

The fitted dose-dependent curve is illustrated in Figure 1B. In this study, we used 10 μM and 0.5 μM for ranolazine and HW-0168, respectively. The above-fitted equations were finally incorporated into the ‘Healthy + CO’ and ‘HF + CO’ cell models.

The ionic current under the action of the drug is calculated by Eq. 8: I i o n D r u g = I i o n ⋅ f i o n D r u g (8) where f i o n D r u g represent the effect of a drug on a certain ionic current.

In addition to ranolazine and HW-0168, we also selected six multi-channel blockers (i.e., amiodarone, verapamil, nifedipine, quinidine, vanoxerine, and bepridil) and four specific I _Kr activators (i.e., KB130015, ICA-105574, NS1643, NS3623) for efficacy simulation and screening of the drugs. A simple pore block theory (Brennan et al., 2009) was used in this study to model the interactions between drugs and ion channels. Based on this theory, the effect of drugs blocking ion channels was fitted by the following formula: θ = 1 1 + ( I C 50 / [ D ] ) n H (9) where θ is the blocking efficiency, [D] is the concentration of a drug, IC ₅₀ is the half-maximal inhibitory concentration, and nH is the Hill coefficient.

The effect of drugs activating ion channels was fitted by Eq. 11: Y = A c t max − 1 1 + ( E C 50 / [ D ] ) n H (10) where Y is the activation efficiency, and Act _max is the maximum activation efficiency, EC ₅₀ is the compound concentration resulting in 50% of the Act _max .

The six multi-channel blockers act on related ion channels in a dose-dependent manner, and the related parameters are listed in Table 2. To evaluate the drug efficacy more objectively, we explored all drugs at three doses based on their C_max, as shown in Table 3. The four specific I _Kr activators activated I _Kr currents in a dose-dependent manner as well, and the relevant parameters are shown in Table 4.

The 1D transmural ventricular strand model, which is a linear syncytium formed by coupling multiple cells, can be calculated by adding a diffusion term to the cell model equation: ∂ V m ∂ t = D ( ∂ 2 V m ∂ x 2 ) − I i o n C m (11) where D is the scalar diffusion coefficient that decides the conduction velocity of APs.

The 1D transmural strand was 15 mm long, which was close to the normal range of the human transmural ventricle width (∼4.0–14.0 mm) (Drouin et al., 1995; Yan et al., 1998). The strand was discretized into 100 interconnected nodes with a spatial precision of 0.15 mm, which was consistent with the reported cell length [i.e., 80–150 μm (Hinrichs et al., 2011)]. The proportions for transmural cell types were set to 25:35:40 for ENDO, MID, and EPI cells, which were identical to that used in previous studies (Zhang and Hancox, 2004; Luo et al., 2017). Such proportions reliably reproduced a positive T wave in the computed pseudo-ECG under control (healthy) conditions. The diffusion coefficient D was set to 0.127 mm²/ms, giving a CV of planar excitation waves of 70 cm/s through the strand, which matched well with the experimental data from human ventricles (Taggart et al., 2000).

To further quantify the critical size of EAD cells for overcoming the source-sink effect and initiating triggers in ventricular tissue, we simulated a 15 mm homogenous ventricular strand consisting of only MID cells for the failing heart, with the center of the strand (Figure 2, red region) containing a variable number of contiguous cells affected by CO. The number of cells in the susceptible region was gradually increased until the synchronously occurred EADs overcame the source-sink effect and trigger a premature beat. The critical cell number was recorded as a metric for measuring the susceptibility to arrhythmias.

The pseudo-ECG was calculated from the constructed 1D strand model by the following equation (Gima and Rudy, 2002): ϕ e ( x ′ ) = a 2 4 ∫ ( − ∇ V m ) ⋅ [ ∇ 1 r ] d x (12) where ϕ e is a unipolar potential generated by the strand, a is the radius of the strand, dx is the spatial resolution, and r is the Euclidean distance from a point x to another point x′.

As shown in Figure 3, the period from the earliest appearance of the QRS complex to the end of the T-wave was defined as the QT interval, measured in milliseconds. The end of the T-wave was defined as the return of the descending limb to the TP baseline.

To demonstrate the robustness of the reported findings, we constructed cell population models with reference to previous studies (Britton et al., 2013; Sutanto and Heijman, 2020). Specifically, the maximum conductance of the nine major ionic currents (I _Na, I _NaL, I _CaL, I _Kr, I _Ks, I _K1, I _to, I _NaCa, and I _NaK) in the original deterministic model was scaled by a group of factors that follow a normal distribution with mean 1.0 and standard deviation 0.2. In this way, 1,000 population model variants were obtained.

The CV dynamic restitution curves were obtained using a dynamic pacing protocol. Specifically, the 1D strand model was paced with a certain basic cycle length (BCL) until reading its steady state upon which the CV value was recorded for that BCL. The initial BCL was set to 3,000 ms and was decreased gradually until the model failed to produce excitation waves. Based on the ‘CV-BCL’ pairs generated by the above protocol, CV restitution curves could be plotted against BCL.

Similar to the 1D strand model, the monodomain equation (Eq. 11) was adopted to describe the propagation of excitation waves in the ventricular slice. Isotropic propagation was assumed, and the diffusion coefficient D was set to 0.154 mm²/ms, to produce a CV of 0.74 m/s (Taggart et al., 2000). The spatial step was set to 0.15 mm to be consistent with that in 1D models. To mimic the physiological characteristics of the Purkinje fibers, a series of supra-threshold stimuli were applied to several pacing sites on the endocardium of the slice.

The influences of the aforementioned drugs were also evaluated under the heart failure condition. Simulated actions of ranolazine on CO-affected cells and tissues of heart failure are presented in Figure 8. Overall, ranolazine exacerbated the CO and heart failure-induced arrhythmias. In detail, the CO-induced 2:1 alternated EADs in MID cells became 1:1 consecutive EADs (Figure 8Aii), resulting in complete repolarization failure. Ranolazine also led to the occurrence of EAD in EPI cells (Figure 8Aiii). Above EAD activities in single cells did not develop into ectopic beats in 1D ventricular strands due to the ‘source-sink’ effect (Xie et al., 2010); however, ranolazine resulted in the 1:1 conduction failure of excitation waves at the pacing frequency of 1.25 Hz (Figure 8Bi). For the pseudo-ECG, ranolazine did not eliminate the CO-induced ECG morphological changes in heart failure tissue and further led to failed depolarization due to the considerably prolonged repolarization phase of the last cycle (Figure 8Bii).

Figure 9 shows the effects of the other six multi-channel blockers on ECG morphology in heart failure conditions. Due to the remodeled transmural gradient of repolarization in the heart failure condition, the T-wave was almost flattened. In terms of the QT-interval, amiodarone (0.0005 μM), nifedipine (0.005 μM), and verapamil (0.03 μM) had almost no effect on the QT interval, and bepridil (0.01 μM) slightly prolonged the QT interval. In addition, quinidine (1 μM) and vanoxerine (0.005 μM) caused depolarization failure. Overall, all six drugs were not effective against CO-induced arrhythmias in heart failure conditions.

The baseline model of HF + CO showed that CO could induce pronounced EAD activities in MID cells, but these EADs did not evolve into ectopic beats in 1-D tissue due to the ‘source-sink’ effect (i.e., the depolarization force of EAD is not able to trigger an excitation due to the limited number of EAD cells) (Xie et al., 2010). Applying ranolazine did not trigger ectopic beats in the tissue either; however, it did diminish the repolarization ability in terms of the cellular action potential (Figure 8A). To give a more intuitive presentation of the increased proarrhythmic risk of ranolazine, we quantified the risk by measuring the critical number for generating the ectopic beat. Specifically, we constructed a 1D model of HF MID cells, with its central segment being set to CO-affected, and the minimum number of affected cells that could overcome the source-sink effect and lead to ectopic beats was recorded as the critical cell number. As shown in Figure 10, simulations suggested that the critical cell number under CO conditions was 68, corresponding to a tissue length of 10.2 mm. In contrast, the critical cell number was only 58 after the addition of ranolazine, which suggested an increased susceptibility to ectopic beats. Action potentials of representative cells within the CO-affected region (marked ‘*’ and ‘**’ in Figure 10) were plotted in the right panels of Figure 10.

In our previous study (Jiang et al., 2022), we have shown that the suppression of I _Kr is the main factor responsible for the CO-induced prolongation of APD and QT interval. Considering the critical role of I _Kr in the pathological pathway and the bad efficacy of multi-channel blockers, we evaluated several specific I _Kr activators in this section. For simplicity, the simulation results of a representative drug HW-0168 (full name: N-(2-(tert-butyl)phenyl)-6-(4-chlorophenyl)-4-(trifluoromethyl) nicotinamide) (Dong et al., 2019) are presented in detail (Figures 11, 12), whereas only the effective dose is recorded for the other activators (Table 4).

On the ‘healthy + CO’ condition, it can be observed that the HW-0168 at a dose of 0.5 μM [therapeutic range suggested in clinical: 0.5–1 μM (Dong et al., 2019)] effectively shortened the APD prolongation caused by CO and reversed the prolonged APD₉₀ to almost the same as the healthy condition. Generated pseudo-ECGs using 1D transmural ventricular strand models showed consistent results——HW-0168 restored the prolonged QT interval to a level that was almost identical to the control condition (Figure 11Bii). In addition, HW-0168 also improved the conduction properties of excitation waves and shortened the conduction wavelength of the tissue (Figure 11Bi).

The efficacy of HW-0168 under heart failure conditions is presented in Figure 12. Simulation results showed that HW-0168 effectively reversed the proarrhythmic effects (i.e., prolonged APDs and EADs) of CO in all three cell types (Figure 12A), and shortened the excitation wavelength in the heart failure tissue (Figure 12Bi). For the ECG, although the drug did not restore the altered T-wave morphology in heart failure, it eliminated the QT interval prolongation effects by CO.

According to the above results, the selective I _Kr activator achieved desired treatment for CO-induced arrhythmias. Therefore, more existent I _Kr activators (i.e., KB130015 (Gessner et al., 2010), ICA-105574 (Asayama et al., 2013), NS1643 (Casis et al., 2006), NS3623 (Hansen et al., 2006)) were tested and the doses of drugs under which the QT-interval was restored were recorded in Table 6. According to our simulation results, ICA-105574 was the most sensitive one, which restored the QT-interval and suppressed EADs (under heart failure conditions) at a dose of only 0.25 μM.

Considering the potential influence of intercellular or intersubject variability on the reported findings, we built cell population models and performed additional simulations based on them. The simulation results are shown in Figure 13. It can be observed that EADs occurred occasionally under the HF condition, with a ratio of only 2.6%. Next, after considering the effects of CO, APDs of cell populations were generally prolonged, and the ratio of cells with EAD increased to 18.5%. The administration of ranolazine aggravated the situation, and the ratio of EAD cells increased dramatically to 58.2% (as shown in panel Aiii). In contrast, the addition of HW effectively alleviated the above arrhythmogenesis at the cellular level, which was evidenced by the complete suppression of EAD activities and the generally shortened APDs.

To avoid the potential difference caused by the simplified model geometry, we conducted simulation experiments for two representative drugs, i.e., ranolazine and HW-0168, using a 2D realistic ventricular slice model. The simulation results are shown in Figure 14. It can be observed obviously that the tissue slice with ranolazine took more time to repolarize than that with HW-0168 (Figure 14A). In terms of the ECG, the 2D-based ECGs are consistent with the 1D-based ones (Figure 14B). For example, ranolazine further prolonged the QT interval based on CO and therefore exacerbated the proarrhythmic effect. On the other hand, HW-0168 still exerted the antiarrhythmic effects of ranolazine by restoring the QT interval.

The severe cardiotoxic consequences of CO urgently require an effective therapeutic strategy to treat them. In this study, we evaluated the efficacy of various multi-channel blockers and specific I _Kr activators against CO-induced ventricular arrhythmias in healthy and failing hearts. The major findings are as follows: 1) The tested existent antiarrhythmic drugs failed to rescue the heart from CO-induced arrhythmias, and most of them even aggravated the arrhythmogenic condition, which was evidenced by the more frequent EAD activities and decreased critical cell numbers for triggering ectopic beats. 2) In contrast, specific I _Kr activators demonstrated good efficacy according to the improved biomarkers at both cellular and tissue levels. All of the tested I _Kr activators restored the prolonged QT intervals in both healthy and heart failure conditions, and the EADs in MID cells were successfully suppressed as well. 3) In-depth case analysis with ranolazine and HW-0168 revealed the critical role of I _Kr in the CO-induced functional changes in cardiac electrophysiology, and neither I _CaL nor I _NaL blockers were able to offset the decreased repolarization forces caused by the CO-induced I _Kr inhibition. 4) Of note, the drug ranolazine was previously suggested as a potential strategy in dealing with CO-induced arrhythmogenesis due to its good efficacy demonstrated in rats, and the failure of ranolazine in the human tissue in this study hinted the crucial role of inter-species variances when determining the pharmacotherapeutic strategy.

Ranolazine was first suggested in Dallas et al.’s study (Dallas et al., 2012) for the treatment of CO-induced arrhythmias. Based on the experimental results obtained from rats, they proposed that CO-induced EADs arouse from the activation of NO synthase, which in turn leads to the NO-mediated nitrosylation of Na_V1.5 and the enhanced I _NaL. Correspondingly, the I _NaL inhibitor ranolazine abolished the EADs and was considered to be effective in dealing with CO-induced arrhythmias. Similarly, Morita et al. also observed the antiarrhythmic effects of ranolazine for its suppression of reentrant and multifocal ventricular fibrillation in rat ventricles (Morita et al., 2011). However, APs in rats are distinctly different from those in humans, and the such discrepancy may lead to species-dependent effects of the same drug. This hypothesis was explored in Al-Owais et al.’s study (Al-Owais et al., 2017), where the effects of ranolazine were examined in guinea pigs—a species with action potentials more closely resembling that of humans. Interestingly, ranolazine failed to abolish CO-induced EAD and even exacerbated such proarrhythmic factors.

Our simulations suggested that ranolazine exerted similar proarrhythmic effects in human hearts. Specifically, ranolazine further prolonged AP durations and QT intervals in healthy human simulations (Figure 4), while in heart failure conditions it led to more pronounced EADs in MID and EPI cells (Figure 6). The above model-dependent effects of ranolazine arose from the differences of I _Kr, a major outward current responsible for the repolarization in human APs but are almost negligible in rat myocytes (Pandit et al., 2001). Although the I _NaL inhibition effects of ranolazine tend to suppress EAD, the drug can also reduce the repolarization force by inhibiting I _Kr. Further assessment using a 1D homogenous ventricular strand consisting of only MID cells found that ranolazine decreased the critical cell number for triggering ectopic beats (from 68 to 58), which also suggested the increased arrhythmogenic risk of the drug. These findings provide new insights into the side effects of ranolazine on the treatment of CO-induced arrhythmias. They also highlighted that the drug effects obtained in rats need to be carefully interpreted in clinical trials due to the species-dependent differences.

In addition to ranolazine, we evaluated more existent antiarrhythmic drugs to find potential drug strategies for CO-induced arrhythmias. Calcium current blockers were focused on in hopes of attenuating the depolarization force in the plateau phase and therefore shortening the action potential and the QT interval. However, none of the six drugs was able to rescue the heart from arrhythmogenesis, and most of them even worsened the conditions, evidenced by the further prolonged QT intervals and more frequently observed EAD activities. By analyzing the separate role of each channel current in the integral effect of multi-channel drugs, we found that blocking I _CaL and I _NaL was not able to offset the reduction of I _Kr by CO; furthermore, most of these multi-channel blockers also inhibited I _Kr with a relatively low affinity. Indeed, the hERG channel that conducts I _Kr is a highly sensitive target and it accounts for the majority of drug withdrawal events in the last 2 decades (Brown, 2004; Stockbridge et al., 2013; Villoutreix and Taboureau, 2015). On the other hand, there are few drugs available in the current antiarrhythmic category exerting I _Kr activating effects (Lei et al., 2018), making it difficult to find a proper drug strategy. We have also tried pinacidil (an I _KATP activator) in the model, but it did not produce any significant efficacy as well (data not shown). This can be attributed to the fact that the K-ATP channel barely opens under normoxic conditions due to its ATP-sensitive characteristic (Dart and Standen, 1995); therefore, the I _KATP would not make obvious differences even a high magnification ratio was used in the model.

In-depth analysis has demonstrated that I _Kr plays a major role in CO-induced arrhythmogenesis (Jiang et al., 2022). Considering that existent multi-channel antiarrhythmic drugs did not achieve idealized efficacy, we turned to evaluate the potential phrenological effect of specific I _Kr activators. In line with expectations, the simulation results showed that I _Kr activators could effectively reverse the proarrhythmic effects of CO. All the tested drugs notwithstanding in different doses restored AP and ECG morphologies almost to their control levels in healthy human simulations, and they also suppressed EADs and ectopic beats in heart failure human simulations. These findings suggest that the I _Kr activator is a promising pharmacotherapy for the treatment of CO-induced arrhythmias.

This study lacks validation of heart failure models. Though we have adopted a well-established cell model under heart failure conditions and replicated several known electrophysiological changes in failing hearts, for example, the prolonged APD (Akar and Rosenbaum, 2003; Lou et al., 2012), the decreased conduction velocity (Akar et al., 2004), the widened QRS complex (Shenkman et al., 2002; Sandhu and Bahler, 2004), and the prolonged QT interval (Davey et al., 2000; Medina-Ravell et al., 2003); however, we did not find enough tissue-level experimental data to validate other observations such as the flattened T-wave.

The above limitations shall not change the main conclusions of this study. Specifically, most of the observations and conclusions in the present study were based on the damaged cellular repolarization and the consequent QT prolongation in failing hearts, which were well-established in biological experiments (Davey et al., 2000; Medina-Ravell et al., 2003; Lou et al., 2012; Ng et al., 2014). In addition, for the EAD phenomenon, we adopted relatively conservative parameters (i.e., no EAD phenomenon occurred in pure heart failure conditions) to avoid exaggeration of the experimental results.

The experimental data on CO effects, drugs, and currents used in this study were obtained from different species, and the CO effects were obtained at room temperature. Interspecies differences and temperature dependence should be taken into account when interpreting and translating the results. The effects on APD in this study were measured in individual isolated ventricular myocytes, and the potential cell-coupling effects on the APD in high-dimensional models were not considered. Besides, the pathological model of CO was constructed based on experimental data obtained from different CORM-2 doses (10–30 μM) (Al-Owais et al., 2021), which should be considered in future studies. As for the drugs, the I _Kr activators proposed in this study for the treatment of CO-induced arrhythmias currently face some disadvantages and unknowns. Specifically, compared with the FDA-approved drugs such as ranolazine and amiodarone, I _Kr activators represented by HW-0168 are currently only used in biological experiments and simulation experiments, and their effective doses have not been clinically verified and side effects are not being disclosed. Moreover, whether these I _Kr activators interact with ion channels other than I _Kr remain unknown. If this is the case, then they must be treated as multiple-channel drugs and the potential offset or synergy effects among the involved ion currents should be considered.

Finally, according to our previous research review (Zhang et al., 2021), CO was also known to affect multiple cellular pathways other than the ion channels in this study. The present study mainly considered arrhythmias caused by changes in ionic currents directly induced by CO, without considering the mitochondrial toxicity of CO and some other complicated electrophysiological remodeling induced by cellular ischemia. Specifically, CO poisoning will increase ROS and RNS (Piantadosi, 2008), which further impair the chondrial energetics and can alter the intracellular calcium handling as well (Hegyi et al., 2021). This alteration will subsequently impact the expression and trafficking of channels (Sutanto et al., 2020). These cellular pathways warrant further investigations in the future.