The equations for I_CaL in the TP06 model (Ten Tusscher and Panfilov, 2006a) were modified to incorporate experimental data on CACNA1C mutation-induced changes. First, we determined the modifications to the original I_CaL model to reproduce the behavior of the mutant I_CaL (Figure 1A) during the same voltage-clamp employed in experiments (Fukuyama et al., 2014). Theoretical steady-state activation and inactivation curves which were used to simulate G1783C, wild-type (WT), P381S, M456I, A582D, and R858H I_CaL, are shown in Figure 1C. Second, based on experimental I-V relationships (Figure 1B), mathematical models of I_CaL were constructed (formulations are listed in Supplementary Material). This was achieved by simulating the experimental voltage-clamp protocol (Fukuyama et al., 2014) and scaling relative current proportions for G1783C, WT, P381S, M456I, A582D, and R858H conditions. Peak inward G1783C, P381S, M456I, A582D and R858H I_CaL density was, respectively, ~0.86-, ~1.04-, ~1.08-, ~1.19-, ~1.54-fold that for WT I_CaL (Figure 1B). The simulated I-V curves (Figure 1D) matched closely with the experimental observation (Figure 1B).

Upon biophysical analysis of the experimental data on CACNA1C mutations (Fukuyama et al., 2014), four major changes to I_CaL were considered. These changes included: I_CaL current densities, steady-state activation curves, voltage-dependent inactivation curves as well as the time constant for the voltage-dependent inactivation. Parameters for I_CaL, including the scaling factor of the I_CaL conductance (CSF), the midpoint voltage of the voltage-activation curve (V_{a, 0.5}), the slope of the voltage-activation curve (S_a), the midpoint voltage of the voltage-inactivation curve (V_{ina, 0.5}), the slope of the voltage-inactivation curve (S_ina) and the scaling factor of voltage-inactivation time constant (TCSF), were modified to reproduce the experimental I–V relationships (Comparison between simulation and experimental results for I_CaL can be found in Table S1).

The TP06 model, based on human experimental data, was developed to reproduce transmural heterogeneity of electrical properties (Ten Tusscher et al., 2004; Ten Tusscher and Panfilov, 2006a) by changing maximum conductivities of the transient outward potassium channel current (I_to) and slow delayed rectifier potassium channel current (I_Ks). The I_CaL formulations were the same for endocardial (ENDO), midmyocardial (MCELL) and epicardial (EPI) cell models. The model (Ten Tusscher et al., 2006) was suggested to be suitable for simulating wave dynamics at the tissue and organ levels (Ten Tusscher and Panfilov, 2006b; Ten Tusscher et al., 2009). In 2013, the calcium-induced-calcium release flux (I_rel) was modeled as the combination of both SR calcium release and SR calcium leak by Lascano et al. (2013). The same modifications were employed in our previous studies (Bai et al., 2016b, 2017; Liu et al., 2016; formulations are listed in Supplementary Material).

In single cell simulations, APs were elicited by pre-pacing the models for 100 cycles to reach a stable steady state. APD was computed as AP duration at 90% repolarization (APD₉₀). Changes in I_CaL, sodium-calcium exchanger current (I_NCX), I_rel, cytoplasmic calcium concentration ([Ca²⁺]_i) and SR calcium concentration ([Ca²⁺]_SR) were used to analyze CACNA1C mutations-induced calcium handling. Differences of electrical properties between endocardial- and epicardial- cells (ENDO-EPI), between endocardial- and midmyocardial- cells (ENDO-M) and between epicardial- and midmyocardial- cells (EPI-M) may contribute to transmural electrical heterogeneities in tissues (Zhang et al., 2008). Transmural electrical heterogeneity caused by CACNA1C mutations was assessed by quantifying [Ca²⁺]_i amplitude ([Ca²⁺]_i(m)), SR calcium content ([Ca²⁺]_SR(m)), I_rel amplitude (I_rel(m)), and APD. The rate dependence of APD was also investigated by using a pacing cycle length (PCL) of 2000, 1000, and 500 ms, respectively.

The maximum slope of the APD restitution (APDR) curve is a main marker for determining whether alternans and spiral breakup will occur (Ten Tusscher and Panfilov, 2006a). We used both the standard S1-S2 and dynamic protocols to determine APDR. For the S1-S2 protocol, 30 S1 stimuli were applied at a PCL of 1,000 ms, and the S2 stimulus was applied at varying diastolic intervals (DI) after the AP evoked by the last S1 stimulus. For the dynamic protocol, a series of 30 stimuli were applied at a PCL of 1,000 ms, after which the PCL was decreased. APDR curves were obtained by plotting APD against DI.

Figure 2 shows multicellular tissue models which consist of a 1D transmural ventricular fiber (Figure 2A), a 1D MCELL homogeneous cable (Figure 2B), a 2D transmural ventricular sheet (Figure 2C), and a 2D MCELL tissue (Figure 2D).

The 1D transmural ventricular fiber of length 24.75 mm (Figure 2A), which contained a 9 mm long ENDO region, a 6.75 mm long MCELL region and a 9 mm long EPI region, was constructed to compute a pseudo-ECG and investigate the tissue vulnerability to unidirectional conduction block. The fiber model consisted of 60 ENDO nodes, 45 MCELL nodes and 60 EPI nodes (O'Hara et al., 2011). Following the application of a sequence of 10 conditioning S1 pulses (with amplitude of −40 μA/cm² for 3 ms) applied to a 0.45 mm long ENDO segment at a PCL of 1,000 ms, repolarization characteristics were evaluated by computing repolarization time (RT), DOR, and maximum spatial gradient of APD (MSG). RT was computed as the largest APD in the transmural cable. DOR was measured as the difference between the largest and smallest APD of cells in the fiber. Spatial gradient of APD (SG) along the transmural cable was calculated as changes of APD per millimeter. In addition, the pseudo-ECG was obtained by using the method of Gima and Rudy (Gima and Rudy, 2002; formulations are listed in Supplementary Material). QT interval, T-wave width and T-wave amplitude were quantified to examine the changes in the ECG due to CACNA1C mutations. QT interval was estimated as the time interval between the Q-wave onset and the T-wave end, T-wave width was computed as the peak and the end of the T-wave, and T-wave amplitude was defined as the peak voltage of the T-wave. The T-wave end was determined by the intersection point of the baseline (y = 0 mV) and the T wave. The inducibility of unidirectional conduction block for each mutant CACNA1C was also quantified by computing the maximum PCL (MPCL) that produced 2:1 block.

A 1D MCELL cable of length 24.75 mm (Figure 2B) was used to measure conduction velocity (CV) by calculating the time Δt for the wavefront to propagate from x − Δx to x + Δx, defining CV = 2 Δx/Δt. To examine whether focal activity in the 1D cable was induced at a PCL of 500 ms, solitary planar waves were initiated by applying S1 pulses (with the same size, strength and duration as the one used for the 1D transmural simulation).

A transmural ventricular sheet (Figure 2C), which was constructed by expanding the 1D transmural fiber into a sheet with a length of 24.75 mm and a width of 150 mm, was developed to examine tissue vulnerability to the initiation of reentrant waves. A S1-S1 stimulation (with the same strength and duration as the one used for the 1D simulations) was applied to a region of 0.45 × 75 mm², located at the center of the left ENDO boundary.

A 2D MCELL tissue model of 375 × 375 mm² (Figure 2D) was developed to investigate effects of CACNA1C mutations on the spatiotemporal behavior of spiral waves. Spiral waves were induced by a standard S1–S2 stimulation. A plane wave was initiated by applying the S1 stimulus to the left side of the domain (0.75 × 375 mm²). Once the plane wave had passed over the first half of the domain, the S2 stimulus was applied to the first quarter of the domain so that a spiral wave was produced. AP was recorded from the representative point (x = 187.5 mm, y = 187.5 mm) and the fundamental frequency was obtained from the power spectra of the AP.

For the 3D model, simulations were performed using an anatomical human ventricular geometry developed in our previous studies (Bai et al., 2015, 2016c). It has a spatial resolution of 0.5 mm with ~6.38 million cells in total. For both left and right ventricles, the tissue was segmented into distinctive ENDO, MCELL, and EPI layers with similar contiguous figurations in the transmural wall as in the 1D transmural ventricular fiber model. The S1-S1 stimulation was applied to a 2.5 mm wide region of the ENDO layer. For the purposes of the present study, the 3D anatomical model was assumed to be electrically homogeneous.

The cell models were incorporated into a parabolic partial differential equation (PDE) to construct mono-domain models of cardiac electrophysiology to describe the reaction-diffusion system in simulating cardiac dynamics (Clayton and Panfilov, 2008). The governing equation is

where C_m = 1 μF/cm² is the capacitance, D is the effective diffusion constants, and I_ion is the total transmembrane current.

We used a forward-Euler method for marching, with a time step (Δt) of 0.02 ms and a space step (Δx = Δy = Δz), to solve the PDEs (Equation 1). The value of D is set to be 0.0385 mm²/ms for simulating excitation waves in the 1D transmural ventricular fiber, the 1D homogeneous cable, the 2D transmural ventricular sheet and the 3D ventricular model. The spatial resolution in 1D, 2D and 3D models are as follows: 1D ventricular model: Δx = 0.15 mm; 2D transmural ventricular model: Δx = Δy = 0.15 mm; 3D ventricular model: Δx = Δy = Δz = 0.5 mm. For the 2D MCELL ventricular model, D is set to be 0.154 mm²/ms and the spatial resolution is chosen to be 0.25 mm to investigate the role of APDR in the occurrence of electrical instability (Ten Tusscher et al., 2004; Vandersickel et al., 2014; Nayak and Pandit, 2015; Zimik et al., 2015; Nayak et al., 2017). We also used Neumann (i.e., no-flux) boundary conditions (Clayton and Panfilov, 2008). Simulations were carried out on a 64-bit Intel core i7-3930K CPU system with 64 GB memory. Efficient parallelization was implemented using GPU acceleration (Bai et al., 2015). Although different time, space and diffusion coefficient were used, solutions to the mono-domain model with isotropic diffusion fulfilled the stability criterion: (i.e., DΔt/ Δx² < 1/2; Clayton and Panfilov, 2008).

I_CaL arising from the Ca_V1.2 mutations altered intracellular calcium handling and prolonged AP as shown in Figure 3A. For ENDO cells, the I_CaL amplitude was increased progressively from 4.12 pA/pF in the G1783C condition (Magenta) to 4.28 pA/pF (WT, Black), 4.35 pA/pF (P381S, Red), 4.44 pA/pF (M456I, Green), 4.71 pA/pF (A582D, Blue), and 5.0 pA/pF (R858H, Cyan), respectively (Figure 3B). Altered intracellular calcium handling resulted from increased I_CaL was shown by the time courses of I_NCX(Figure 3C), I_rel (Figure 3D), [Ca²⁺]_i (Figure 3E) and [Ca²⁺]_SR(Figure 3F). Changes in [Ca²⁺]_i(m), [Ca²⁺]_SR(m), and I_rel(m) were related to an increase in I_CaL (changes are summarized in Table S2). Increased I_CaL during the AP plateau triggered a large I_rel, leading to cytoplasmic calcium overload and APD prolongation. The measured APD was 264.4 ms (G1783C), 267.8 ms (WT), 271.8 ms (P381S), 277 ms (M456I), 282.6 ms (A582D), and 294 ms (R858H), respectively.

Simulations of transmural electrical heterogeneity were performed for each mutation, producing APs of ENDO, MCELL, and EPI cells. [Ca²⁺]_i(m), [Ca²⁺]_SR(m), I_rel(m), and APD were used to quantify electrical properties of the different cell types (listed in Table S2). The electrical differences associated with each mutant CACNA1C for ENDO-EPI, ENDO-M as well as EPI-M are summarized in Table S3. There were small electrical differences for ENDO-EPI, whilst large electrical heterogeneities were found for EPI-M and ENDO-M. For one “mutant” I_CaL tested (R858H), [Ca²⁺]_SR(m) of ENDO-EPI (~0.1 mM) was smaller than that of EPI-M (~0.64 mM) as well as ENDO-M (~0.74 mM). Among in these mutations, the Ca_V1.2 R858H mutation with the largest effect on electrophysiological heterogeneity was present. For the [Ca²⁺]_SR(m), the ENDO-M difference for the R858H condition was larger (~0.74 mM) than the one (~0.64 mM) produced by the presence of the WT I_CaL, while other mutations caused changes ranging from ~0.59 mM to ~0.69 mM. Other characteristics of transmural heterogeneity are listed in Table S3, which demonstrate that the largest electrophysiological heterogeneity is with the R858H condition.

The effects of I_CaL associated with CACNA1C mutations on AP shape are shown in Figure 4. Cell simulations were conducted by increasing the pacing frequency from 0.5 to 2 Hz [corresponding to PCL of 2,000 ms (Green), 1,000 ms (Red) as well as 500 ms (Black)] and the measured APDs for G1783C, WT, P381S, M456I, A582D, and R858H cells are listed in Table S4. For the ENDO (Figure 4, left column), MCELL (Figure 4, middle column) and EPI (Figure 4, right column) cells, APD₉₀ was abbreviated with PCL and no afterdepolarizations were triggered under the G1783C, WT, P381S, M456I, and A582D conditions. However, under the R858H condition, a DAD in the MCELL cells was induced, but no ectopic beats in the EPI and ENDO cells were observed, when the PCL was 500 ms.

To illustrate afterdepolarizations generating events, the time courses of AP, I_CaL, I_NCX, I_rel, [Ca²⁺]_i, and [Ca²⁺]_SR for the ENDO (Figure 5, left column), MCELL (Figure 5, middle column) and EPI (Figure 5, right column) cells at PCL of 500 ms (Red) and 1,000 ms (Black) are shown. As can be seen, I_CaL arising from the R858H mutation contributed to an increase in [Ca²⁺]_SR, which, consequently, enhanced spontaneous I_rel, accompanied by an inward I_NCX that depolarized the cell. Subsequently, EADs (marked with a blue circle) and DADs (marked with red circles) were triggered in R858H-MCELL cells. Compared with R858H-ENDO and R858H-EPI cells, spontaneous I_rel frequency was higher in R858H-MCELL cells. Additionally, spontaneous I_rel occurred during repolarization and therefore led to the EADs, however, DADs were induced by spontaneous I_rel after repolarization. In addition, the APD (≈421.2 ms) of the DAD-induced AP was larger than that (348.2 ms) of the R858H-MCELL cells at a PCL of 500 ms, indicating QT interval prolongation at rapid heart rates.

The effects of the mutant I_CaL on ventricular APDR are shown in Figure 6. The APD reduction was rate-dependent for ENDO (Figures 6A,D), EPI (Figures 6B,E) and MCELL (Figures 6C,F) cells. Across the range of DIs tested, the measured APD was larger for the R858H (Cyan) condition than for the G1783C (Magenta), WT (Black), P381S (Red), M456I (Green), and A582D (Blue) conditions. For the R858H settings, the maximum restitution slope of 2 for MCELL cells (Figure 6E) was larger than those for ENDO (Figure 6D, 1.35) as well as EPI cells (Figure 6F, 1.5). For MCELL cells, the computed APDR slopes of MCELL cells (Figure 6E) were 2 (G1783C), 1.8 (WT), 1.7 (P381S), 1.7 (M456I), 1.7 (A582D), and 2 (R858H), respectively. Details of maximum APDR slopes for three cell types can be found in Table S5.

To examine the manifestation of the temporal and spatial dispersion of AP in the pseudo-ECG, AP propagation (Figure 7A), spatial distribution of APD (Figure 7B), spatial gradient of APD (Figure 7C) and the computed pseudo-ECG (Figure 7D) of 1D transmural strand at a PCL of 1,000 ms for each “mutant” I_CaL are shown. Repolarization time (Figure 7E), DOR (Figure 7F), the maximum spatial gradient of APD at the EPI-M junction (Figure 7G), QT interval (Figure 7H), and T wave width (Figure 7I) were computed for the G1783C, WT, P381S, M456I, A582D, and R858H conditions. As can be seen in Figure 7D, the QT interval was prolonged from 397.1 ms (G1783C) to 398.5 ms (WT), 401.2 ms (P381S), 408.2 ms (M456I), 415.9 ms (A582D), and 425.7 ms (R858H), respectively (Figure 7H). The QT interval prolongation underlying increased APD in single cells has been linked to repolarization time (RT; Gima and Rudy, 2002). Therefore, the effects of the mutant I_CaL on RT were examined. The computed RT was 360.48 ms (G1783C), 363.6 ms (WT), 367.34 ms (P381S), 374.84 ms (M456I), 381.14 ms (A582D), and 393.32 ms (R858H), respectively (Figure 7E). The QT prolongation was consistent with the prediction of RT. Moreover, the APD of MCELL cells with the longest repolarization time was increased from 359 ms (G1783C) to 364.2 ms (WT), 370.4 ms (P381S), 379.6 ms (M456I), 389.4 ms (A582D), and 410.8 ms (R858H), respectively (listed in Table S6). The QT interval prolongation can be attributed to I_CaL due to CACNA1C mutations that influence APD in single cells and ventricular RT in ventricular tissues.

In addition, changes in the T wave were examined. As is shown in Figure 7H, the T-wave width was 60 ms (G1783C), 59.9 ms (WT), 59.9 ms (P381S), 61.7 ms (M456I), 63.9 ms (A582D), and 64.5 ms (R858H), respectively (Figure 7I). The augmented T-wave width has been attributed to an increase in DOR (Antzelevitch et al., 1998), therefore the effects of the mutant I_CaL on APD dispersion in the transmural strand were also investigated. The computed DOR was enhanced from 60.94 ms (G1783C) to 61.3 ms (WT), 62.2 ms (P381S), 64.12 ms (M456I), 65.64 ms (A582D), and 67.96 ms (R858H; Figure 7F). In addition, it has been suggested that DOR is a marker of electrical heterogeneity in APD (Antzelevitch, 2001). Comparatively, the maximal APD difference for EPI-M was 93 ms (G1783C), 94.8 ms (WT), 97 ms (P381S), 101 ms (M456I), 104.8 ms (A582D), and 115.2 ms (R858H), respectively (summarized in Table S3). The APD heterogeneity induced by CACNA1C mutations may contribute to DOR and thereby to T-wave width.

Similarly, T-wave amplitudes were calculated. The T-wave amplitude was 1.051 mV (G1783C), 1.053 mV (WT), 1.056 mV (P381S), 1.075mV (M456I), 1.087 mV (A582D), and 1.104 mV (R858H), respectively (listed in Table S6). Compared with the WT condition, there was no evident difference in T-wave amplitude. Changes of T-wave amplitude can be attributed to altered temporal and spatial gradients in membrane potential (Gima and Rudy, 2002; Zhang et al., 2008). Therefore, the APD gradient (Figure 7C) and the membrane potential gradient (δ, Figure 8) were computed. Figure 8A shows simulated ENDO (Black), EPI (Red), MCELL (Green) APs for each mutant condition whilst Figure 8B shows corresponding time-course plots of δ for EPI-ENDO (Light Gray), EPI-M (Gray), and ENDO-M (Dark Gray). The maximal δs for EPI-M (Figure 8D) and ENDO-M (Figure 8E) were greater than that for EPI-ENDO (Figure 8C). This was consistent with the prediction of electrical heterogeneities of AP in single cells (summarized in Table S3). There was no significant change in the maximum EPI-ENDO δ between these mutations (Figure 8F). However, the maximum EPI-M δ was 73.8 mV (G1783C), 74.2 mV (WT), 75.2 mV (P381S), 75.6 mV (M456I), 77.1 mV (A582D), and 79.4 mV (R858H), respectively (Figure 8G), and the maximum ENDO-M δ was 91.59 mV (G1783C), 92.06 mV (WT), 93.06 mV (P381S), 93.14 mV (M456I), 94.69 mV (A582D), and 96.27 mV (R858H), respectively (Figure 8H). Also, the maximum spatial gradient (MSG) of APD at the EPI-M junction was 18 ms/mm (G1783C), 18.1 ms/mm (WT), 18.5 ms/mm (P381S), 18.9 ms/mm (M456I), 19.2 ms/mm (A582D), and 19.5 ms/mm (R858H), respectively (Figure 7G). Changes in MSG at the EPI-M junction and maximal EPI-M δ were consistent with the altered T-wave amplitude.

Taken together, the ECG phenotypes were dependent on changes of AP caused by altered I_CaL due to CACNA1C mutations in single cells and the longest QT interval was obtained under the R858H condition (listed in Table S6).

To investigate the cellular level conditions required for afterdepolarizations to trigger a premature ventricular complex (PVC) at the multicellular tissue level, a 24.75-mm long MCELL strand model was constructed and space-time plots of AP propagation in the 1D homogeneous cable at the PCL of 500 ms are shown in Figure 9. For the G1783C, WT, P381S, M456I and A582D conditions, no focal activity was observed whilst a R858H-mediated focal activity (marked with a white rectangle) was triggered at Time = 35,000 ms. Although a DAD mediated-AP was induced in MCELL cells (Figure 5, middle column), such simulations may not necessarily reflect the situation for intact tissue, in which electrical coupling occurs between cells and may smooth out these electrical differences. Thus, the propagation of the DAD mediated-AP was not observed (Figure S1) and thereby no PVC occurred under the R858H condition.

To identify the electrophysiological substrates that promote arrhythmogenesis, tissue vulnerability to unidirectional conduction block necessary to the genesis of reentry was investigated, and space-time plots of AP propagation in the 1D transmural cable at different PCLs are shown in Figure 10. As can be seen, the tissue susceptibility to unidirectional conduction block showed a RT dependency and increased with a decrease in PCL. For instance, the MPCL (listed in Table S6) increased from 350 ms (G1783C) to 352 ms (WT), 354 ms (P381S), 356 ms (M456I), 365 ms (A582D), and 370 ms (R858H), respectively, with RT prolongation from 360.48 ms (G1783C) to 363.6 ms (WT), 367.34 ms (P381S), 374.84 ms (M456I), 381.14 ms (A582D), and 393.32 ms (R858H), respectively. As for R858H conditions, no unidirectional conduction block was observed at PCL = 1,000 ms, while unidirectional conduction block occurred at PCL = 370 ms. Bidirectional conduction was obtained in G1783C, WT, P381S, M456I, and A582D conditions at PCL = 370 ms. Therefore, R858H tissue was more vulnerable to unidirectional conduction block than other type tissues.

To examine if mutations-induced changes in APDR curves promote the breakup of a spiral wave, an idealized homogeneous tissue model was constructed and spiral waves were induced by S1-S2 stimulation. Figure 11 shows snapshots of spiral waves (Time = 5,000 ms, left column), time series of action potentials (middle column) from the point indicated by an asterisk (left column) and power spectra (right column) for the G1783C, WT, P381S, M456I, A582D, and R858H conditions. Reentrant waves were stable and persistent. The fundamental frequencies obtained from these power spectra of action potentials were 2.8 Hz. These ordered plots of action potential provided evidence for the rotating spiral state.

To determine if R858H-induced changes were necessary to induce ventricular arrhythmias, an idealized transmural tissue model was developed to investigate the initiation of reentry with a S1-S1 stimulation. As is shown in Figure 12, when the PCL was 370 ms, the first S1 stimulus (t = 0 ms) produced a wave to propagate from the ENDO layer to the EPI layer. The second S1 stimulus (t = 370 ms) produced a conditioning wave and no unidirectional conduction block occurred under G1783C, WT, P381S, M456I, and A582D conditions (t = 750 ms). A unidirectional conduction block was initiated under the R858H condition, leading to the genesis of spiral waves (t = 750 ms) and fibrillation-like activity. The reentry wave for the R858H conditions was unstable and promoted self-termination when it collided with its own prolonged repolarization tail or tissue borders (t = 1,000 ms). These results support the previous notion that the MPCL of 2:1 block for the R858H mutation is 370 ms (shown in Figure 10) and unidirectional conduction block is responsible for the genesis of reentry. These 2D results concur with the 1D simulations data, further illustrating the pro-arrhythmic effects of the R858H mutation. In addition, reentry can be induced under other conditions (Figure S2) if the S1-S1 interval was decreased. For instance, the S1-S1 interval for initiating spiral waves was 365 ms (R858H), 340 ms (A582D), 335 ms (M456I), 320 ms (P381S), 315 ms (WT), and 310 ms (G1783C), respectively.

Further simulations were performed in a 2D slice model. Figure 13 shows that the R858H I_CaL induced reentrant spiral waves, which lead to sustained multiple reentrant wavelets in a 2D ventricular slice. Snapshots of excitation waves at different time points (t = 10, 350, 440, 520, and 990 ms) are shown for the G1783C, WT, P381S, M456I, A582D, and R858H conditions. As can be seen, a S1 stimulus was applied to the four pacing sites on the ENDO layer (marked with white asterisks, t = 10 ms) and excitation waves were initiated. At t = 350 ms, another S1 stimulus was used to induce spiral waves. At t = 440 ms, excitation waves propagated to the whole sheet under the G1783C and WT conditions, local conduction block was observed under the P381S, M456I, and A582D conditions, and unidirectional conduction block occurred under the R858H condition. Several spiral waves were produced at t = 520 ms and sustained reentry was observed at t = 990 ms. It was found that spiral waves can be induced under other conditions by decreasing S1-S1 interval. The S1-S1 interval for initiating spiral waves was 310 ms (G1783C), 315 ms (WT), 320 ms (P381S), 325 ms (M456I), 330 ms (A582D), and 345 ms (R858H), respectively (Figure S3). Detailed movies of spiral wave initiation can be found in Videos S1–S6.

To examine if R858H-induced changes promote ventricular arrhythmias in human ventricles, further simulations were performed using a 3D human heart geometry. As is shown in Figure 14, the first S1 stimulus (t = 10 ms) produced a wave to propagate from the ENDO layer to the EPI layer. Ventricular repolarization (marked with a black circle or a blue circle) under the WT condition occurred much earlier than that for the R858H settings (t = 340 ms). At t = 350 ms, another excitation wave initiated by the second S1 stimulus conducted bidirectionally (marked with a bidirectional arrow) under the WT condition whilst the wave was locally blocked (marked with unidirectional arrows) by unrecovered tissues (marked with a red circle). Therefore, spiral waves (marked with unidirectional arrows) were induced under the R858H condition (t = 700 ms and t = 800 ms), and re-entrant waves were persistent throughout the simulation (t = 990 ms). Moreover, when the S1-S1 interval was gradually decreased, spiral waves can also be induced (Figure S4). In details, the S1-S1 interval for initiating spiral waves was 250 ms (G1783C), 300 ms (WT), 309 ms (P381S), 340 ms (M456I), 345 ms (A582D), and 348 ms (R858H), respectively. Detailed movies of the spiral wave initiation can be found in Videos S7–S12.

Taken together, these data demonstrated that the R858H mutation facilitates initiation of reentrant excitation waves and suggested that the R858H tissue was more vulnerable to the initiation of reentrant excitation waves than other type tissues. A summary of the effects of CACNA1C mutated I_CaL on simulated human ventricular electrical activity is listed in Table 1.

To our knowledge, this is the first study to investigate mechanisms underlying the genesis of VF in patients with the CACNA1C R858H mutation. The study found five major findings, (i) due to the R858H mutation I_CaL caused intracellular calcium overload, resulting in afterdepolarizations in single cells and focal activity in 1D fiber tissues; (ii) The R858H mutation induced an increase in I_CaL prolonged APD and augmented RT, leading to QT interval prolongation; (iii) The R858H mutation-induced electrical differences between cells augmented electrical heterogeneity, causing repolarization dispersion and thereby increasing T-wave width; (iv) Although, the R858H mutation steepened the APDR relationships in single cells, a stable spiral wave remained in homogeneous tissues; (v) Changes in cellular electrophysiology modulated wave conduction at tissue level facilitating unidirectional conduction block and thereby increasing tissue susceptibility to VF genesis in transmural ventricular tissues.

These simulation data in this study constitute novel evidence that the pro-arrhythmic effects of I_CaL associated with Ca_V1.2 R858H mutation involve both increased cell susceptibility to afterdepolarizations and tissue vulnerability to the reentry. The effects of the R858H mutation were investigated at cellular, 1D strand and 2D tissue and 3D organ levels, showing not only alterations in calcium handling and QT prolongation, but also repolarization dispersion and reentry.

In this study, the electrophysiological consequences of CACNA1C mutations were investigated by modeling the I_CaL, as in previous studies (Faber et al., 2007; Zhu and Clancy, 2007; Sung et al., 2010; Morotti et al., 2012; Boczek et al., 2015a; Bai et al., 2016c). Based on experimental data on I_CaL current densities, steady-state activation curves, voltage-dependent inactivation curves as well as the time constant for the voltage-dependent inactivation (Fukuyama et al., 2014), the I_CaL model was developed by changing the I_CaL conductance, the half-activation voltage, the slope of the activation curve, the half-inactivation voltage, the slope of the inactivation curve and the scaling factor of inactivation time constant. The “mutant” I_CaL models successfully reproduced the I-V relationships obtained from experimental studies (Fukuyama et al., 2014). These I_CaL models were incorporated into a human ventricular cell model, allowing us to relate changes in I_CaL to AP at the cellular level, in agreement with previous studies (Splawski et al., 2005; Sung et al., 2010; Wemhöner et al., 2015; Bai et al., 2016c). According to transmural ventricular wedge preparation models (Gima and Rudy, 2002; O'Hara et al., 2011), the 1D transmural strand model, integrating human ENDO-, MCELL- and EPI- cells, was constructed to compute unipolar pseudo-ECGs. In the cable, the proportion of each region composed of each distinct cell type is consistent with that used in other studies (O'Hara et al., 2011). In agreement with clinical findings, R858H-induced I_CaL led to a QT interval of 425.7 ms, which is within the QT range (420–476 ms) of R858H-porbands (Fukuyama et al., 2014). These models could be considered as a first step toward the validation of electrophysiological model. Thus, cardiac models provide a powerful tool for the study of mechanisms underlying ventricular arrhythmias caused by CACNA1C mutation effects.

Previous studies have demonstrated that distinct mutations can have variable effects on current morphology and lead to varying degrees of electrophysiological consequences, depending on kinetic changes induced by the CACNA1C mutation (Sung et al., 2010; Wemhöner et al., 2015; Sutphin et al., 2016). The “mutant” I_CaL represents state-specific kinetic properties of ion channels, allowing us to relate functional changes in I_CaL to AP at the cellular level, the electrocardiogram characteristics at the fiber tissue level and the spatiotemporal behavior of excitation waves at the sheet tissue level. In the present study, although G1783C, P381S, M456I, A582D, and R858H mutations changed Ca_V1.2 in activation and inactivation, VF was only observed in R858H patients (Fukuyama et al., 2013, 2014). Therefore, we have focused on the R858H mutation marked by an increase in the I_CaL and have investigated mechanisms underlying the genesis of reentry.

These analyses of R858H revealed an unexpected, distinct electrophysiological phenotype from the classical TS mutations (G406R/G402S; Splawski et al., 2004, 2005). Electrophysiological studies of G406R in exon 8A and G406R/G402S in exon8 showed almost complete loss of inactivation of Ca_V1.2 (Splawski et al., 2004, 2005). In contrast, R858H showed a significant gain of current density, ~2 mV negative shift of activation and ~2 mV positive shift of inactivation, resulting in an increased window current (Fukuyama et al., 2014). Due to the different electrophysiological phenotypes between the canonical G406R/G402S mutations and R858H, we performed modeling studies to better understand how the increased window current may affect the ventricular AP. Like the modeling studies of G406R/G402S (Splawski et al., 2004; Zhu and Clancy, 2007), our results showed prolongation of APD in TP06 and O'Hara-Rudy dynamic (ORd) human ventricular cell models (Figure S5). At fast pacing rates, R858H-induced SR calcium loading, which may lead to spontaneous afterdepolarizations, was also observed (Figure S6). The result of R858H-induced afterdepolarizations is consistent with other reports that pro-arrhythmic actions of CACNA1C mutations are due to SR calcium leak, thereby enabling an I_NCX that triggers afterdepolarizations (Splawski et al., 2005). In agreement with experimental findings (Priori and Corr, 1990), the R858H induced EADs and DADs share the same mechanism, where NCX-mediated inward current is involved.

These afterdepolarizations may mainly result from an increase I_CaL arising from the R858H mutation, because the VF in the R858H patients was completely prevented by atenolol, a beta-blocker (Fukuyama et al., 2014). Thus, our results support the concept that the R858H-induced defect is sufficient to cause significant action potential prolongation and afterdepolarizations in ventricular myocytes. Our data highlight how small changes in cellular calcium entry through Ca_V1.2 can lead to far-reaching changes in action potentials and intracellular calcium handling.

Computational modeling has shed valuable light on the cellular mechanisms underlying ventricular arrhythmia. Splawski et al. showed that a G406R mutation prolonged action potentials, altered intracellular calcium handling, and triggered DADs (Splawski et al., 2005). Sung et al. showed that G406R-induced EADs and DADs shared the same mechanism, with NCX-mediated inward current (Sung et al., 2010). However, how the G406R mutation is linked to ventricular arrhythmias is not yet fully understood. Thiel et al. showed that a G406R mutation altered intracellular calcium handling and calcium signaling, which could in turn contribute to DADs due to calcium signaling-induced SR calcium overload (Thiel et al., 2008). Recently, Dick et al. suggested that TS mutations (G406R/ G402S) disrupted calcium-dependent inactivation, which played an important role in APD prolongation and development of EADs (Dick et al., 2016). Different from the mechanism underlying G406R-induced arrhythmias, an I1166T mutation led to an EAD due to I_CaL reactivation (Boczek et al., 2015a). Other CACNA1C mutations also caused QT prolongation and cardiac arrhythmias (Gillis et al., 2012; Boczek et al., 2013, 2015b; Fukuyama et al., 2013, 2014; Hennessey et al., 2014; Wemhöner et al., 2015; Landstrom et al., 2016; Sutphin et al., 2016), but their cellular mechanisms underlying cardiac arrhythmias remain relatively unexplored, and the causal link between the CACNA1C mutations and ventricular arrhythmias is not addressed directly. In this study, the cellular mechanism that I_CaL due to the R858H mutation prolongs APD and increases SR calcium load, resulting in spontaneous I_rel and afterdepolarizations, is proposed. These simulation data of the present study add to the growing weight of evidence implicating cellular mechanisms of CACNA1C mutations.

The LQT8 is associated with malignant ventricular arrhythmias and a patient with a R858H mutation experienced an episode of VF (Fukuyama et al., 2013, 2014). Our data showed that R858H-induced electrical heterogeneity increased the transmural DOR and favored unidirectional conduction block necessary to the development of VF (Zhu and Clancy, 2007; Sung et al., 2010). The simulations may provide evidence for the pro-arrhythmic effects of the R858H mutation in facilitating reentrant excitation waves. They support the notion that DOR is a marker of tissue electrical heterogeneity and a substrate for reentrant arrhythmias (Antzelevitch and Fish, 2001; Chauhan et al., 2006). In agreement with previous studies (Zhang et al., 2008), cell simulation results indicate that membrane potential differences and inhomogeneous APDs of ventricular cardiomyocytes contribute to the tissue electrical heterogeneity. As for membrane potential differences, afterdepolarizations became inducible under the R858H condition when the SR calcium load was progressively increased in the MCELL cells. The maximal membrane potential differences both ENDO-M and EPI-M APs were significantly increased. As for inhomogeneous APD, the R858H mutation produced inhomogeneous APD prolongation in EPI, MCELL and ENDO cells, leading to the lager spatial gradient of APD along the transmural ventricular wall. Therefore, these electrical differences due to the R858H mutation at the cellular level contributed to an augmented transmural DOR, increasing tissue vulnerability to unidirectional conduction block.

Linking APD prolongation to unidirectional conduction block can provide insights into the risk of cardiac arrhythmias. The APD corresponding to the R858H model was longer than that of other models at S1-S1 = 1000 ms. The RT prolongation was consistent with the prediction of APD. At S1-S1 = 370 ms, the AP maps showed 2:1 block under the R858H condition. A similar 2:1 block was also observed under rapid conditions in the other heart models. Moreover, in the R858H heart, a local conduction block corresponded to unrecovered tissues and spiral waves were induced (Figure 14, the R858H heart model at S1-S1 = 348 ms). As the S1-S1 interval was further decreased (at PCL = 300 ms), we observed unidirectional conduction block and reentrant waves in the WT ventricles (Figure S4). This suggested that the prolonged APD tissues formed excitable obstacles which can anchor spiral waves.

Our simulations indicate that the pro-arrhythmic effect of the R858H mutation is reflected by the increased MPCL of 2:1 block. This is compatible with observations from previous studies (Zhu and Clancy, 2007), in which a long PCL led to 2:1 block under the G406R/G402S mutation condition. Compared with other mutations, unidirectional conduction block under the R858H condition occurred at a longer S1-S1 interval, implying higher risk of cardiac arrhythmias in R858H patients at relatively slow heart rates. Furthermore, the 2D and 3D simulated results support the notion that the R858H mutation increases the proarrhythmic risk. Although the spiral waves can be induced under all conditions, a relatively long S1-S1 interval is required for the initiation of reentrant waves under the R858H condition.

VF is the most common cause of sudden cardiac death. Fibrillation results when an electrical wave break induces reentry and triggers a cascade of new wave breaks (Weiss et al., 2005). In our simulations, three potential mechanisms underlying VF caused by the R858H mutation were considered, (a) the restitution hypothesis; (b) afterdepolarizations-mediated fibrillation hypothesis; and (c) dispersion of refractoriness hypothesis.

According to the restitution hypothesis, the chaotic excitation dynamics during VF are the result of dynamical instabilities in APD, the occurrence of which requires that the slope of the APD restitution curve exceeds 1 (Fenton et al., 2002). Although the APDR relationships were steepened by the R858H mutation, stable spiral wave dynamics occurred in homogeneous ventricular tissues. The measured maximum slope of APDR for the EPI cell was 1.5 which didn't exceed the minimum APDR slope (1.5) for the occurrence of spiral breakup in the TP06 model (Ten Tusscher and Panfilov, 2006a). Also, a spiral wave remained in the MCELL tissue with the APDR slope of 2 which was within the range of experimentally measured slopes (1.28–3.78; Pak et al., 2004). In addition, the slopes measured in the ORd model were within 1 under the R858H condition (Figure S7). Thus, our study confirmed that the effects of R858H-induced I_CaL are model independent. These results indicated that steep APDR is not the mechanism of VF caused by the R858H mutation.

According to afterdepolarizations-mediated fibrillation hypothesis, abnormal electrical excitations caused by afterdepolarizations can disrupt the normal propagation of electrical waves and cause life-threatening arrhythmias like VF. The mechanisms of afterdepolarizations generation include the reactivation of I_CaL and the enhancement of I_NCX arising from spontaneous calcium release. Simulation studies showed that three types of spiral fibrillation emerged when EAD was mainly induced by reactivation of I_CaL (Vandersickel et al., 2014; Zimik et al., 2015). For the R858H condition, afterdepolarizations were caused by the enhancement of the I_NCX arising from spontaneous calcium release. Although, focal activity occurred in the 1D fiber tissue, electrical coupling between cells suppressed the formation of PVC under the R858H condition. Experimental and computational studies were done to investigate how clumps of cells, eliciting afterdepolarizations in synchrony, give rise to triggered activities, which can disturb any prevailing course of wave propagation and induce electrical-wave turbulence (Sato et al., 2009; de Lange et al., 2012; Myles et al., 2012). Local synchronization of afterdepolarizations cells was suggested the mechanism underlying the formation of PVCs. Of course, our previous study has demonstrated that afterdepolarizations induced by SR calcium overload are synchronized to overcome the source-sink mismatch and produce PVCs in cardiac tissues (Bai et al., 2017). These studies support the notion that R858H-induced afterdepolarizations in MCELL cells are implicated in the induction of electrical-wave turbulence.

According to the dispersion of refractoriness hypothesis, wave break is produced by electrophysiological and anatomic heterogeneities in the tissue (Xie et al., 2001). For the R858H condition, intrinsic electrophysiological heterogeneity was sufficiently increased by the altered I_CaL due to a R858H mutation. Electrophysiological heterogeneity produced DOR responsible for 2:1 conduction block necessary for the fibrillation-like activity. Simulation studies also demonstrated that this breakup mechanism did not require steep APDR (Xie et al., 2001; Fenton et al., 2002). Thus, these data indicated that DOR was a breakup mechanism for the R858H condition.

Taken together, steep APDR is not the mechanism responsible for VF in the current work. Our results suggest that R858H induced afterdepolarizations and transmural APD dispersion are responsible for the potential mechanisms underlying the formation of VF. For afterdepolarizations-mediated fibrillation, it must be synchronized across many cells (Xie et al., 2010; Myles et al., 2012). And DOR caused by the R858H mutation increased the susceptibility to dynamic instabilities.

The presented results showed a QT interval prolongation from 397.1 ms (WT) to 425.7 ms (R858H), which was within the range of the QTc interval (420–476 ms) in R858H-probands (Fukuyama et al., 2014). Thus, a bridge between CACNA1C mutations and pro-arrhythmic phenotypes was built by the cardiac model. Also, the model can be used to investigate the effect of therapies targeting selected ion channels. More importantly, defining a classification based on a reliable VF risk factor would be very useful in guiding the selection of LQT8 patients who are in need of a ventricular defibrillator. In the current work, we have shown that, DAD-mediated AP at fast pacing rates isn't shortened but further prolonged, indicating QT interval prolongation at fast pacing (Viskin et al., 2010; Adler et al., 2012; Kaye et al., 2013). The QT abnormality may become a clinical marker to determine the VF risk in LQT8 patients.

The TP06 model was used to investigate the pro-arrhythmic effects of CACNA1C mutations and its limitations were discussed elsewhere (Zhang et al., 2008; Bai et al., 2016c). Here, we address several limitations in this study. (i) Due to the lack of a precise model of complex calcium cycling, this study does not consider the effect of a calcium spark on the genesis of afterdepolarizations and further refinement of the model is required. (ii) Our simulations show that R858H-induced defects can lead to the cellular arrhythmogenesis in the absence of neural influences. However, clinically observed arrhythmias typically occur in the setting of increased sympathetic tone (Best and Kamp, 2010), such as crying, and so future efforts will expand the study to investigate the effects of the autonomic nervous system on the development of cardiac arrhythmias. (iii) To maintain structures of the TP06 (Ten Tusscher and Panfilov, 2006a) and ORd (O'Hara et al., 2011) models, the original I_CaL conductance was used for the WT condition and the maximum conductance of each “mutant” I_CaL was determined by scaling relative current proportions. In agreement with the TP06 model (Ten Tusscher et al., 2004; Ten Tusscher and Panfilov, 2006a), the I_CaL was same for ENDO, MCELL and EPI cells in the ORd model (O'Hara et al., 2011) and electrophysiological heterogeneity caused by different I_CaL densities for three cell types was not considered in this study. Thus, special attention must be paid to explain these simulated data. (iv) In the multicellular tissue model, due to a lack of detailed experimental data on cell properties and well-delineated regions, the proportion of each region composed of each distinct cell type was chosen to produce a positive T-wave in the 1D transmural cable, similar to those used in other studies (Gima and Rudy, 2002; O'Hara et al., 2011). However, simulated ECG, which might be influenced by torso effects and the dimensions of the regions, wasn't completely consistent with clinical findings (Fukuyama et al., 2014). (v) Fibrillation at the tissue level resulted from afterdepolarizations at a single-cell level was not observed, but modeling and experimental studies have suggested that spatial-temporal synchronization of SR calcium overload and release can overcome the source-sink mismatch and produce focal arrhythmia in the heart (Xie et al., 2010; Myles et al., 2012), and so future efforts will expand the study to investigate that how the heart produces spatial-temporal synchronization of afterdepolarizations and the R858H mutation triggers VF in the heart. (vi) In the 1D and 2D tissue models, we assumed isotropic cell-to-cell electrical coupling. On the one hand, possible anisotropic intercellular electrical coupling may play important roles in the initiation and perpetuation of reentry. On the other hand, omitting anisotropy can be useful, in that changes to tissue behavior observed in the present study can be attributed with confidence to the implemented modifications to I_CaL. (vii) The space step for the 2D homogeneous ventricular model is different from that used for 2D transmural ventricular sheet. The space step is chosen to produce a stable spiral wave close to those in other studies (Ten Tusscher and Panfilov, 2006a; Vandersickel et al., 2014; Nayak and Pandit, 2015; Zimik et al., 2015; Nayak et al., 2017), in that rotating spiral can be attributed with confidence to the APDR slope which is below the minimum value for the occurrence of electrical instability. (viii) Although the 2D homogeneous tissue was designed to examine if mutations-induced changes in APDR curves promote the breakup of a spiral wave and we expected the simulated results to hold in tissues of reasonable size, a reasonable domain cannot avoid the frequent termination of spiral waves by collisions with the boundaries of the simulation domain. Therefore, the size of the 2D homogeneous tissue chosen is similar to that used in other studies (Ten Tusscher and Panfilov, 2006a; Vandersickel et al., 2014; Nayak and Pandit, 2015; Zimik et al., 2015; Nayak et al., 2017), but it is larger than that of a typical human heart. Special attention must be paid to explain these simulated data. (ix) Due to a lack of fiber orientation data, fiber orientation was not included in the 3D model and its effect on reentry was not studied. Nevertheless, whilst it is important to make explicit the potential limitations of the approaches adopted in the present study, these potential limitations are not expected to influence fundamental conclusions that can be drawn from our data on the mechanisms by which increased I_CaL due to the R858H mutation can increase risk of ventricular arrhythmias.

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.