Protein-Mediated Electroporation through a Cardiac Voltage Sensing Domain Due to a nsPEF Stimulus

1. Introduction

nsPEF, a technology that emerged in 1995 [1], has seen a significant surge in research interest since 2005 [2]. nsPEF’s ability to elicit specific cellular effects [3] has resulted in an impressive range of applications, including the activation of neurons [4,5,6,7,8,9] and myocytes [10,11,12,13], wound healing [14,15,16,17], manipulation of phenotype [18], modulation of gene expression [19,20,21,22,23,24], antiparasitic effects [25,26,27], enhancement of immune response [28,29,30,31,32,33], cell proliferation [18,34,35,36], improved fermentation [37,38], sterilization for the food industry [39,40,41], seed germination [42,43,44], and most notably, the development of novel cancer therapies [2]. Recently nsPEF have been proposed for virus inactivation [45].

The nsPEF technique involves delivering high electric field pulses (∼1-300 kV/cm) in nanoseconds or even picoseconds into biological tissues or cells. Molecular Dynamics studies have shown that nsPEF is theoretically capable of forming membrane nanopores [46,47]. Experimental results have also strongly supported the formation of nanopores due to nsPEF stimulus [5,48,49]. While direct evidence of nanopore formation is lacking due to their small size (~2nm) and transient nature, the scientific community largely supports this phenomenon as the primary effect of nsPEF on cells. Further experimental results have shown that the activation of voltage-gated channels (VGCs) is another primary effect of nsPEF [4,8,50,51,52,53]. Therefore, the effects of nsPEF on cells can be broadly divided into two categories: formation of membrane nanopores and activation of VGCs.

However, recent research using molecular dynamics (MD) simulations has suggested an additional effect of nsPEFs: the creation of pores in transmembrane proteins. In 2018, the first study to investigate pore formation in transmembrane proteins (specifically, human aquaporin) using MD simulations was published [54]. Later, in 2020, Rems et al. [55] also used MD simulations to investigate nsPEF-induced pore formation in three distinct voltage-gated channels (a bacterial VGNC, an eukaryotic VGNC, and a human hyperpolarization-activated cyclic nucleotide-gated channel), observing both simple and complex pore formation in the VSD. Complex pores were proposed to be the source of theorized lipidic pores and may be stabilized by the presence of ions and other channel components, such as TMHs [56,57,58]. Recently, our group also observed the formation of complex pores in a VSD of a human VGCC using an MD simulation, under the application of an external electric field of 0.2 V/nm lasting for 50 ns to mimic a nsPEF stimulus [59]. However, there is currently no direct experimental evidence for the formation of nanopores or transmembrane pores due to nsPEF stimulation. Therefore, we strongly encourage further experimental investigation in this area. The potential biotechnological applications for such an nsPEF device capable of creating protein pores are highly intriguing.

Ubiquitous and structurally similar integral membrane proteins, known as voltage-gated channels (VGCs), include Na${}^{+}$, Ca${}^{2+}$, and K${}^{+}$ channels [60]. While voltage-activated Na${}^{+}$ and Ca${}^{2+}$ channels consist of monomers with varying auxiliary subunits, voltage-activated K${}^{+}$ channels are tetramers. All of these channels have four repeated structures, each formed by six $\alpha$-helix, as shown in Figure 1. The VSD is composed of helix S1 to S4, while the pore domain is formed by helix 5 and 6 [61]. Notably, the four VSD structures are present in all Na${}^{+}$, Ca${}^{2+}$ and K${}^{+}$ [62].

In this article, we investigate the formation of simple and complex pores in the VSD of a Nav1.5-E1784K voltage-gated sodium channel (VGNC). Nav1.5 is the primary VGNC channel in the heart, and mutations in Nav1.5 have been linked to various cardiac disorders, including type 3 long QT syndrome (LQT3) and Brugada syndrome (BrS). The E1784K mutation is commonly found in patients with both LQT3 and BrS [63]. The Nav1.5 VGNC used in our study is in its native conformation and is not bound to any external molecule. This mutated VGNC have a structure nearly identical to that of the wild-type human Nav1.5. Interestingly Nav1.5, present S4 helix with three arginines instead of four, as opposed to the four arginines found in most VGCs. These arginines play a crucial role in the activation of VGCs [64,65,66,67,68,69], as they provide a high charge density that can respond to changes in membrane potential. Therefore, investigating the formation of VSD pores with only three arginines in the S4 helix, which is expected to be the driving force for pore formation, is of particular interest. The two existing studies that have explored pore formation due to a nsPEF stimulus have used VSDs with four arginines in their S4 helix [55,59].

Recently, nsPEF stimulation has been proposed as a new defibrillation method to achieve a higher efficiency of defibrillation on the first shock, reduce shock energy, minimize side effects, and lower the probability of reinduction of arrhythmias [70]. According to Azarov et al. [12], applying nsPEF stimulus to cardiomyocytes results in the formation of membrane nanopores, which are less harmful to transmembrane leaks compared to conventional defibrillation using millisecond pulses. In addition, transient inhibition of Na${}^{+}$ and Ca${}^{2+}$ VGCs [71,72] may aid the antiarrhythmic effect of nsPEF defibrillation. Nevertheless, the mechanism behind nsPEF-induced action potentials generation remains poorly understood. Understanding the formation of VSD pores in cardiomyocytes due to nsPEF stimulus may contribute to comprehending the generation of action potential in these cells. Therefore, studying the formation of VSDs pores in cardiomyocytes may shed light on the biophysics of using nsPEF as a defibrillator.

To further investigate the formation of VSD pores, we created a simulation box for molecular dynamics (MD) using a pure 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) bilayer, a common phospholipid in MD simulations. The VSD was immersed within the POPC bilayer for subsequent analysis. To determine the external electric field that would best induce a change in VSD conformation to express the most VSD pore-rich amalgam. The simulation box was subjected to various electric fields. The optimal electric field, found to be 0.13 V/nm, and 200 replicates at this electricfield was performed. The resulting structures were clustered with a machine learning tool to identify the most representative ones. This analysis will facilitate further studies of VSD pore formation.

2. MD simulations methods

The simulations were conducted using the GROMACS molecular dynamics package 2020.6 [73], and the CHARMM36 force field was utilized to account for parameters including bond length, angle bending, angle torsions, and non-bonded interactions of the molecular systems [74]. A rectangular simulation box was constructed with a VSD at the center of the bilayer, containing 166 POPC molecules in a bilayer arrangement, 9,713 H${}_2$O molecules, and 1 K${}^{+}$ ion to neutralize the system. Prior to the molecular dynamics simulation, a steepest descent minimization algorithm was applied to relax the molecular system, followed by 6 steps of equilibration, as detailed in SI section 1. Box 1 was equilibrated with a size of 7.79 nm x 7.79 nm x 8.50 nm, which ensured that the VSD did not see its periodic images. The leap-frog integrator algorithm was utilized with a time step of 2 fs to integrate Newton’s equations of motion. The Lennard Jones potential (LJ) was used to simulate Van der Waals interactions, and the Particle Mesh Ewald (PME) method was used to calculate long-range electrostatic interactions [75].

The interactions were truncated using a cut-off of 1.2 nm for both the LJ and electrostatic forces. The bond lengths were constrained using the linear constraint solver (LINCS) [76]. The simulations were performed in an NPT ensemble at a temperature of 310 K and a pressure of 1 bar, with coupling to a Nose-Hoover thermostat [77] and a Parrinello-Rahman barostat [78]. Time constants of 1 ps and 5 ps were used for temperature and pressure, respectively. Periodic boundary conditions were applied in all directions. The TIP3P (transferable intermolecular potential with 3 points) water model [79], which assigns charges and Lennard-Jones parameters to each of the 3 atoms in an explicit 3-atom rigid water molecule, was used for all simulations.

The Membrane Builder module of CHARMM-GUI [80] was utilized to embed the VSD in the bilayer. The equilibrated box underwent a 200 ns simulation Figure 2, after which different values of $\overrightarrow{E}$ were explored to determine the optimum one capable of creating simple and complex pores in the VSD. Eight replicas of each $\overrightarrow{E}$ were simulated, with a range of 0.1 V/nm, 0.12 V/nm, 0.14 V/nm, 0.16 V/nm, 0.18 V/nm, and E0.2 V/nm. The Root Mean Square Deviation (RMSD) of the VSD was analyzed for each $\overrightarrow{E}$, and a $\overrightarrow{E}$ of 0.13 V/nm was selected. A continuation of 200 replicas was then simulated at E=0.13 V/nm for 50 ns to explore all phases of the VSD structural changes. The direction of $\overrightarrow{E}$ was consistent across all simulations, with the vectorial force aligned antiparallel to the z-axis and perpendicular to the bilayer. The VSD used in the study was obtained from electron microscopy of a human Nav1.5-E1784K channel (PDB ID: 7DTC), with a resolution of 3.30 Å.

3. Simulation Results

3.1. Exploring different external Electric Fields

First, we performed an exploration by applying different $\overrightarrow{E}$ strengths along the z-axis to a Box containing an equilibrated cardiac VSD. Our goal was to determine the optimal $\overrightarrow{E}$ that would induce the most significant structural changes in the VSD, leading to the formation of both protein pores and complex-pores. The $\overrightarrow{E}$ values used ranged from 0.1 V/nm to 0.2 V/nm in increments of 0.02 V/mm. Previous studies have employed voltage-sensitive dyes [81,82] and microelectrodes [83] to determine that a membrane potential of 1.5 V represents the largest potential that the membrane can support before discharging, others studies support this voltage value [82,84,85,86,87,88]. Gumbart et al. [89] calculated that the best approximation for determining membrane potential in molecular dynamics simulations along the z-axis is $\overrightarrow{E_z}$$L_z$,, where $L_z$ represents the length of the simulation box of the z-axis. We therefore explored $\overrightarrow{E}$ values up to 0.2 V/nm, given that the membrane potential for an $L_z$ value of 8.5 is 0.17 V, which is above the limit of 1.5 V. We quantified the structural changes in the VSD induced by each $\overrightarrow{E}$ value by measuring the RMSD of eight replicas per $\overrightarrow{E}$ value.

As evidenced in Figure 3A and Figure 3B, at an electric field strength of $\overrightarrow{E_z}$=0.10 V/nm and $\overrightarrow{E_z}$=0.12 V/mn, respectively, the VSD displayed in Figure 3 exhibits negligible changes in RMSD across each of the replicas. This suggests that the internal forces governing the structural integrity of the VSD remain resilient to external stimuli up to a certain threshold. Conversely, when exposed to electric fields of $\overrightarrow{E_z}$=0.14 V/nm or greater, the majority of replicas displayed higher RMSD values, with a final RMSD of 2.9001150 nm at 50 ns, resulting in a 90°rotation of the membrane and yielding an output structure that lacks biological significance (Figure 3C). This phenomenon was observed across most of the replicas subjected to electric fields beyond $\overrightarrow{E_z}$=0.14 V/nm (Figure 3D, E, F). In Figure S1 of the Supplementary Information, we illustrate the bilayer’s turning point in response to external electric fields ranging from $\overrightarrow{E_z}$=0.14 up to 0.20 V/nm. For this reason, we opted to employ clusterization in machine learning analysis of VSD structures forming pores under external electric fields of $\overrightarrow{E_z}$=0.13 V/nm for 50 ns. We conducted 200 replicas at $\overrightarrow{E_z}$=0.13 V/nm to ensure complete phase space coverage and obtain the difference in free energy ($\Delta G$) between clusters. The details of the methodology used to compute $\Delta G$ are elucidated in the subsequent section. Replicas subjected to an external electric field of $\overrightarrow{E_z}$=0.13 V/nm and having the highest final RMSD of 2.5705926 nm exhibited a POPC membrane that retained its bilayer arrangement. Therefore, is expected that lowers final RMSD of the rest of the replicas under an external $\overrightarrow{E_z}$=0.13 V/nm do not turn over it bilayer being consisting for further analysis. Despite replicas under an external $\overrightarrow{E_z}$=0.13 V/nm ending in higher RMSD not turn over they membrane complete, the encompass a evident membrane rearrangement, as can me observed in Figure 4.

3.2. Clusterization

To identify the most representative conformational changes of the voltage-sensing domain (VSD), we employed a clustering approach based on different replicas. This analysis involved an evaluation of the fraction of native contacts (FNC) and the RMSD of the VSD. Initially, we discarded unstable replicas that did not attain a stable structure. To do so, we assessed the change in the first derivative of the RMSD over the last 10 nanoseconds of the simulation, and only those replicas with a change of less than 15% were retained. We selected 162 stable replicas that met this criterion. We further filtered these replicas to exclude those that ended with an RMSD of less than 0.35 nm and did not form a conclusive VSD pore (as depicted in Figure 5). This filtering step yielded a final set of 142 replicas that were suitable for clusterization. This filtering step improved the quality of the clusterization by ensuring that the resulting structures had full VSD pores. To gain a better understanding of how VSD pores change with the final RMSD, refer to Figure 6. For additional examples of formation of pores see Figure S3. Too have an idea of the kinetic formation of pores, we construct videos of the water density along z axis in the 50 ns of simulation for different replicas with different final RMSD. To see these videos go to $https://github.com/DLab/article-nspef-pores.$

To clusterize the data, we constructed a 3D cube with tuples that contained 3 coordinates, which were represented in a 3D graphic. The first coordinate was defined as the RMSD between two frames. In this context, frames refer to each VSD structure presented during the last 10 ns of simulation of the 142 replicas that were selected.

The second coordinate was defined as the FNC between the initial frame of the replica containing the first frame selected to obtain the RMSD and the first frame used to obtain the RMSD. The third coordinate was defined as the FNC between the initial frame of the replica containing the second frame selected to obtain the RMSD and the second frame used to obtain the RMSD.

A total of 142,000 frames were obtained from the 142 replicas, given that the total time of simulation was 10 ns and frames were extracted every 10 ps. This yielded a total of 142,000${}^3$=2.86x10${}^{15}$ coordinates, which was computationally unmanageable. To address this issue, we reduced the number of frames by obtaining frames every nanosecond, resulting in 10 frames every 10 ns.

Then, we constructed a vector with 1,420 RMSD values and two additional vectors with 1,420 FNC values. We performed combinatorics on the components of each vector to construct the tuples, yielding a total of 2,863,288,000 tuples. To work with this large quantity of tuples, we normalized them, resulting in a reduction to 5,616,518 tuples, each of which had a corresponding weight.

The weight of each tuple was different since tuples could share the same coordinates. We plotted a reasonable quantity of 5,616,518 tuples, Figure S4, which was manageable with our computing power.

. One technique of clusterization involves the unsupervised K-means algorithm, is commonly employed to segregate data into distinct groupings based on shared characteristics. For a comprehensive explanation of the K-means methodology, refer to reference [90]. In essence, K-means is a machine learning algorithm that assigns objects to k clusters based on their features. This clustering procedure involves minimizing the aggregate distance between each object and the centroid of its designated cluster.

The elbow method [91] is a popular technique used to determine the optimal number of clusters in a dataset for unsupervised learning algorithms such as K-means clustering. The method involves plotting the number of clusters against the corresponding values of the Within Cluster Sum of Squares (WCSS), which measures the total distance between data points within each cluster.As the number of clusters increases, the WCSS decreases, indicating a better fit to the data. However, there is a point at which the rate of improvement decreases, resulting in a gradual change in the slope of the plot. This point is called the elbow point, and it represents the optimal number of clusters. The elbow point can be visually identified as the point where the curve starts to flatten out, indicating that the incremental gain in clustering performance is not worth the additional complexity of adding more clusters. The elbow method is a simple yet effective tool for determining the appropriate number of clusters for a given dataset.

The clusterization method yielded 3 clusters, named from 1 to 3 in descending order, with respect to their size. Cluster 1 contains 2,269,895 tuples, cluster 2 contains 2,200,065 tuples and cluster 3 contains 1,146,558 tuples. The cloud of the 5,616,518 tuples was plotted in a 3D graph, and each resulting cluster was coloured differently, additionally each cluster was plotted separated, Figure 7.

When considering the number of executed independent replicas and produced frames, resulting in a 3D data cube with over 5 million points, it can be argued that the phase space describing the conformational changes of the VSD under the influence of external stimuli mimicking nsPEF has been thoroughly sampled. The total set of frames collectively represents the kinetically accessible set of possible conformations of the molecular system. Therefore, it can be assumed that all possible conformations have been sampled, resulting in an equilibrium between clusters. The free energy between clusters can now be computed by relying on established methods. To see how the tuples cover the phase space go to Figure S4.

$$\Delta G^0=-RTln\left(\frac{P1}{P2}\right),$$

(1)

where R is the molar gas constant, T is the temperature, and in our case $P1$ and $P2$ represent the number of tuples on each cluster. The calculated $\Delta$G between all clusters appears in Table 1.

A 2D representation depicting the density of the Tuples for each cluster, considering the RMSD and the FNC, can be seen in Figure 8

From the density profile of each cluster, the closest tuple to the maximum density was obtained using the least squared method. This method was used over a kernel density estimation (KDE) constructed with all tuples using the quadratic error to find the tuple closest to the maximum density point in the KDE. Each of these tuples has an associated frame from the simulations that correspond to the most representative conformation for each cluster. In Figure 9 are represented each representative VSD of each cluster in new-cartoon, seen above the x-y plane and above the z-y plane, and finally are represented each VSD new cartoon with water in van der waals spheres.

4. Discussion

The activation of ionic channels occurs due to the movement or torsion of the S4 helix of voltage-sensing domains (VSDs), which is triggered by a change in the membrane potential. This phenomenon is attributed to the highly charged nature of the S4 helix, which contains four charged arginine residues in its structure. Various models have been proposed to explain the transfer of charge during the activation of ionic channels, and all of them are linked to the movement of the S4 helix [92,93]. These models include the helical screw-sliding model [64,65], the kinetics model [66], and the paddle model [67].

Thus, it is expected that the S4 helix would be the primary factor driving the formation of pores and complex pores in response to a nsPEF stimulus. However, this does not appear to be the case. This study involving a single VSD from the cardiac sodium channel Nav1.5-E1784K, which only contains three arginine residues in the S4 helix (most VSDs have four arginine residues), showed that pores and complex pores were only generated at a membrane potential of $\overrightarrow{E_z}$$L_z$=1.1 V. In contrast, other studies using VSDs with four arginine residues required 1.5 V to generate pores and complex pores in the same POPC bilayer [55], and 1.75 V in a more resistant membrane with a cholesterol:POPC ratio of 3:1 [59]. These results strongly suggest that the S4 helix is not the main driver of pore and complex pore formation. This observation raises the question of which VSD properties make them more prone to pore formation. While it is possible that a lower number of arginine residues in the S4 helix may increase the susceptibility to pore or complex pore formation, more studies are required to confirm this conclusion.

One important observation is that the structural changes in the VSD in response to an external electric field are accompanied by a rearrangement of the surrounding membrane, see Figure 4. Considering that membranes can act as allosteric regulators of protein structure and function [94,95], it is hypothesized that the formation of pores and complex pores involves a cooperative rearrangement between the VSD and the membrane. This is supported by a previous study that found pores and complex pores forming in a cholesterol:POPC 3:1 membrane at an electric field strength of $\overrightarrow{E_z}$=0.2 V/n, but not in a pure POPC membrane [59].

Regarding the density map in Figure 8, cluster 3 appears to be the largest, but this is actually due to the data being more dispersed, as can be observed in Figure 7D.

The clusterization throws 3 structures, two very similar between them that present simple pores. VSD of cluster 2 present a slightly more open VSD that the VSD of cluster 1 letting the pass of more water. The VSD of cluster 3 present a more open pore that form a complex pore.

5. Conclusions

The primary determinant governing the formation of pores in the VSD due to nsPEF stimulation is not the S4 helix, as one might expect intuitively. This specific cardiac VSD, containing only three arginines in the S4 helix, is less resistant to external electric fields compared to similar studies [55,59]. $\overrightarrow{E_z}$ values above 0.14 V/nm, resulting in replicas with membranes rotated by 90°, which are not biologically significant. The final RMSD of the simulated replicas that differentiate between pores and complex pores is approximately 0.8 nm, and above this value, the complex pores are similar, as shown in Figure 6. Full VSD pores start forming at an RMSD greater than 0.35 nm.

The clusterization present 3 cluster, and the most representative VSDs of each cluster are no very different beteewn them, but the water that pass through the is very different. 2 VSDs form a simple pore, and VSD of cluster 3 present a complex pore. The thermal noise of 2.47 kJ/mol indicates that all representative structures of the VSD obtained through clusterization are oscillating in thermal equilibrium with each other

The observation of protein-mediated electroporation through the formation of pores and complex pores in VSDs by nsPEF stimulation, as demonstrated in this MD simulation study, must be considered for a better understanding of the biophysics underlying the effects of nsPEF applications, especially in the case of cardiac VSDs. This newly discovered phenomenon of protein-mediated electroporation due to a nsPEF stimulus is not well-described and should be considered an additional effect of nsPEF in cells. For more information about protein-mediated electroporation go to refer [3] The application of this phenomenon could have significant impacts in various scientific fields. We strongly support the exhortation made by Rems et al. [55] for experimentalists to conduct further investigations on this issue.