In order to simulate triggered and ectopic activity underlain by spontaneous Ca²⁺ release events, a spatial model of intracellular Ca²⁺ 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 Ca²⁺ 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 Ca²⁺ concentrations: the intracellular spaces of the dyadic cleft space ([Ca²⁺]_ds), subspace ([Ca²⁺]_SS) and bulk-cytosolic space ([Ca²⁺]_i), and the network and junctional SR ([Ca²⁺]_nSR, ([Ca²⁺]_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:

where transfer between compartments is given by:

where v_x 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.

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 (I_Na, I_Kr, I_Ks, I_K1, I_to, I_NaK) were included, without the further details of phosphorylation included in the original study; Ca²⁺ currents (I_CaL, I_NaCa, I_Cap, I_Cab) are described by the Ca²⁺ handling model. These simplifications were implemented to improve computational efficiency and for integration with the general stochastic Ca²⁺ 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 Ca²⁺ 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 Ca²⁺ transients, for the three transmural cell types of the simplified ORd model with stochastic Ca²⁺ handling, are shown in Figure 1B. The updated cell model exhibits action potential and Ca²⁺ 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 Ca²⁺ 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 Ca²⁺ handling system and thus the Ca²⁺ 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.

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 (I_Ks) maximal conductance by 50%, similar to LQTS1, which we term “prolonged QT variant a” (PQTa) in the remainder of the manuscript; (ii) downregulation of I_Kr maximal conductance by 50% (PQTb), similar to LQTS2; and (iii) upregulation of I_Ca,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 I_Kr were modeled by modifying the I_Kr 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 I_Kr 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 I_Na (which are not as well characterized as its effects on I_Kr; Guo et al., 2014) were simulated by reducing the maximal conductance of I_Na by 0, 40, and 80%.

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 mm²ms⁻¹, 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⁻¹ 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.

Performing tissue level simulations using the fully detailed spatial Ca²⁺ 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 Ca²⁺ release in single cell, but rather to understand the considerations determining the manifestation of sub-cellular Ca²⁺ release as triggered action potentials and propagating electrical excitation in the presence of prolonged APD and its pharmacological modulation.

We therefore implemented a “non-spatial” simplification of the Ca²⁺ handling model (described in Isolated Cell Models–Intracellular Ca²⁺ Handling, above) for use in tissue-level simulations, to capture spontaneous Ca²⁺ release events at significantly reduced computational cost and with complete controllability. The non-spatial model consists of a single CRU (with RyR and LTCC dynamics solved deterministically) with additional analytical functions which describe the RyR waveform associated with whole-cell spontaneous release events, derived from analysis of the fully detailed spatial cell model, similar to the approach used in Campos et al. (2015) and Colman et al. (2015). For a simple transient-spike morphology (Figure 2) this function has the form:

where t_i is the initiation time of the spontaneous release and t_up, t_decay and N_{RyR_peak} describe the shape of the waveform and are all determined from the duration (Figures 2, 3 and described below). The function for the plateau-like waveform (corresponding to durations longer than 250 ms) is derived from the same parameters (note the peak time, t_p, is t_i+t_up):

Thus, the waveform is completely described by the initiation time and the duration. The stochastic nature of spontaneous Ca²⁺ release is captured by randomly selecting these parameters from functions describing their physiological distributions. The distributions describing initiation time of whole-cell spontaneous release events in the spatial cell model are typically skewed (Figure 3Aa) and well approximated by two sigmoidal functions, split around the cumulative probability of 0.25 (corresponding to a specific initiation time, t_{i_sep}; Figure 3Ab). The initiation time can therefore be determined by passing a random number into the inverse of the two sigmoidal functions:

where the gradient parameters of the two sigmoidal functions (k_F1, k_F2) and the t_i at the cumulative probability of 0.25, t_{i_sep}, completely control the resulting distribution. The duration, D, of the RyR waveform can be determined from distributions in an analogous manner:

where MD refers to the median and k_{1_MD} and k_{2_MD} are functions of the median (in conditions where longer waveforms are observed, the variability in waveform duration between simulations is larger; Figure 3Ba–c):

Thus, the distribution is entirely described by the median. Finally, t_up and N_{RyR_peak} can be determined from the given duration (Figure 3C):

Where MRYR refers to the median and is given by:

And if duration > 250 ms, it is also necessary to compute N_{RyR_plateau}:

These formulations therefore allow complete control over spontaneous release dynamics through just four parameters (k_F1, k_F2, t_{i_sep}, median duration). We fix t_{i_sep} in simulations such that only the width of the initiation time distribution (k_F1, k_F2) and the median of the duration distribution are varied: t_i distribution widths of 350, 550, and 1,000 ms (Figure 3Ac) were used in tissue simulations; duration medians of 150, 200, 250, 300, and 350 ms (Figure 3Bd) were used in single-cell and tissue simulations. Derivation of these equations from the fully detailed spatial cell model ensures self-consistency and physiological validity of the resulting waveforms.

Models were coded in C/C++ and run on a Linux desktop machine, or using the University of Leeds ARC2 High Performance Computing facilities. Equations for isolated cell models were solved using a forward Euler method with a time step of Δt = 0.05 ms; ion channel gating equations were solved using the Rush-Larsen scheme (Rush and Larsen, 1978). For tissue models, the monodomain equation was solved using a Forward Time Centred Space method: Space steps of Δx = 0.2 mm were used in the 1D strand, Δx = Δy = 0.2 mm in the 2D tissue, and Δx = Δy = 0.425 mm and Δz = 0.5 mm in the 3D wedge model (as defined by the diffusion tensor MRI dataset). Parallelisation was implemented with OpenMP. Cell APD was measured from the time the membrane potential crossed −80 mV during the upstroke of the action potential, to when the membrane potential crossed back over −80 mV during the repolarisation phase.

The single cell model was paced at a cycle length of 1,000 ms for 100 beats under control (WT), PQT and PQT + MC-II-157c conditions to evaluate the efficacy of MC-II-157c on reversing PQT induced APD prolongation (Figure 4). Figure 4A shows effects on AP morphology of the PQTa, PQTb, and PQTc conditions in endocardial, midmyocardial, and epicardial cells, with WT action potentials shown as a reference. The PQTb condition (downregulation of I_Kr) has the largest effect on prolonging APD (from 360 to 649 ms in midmyocardial cells, an increase of 80%), due to the larger I_Kr conductance in human cells relative to I_Ks conductance, and its primary role in AP repolarisation. All three PQT conditions increase the transmural difference in APD (the difference between the longest and shortest APDs in the three cell types), from 88 ms in WT to 97, 279 and 122 ms with PQTa, PQTb, and PQTc, respectively.

Figure 4B shows how the hERG activator MC-II-157c reduces APD back towards control levels in all cell types and with all PQT conditions. The drug reduced APD back to control levels in PQTa (downregulation of I_Ks) due to the minimal effect this condition has on initially prolonging APD, but the largest decrease was seen with PQTb in midmyocardial cells, where APD was reduced from 649 to 540 ms, a reduction of 17% (Figure 4C). It should be noted, however, that MC-II-157c has transmurally heterogeneous effects; that is, the degree of APD reduction seen in the three different cell types is not identical, with the drug under PQTb conditions (for example) giving a 5% decrease in endocardial cells, 17% in midmyocardial cells and 4% in epicardial cells. Consequently, the maximal transmural difference in APD in PQTb reduces from 279 ms with the PQT condition alone to 184 ms with PQTb plus MC-II-157c, but does not reduce transmural difference in APD back down to WT levels (184 vs. 88 ms).

ECGs were computed for the different conditions using the 3D ventricular wedge model under normal pacing (BCL = 1,000 ms). QT-prolongation was observed for all three remodeling types (Figure 4C), with PQTb exhibiting the longest QT-interval, congruent with single cell results (QT = 310 ms in WT compared to 326, 418, and 329 ms in PQTa-c, respectively). MC-II-157c resulted in a delay in the QRS peak as well as earlier absolute repolarisation time and consequent shortening of the QT interval (QT = 309 ms, 403 and 310 ms in PQTa-c + MC-II-157c; Figure 4C); in PQTa and PQTc, MC-II-157c fully reverses QT prolongation.

The mechanism by which MC-II-157c shortens APD, despire reducing the maximal conductance of I_Kr, is illustrated in Figure 5: the shifts in the activation and inactivation curves and the slowing of deactivation kinetics result in an increased current during both voltage clamp and AP clamp experiments. These results are generally congruent with the original study of Guo et al. (2014), although the current traces do differ in morphology and extent of effect of MC-II-157c (see Limitations).

The vulnerability to the emergence of whole-cell spontaneous Ca²⁺ release events (such as intracellular Ca²⁺ waves, Figure 6A) is primarily controlled by the dynamics of the intracellular Ca²⁺ handling system and the SR Ca²⁺ load, wherein cells typically exhibit an SR load threshold above which the probability of spontaneous release events significantly increases (Wagner et al., 2015). The focus of this study is not on the mechanisms of spontaneous Ca²⁺ release, and therefore investigation of Ca²⁺ handling remodeling is beyond its scope. However, APD prolongation associated with PQT, and its subsequent reversal by MC-II-157c, may influence SR loading and consequently the vulnerability to the emergence of whole-cell spontaneous Ca²⁺ release.

An SR-loading protocol was used to analyse this behavior, by pacing the spatial cell model at a rapid rate (cycle length of 400 ms). The maximal flux rate of intracellular uptake was increased, simulating the effect of sympathetic stimulation, in order to promote SR loading. An increase by a factor of two was chosen as this loaded the SR-Ca²⁺ in the WT model to just above the spontaneous release threshold; the most suitable region to reveal the consequence of APD modulation on spontaneous activity.

The time series of SR-Ca²⁺ in WT and the PQTc condition in isolated epicardial cells illustrates the effect on SR-Ca²⁺ loading under rapid pacing and highlights that APD prolongation associated with PQT can promote loading (Figure 6B). The peak of the SR Ca²⁺ concentration during this time provides a measure of SR Ca²⁺ loading; these data are shown in Figure 6C for WT and all PQT conditions with and without MC-II-157c. The PQTa condition (downregulation of I_Ks) does not promote SR loading compared to WT (peak SR Ca²⁺ load of 1.07 mM in WT and 1.06 mM in PQTa), but the PQTb (downregulation of I_Kr) and PQTc (upregulation of I_Ca,L) conditions do (peak SR Ca²⁺ load of 1.106 and 1.13 mM, respectively), with PQTc having the largest effect due to the increase in the transmembrane Ca²⁺ current in this condition. MC-II-157c reverses SR Ca²⁺ loading in PQTb, reducing peak SR Ca²⁺ load at periodic steady-state to 1.05 mM, which is below WT levels. However, the drug has only a very small effect in PQTc (down to 1.108 mM) due to the minimal effect MC-II-157c has on APD in this condition, and the subsequently small change in the time course of I_Ca,L.

These relatively small changes in SR-Ca²⁺ can manifest as significant differences in the vulnerability to spontaneous Ca²⁺ release due to the non-linear threshold dependence on SR Ca²⁺ concentration (Figures 6D,E). The probability of whole-cell spontaneous release events and the probability distributions describing the initiation time can be computed from a large set of simulations (N = 1,000 per condition). Example distributions are shown in Figure 6E for the WT and PQTc with and without MC-II-157c, highlighting that even the small change in SR Ca²⁺ as a result of MC-II-157c significantly reduces the probability of spontaneous release (~50% in PQTc with MC-II-157c compared to 98% in PQTc alone) as well as widening the distribution (although not fully reversed to WT). Due to the choice of loading parameters giving the WT close to threshold, no spontaneous release occurs for either PQTa or PQTb under the application of MC-II-157c as in these conditions threshold SR Ca²⁺ is not reached.

The effect of MC-II-157c on the probability of DADs manifesting as full triggered action potentials was investigated using the simplified, non-spatial model such that spontaneous release waveforms could be directly controlled. The initiation time was set to 1,000 ms and behavior of the cell models for WT, PQT and PQT + MC-II-157c was compared over multiple simulations (N = 1,000 per condition) for four different RyR waveform duration distributions (see details in Simulating Spontaneous Release Events in Tissue Models; the distribution determines the values of duration which can be selected from a random number input and thus the actual value of the duration will vary randomly within the 1,000 simulations according to the given distribution).

Examples of 100 simulations for WT and PQTb, with and without MC-II-157c, are shown in Figure 7A. In PQTb, more DADs turn into triggered activity compared to in WT (top two panels), while simulated application of MC-II-157c reduces the occurances of triggered activity (bottom two panels). These data are summarized for all conditions and two duration distributions (median 300 and 350 ms, see Simulating Spontaneous Release Events in Tissue Models) in Figure 7B, wherein the degree of block of I_Na associated with MC-II-157c is varied (0, 40, and 80% block). In all cases, the I_Kr modification reduces the probability of triggered activity (defined as the number of simulations in which triggered activity occurred as a proportion of the total), but the role of I_Na is less clear; blocking I_Na is important, but the degree of I_Na block has different effects depending on the condition. The mechanism of the drug's action is shown in Figure 7C: principally, an increase in repolarising I_Kr acts during the DAD to keep the cell's membrane potential below the threshold for triggered activity; reduced I_Na also pays a role (although not as great as that of I_Kr) in the drug's mechanism of action as it reduces excitability of the cell, again reducing the ease with which the cell's membrane potential can reach threshold.

The vulnerability to the development of ectopic activity (i.e., propagation of triggered activity in all directions through the tissue) was assessed in a 2D homogeneous sheet to provide a medium for the synchronization of independent stochastic release events. Ten simulations were performed for each condition (WT and PQT variants ± MC-II-157c for the combinations of t_i and duration median distributions; Epicardial cell model) wherein the tissue model was paced to steady state and then left quiescent for two simulation seconds within which time the spontaneous release occurs.

The occurrence of ectopic activity exhibits a largely “all-or-nothing” response, where most conditions lead to either 0 or 100% of simulations resulting in a premature excitation (Table 1). Ectopic activity was promoted by narrow distributions of initiation time (i.e., tight synchronization) and short RyR waveforms (i.e., large spontaneous Ca²⁺ transients), and conversely inhibited by wide distributions of initiation time (i.e., lose synchronization) and long RyR waveforms (i.e., small spontaneous Ca²⁺ transients). For example, no ectopic activity was observed for any condition with t_i distribution widths of 1,000 ms or duration medians of 300 ms or longer, whereas a t_i distribution width of 350 ms combined with median durations of 150 or 200 ms resulted in ectopic activity occurring in 100% of simulations (Table 1).

PQT variants were more susceptible to the development of ectopic activity in tissue than WT, congruent with the single-cell results (See MC-II-157c Is Effective in Inhibiting DADs Turning into Triggered Activity). Similarly in-line with single-cell results, MC-II-157c can inhibit ectopic activity (Table 1): for example, it reduces or entirely inhibits the occurrence of premature excitation in four conditions: PQTa with t_i width of 550 ms and duration median 150 ms; PQTc with t_i width of 550 ms and duration median 200 ms; PQTc with t_i width of 550 ms and duration median 150 ms; and PQTc with t_i width of 350 ms and duration median 250 ms. The effect of I_Na block is also congruent with single cell results: it contributes to the inhibition of ectopic activity, but to a smaller extent than I_Kr modification.

Figure 8A shows an example of synchronized triggered activity initiating in a region of the 2D tissue (shown by the clustered peaks the images), before this triggered activity spreads throughout the tissue (shown by the yellow region) as an ectopic propagation. When the same situation is simulated with MC-II-157c (Figure 8B), triggered activity, and therefore ectopic propagation, is inhibited.

An S1-S2 pacing protocol was applied to the 1D strand model in order to compute the vulnerability window: S2 stimuli were applied across a range of time intervals, centered on one in every five cells of the 100 comprising the model, from the tenth to the 90th.

Examples of propagation applied during the repolarisation phase in 1D tissue simulations are shown in Figure 9A: if the triggered activity occurs early (at 311 ms in this example), the resultant excitation is surrounded by refractory tissue and the triggered activity dies out without propagating; If the triggered activity occurs slightly later (e.g., 316 ms) then refractory tissue is encountered at only one side of the triggered activity site and unidirectional block (or unidirectional propagation) occurs, in the retrograde direction (back toward the endocardium) in this example; If triggered activity occurs later than this (321 ms in this example) then all surrounding tissue has recovered and the triggered activity propagates in both directions along the strand (i.e., ectopic propagation, analogous to the situation shown in 2D tissue in Figure 8A). It is the unidirectional block situation (i.e., when the triggered activity occurs in the VW) that can lead to re-entrant arrhythmias if this situation occurred in 2D or 3D tissue. The VW identifies occurrences of trigger and substrate interaction that may lead to arrhythmias, and so quantifying the size VW is a convenient method to examine effects of disease conditions and drugs on trigger and substrate interaction. The baseline VWs for WT and the PQTa condition are mapped out in Figure 9B, as well as VWs in these two situations with simulated addition of MC-II-157c (I_Kr modification plus 50% I_Na block), as well as with only the MC-II-157c I_Kr modification. These VWs are quantified by length (over which they occur in the 1D strand) and area (length × temporal width). In both WT and PQTa conditions, the drug increases both the length and the area of the VW (for example, length increases from 9 to 16 mm in the PQTa condition, and area increases from 78 to 101 mm.ms). Again, the influence of I_Na is different depending on the condition: including the effects of I_Na (i.e., the full MC-II-157c simulations compared to the I_Kr only simulations) increases the length of tissue over which the VW occurs in all conditions (e.g., from 12 to 16 mm in PQTa), but reduces the overall area of the VW (e.g., from 126 to 101 mm.ms).

Our key findings are that: (i) Despite the hERG activator MC-II-157c reducing the maximal conductance of I_Kr by 12% compared to WT, the drug's modifications to I_Kr activation, inactivation and deactivation kinetics result in an overall increase in I_Kr during the action potential, and a concomitantly reduced APD under all PQT conditions in all transmural cell types; (ii) Although MC-II-157c acts on membrane ionic currents carrying K⁺ and Na⁺, the drug has indirect effects on intracellular Ca²⁺ handling, particularly SR Ca²⁺ loading and related spontaneous SR Ca²⁺ release events and subsequent DADs, through modulation of the AP; (iii) Increased I_Kr (through its repolarising effects) and reduced I_Na (by decreasing cell excitability) act to reduce the probability of a DAD reaching threshold and developing into triggered activity at both cellular and tissue scales; and (iv) Despite MC-II-157c reducing triggered activity at the cell level, the drug can increase both the spatial region of tissue over which a VW for unidirectional conduction block occurs, as well as the temporal width of the VW at all points along the tissue, and in doing so increases the total spatiotemporal size of the VW. These results highlight the complex considerations which underlie overall vulnerability to arrhythmia at multiple scales.

At the isolated cell level, MC-II-157c reduces SR Ca²⁺ loading, reduces the occurrences of spontaneous SR Ca²⁺ release events, reduces DAD occurrences, and reduces the probability of a DAD developing into a triggered action potential. At the tissue level, the combination of these factors leads to a significantly reduced vulnerability to the development of ectopic beats; MC-II-157c reduces the probability of an ectopic beat development at given spontaneous release function distributions (i.e., synchronization degree) as well as inhibiting SR Ca²⁺ loading and thus reducing synchronization (which, itself, reduces the probability of ectopic beats). In relation to the development of Ca²⁺ induced triggers at single cell and tissue levels, therefore, our simulation results suggest that the drug has efficacy as an anti-arrhythmic.

However, even though arrhythmic triggered activity is reduced with the drug, the spatiotemporal area of the VW increased, and so the anti-arrhythmic effects of the reduced probability of an arrhythmia trigger occurring is opposed by the pro-arrhythmic effects of an increased arrhythmia substrate. This highlights the importance of examining the effects of pharmacological compounds on both the triggers and the substrates that underlie arrhythmia initiation: there is a delicate balance between increased/decreased trigger and increased/decreased substrate (i.e., trigger-substrate interaction) that determines whether any given trigger stimulus (e.g., a spontaneous SR Ca²⁺ release event) will result in unidirectional propagation and, potentially, initiation of a re-entrant arrhythmia. This is one potential reason that single cell studies showing efficacious effects of putative anti-arrhythmic drugs may not translate to the clinic.

The effects of MC-II-157c on the VW occur through two mechanisms due to the drug's action on both I_Kr and I_Na, both of which modify transmural dispersion of repolarisation. Activation of I_Kr causes heterogeneous changes to APD at the tissue level, the same mechanism as in our previous studies examining the effects of NS1643 (Peitersen et al., 2008). Block of I_Na results in slowed transmural conduction (Kléber and Rudy, 2004) and therefore delayed activation of epicardial (but not endocardial) tissue: this in turn modifies transmural dispersion of repolarisation, even though the change to APD is minimal with I_Na block. Block of I_Na also reduces excitability of the tissue, necessitating a larger trigger to initiate propagation, which also contributes to the change in the VW.

The role of I_Na loss of function in arrhythmogenesis has been examined in detail previously (see Clancy et al., 2015 for a review); however, one intriguing finding from our cell and tissue simulations was the varying effects that different magnitudes of I_Na block had on trigger development and substrate size. At the cell level, the probability of triggered activity developing from DADs did not follow a simple monotonic change with increasing I_Na block in all cases. Take, for example, the 350 ms Ca²⁺ release duration distribution results shown in the top panels of Figure 7B: In WT and the PQTb condition, I_Na block (in addition to the I_Kr modification) reduces the probability of DADs developing into triggered activity (relative to the I_Kr modification alone), with more block reducing this probability; that is to say, the greater the I_Na block, the more anti-arrhythmic (in terms of reducing triggered activity) the effects. However, in the PQTa condition, while 40% block of I_Na reduced the probability of triggered activity occurring, increasing block of the current to 80% slightly increased the probability of triggered activity; in this condition, a small amount of I_Na block has anti-arrhythmic effects, but increasing this small level of block is pro-arrhythmic. In the PQTc condition, I_Na block of any magnitude increased the probability of triggered activity (i.e., I_Na block is pro-arrhythmic), with the probability of triggered activity occurring increasing as I_Na block is increased. Similarly varied results were found for the 300 ms Ca²⁺ release duration distribution (lower panels in Figure 7B), although the anti-/pro-arrhythmic effects did not necessarily match those seen with the 350 ms distribution.

One further note of caution with regards to I_Na block comes from our tissue-level VW results in Figure 9B. Quantification of the length of tissue over which the VW occurs shows that 50% I_Na block (compare the full MC-II-157c effects to the effects with the I_Kr modification alone) increases this length, potentially due to an increase in transmural dispersion of repolarisation with I_Na block; analyses of these results alone would conclude that I_Na block is pro-arrhythmic. However, despite the spatial width of the VW increasing, block of I_Na results in the total spatiotemporal area of the VW decreasing (again, compare the full MC-II-157c effects to the effects with the I_Kr modification alone), i.e., an anti-arrhythmic result, likely due to the reduced excitability that results from I_Na block reducing the likelihood that any triggered activity would propagate.

Our tissue-level findings therefore indicate that I_Na block per se is an effective antiarrhythmic strategy (as seen with class I antiarrhythmic drugs; Camm, 2012), but our cell-level findings highlight that the magnitude of I_Na determines cell (and by extension, tissue) electrophysiological consequences in a manner that is not intuitive. The mechanisms underlying these varying effects of I_Na block on arrhythmogenesis, particularly on DAD initiation and their development into triggered activity, remain to be elucidated.

Although this study did not focus on the mechanisms of arrhythmia triggers (in that we prescribed SR Ca²⁺ release events under certain conditions), one interesting finding did emerge in relation to arrhythmia triggers that can result from abnormal intracellular Ca²⁺ handling: Modification of membrane ion channels carrying ions other than Ca²⁺ (K⁺ and Na⁺ in this case) can have beneficial effects in terms of restoring abnormal intracellular Ca²⁺ handling, through their actions in shortening APD. This was shown in Figure 6, where the increased SR Ca²⁺ loading (and resultant spontaneous SR Ca²⁺ release events) seen under PQT conditions was reversed by upregulating I_Kr and downregulating I_Na, which in turn shortened APD. This reduces the duration over which I_Ca,L is activated, reduces the amount of Ca²⁺ crossing the cell membrane and entering the cell via that current, and thus reduces SR Ca²⁺ load. Furthermore, the combined action of both of these current modifications reduced the probability of DADs manifesting as triggered action potentials in single cell as well as triggered action potentials manifesting as fully propagating ectopic beats at the tissue scale. Targeting cell membrane ion channels carrying ions other than Ca²⁺ in order to restore abnormal intracellular Ca²⁺ handling (as seen in heart failure, for example; Lou et al., 2012) may be beneficial in cases where up/downregulation of Ca²⁺-specific drug targets (ryanodine receptors or SERCA, for example) will alter the delicate homeostasis of an already-compromised system, yielding negative results (Ratner, 2015).

One-dimensional models (other than 1D rings; e.g., Vinet and Roberge, 1994) cannot simulate re-entrant activity, and so it is necessary to use 2D and 3D models to examine how unidirectional propagation develops into re-entry. Although the quantitative characteristics of the VW examined in 1D tissues may change in 2D and 3D depending on the spatial locations of cell types (e.g., “base-apex” as well as transmural distributions), based on our previous work we would expect the qualitative characteristics to remain similar (Benson et al., 2007, 2008, 2011a). We show examples of 2D and 3D modeling in Figure 10: In Figure 10A, a trigger occurring in the VW (resulting in unidirectional conduction block) develops into re-entry in a simple 2D model; In Figure 10B, a trigger occurring outside the VW develops into ectopic propagation (i.e., not re-entrant) in a detailed 3D model of a human left ventricular wall slab, which manifests as significant differences in the body surface potential activation maps (Figure 10C). The advantages of using detailed 3D models (in this case, where the geometry is obtained from diffusion tensor MRI) lie in their ability to reproduce the orthotropic conduction velocities resulting from complicated tissue architecture (i.e., fiber and sheet structure; Benson et al., 2007; Smaill et al., 2013) and the boundary and curvature effects that can modulate electrotonic coupling and propagation (Walton et al., 2013; Campos et al., 2015): these effects are crucial in understanding the complex and chaotic propagation patterns underlying cardiac arrhythmias. Nevertheless, the 1D models used in this study allow us to examine arrhythmia trigger-substrate interaction in a simple and methodological manner.

In this study, a general model of intracellular Ca²⁺ handling was integrated with a simplified formulation of a human ventricular AP model. Due to the general nature of the Ca²⁺ handling model, details of heterogeneity in Ca²⁺ handling were not included in order to avoid introducing artifacts. For this reason, and due to the multi-scale focus of the study, detailed investigation of the mechanisms of spontaneous release in single cell, and their regional dependencies, was not performed, and investigation was instead limited to the potential effect of MC-II-157c on reversing SR Ca²⁺ loading. There are many possible mechanisms of spontaneous release associated with diseases linked to prolonged QT interval (e.g., hyper phosphorylation of the RyRs; upregulation of SERCA; detubulation) which were not considered in the present study. However, the aim of this study was to investigate the multi-scale interaction between trigger and substrate, and the use of the simplified spontaneous release functions in single-cell and tissue simulations allowed this to be analyzed in a general manner and across a large range of conditions, independent of spontaneous release mechanism. In future, combining detailed single cell studies of the mechanisms of spontaneous release in disease conditions with tissue simulations of the same conditions would provide further mechanistic insight; a dynamic simplified spontaneous release model would furthermore allow the study of long-term interactions between trigger and substrate in tissue e.g., during re-entry.

The factors underlying the synchronization and propagation of ectopic activity are highly complex and it is therefore worth making explicit that our simplified approach (which does not consider, for example, heterogeneity in the distributions describing Ca²⁺ release) is primarily suitable for interpretation of general trends only, rather than as a quantitative analysis of the efficacy of MC-II-157c in modulating Ca²⁺ release dependent triggers; it is encouraging to note that our results are consistent with those of a previous study (Campos et al., 2015), with the emergence of premature excitation from independent stochastic events overcoming electrotonic load, and the steep, “all-or-nothing” relationship observed in tissue. This steep relationship also likely accounts for the lack of effect of MC-II-157c on PQTb (the condition exhibiting the highest vulnerability to ectopic activity), wherein the distributions selected were not close enough to the threshold region to reveal an effect.

The relative contribution of I_Kr to repolarisation is different in different models of the human ventricular action potential (Mirams et al., 2014): for example, I_Kr contributes more repolarising current in the ORd model (O'Hara et al., 2011) than in the ten Tusscher & Panfilov model (ten Tusscher and Panfilov, 2006), despite both models being validated against experimental data. It is therefore possible that our results overestimate the effects of the hERG activator MC-II-157c, although the identified pro- and anti-arrhythmic mechanisms will still be relevant.

The limitations of using 1D and 2D simplifications of 3D cardiac tissue have been discussed in detail previously (Clayton et al., 2011). Here we only note that the 1D strand model of the left ventricular wall allows us to examine mechanisms underlying how arrhythmia triggers and substrates interact to initiate unidirectional propagation, without the added complicating effects that geometrical (shape) and architectural (fiber, sheet etc.) considerations would bring. Nevertheless, it is anticipated that these geometrical and architectural effects will play a role not only in the transition from unidirectional propagation to re-entry, but on the initiation of the unidirectional propagation itself through, for example, electrotonic effects (Benson et al., 2007; Walton et al., 2013). Similarly, our 2D and 3D tissues were isotropic, i.e., no fiber or sheet structure, and so any conclusions drawn from these simulations should be interpreted with this in mind. Elucidating the roles that tissue geometry and architecture play in arrhythmogenesis (Smaill et al., 2013) is an important next step in fully understanding arrhythmia trigger and substrate interaction.

The spatial distributions of endocardial, midmyocardial, and epicardial cell types across the human ventricular wall has still not been confirmed: some studies suggest that midmyocardial cells are found predominantly in isolated regions of the subendocardium (Glukhov et al., 2010), while others suggest a continuous population of midmyocardial cells in the subepicardial region (Drouin et al., 1995); these distributions may be dependent on species, location in the ventricular wall, and disease state (Antzelevitch, 2010; Strom et al., 2010). Because of this uncertainty, we set the spatial distribution of the three cell types to be equal in the transmural direction, but the effects of the electrotonic interactions of different regions of cell types are likely to be qualitatively similar if these distributions are altered.

The simulated model of MC-II-157c reproduced qualitatively the key features of the effect of the compound on I_Kr (i.e., an increased activity during a depolarising pulse), but it is important to note that there were significant differences between the simulation and experimental data (Figure 5): firstly, the formulation of I_Kr implemented does not have a time-dependent inactivation, which is observed in the experimental trace; secondly, the simulated data exhibited a larger difference in the magnitude of the current during the depolarisation step between baseline and MC-II-157c than observed experimentally. A more detailed model of I_Kr in both basal and MC-II-157c conditions would be essential for future and more detailed analysis of the compound specifically.

Furthermore, recent work, (e.g., Li et al., 2017), has highlighted that simple modulation of Hodgkin-Huxley current formulations, as used to simulate I_Kr in the ORd model, may not sufficiently capture complex drug-channel interaction dynamics, and that more complex Markov model formulations may be necessary to simulate such dynamics. However, until a full experimental characterisation of the dynamics of MC-II-157c effects on I_Kr under a range of conditions is carried out, allowing a validated Markov model of drug-channel interactions to be developed, we make use of the available data (Guo et al., 2014) to modify the Hodgkin-Huxley formulation used in the ORd model. The effects of MC-II-157c on I_Na are not as well characterized as its effects on I_Kr (Guo et al., 2014). We simulated the drug's action on I_Na by a simple reduction in the maximal conductance of the current. It is possible, however, that MC-II-157c also modifies the current's kinetics (i.e., shifts to the current's activation and inactivation curves, and/or a change to the time constants associated with these processes) in a similar manner to the way in which the kinetics of I_Kr are modified. Thus, further experimental characterisation of the drug's effects on both I_Kr and I_Na are required.

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 KF and handling Editor declared their shared affiliation.