Ranolazine as an Alternative Therapy to Flecainide for SCN5A V411M Long QT Syndrome Type 3 Patients

The prolongation of the QT interval represents the main feature of the long QT syndrome (LQTS), a life-threatening genetic disease. The heterozygous SCN5A V411M mutation of the human sodium channel leads to a LQTS type 3 with severe proarrhythmic effects due to an increase in the late component of the sodium current (INaL). The two sodium blockers flecainide and ranolazine are equally recommended by the current 2015 ESC guidelines to treat patients with LQTS type 3 and persistently prolonged QT intervals. However, awareness of pro-arrhythmic effects of flecainide in LQTS type 3 patients arose upon the study of the SCN5A E1784K mutation. Regarding SCN5A V411M individuals, flecainide showed good results albeit in a reduced number of patients and no evidence supporting the use of ranolazine has ever been released. Therefore, we ought to compare the effect of ranolazine and flecainide in a SCN5A V411M model using an in-silico modeling and simulation approach. We collected clinical data of four patients. Then, we fitted four Markovian models of the human sodium current (INa) to experimental and clinical data. Two of them correspond to the wild type and the heterozygous SCN5A V411M scenarios, and the other two mimic the effects of flecainide and ranolazine on INa. Next, we inserted them into three isolated cell action potential (AP) models for endocardial, midmyocardial and epicardial cells and in a one-dimensional tissue model. The SCN5A V411M mutation produced a 15.9% APD90 prolongation in the isolated endocardial cell model, which corresponded to a 14.3% of the QT interval prolongation in a one-dimensional strand model, in keeping with clinical observations. Although with different underlying mechanisms, flecainide and ranolazine partially countered this prolongation at the isolated endocardial model by reducing the APD90 by 8.7 and 4.3%, and the QT interval by 7.2 and 3.2%, respectively. While flecainide specifically targeted the mutation-induced increase in peak INaL, ranolazine reduced it during the entire AP. Our simulations also suggest that ranolazine could prevent early afterdepolarizations triggered by the SCN5A V411M mutation during bradycardia, as flecainide. We conclude that ranolazine could be used to treat SCN5A V411M patients, specifically when flecainide is contraindicated.


1
Expanded materials and methods

In silico patch-clamp protocols
To study the sodium channel dynamics and reproduce the experimental results, we reproduced the described conditions (pH, temperature and sodium concentrations) while integrating the sodium current over the described time and voltage protocols. Temperature derived changes were applied with a scaling factor calculated with a Q10 of 3 as follows.
Patch clamp simulations were usually performed at room temperature (according to the reference experimental protocol), while action potential simulations were always carried out at 37ºC.
Both extracellular and intracellular sodium concentrations were fixed during in silico patch-clamp experiments, while intracellular sodium concentrations were dynamically updated in action potential cellular models. Results were also extracted following the methods described in their respective references. These include indications concerning when and where to measure the current as well as how to normalize the data if needed. Experimental data points were manually digitized from their original figures.
Patch-clamp protocols were simulated using custom Matlab2014b (The Mathworks, inc) code including an ordinary differential equation (ODE) solver with variable time step limited to a maximum of 0.5 ms.

General Optimization Protocols
We used the function "fminsrchbnd" (Moreno et al., 2011(Moreno et al., , 2013(Moreno et al., , 2016, which is an implementation of the Nelder-Mead simplex that allows boundaries to be set to the parameters that ought to be fitted. Termination conditions were met when a maximum number of iterations -or cost function evaluations -was reached (300). This could also happen when error variation between consecutive iterations stayed below a tolerance level of 0.01, but that rarely was the case due to the limited number of iterations.
Scores for each test were computed as the sum of squared differences between test results and reference data. Dividing the score by the test's number of samples was done to normalize the contribution of each test. Nonetheless, to compare scores from very different units (i.e.: milliseconds and normalized values), custom scalar factors were added to each test's score. The final error score of a cost function was determined by the sum of the scores of every single test in that function.
Optimizations were prone to being unable to surpass some local minima rather than finding the best solution (Moreno et al., 2016). Specifically, due to the models being already fit to several of the tests, the models were not able to "exit" the current local minimum. Therefore, we gave each optimization a randomly generated vector of initial variables (10% variability) and launched several optimizations for every model. To ensure reproducibility, we controlled the initial vectors by giving the random number generator a seed to start with.

Optimization of the Wild-Type Sodium Current Model
The optimization of the Wild-Type sodium channel contains the following tests of which six were extracted from Moreno et al. 2011(Moreno et al., 2011 including steady state availability (SSA), activation (ACT), recovery from inactivation (RFI), recovery from use-dependent block (RUDB) and time to 50% activation (Tau50). Three new tests were added running cellular simulations to build the APD90 Restitution Curve (REST), optimize the slow sodium current IV (INaL IV) and time course (INaL) according to AP clamp recordings (Horvath et al., 2013;Hegyi et al., 2018). The model was also constrained by maximum upstroke velocity (max dV/dt) and channel Mean Opening Time (MOT).
Before starting the optimization procedures, we stabilized the cellular model by running a 300second train of pulses at several Basic Cycle Lengths or BCLs (300, 400, 500, 1000, 1500 and 2000 ms). These simulations served as steady state starting points for INaL and REST tests.
Every iteration, we first stabilized the INa Markovian model by applying a single pulse from -100 mV to -10 mV during 200 ms, followed by a 5-s, -100 mV, resting membrane potential phase. We used the final states from the latter as a starting point for every test involving patch-clamp protocols.
INaL time course: We applied a train of 40 stimuli at 1 Hz to the modified isolated endocardial model and saved the last beat for analysis. We extracted three parameters from the INaL time course as can be observed in Supplementary Figure S1. Firstly, we defined the "dome" as the maximum late sodium current (registered after the fast sodium current had inactivated). Secondly, we defined the "valley" as the minimum current between peak INa fast and the dome. Finally, we defined the halftime taken by the current to reach the dome from the valley as 1 2 ⁄ . Both valley and 1 2 ⁄ were normalized to the time and current of the dome. We set reference dome current to -0.34 pA/pF (Dutta et al., 2017) , a dome-to-valley value to 0.59 and 1 2 ⁄ value to 0.63 (Horvath et al., 2013).
Steady state availability: We applied a single 25 ms pulse to -10 mV from a 5-second variable potential test pulse from -120 mV to -40 mV (5 mV intervals), extracted peak currents elicited by the second pulse, normalized them to tonic block at -120 mV and plotted them against pulse voltage.
Activation: We applied a 25 ms test pulse to variable potentials from -80 mV to 20 mV from a resting potential of -100mV. Then, we calculated the resulting channel conductance from the elicited peak currents, normalized the results to maximum conductance and plotted them against pulse potential. Conductance values for test potentials above the potential of maximum conductance were set to 1.
Recovery from inactivation: We evaluated recovery from inactivation with a standard double-pulse protocol from -100mV to -10mV. The second pulse was delayed with increasingly higher time intervals ranging from 0.1ms to 6s. We extracted peak currents elicited by the second pulse and normalized them to peak currents elicited during the first pulse, then represented the data against time intervals.
Recovery from use-dependent block: First, we simulated a train of 300, 25-ms pulses, from -100 mV to -10 mV, at a pacing rate of 25Hz. Potential was set back to -100mV before a last pulse to -10 mV was applied after a variable delay ranging from 0.5 ms to 9 s. We extracted maximum current peaks elicited during the second pulse and normalized them to their maximum value. We represented the data against time intervals.
Tau 50% activation: Following the Activation protocols, we measured Tau50 as the time the current took to reach 50% of the maximum peak current at the beginning of each pulse.
APD90 restitution curve: We performed 40-second simulations with the modified endocardial model at 300, 400, 500, 1000, 1500 and 2000 ms BCLs. We carried every simulation in parallel and saved the last beat for further analysis. We extracted the APD90 of saved beats by calculating the interval between de time of max dV/dt and the time of 90% repolarization. We plotted the data against BCL and compared the results to the ORd model's reference (O'Hara et al., 2011).
Maximum upstroke velocity was extracted during BCL 1000 ms simulations. Target value was set to 250 V/s.
MOT was calculated as in Moreno et al. 2011(Moreno et al., 2011.

Optimization of the SCN5A V411M Mutation Model
Optimization of the SCN5A V411M mutation model included protocols testing the following dynamics, namely, activation (ACT), inactivation (INACT), inactivation time constants (TauINACT), current-voltage relationship (IV) and prolongation (Prol). Except the latter two, all protocols follow the same methods as in Horne et al. 2011.
Before the optimization, we paced the single cell model with a train of 300 square stimulus at 1 Hz to reach steady state.
We ran one iteration of the cost function before starting the optimization in order to create a starting set of curves corresponding to the wild type model. Target values were generated and applying the relative wild-type-to-mutation changes in current dynamics observed by Horne et al. (Horne et al., 2011).
Activation: We applied a 200-ms pulse to a variable potential from -80 mV to 30 mV and extracted the resulting maximum conductances from the elicited peak currents, which we normalized to the highest value. We represented the data against pulse potential and fitted it to the following Boltzmann equation: Where Vh is the half-maximal voltage and s is the slope of the curve. Both parameters were compared to their respective wild type model values to check for voltage shift (Vh -Vh_wt) and slope change (s -s_wt).
Inactivation: We simulated a variable voltage pre-pulse of 300 ms from -160 mV to 0 mV from a resting potential of -110 mV, followed by a 20-ms pulse to -20 mV. Then, we extracted peak currents elicited during the second pulse, normalized the data to their maximum value and plotted them against pre-pulse potential. Finally, we fitted the curve to a Boltzmann equation as in ACT and extracted Vh and slope shifts.
Mean inactivation time constants: Inactivation current time courses from the ACT protocol were isolated (from maximum peak current to the end of the pulse) and fitted to a single exponential function such as: Where is the current, is the time and and are constants. The resulting τ values were plotted against pulse potential.
Current-voltage relationship: We applied a variable voltage 250-ms pulse from -80 mV to 40 mV from a resting potential of -110 mV. Elicited peak currents were extracted, normalized to their maximum value and plotted against pulse potential.
APD90 Prolongation: We ran 40-beat simulations with the heterozygous SCN5A V411M isolated endocardial model (50% wild type and 50% mutated currents), saving the last beat for further analysis. Prolongation of the APD90 was measured relative to the wild type model. We used a target prolongation of 16% as a surrogate of the QTc prolongation.
qNaf: the total charge carried by INaf was taken into account during the APD90 prolongation simulations by comparing it to the wild type value. This prevented the model from not depolarizing.

Optimization of Flecainide
Optimization Before starting the optimization, the isolated endocardial cell model was simulated for 300s at 1Hz pacing rate in drug-free conditions to reach steady state.
Flecainide enhances the inward rectifier current (IK1) as found by Caballero et al. 2010(Caballero et al., 2010, and it is estimated that 1.5 µM therapeutic flecainide (DailyMed) should increase its conductance to a 151%. Flecainide also blocks the rapid potassium delayed rectifier current (IKr) with an IC50 of 3.91 µM (Paul et al., 2002), which reduces the its conductance to a 72.27%.
Simulations trying to assess flecainide's INaL IC50 gave very high (around 90 µM) values (not shown) before optimization, revealing a lack of INaL block at therapeutic concentrations. Therefore, normal and bursting state affinities of the charged drug for the channel, but not diffusion nor neutral drug affinities, were allowed to change during the optimization to enable fitting of the fast and late currents block dynamics.
Steady-State Availability: We used the same protocol as in the wild-type model to obtain a similar curve under the effects of 10 µM flecainide.
Recovery from use-dependent block: First, we simulated a train of 100 pulses at 25 Hz from -100mV to -10mV. Then, a second identical pulse was applied after a variable time from 0.5 to 9s at -100mV. Concentrations of flecainide were fixed at 10µM. Peak currents elicited from the delayed pulse were extracted, normalized to tonic block and plotted against time intervals.
INaf block curves: The standard protocol consisted of a 30-ms pulse from -100 mV to -20 mV. This protocol was applied 40 times at a rate of 0.2Hz or 60 times at a rate of 1 or 3 Hz. The maximum peak current elicited by the last pulse was extracted for increasing flecainide concentrations and normalized to drug-free conditions.
INaL IC50: From a -120 mV 200 ms pulse, we applied a 40-ms pre-pulse to -15 mV followed by a 200ms pulse to 40 mV. We measured INaL as the maximum elicited current during a ramp from 40 to -95 mV (-1.35 V/s) following the beforementioned pulses. The protocol was repeated 3 times for each of several drug concentrations until 50% block of the control current (drug-free) was reached. The corresponding concentration was retrieved as the INaL IC50.
APD90 Prolongation: We performed simulations with the single cell model to test the effects of flecainide every iteration. We applied a train of 40 beats at 1 Hz with 1.5 µM therapeutic flecainide. APD90 from the last beat was normalized to drug free conditions to obtain drug-induced prolongation. Target prolongation was adjusted to be as low as possible.
Because flecainide has a pKa of 9.3 (Moreno et al., 2011) it is 99% charged under physiological conditions, before optimizing the complete drug modelin a first phasewe optimized only neutral flecainide for DDUDB (see below) at 10 Hz and RUDB for time intervals ranging from 0.5s to 7s. Then, neutral drug parameters were held constant while the remaining parameters were optimized using all protocols in two additional phases. The second phase included the complete drug versions of the abovementioned protocols as well as SSA, RUDB and the INaf02, INaf1 and INaf3 set. Finally, the third phase added INaL and APD90 prolongation in wild-type to the test batch.
Dose-dependent use-dependent block for neutral flecainide (DDUDB): We simulated a train of 300, 25-millisecond, pulses at 10Hz from -100mV to -10mV. We normalized the current elicited by the last pulse under increasing neutral flecainide concentrations and normalized the value to the first pulse. We plotted the results against drug concentration.
Flecainide optimizations were terminated as usual when the conditions were met. Nonetheless, to improve the INaL IC50, we searched for best value in the parameter history giving similar fitting results in the other tests.

Optimization of ranolazine
Optimization of ranolazine consisted of the following test protocols: steady-state availability (SSA), tonic block of peak and late sodium current (TB), use-dependent block (UDB), recovery from usedependent block (RUDB) and frequency-dependent use-dependent block (FDUDB). Diffusion and affinities for the normal and bursting states were set as in the original model and prevented from being modified.
Steady-state availability: we applied a 100 ms pulse to -10 mV from a variable conditioning potential of -120 mV to -40 mV and extracted the resulting peak current under exposure to 10 µM ranolazine.
Values were normalized to drug-free conditions and plotted against conditioning pulse potential.
Tonic block: A single 500-millisecond pulse to -10 mV from a resting potential of -100 mV was applied to evaluate the effect of increasing ranolazine concentrations. We extracted INaf as the elicited peak current, while INaL was measured as the remaining current at the end of the pulse. Both values were normalized to drug-free conditions and plotted against drug concentrations.
Use-dependent block: We applied 300, 25 ms, pulses to -10 mV from a resting potential of -100 mV at a rate of 5 Hz and under exposure to increasing ranolazine concentrations. We saved the peak current elicited by the last pulse, normalized it to drug-free conditions and plotted the resulting values against drug concentrations.
Recovery from use-dependent block: This protocol is similar to flecainide's RUDB protocol. Intervals ranged from 0.1 s to 10 s and ranolazine concentrations were set to 10 µM.
Frequency-dependent use-dependent block: A train of 300 square 25-ms pulses from -100 mV to -10 mV were applied at 1 Hz, 2 Hz, 5 Hz and 10 Hz rates under exposure to 100 µM ranolazine. The peak sodium currents elicited by the last pulse were extracted and normalized to drug-free conditions, then plotted against pacing rate.
Supplementary Figure S2. Calculation of the three parameters we used to evaluate the time course of INaL (dome, valley and t_1/2). Valley was normalized to the dome value, while t_1/2 was normalized to the total time between dome and valley.
Supplementary Figure S4. Comparison of the simulated action potential time courses for isolated endocardial cells after a train of 40 pulses (continuous lines) and at the steady-state (300 pulses, dashed lines). Wild type (black lines) and V411M mutated (red lines) cells in control and under exposure to therapeutic concentrations of flecainide (green lines) and ranolazine (blue lines). The very small differences between the simulated action potential time courses obtained with 40 and 300 pulses corroborate the validity of applying 40 stimuli to the isolated endocardial models in order to reduce the computational cost during the optimization procedure of the wild type, flecainide and mutation INa models.  Table T2. Equations of the transition rates from the wild type and V411M INa models. P1 to p13 are the parameters fitted during the optimizations and they were restricted to positive values. is the temperature factor and is the membrane voltage.  Table T3. Transition rates and affinities of the flecainide model. P1 to p16 are the parameters that were fitted during the optimizations and they were restricted to positive values.

Supplementary
Transition rates (ms -1 ) equations α11+ and α11n 11 α12+ and α12n 12 β11+ and β11n 11 β12+ and β12n 12  Table T4. Transition rates and affinities of the ranolazine model. P1 to p10 are the parameters fitted during the optimizations and they were restricted to positive values. P11 and p12 were incorporated only for explaining the supplementary Figure S5, but they were not modified during the optimization of the ranolazine model.
Transition rates (ms -1 ) equations α11+ and α11n 11 α12+ and α12n 12 β11+ and β11n 11 β12+ and β12n 12  Supplementary Table T6. INa V411M mutation model optimization error calculations and weights in the cost function. The words containing "data" and "reference" indicate vectors containing the results of our tests (red lines in Figure 4 of the main text) and the reference target values (open squares in Figure 4 of the main text), respectively. N is the number of samples.  Supplementary Table T7. Flecainide model optimization error calculations and weights in the cost function. The words containing "data" and "reference" indicate vectors containing the results of our tests (dots and black lines in Figure 5 of the main text) and the reference experimental values (open squares in Figure 5 of Supplementary Table T8. Ranolazine model optimization error calculations and weights in the cost function. The words containing "data" and "reference" indicate vectors containing the results of our tests (black lines in Figure 6 of the main text) and the reference experimental values (open squares in Figure 6 of the main text), respectively. N is the number of samples.