# Ionic transport through a protein nanopore: a Coarse-Grained Molecular Dynamics Study

Nathalie Basdevant, Delphine Dessaux, Rosa Ramirez

► To cite this version:

Nathalie Basdevant, Delphine Dessaux, Rosa Ramirez. Ionic transport through a protein nanopore: a Coarse-Grained Molecular Dynamics Study. Scientific Reports, 2019, 9 (1), 10.1038/s41598-019-51942-y. hal-02395220

HAL Id: hal-02395220

https://univ-evry.hal.science/hal-02395220v1

Submitted on 5 Dec 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
OPEN

# Ionic transport through a protein nanopore: a Coarse-Grained Molecular Dynamics Study

Nathalie Basdevant<sup>1,†</sup>, Delphine Dessaux<sup>1,†</sup> & Rosa Ramirez<sup>∗</sup>

The MARTINI coarse-grained (CG) force field is used to test the ability of CG models to simulate ionic transport through protein nanopores. The ionic conductivity of CG ions in solution was computed and compared with experimental results. Next, we studied the electrostatic behavior of a solvated CG lipid bilayer in salt solution under an external electric field. We showed this approach correctly describes the experimental conditions under a potential bias. Finally, we performed CG molecular dynamics simulations of the ionic transport through a protein nanopore (α-hemolysin) inserted in a lipid bilayer, under different electric fields, for 2–3 microseconds. The resulting I – V curve is qualitatively consistent with experiments, although the computed current is one order of magnitude smaller. Current saturation was observed for potential biases over ±350 mV. We also discuss the time to reach a stationary regime and the role of the protein flexibility in our CG simulations.

Single-molecule experimental techniques are widely used to understand the kinetics and energetics involved into biological phenomena. Among them, because of its simplicity and amenability to parallelization, nanopore force spectroscopy is an increasingly used technique for ultra-sensitive detection of macromolecules. The method consists in applying an electric potential difference to guide a charged molecule through a nanopore (artificial or biological) inserted in a solid or lipid membrane, in presence of a salt solution<sup>1–5</sup>. A macromolecule going through the pore, partially blocking it, creates a variation in the electrical current according to the nature of the molecule and the pore characteristics: size, conformation, and structure. Ionic current measurements are a probe for the detection. Amongst the various biological nanopore, α-hemolysin is the most widely employed in nano-pore analytics. It has been commonly used as a nanopore for DNA or RNA translocation<sup>1,3</sup>, DNA sequencing or biosensing<sup>6</sup>, DNA unzipping of individual hairpins<sup>7–9</sup>, and the study of protein translocation<sup>4</sup>.

The understanding of the physical mechanisms governing the translocation of ions or molecules through a nanopore is not only helpful to reveal relevant biological phenomena but also for the design of new devices based on such nanopores. The microscopic details of these processes cannot be deduced straightforwardly from the experiments. In this sense, an increasing effort to provide theoretical and computational approaches in order to understand such processes is being done<sup>10</sup>. In general, these approaches require a choice between the amount of computational resources and the microscopic level of accuracy. Among continuum approaches, the Poisson-Nernst Planck theory (PNP)<sup>11–16</sup> is useful to study ion transport at a relatively small computational cost, although not very accurate for narrow pores. At the same level, Brownian or Langevin dynamics methods (BD)<sup>12,17–19</sup> can be an alternative. In this case, the main challenge is the representation of the protein-continuum solvent interface.

In All-Atom Molecular Dynamics (AA-MD) simulations, the individual motion of ions, water, protein and membrane atoms are fully computed based on a classical atomic force field. This approach will be limited by the classical force field approximation used, but AA-MD simulations are undoubtedly excellent to provide microscopic details of biological translocation processes “at short times”. The huge computational cost is the main disadvantage of this method and therefore, there are several ways to accelerate the transport of molecules through the α-hemolysin pore: using very high electric potentials (1 V compared to hundreds of mV for experiments)<sup>8,20,21</sup> or forcing the transport of molecules through the pore using velocity-pulling methods<sup>22</sup>. Although some results of these AA-MD are in accordance with experiments, long time phenomena cannot be observed in such simulations, generally limited to a few hundred nanoseconds. The high-potential/pulling methods can also induce artifacts in the results.

Coarse-Grained (CG) force fields, representing several atoms by one site, for lipid membrane and biomolecules, are a very promising alternative to All-Atoms models, as they enable simulations of time-lengths closer

LAMBE, Univ Evry, CNRS, CEA, Université Paris-Saclay, 91025, Evry, France. *email: rosa.ramirez@univ-evry.fr
to experiments. In this sense, our article focuses on MD simulations where the AA representation is replaced by an “atomistic” CG force field, in our case the MARTINI force field. The CG MARTINI force field<sup>23</sup>, originally developed for lipids, is widely used to simulate membrane systems for time-lengths up to hundreds of microseconds<sup>24–26</sup>. Furthermore, a polarizable CG water model is needed to correctly describe the electrostatic behavior of electrolytes. Among the several CG polarizable solvent models available<sup>27–29</sup>, the Polarizable Water (PW) model<sup>29</sup> has been parametrized to be used together with the standard MARTINI force field. Whereas classical AA simulations of ionic current through α-hemolysin have been done, to our knowledge, such a study using a CG model has not been yet performed.

In this work, our MARTINI CG MD simulations for determining ion currents through α-hemolysin are detailed and discussed. We begin by studying the electrostatic behavior of solvent, ions and membrane using this force field. The ionic conductivity of CG MARTINI ions in PW solution as a function of concentration is discussed. Secondly, we study the electrostatic behavior of a solvated DPPC (dipalmitoylphosphatidylcholine) lipid bilayer with excess ions in solution. We will confirm that applying an electric field in a direction perpendicular to the membrane plan is accurate to represent the experimentally applied potential difference. Finally, we will study the current through the protein nanopore at various electric potentials on microseconds-long MD.

## Results

All our MD simulations were performed using the Gromacs software package<sup>30</sup> (version 5.0.2) and the MARTINI 2.2 Coarse-Grained (CG) model for ions, lipids<sup>23</sup>, proteins<sup>31,32</sup> and the PW (Polarizable Water) model<sup>29</sup> for water. The systems were simulated either in the NVT or the NPT ensemble at a temperature of T = 320 K, where PW water presents a density of 1020 kg/m<sup>3</sup><sup>29</sup>. MARTINI CG ions consist in particles of the same size and mass as a PW particle, carrying positive or negative charge. Although they are called “Na<sup>+</sup>” and “Cl<sup>-</sup>” in the force field parameters files, they do not represent a realistic sodium or chloride ion<sup>33</sup>. Since the MARTINI salt is perfectly symmetrical, we considered that it represents the KCl salt instead of NaCl. Methodological details of the simulations can be found in the Methods. The details of the ionic molar concentration setting for NVT simulations is discussed in the Supplementary information.

**Ionic conductivity.** First of all, we focused on the ionic conductivity of MARTINI ions in bulk PW water. The correct interplay between ions and water force fields is of critical importance to simulate conduction processes. In the case of MARTINI PW water, bulk properties such as density and dielectric constant have been optimized to construct the PW force field, while for ions the equilibrium bulk correlation function ION-PW water has been benchmarked to determine ions parameters<sup>29</sup>. We would like to address the ability of the combination of these ion and water force fields to yield reasonable results when used to simulate transport processes.

Our simulations consisted of a cubic box containing PW water particles, cations (K<sup>+</sup>) and anions (Cl<sup>-</sup>) with concentrations from 0.1 M to 2 M (see Methods for details). For each concentration, three different conditions were simulated: no external electric field, and with electric fields (E) of 0.02 V/nm and 0.04 V/nm in the z direction.

To calculate the ionic conductivity κ of CG ions as a function of concentration, we computed the average total ionic current <I> in the z direction as follows: the z-component of the velocity of the center of mass of all K<sup>+</sup> ions, v<sub>z</sub><sup>+</sup>, and the one of all Cl<sup>-</sup> ions, v<sub>z</sub><sup>-</sup>, were calculated every 500 ps, and the instantaneous current in the z direction was then computed using the formula:

$$I(t) = \frac{1}{L_z} N_{\text{ions}} e(v_z^+ - v_z^-), \quad (1)$$

where N<sub>ions</sub> is the total number of ions of each category, L<sub>z</sub> is the box dimension in z direction and e is the elementary charge. The cumulative sum of the current over the length of the simulation (1 µs) was then calculated and fitted with a linear regression, the resulting slope of this fit yielding the average current <I>. The conductivity κ is defined by κ = <I>/EA, where A is the cross-sectional area, which is here: A = L<sub>x</sub>L<sub>y</sub> (L<sub>x</sub> and L<sub>y</sub> are the boxes dimensions in x and y directions). Since two distinct electric fields were applied in the simulations at each concentration, we computed the average of the two calculated conductivities.

The concentration dependence of the conductivity computed at a fixed temperature of 320 K is shown on Fig. 1, as well as the experimental results and AA simulations using the CHARMM and AMBER force field performed by Pezeshki et al.<sup>34,35</sup> for KCl (numerical values from the figure of ref.<sup>34</sup> were kindly provided by Ulrich Kleinekathöfer). The CG conductivity increases with the ionic concentration for concentrations below 1.15 M and is slightly lower than the reference values. For concentrations above this value, the conductivity decreases with the concentration. This effect is due to the size of the CG particles and principally PW particles: at concentrations larger than 1.15 M, the CG model cannot correctly represent the water Coulomb screening, and non-linear effects appear dramatically. As an example, at 2.3 M, there are only twice more PW particles than CG ions (see Supplementary Information for further discussion on the concentration calculation for the ions CG model).

Using AA models, a saturation of the conductivity was also observed above 1 M, compared to experiments, but there is no decrease. A saturation effect is also measured experimentally for conductivity in the bulk for certain types of ions with a large hydration shell<sup>36</sup>.

The CG ionic conductivities values with PW at 320 K are in very good agreement with experimental and all-atom theoretical data for concentrations of KCl below 1.15 M. This result is quite different from<sup>37</sup> where the authors, using a Nernst-Einstein method, show that conductivity of CG MARTINI ions is close to experiments only if an effective concentration is used.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Ionic conductivity of KCl as a function of molarity. Experimental results (in red) and all-atom simulations using the CHARMM (blue) and AMBER (green) force field performed by Pezeshki et al.<sup>34</sup> for KCl at 320 K, and CG results (black) computed with our CG MARTINI MD simulations.</figcaption></figure>

**Charge distribution and electric potential of a DPPC bilayer.** The electric potential across lipid membranes is essential for various biological processes in the cell, such as signaling, ion transport through membrane channels, or translocation of molecules through pores. *In vivo*, such potential difference across the membrane is caused by ionic concentration gradients.

In nanopore experiments, this difference is established and monitored using two electrodes placed across the membrane¹, while there are several options to simulate the trans-membrane potential *in silico*. In this work, a uniform electric field *E* in the *z* direction perpendicular to the membrane surface is applied to every charged atom of the system³⁸. This method has been previously used by Aksimentiev and Schulten²¹ on the α-hemolysin pore, and successfully applied to various biological problems⁸,¹⁰,³⁴,³⁹. An extensive MD study on several membrane systems showed that this method accurately represents an electric potential difference of $\Delta V = EL_z$, where $L_z$ is the system size in the *z* direction, independently of the membrane thickness⁴⁰.

In the original paper on the PW model²⁹, the electric potential around a DPPC membrane was computed without any ions and external electric field. Therefore, a study of the electrostatic properties of a system composed of MARTINI {DPPC + PW + IONS} under an external electric field is needed to assure the validity of the model for describing the ionic transport. After equilibration of a system composed of DPPC lipids, PW water molecules and ions (see Methods), the MD production was performed in the NPT ensemble for two microseconds at 320 K and 1 bar. Several ionic concentrations were simulated: pure water, 0.2 M, 0.4 M and 0.8 M. Moreover, at each concentration, three different electric fields were applied: no external field, 0.01 V/nm and 0.02 V/nm in the *z* direction. These electric fields correspond to a theoretical transmembrane potential difference of 0 mV, 160 mV and 320 mV, respectively. It should be noted that, in the NPT ensemble, the box dimension in the *z* direction increased at the beginning of the dynamics, and the $L_z$ length was around 16 nm on average during all the membrane MD (from 16.24 nm without ions to 16.41 nm at 0.8 M).

The averaged charge density, on the last 500 ns of the 2µs-simulation along the *z* direction for each species, DPPC, water and ions, was computed for each simulation. Figure 2a shows the total linear charge density and its decomposition (b–d), without an electric field, as a function of the *z* coordinate, without ions and at 0.4 M and 0.8 M concentrations. The total charge density is globally similar to previous AA simulations⁴¹, showing four peaks from either side of the membrane. The decomposition of the charge density is also consistent with the AA one, so that water and ions, when present, rearrange in reaction of the membrane charge density. The DPPC charge density (Fig. 2b) is slightly different when ions are present in the system, showing that the lipids rearrange a little in the presence of salt. Compared to AA simulations⁴¹, the effect of the salt concentration on the membrane potential is weaker since it is limited to the charged head group. Moreover, the influence of the ionic concentration on the area per lipid is much smaller than for AA simulations, decreasing from 0.635 nm² with no ion to 0.617 nm² at 0.8 M, instead of 0.62 nm² with no ions to 0.5 nm² at 0.8 M. The solvent charge density is very different with and without the presence of ions (Fig. 2c), since the dipolar moments of water molecules react to the ionic charge distribution. Figure 2e,f show the particle density for each ions species ($K^+$ (e) and $Cl^-$ (f)). From these figures, we can infer that there is a layer of cations close to the negative charges of the phosphate groups, corresponding to the two peaks around $|z| = 2$ nm on Fig. 2e. There is also an overlay of anions close to the positive beads of the lipid heads, corresponding to the two peaks around $|z| = 2.5$ nm on Fig. 2f. This ionic structure on the membrane surface is also consistent with AA simulations⁴¹.

The electrostatic potential along the *z* direction for each was then calculated by the double numerical integration of the charge density, according to⁴²,⁴³. In the absence of an external electric field, the system, centered on the bilayer, is theoretically symmetrical in the transformation (*x*, *y*, *z*) to (*x*, *y*, −*z*). Therefore, for the systems without electric field, the charge density $\rho(z)$ was symmetrized before calculation of the potential. Figure 3a shows the decomposition of the electric potential without ions. When no ions are present, the PW particles are oriented to compensate the electric potential created by the lipids. However, the total electric potential inside the membrane is negative, contrary to AA⁴¹ for which the total electric potential is slightly positive inside the membrane. The PW solvent electric potential does not compensate enough the potential created by the DPPC membrane, as it was previously shown in the original PW paper²⁹.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** (a) Total linear charge density along the *z* axis around a CG DPPC bilayer with (0.4 M in red and 0.8 M in green) and without ions (in black), and its decomposition: (b) DPPC charge density, (c) solvent, (d) ions. (e–f): Particle density for K⁺ (e) and Cl⁻ (f) ions. The *z* axis is centered at the middle of the lipid bilayer.</figcaption></figure>

On Fig. 3b, the decomposition of the electric potential at 0.8 M ionic concentration is shown. The total electric potential inside the membrane is less negative than without ions. The potential created by the membrane increases with the ion concentration, as in AA simulations⁴¹, proving that the bilayer is sensitive to the ions presence. When ions are present, the solvent potential is much lower and the rearrangement of ions is compensating the lipids potential. This effect is consistent with the AA study by Rodriguez and Garcia⁴¹, that showed that in the presence of ions, water plays a significantly smaller contribution to the total electric potential. Therefore, although the total electric potential inside the membrane is negative instead of positive, CG ions and water adopt the same behavior around a lipid membrane as an AA force-field.

Figure 3c shows the same decomposition for simulations with an external electric field. The total potential thus also contains the applied external potential (in dashed line). The figure shows that the ionic polarization potential is strongly altered by the external potential. This is not the case for other species whose charge density distribution is much less modified by the external field. On Fig. 3d, the total electrostatic potentials for simulations with pure water and with 0.8 M ions are represented. In the bulk, the reaction polarization potential is opposed to the external imposed one, and the resulting total electric potential difference is therefore only effective over the membrane. This result is consistent with previous AA-MD⁴⁰ studies which attest that the electric potential does not depend on the thickness of the membrane. The authors showed that the reaction potential of the solution, composed of solvent and ions, to the external field, combined to the external electric field, results in an electric potential difference across the membrane only, with a zero-potential difference within the two bulk regions, proving that using the relation $\Delta V = EL_z$ is accurate. Our results show that this assumption is also accurate with the MARTINI CG model.

We showed that the MARTINI force field with the PW solvent reproduces accurately the ionic conductivity in the bulk and the electrostatics properties around a lipid membrane in the presence of an external electric field. We can therefore now try to simulate ionic transport through a nanopore using this methodology.

**Ionic current through the nanopore.** In this section, we present the results of the simulations for the system composed by the α-hemolysin (PDB entry: 7AHL⁴⁴) inserted in a DPPC bilayer and surrounded by PW water, ions and counter-ions in the presence of an external electric field. A value of 1.1 M ± 0.1 M ionic concentration was measured at the bulk region. Further details of these procedures and the parameters used can be found in Methods.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** (a–c) Total electric potential along the *z* axis of a CG DPPC membrane (black) and its decomposition, without ions (a), with 0.8 M ionic concentration (b), with no external electric field, and in the presence of an external electric field of 0.02 V/nm (c). (d) Total electrostatic potentials with no electric field (black), electric fields of 0.01 V/nm (red) and 0.02 V/nm (green) without ions (solid line) and with 0.4 M and 0.8 M concentrations (dashed lines).</figcaption></figure>

The initial condition for the ensemble α-hemolysin suspended in a lipid bilayer membrane was taken from final state of Prof. Sansom simulations⁴⁵, where the lipid bilayer grows spontaneously around the protein. We observe a curvature of the membrane in the pore vicinity for our system (see Fig. 4). This distortion could already be measured in the original membrane arrangement⁴⁵. Moreover, this kind of membrane curvature was also observed in AA simulations of α-hemolysin in other lipids bilayers and it appears to be characteristic of the lipids involved and of the way they interact with the protein⁴⁶.

MD simulations of α-hemolysin in DMPC lipid membrane using MARTINI force field were recently studied⁴⁷ and showed a considerable distortion of the protein stem. In our simulations, the ELNeDyn network⁴⁸ was applied inside and between the different protomers, and, as a result, the protein stem remains stable during the whole simulation.

After equilibration of the system, production simulations were carried out in the NPT ensemble ($T = 320$ K and $P = 1$ atm), with an external electric field applied perpendicular to the lipid layer, ranging from −0.04 V/nm to +0.04 V/nm. The potential bias $\Delta V$ value for a given electric field is calculated using an average $L_z = 18.5$ nm, $\Delta V = EL_z$. This $L_z$ value varies less than 5% during the simulation time. All simulations, including the no-field simulation, start from the same initial condition (see Methods). The complete list of simulations performed with their corresponding electric fields and simulation times are shown in Table 1.

To provide reasonable statistics, we performed simulations of 2 and 3 µs, depending on the applied external electric field. We do not present data for an electric field $E = -0.04$ V/nm, since the membrane broke early during the simulation. This is also the case in experiments where the allowed potential differences are $|\Delta V| < 350$ mV.

The instant current was calculated computing the number of cations and anions ($n_c(t)$ and $n_a(t)$) crossing a *xy* plane perpendicular to the *z* axis, with *z* abscissa equal to the position of the center of mass of the groups in the *trans* end of the pore, between *t* and $t + \Delta t$ ($\Delta t = 100$ ps) (see Supplementary Information). The sign represents the direction of the velocity in the *z* axis. These values are normally 0 or 1, so we define the number of ions passing through the plane in the interval $[0, t]$ as $N_{c,a}(t) = \int_0^t n_{c,a}(t')dt'$, and the total charge going through the pore, which is the cumulative function of the current, is: $Q(t) = e(N_c(t) - N_a(t))$. The instant current is defined then as the slope of the function $Q(t): I(t) = \frac{dQ(t)}{dt}$.

The number of ions crossing the pore during the whole simulation time, $N_a(t_s)$ and $N_c(t_s)$ are represented in Table 1, together with the total charge variation during the whole simulation $Q(t_s)$ in *e* units. These values are
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Representation of simulation setup. Water has been removed for clearness purposes. Section of the pore (orange) along the z axis. The DPPC membrane (blue) is curved. K⁺ (green) and Cl⁻ (magenta) ions are going through the pore channel.</figcaption></figure>

<table><caption>**Table 1.** Simulation parameters and number of ion passages during simulations. Electric field E, computed potential bias ΔV, total simulation time t<sub>s</sub>, number of cations and anions traversing during the total simulation time N<sub>c</sub>(t<sub>s</sub>), N<sub>a</sub>(t<sub>s</sub>) and total charge traversing during the simulation time Q(t<sub>s</sub>).</caption><thead><tr><th rowspan="2">E (V/nm)</th><th colspan="5">Positive Bias</th><th colspan="5">Negative Bias</th></tr><tr><th>ΔV(mV)</th><th>t<sub>i</sub>(μs)</th><th>N<sub>c</sub>(t<sub>i</sub>)</th><th>N<sub>a</sub>(t<sub>i</sub>)</th><th>Q<sub>in</sub>(e)</th><th>ΔV(mV)</th><th>t<sub>i</sub>(μs)</th><th>N<sub>c</sub>(t<sub>i</sub>)</th><th>N<sub>a</sub>(t<sub>i</sub>)</th><th>Q<sub>in</sub>(e)</th></tr></thead><tbody><tr><td>±0.00</td><td>0</td><td>2</td><td>-4</td><td>6</td><td>-10</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>±0.005</td><td>92</td><td>3</td><td>45</td><td>-108</td><td>153</td><td>-92</td><td>3</td><td>-21</td><td>103</td><td>-124</td></tr><tr><td>±0.01</td><td>185</td><td>3</td><td>85</td><td>-180</td><td>265</td><td>-185</td><td>3</td><td>-33</td><td>209</td><td>-242</td></tr><tr><td>±0.015</td><td>277</td><td>3</td><td>118</td><td>-212</td><td>330</td><td>-277</td><td>3</td><td>-31</td><td>269</td><td>-300</td></tr><tr><td>±0.02</td><td>370</td><td>2</td><td>144</td><td>-228</td><td>372</td><td>-370</td><td>2</td><td>-28</td><td>218</td><td>-246</td></tr><tr><td>±0.025</td><td>462</td><td>2</td><td>80</td><td>-130</td><td>210</td><td>-463</td><td>2</td><td>-26</td><td>178</td><td>-204</td></tr><tr><td>±0.03</td><td>555</td><td>2</td><td>109</td><td>-113</td><td>222</td><td>-555</td><td>2</td><td>-32</td><td>194</td><td>-226</td></tr><tr><td>±0.04</td><td>740</td><td>2</td><td>241</td><td>229</td><td>470</td><td>-740</td><td>X</td><td>X</td><td>X</td><td>X</td></tr></tbody></table>

in the limit of “reasonable” statistics although, as we will see below, they evolve with time. The rectification and selectivity of the pore can be also deduced from this global values. The instant cumulative current Q(t) is represented in Fig. 5 for potentials under (side a) and over (side b) 300 mV. Since the instant current I(t) is the slope of this function at each time t, the linear increase of the cumulative current with time would indicate a stationary current. This seems to be the case for potential bias under 300 mV, although it should be pointed out that the stationary current is established only after some hundreds of nanoseconds, especially in the case of positive potential bias. The fact that the slope of the function increases with increasing potential is the reflect of a linear-response regime.

This is not the case for potentials over 300 mV, as shown in Fig. 5b). The cumulative currents are not proportional to the potential bias at the end of the simulation, implying that the simulated systems are not in the linear response regime. Moreover, we cannot be confident that the stationary state has been reached at the end of some simulations as the slope of the cumulative current is still varying.

Figure 6 shows the evolution in time of the I-V curve of the system averaged at a fix number of crossings (N=40 events). For every simulation, an initial time t<sub>ini</sub> is defined as Q(t<sub>ini</sub>) = 5 e, which is the time when five charges have already crossed the pore. From this origin t<sub>ini</sub>, we computed the total time t<sub>end</sub> such as ΔQ = Q(t<sub>end</sub>) - Q(t<sub>ini</sub>) = 40 e. In this interval, the current Δt = t<sub>end</sub> - t<sub>ini</sub>, is estimated as $\bar{I} = \frac{\Delta Q}{\Delta t}$. The error bars were calculated with the standard deviation $\sigma = \frac{1}{N} \sum_{i=1..N} |\Delta Q_i - I\Delta t_i|$ so that $\delta I = \frac{2\sigma}{\Delta t}$.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Cumulative current as a function of time, for potentials (a) under or (b) over 300 mV computed during 3 or 2 µs, respectively. (a): ±92 mV (green), ±185 mV (red) and ±277 mV (blue). (b): ±370 mV (green), ±462 mV (red), ±555 mV (blue) and ±740 mV (violet).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Average current during different passages windows as a function of potential bias. Note the difference between the beginning of the simulation (5–45 crossings) and the end of the simulation (195–235 crossings). The black solid line corresponds to $I^{\text{ref}}/10$.</figcaption></figure>

In order to compare our results to experimental values, we have defined a reference current in the linear regime, as $I^{\text{ref}} = R^{\text{ref}}\Delta V$, where $R^{\text{ref}} = 1 \text{ G}\Omega$ is a typical resistance observed experimentally for these systems³⁶. For negative electric fields, $I^{\text{ref}}$ value is multiplied by 0.7 to take into account the observed rectification of 30% in experimental setups. The black continuous line observed in the figure is the function $I^{\text{ref}}/10$. Thus, in our simulations, we observe a $Q_s$ value ten times smaller than the expected one. This is a consequence of the CG model: the resistance of the pore increases since the mobility of the ions inside the channel is highly reduced because of the larger size of PW water molecules. This ratio increases for higher potential differences where the linear regime is not valid anymore.

The measured current is clearly evolving during our simulations. In spite of this, a linear response regime can be established for potential biases below ~350 mV. Details of currents for potential biases between –370 and +370 mV and their anionic and cationic contributions are detailed in Table 2. In this regime, the final current appears slightly smaller than the one measured at the beginning of the simulation and, as mentioned, it is one order of magnitude smaller than the experimental one. For potentials above this value, a limiting current is observed. We find some asymmetry between positive and negative biases: while, for negative biases, the limiting current seems to converge to a well defined value, for the positive ones, this value appears much more noisy.

The origin of the non-linear ohmic behavior of the current-voltage curve has been explained as due to the crowding or depletion regions at the *trans* and *cis* ends of the pore in presence of an applied external potential⁴⁹. Both limiting and over-limiting currents in nanopores have been discussed in van Oeffelen *et al.*⁵⁰. In our simulations, we observe a limiting current at a quite low voltage ~350 mV. This value is probably higher in AA simulations, although it is difficult to observe as the simulated time is not enough to get to the steady state²¹. In ref.⁵¹, the authors compare linear-response theory with AA-MD simulations to determine the departure from the linear regime for two pore models, one of them being a modified α-hemolysin pore. Their conclusion, and also ours as regard of our results, is that α-hemolysin has a small linear response range and that currents at low voltages
<table><caption>**Table 2.** Currents in the linear response regime (±400 mV). The total current values ($I_{tot}$) are given for 195–235 crossings and 5–45 crossings between parenthesis. NC stands for systems where we cannot get to 195–235 crossings. The current presented for anions (Cl⁻) and cations (K⁺) are computed on average during the whole simulation. The same applies for the total current with no bias.</caption><thead><tr><th>V (mV)</th><th>370</th><th>277.5</th><th>185</th><th>92.5</th><th>0.0</th><th>-92.5</th><th>-185</th><th>-277.5</th><th>-370</th></tr></thead><tbody><tr><th>I<sub>tot</sub> (pA)</th><td>33.6 (30.6)</td><td>24.6 (31.2)</td><td>14.8 (15.9)</td><td>NC (8.9)</td><td>-0.8</td><td>NC (-5.2)</td><td>-12.5 (-13.6)</td><td>-15.9 (-19.2)</td><td>-17.8 (-20.2)</td></tr><tr><th>I<sub>K+</sub> (pA)</th><td>11.5</td><td>6.3</td><td>4.5</td><td>2.4</td><td>-0.3</td><td>-1.1</td><td>-1.8</td><td>-1.7</td><td>-2.2</td></tr><tr><th>I<sub>Cl-</sub> (pA)</th><td>-18.2</td><td>-11.3</td><td>-9.6</td><td>-5.8</td><td>0.5</td><td>5.5</td><td>11.1</td><td>14.3</td><td>17.4</td></tr></tbody></table>

cannot be safely extrapolated from the values obtained at higher potential biases. The limiting current cannot be experimentally observed since, as mentioned, the membrane breaks over ~300 mV.

It is not clear from our results if an over-limiting current appears over +700 mV⁵⁰. This could be the case but it cannot be excluded that the stationary state has not been reached, that the effective protein charge (protein plus decorating ions) evolved, or that the periodic boundary conditions are inducing such phenomenon.

As observed, the average current cannot be accurately computed before reaching the stationary state and therefore, an interesting issue is which simulation time is needed to get to this regime. In van Oeffelen et al.⁵⁰, it is pointed out that, due to slow co-diffusion of ions, the time for the charge-current to decay to a stationary state value is determined by the diffusion of particles. This time is estimated as a typical diffusion time $\tau = L^2/D$, where L is a typical distance and $D = 2(D_+D_-)/(D_+ + D_-)$ is an averaged diffusion coefficient for cations and anions. In our simulations, this time is estimated to ~140 ns, where we have used L the length of the system occupied by the solvent, $D_+ = 2.54 \times 10^{-9} \text{ m}^2/\text{s}$ and $D_- = 2.45 \times 10^{-9} \text{ m}^2/\text{s}$ measured during our conductivity calculations. We have computed the evolution of the number of ions inside the stem, the region of the pore going from the trans side to 4 nm above (see Supplementary Information). We observed that the time to get to a stationary number of ions goes from ~100 to 300 ns, being shorter for negative voltage biases. The asymmetry for positive and negative biases is in agreement with AA simulations, as well as the number of ions in the stem computed in previous AA and PNP studies¹²,⁵².

The time to get to a stationary regime is longer than this estimation according to our simulations. The relatively quick convergence of the number of ions inside the pore is not a guarantee for the stationary state, where anions and cations currents, or conversely charge-current and particle-current, should be, in average, time independent. This is not the case for all our simulations. We observe the convergence of the current for voltage biases under +350 mV as well as for negative higher biases. Positive high biases would require longer simulations to get confident results. We can conclude that a 1–1.5 µs simulation should be enough to get a stationary state in the linear regime.

A final pertinent question about CG-MD simulations of such complex systems is if the pore dynamics play a determining role on the current. If this was not the case, simulation time could be saved by tabulating the pore and membrane potential. To test this possibility, we have performed CG-MD simulations where the backbone of the α-hemolysin was constrained to a fixed configuration. The external electric fields imposed were ±0.01 V/nm, ±0.02 V/nm and ±0.03 V/nm. The measured current in this case is much smaller than the corresponding one for the flexible protein and we found clogging periods during the simulation, going up to 200 ns. The current saturates for small potential biases going to a value close to zero (see the Supplementary Information). This result proves that the flexibility of the protein is essential to simulate the translocation process.

## Discussion

The use of MARTINI force field to simulate the transport of ions through the α-hemolysin in a DPPC bilayer is studied for the first time, to our knowledge, in this article.

In first place, we have investigated the transport properties of the original Martini hydrated ions as a function of concentration. This leads us to the problem of comparing CG results to AA or experimental ones. Since the ratio ion/water valid for AA systems is no longer a suitable parameter to set a molar concentration for the CG simulations, we have decided to use the number of ions per volume as concentration. We are aware that this definition may not represent exactly the same concentration in NVT and NPT ensembles, and definitely cannot be directly compared to AA results. In spite of this incertitude, we can conclude that MARTINI ions in PW water present a conductivity that could represent a KCl solution below 1.1 M concentration.

Secondly, we have investigated the properties of the {DPPC + IONS + PW} system in the absence and under different applied electric fields. It is already known that the MARTINI force field leads to a unrealistic positive potential inside the membrane in the absence of an applied electric field²⁹. In the presence of an external field, we found that the expected potential differences through the membrane are recovered, such as in previous AA simulations⁴⁰. We have also pointed out the role of the polarization of water to describe the correct electrostatics, proving that it is necessary to use a polarizable model of water.

Finally, we have studied the whole system {α-hemolysin + DPPC + PW + IONS} under different applied external fields focusing on the ionic transport through the pore. Principally due to the size of MARTINI PW water, the pore presents a smaller conductance than the experimental system. We have also performed simulations with a fixed protein which shows that the channel flexibility plays also an important role on the reduced CG conductance since we observed almost no current. This deviation in conductance seems to be systematic for low applied external fields. The qualitative behavior of the system nevertheless appears to be correct. Our I–V curves show a small range ±350 mV where the linear response regime is established. For higher potential biases,
a saturation of the current is observed. These limits cannot be compared to experiments, for practical reasons, or to AA-MD simulations, since the simulated time is not enough as to get to a stationary state.

PNP solutions have provided to be extremely useful to describe the stationary state of ionic currents in many systems, but a precise description of this particular system would require excluded volume effects, water-ion correlations, hydrodynamics, finite system size correlations and the possibility to take into account the flexibility of the protein. On the other side, AA-MD simulations provide good current statistics at reasonable simulation times for high voltages but a) reaching the steady state seems to take a long time, if it is possible due to the finite system size, and b) the system is at a non-linear regime. Although the non-linearity will appear at higher potentials in AA simulations than in CG simulations, in ref.⁵¹, they estimate this limit to be below 500 mV in a AA-MD study of a modified α-hemolysin, interpolating low voltage values from high voltage ones can be somehow risky. At low voltages, direct sampling of net ion fluxes with AA simulations would require very long simulation times to acquire sufficient statistics.

Our conclusion, as regard to our results, is that CG-MD simulations could result a good alternative to study the details, such as mutation effects, of translocation through protein nanopores, saving some of the simulation time required for the equivalent AA simulations for these kind of problems.

## Methods

**Simulation details.** We performed Molecular Dynamics (MD) using the Gromacs software package³⁰ (version 5.0.2) and the MARTINI 2.2 CG model for ions, lipids²³, proteins³¹,³² and the PW (Polarizable Water) model²⁹ for water. Conductivity simulations systems contained only water and ions at different concentrations, whereas membrane simulations systems contained a DPPC bilayer, water, and ions at different concentrations. In the nanopore simulations, the protein nanopore was inserted in the lipid membrane.

A timestep of 20 fs was used for all MD simulations, which is the recommended timestep for CG MD simulations using MARTINI force field. Periodic boundary conditions in the three dimensions were applied. All MD were performed using the PME method for electrostatic interactions, with a 2 Å spacing for the Fourier grid and a direct space cut-off radius of 13 Å. A relative dielectric constant of $ε_r = 2.5$ was used for PW as recommended²⁹. Van der Waals interactions were shifted at 13 Å. A Nosé-Hoover thermostat⁵³,⁵⁴ was used to maintain the temperature at 320 K.

**Conductivity simulations.** For conductivity simulations, the system was originally built as a cubic box of 15.8 nm³, containing 32,756 PW water particles. The system was minimized and slightly equilibrated (2 ns at NVT) before the introduction of ions at different concentrations. The *genion* tool of Gromacs was used to replace random water particles by positive (the MARTINI $Na^+$ ion) and negative ($Cl^-$) ions in the same quantity, in order to keep the system neutral. The system containing water and ions was then again minimized and equilibrated with a short MD (10 ns) in the NVT ensemble, using the Berendsen thermostat. Then, a one-microsecond-long MD was performed in the NVT ensemble at 320 K, using the Nosé-Hoover thermostat, with a characteristic time interval $τ_T = 1$ ps for both groups (ions and water). For each concentration, ranging from 0.1 M to 2 M, three different simulations were performed: one with no external electric field, one with an electric field of 0.02 V/nm in the *z* direction, and one with an electric field of 0.04 V/nm in the same direction.

**Membrane simulations.** As for the membrane simulations, a box containing 128 DPPC CG lipids organized in a bilayer was downloaded from the MARTINI website⁵⁵. 4,232 water molecules (PW) were added around the membrane in a box of originally 6.7 × 6.1 × 12.0 nm³, the *z* axis oriented perpendicularly to the membrane plan. The system was minimized and shortly equilibrated before addition of ions at different concentrations using the same procedure as the conductivity simulations. Simulations without ions were also performed using the same methodology. After another minimization of the system, a NVT MD simulation was performed for 2 ns using the Berendsen thermostat, followed by 120 ns of NPT MD using the Berendsen thermostat and Berendsen semi-isotropic barostat at 1 bar. The MD production was performed for two microseconds using the Nosé-Hoover thermostat at 320 K (with a characteristic time interval $τ_T = 1$ ps for each group) and the Parrinello-Rahman barostat to keep the pressure at 1 bar (with a characteristic time interval of $τ_P = 5$ ps and a compressibility of $4.5 \times 10^{-5}$ bar⁻¹). It should be noted that, in the NPT ensemble, the box dimension in the *z* direction increased at the beginning of the dynamics, and the $L_z$ length was around 16 nm on average during all the membrane MD. In addition to simulations without ions, three different KCl concentrations were simulated: 0.2 M (60 ions of each species), 0.4 M (124 ions) and 0.8 M (248 ions). For each concentration, three different simulations were performed: one with no external electric field, one with an electric field of 0.01 V/nm in the *z* direction, and one with an electric field of 0.02 V/nm. The averaged charge density was calculated on the last 500 ns of the 2 µs-simulation, using the *gahl potential* tool of Gromacs, decomposing the system into 100 slices in the *z* direction.

**Nanopore simulations.** The crystallographic structure of α-hemolysin was taken from the PDB (entry 7AHL⁴⁰). The coordinates of the missing atoms were added using the *pdb2pqr* software and the crystallographic water molecules were removed. The protein was reduced in coarse grains using the *martinize.py* script. The ElneDyn MARTINI force field was applied on the protein⁴⁸. It is composed of an elastic network model that maintains the secondary structures during the dynamics using springs (500 kJ.mol⁻¹.nm⁻² elastic bond strength and a 0.9 nm cutoff). After a short energy minimization, the resulting CG protein structure was inserted in a DPPC bilayer using the following procedure: the PDB structure of a CG α-hemolysin inserted in a DPPC bilayer of 756 lipid molecules obtained by Prof. Sansom⁴⁵ was downloaded from CGDB (coarse grain database, a database of inserted membrane proteins inside DPPC bilayers⁵⁶). Since the protein CG model used for CGDB is different from the
MARTINI one, a least-squared fit of the backbone of our CG protein on the CGDB one was performed, then the initial protein from CGDB was deleted to keep only our MARTINI CG protein inserted in the DPPC bilayer. The size of the simulation box was initially set at 15.85 × 15.85 × 20 nm. The system was solvated with 28,758 standard W MARTINI CG water molecules using the *genbox* Gromacs command. Another minimization was carried out followed by two equilibration steps: 1 ns in the NVT ensemble (using a 10 fs timestep) followed by 5 ns in the NPT ensemble (using a 20 fs timestep).

The ionic concentration was set around 1 M, corresponding to 2,054 ions of each species (calculated according to the number of initial water molecules and using the ions/PW molecules ratio from conductivity simulations concentrations). Under normal conditions, the nanopore carries a global charge of +7e. Therefore, in order to neutralize the system, 2,050 K⁺ ions and 2,057 Cl⁻ ions were added to the system. A short dynamics of 5 ns was then performed at NPT. After transformation of MARTINI W water into polarizable PW water, a short minimization was carried out, followed by two additional equilibration dynamics of 1 ns at NVT and 10 ns at NPT, both using position restraints on the protein backbone (with a 20 fs integration timestep). A last dynamics of 1 ns at NPT without position restraints was performed using the PME method for electrostatic interactions. The resulting system was used as initial state for all the simulations of the pore as well as the dynamics with the pore under position restraints on the protein backbone.

Several electric fields (E) were applied in the z direction of the system, perpendicular to the membrane plane, both positive and negative: no electric field, +/-0.005 V/nm, +/-0.01 V/nm, +/-0.015 V/nm, +/-0.02 V/nm, +/-0.025 V/nm, +/-0.03 V/nm and +/-0.04 V/nm. Molecular dynamics were carried out for 2 to 3 µs (+/-0.005 V/nm, +/-0.01 V/nm, +/-0.015 V/nm). See Table 1 for further details.

Received: 7 November 2018; Accepted: 17 July 2019;

Published online: 31 October 2019

## References

1. Kasianowicz, J. J., Brandin, E., Branton, D. & Deamer, D. W. Characterization of individual polynucleotide molecules using a membrane channel. *Proc. Natl. Acad. Sci.* **93**, 13770-13773, https://doi.org/10.1073/pnas.93.24.13770 (1996).

2. Dudko, O. K., Mathé, J. & Meller, A. Nanopore force spectroscopy tools for analyzing single biomolecular complexes. *Methods Enzymol.* **475**, 565-589, https://doi.org/10.1016/S0076-6879(10)75021-7 (2010).

3. Meller, A., Nivon, L., Brandin, E., Golovchenko, J. & Branton, D. Rapid nanopore discrimination between single polynucleotide molecules. *Proc. Natl. Acad. Sci.* **97**, 1079-1084, https://doi.org/10.1073/pnas.97.3.1079 (2000).

4. Oukhaled, G. et al. Unfolding of Proteins and Long Transient Conformations Detected by Single Nanopore Recording. *Phys. Rev. Lett.* **98**, 158101, https://doi.org/10.1103/PhysRevLett.98.158101 (2007).

5. Oukhaled, G., Bacri, L., Mathé, J., Pelta, J. & Auvray, L. Effect of screening on the transport of polyelectrolytes through nanopores. *Europhys. Lett.* **82**, 48003, https://doi.org/10.1209/0295-5075/82/48003 (2008).

6. Deamer, D., Akeson, M. & Branton, D. Three decades of nanopore sequencing. *Nat. Biotechnol.* **34**, 518-524, https://doi.org/10.1038/nbt.3423 (2016).

7. Mathé, J., Visram, H., Viasnoff, V., Rabin, Y. & Meller, A. Nanopore Unzipping of Individual DNA Hairpin Molecules. *Biophys. J.* **87**, 3205-3212, https://doi.org/10.1529/biophysj.104.047274 (2004).

8. Mathé, J., Aksimentiev, A., Nelson, D. R., Schulten, K. & Meller, A. Orientation discrimination of single-stranded DNA inside the α-hemolysin membrane channel. *Proc. Natl. Acad. Sci.* **102**, 12377-82, https://doi.org/10.1073/pnas.0502947102 (2005).

9. Muzard, J., Martinho, M., Mathé, J., Bockelmann, U. & Viasnoff, V. DNA Translocation and Unzipping through a Nanopore: Some Geometrical Effects. *Biophys. J.* **98**, 2170-2178, https://doi.org/10.1016/j.bpj.2010.01.041 (2010).

10. Modi, N., Winterhalter, M. & Kleinekathöfer, U. Computational modeling of ion transport through nanopores. *Nanoscale* **4**, 6166-6180, https://doi.org/10.1039/c2nr31024d (2012).

11. Roux, B., Allen, T., Bernéche, S. & Im, W. Theoretical and computational models of biological ion channels. *Quaterly Rev. Biophys.* **37**, 15-103, https://doi.org/10.1017/S0033583504003968 (2004).

12. Noskov, S. Y., Im, W. & Roux, B. Ion permeation through the α-hemolysin channel: theoretical studies based on Brownian dynamics and Poisson-Nernst-Planck electrodiffusion theory. *Biophys. J.* **87**, 2299-2309, https://doi.org/10.1529/biophysj.104.044008 (2004).

13. Cervera, J., Schiedt, B. & Ramírez, P. A Poisson/Nernst-Planck model for ionic transport through synthetic conical nanopores. *Europhys. Lett.* **71**, 35-41, https://doi.org/10.1209/epl/12005-10054-x (2005).

14. Constantin, D. M. C. & Siwy, Z. S. Poisson-Nernst-Planck model of ion current rectification through a nanofluidic diode. *Phys. Rev. E* **76**, 041202, https://doi.org/10.1103/PhysRevE.76.041202 (2007).

15. Cozmuta, I., O'Keeffe, J. T., Bose, D. & Stolc, V. Hybrid md-nernst planck model of α-hemolysin conductance properties. *Mol. Simul.* **31**, 79-93, https://doi.org/10.1080/0892702041233108467 (2005).

16. Birlea, M. & Birlea, S. I. The current-voltage relation of a pore and its asymptotic behavior in a Nernst-Planck model. *J. Electr. Bioimpedance* **3**, 36-41, https://doi.org/10.5617/jeb.296 (2012).

17. De Biase, P. M., Markosyan, S. & Noskov, S. Y. Microsecond simulations of DNA and ion transport in nanopores with novel ion-ion and ion-nucleotides effective potentials. *J. Comput. Chem.* **35**, 711-721, https://doi.org/10.1002/jcc.23544 (2014).

18. De Biase, P. M. et al. What controls open-pore and residual currents in the first sensing zone of α-hemolysin nanopore? Combined experimental and theoretical study. *Nanoscale* **8**, 11571-11579, https://doi.org/10.1039/c6nr00164e (2016).

19. Millar, C. et al. Brownian simulation of charge transport in α-haemolysin. *J. Comput. Electron.* **7**, 28-33, https://doi.org/10.1007/s10825-008-0230-6 (2008).

20. Aksimentiev, A., Heng, J. B., Timp, G. & Schulten, K. Microscopic Kinetics of DNA Translocation through synthetic nanopores. *Biophys. J.* **87**, 2086-2097, https://doi.org/10.1529/biophysj.104.042960 (2004).

21. Aksimentiev, A. & Schulten, K. Imaging α-hemolysin with molecular dynamics: ionic conductance, osmotic permeability, and the electrostatic potential map. *Biophys. J.* **88**, 3745-3761, https://doi.org/10.1529/biophysj.104.058727 (2005).

22. Di Marino, D., Bonome, E. L., Tramontano, A. & Chinappi, M. All-Atom Molecular Dynamics Simulation of Protein Translocation through an α-Hemolysin Nanopore. *J. Phys. Chem. Lett.* **6**, 2963-2668, https://doi.org/10.1021/acs.jpclett.5b01077 (2015).

23. Marrink, S. J., Risselada, H. J., Yefimov, S., Tieleman, D. P. & de Vries, A. H. The MARTINI force field: coarse grained model for biomolecular simulations. *J. Phys. Chem. B* **111**, 7812-7824, https://doi.org/10.1021/jp071097f (2007).

24. Ingólfsson, H. I. et al. Lipid organization of the plasma membrane. *J. Am. Chem. Soc.* **136**, 14554-14559, https://doi.org/10.1021/ja507832e (2014).

25. Pinot, M. et al. Polyunsaturated phospholipids facilitate membrane deformation and fission by endocytic proteins. *Science* **345**, 693-697, https://doi.org/10.1126/science.1255288 (2014).
26. van Eerden, F. J., de Jong, D. H., de Vries, A. H., Wassenaar, T. A. & Marrink, S. J. Characterization of thylakoid lipid membranes from cyanobacteria and higher plants by molecular dynamics simulations. *BBA-Biomembranes* **1848**, 1319–1330, https://doi.org/10.1016/j.bbamem.2015.02.025 (2015).

27. Ha-Duong, T., Basdevant, N. & Borgis, D. A polarizable coarse-grained water model for coarse-grained proteins simulations. *Chem. Phys. Lett.* **468**, 79–82, https://doi.org/10.1016/j.cplett.2008.11.092 (2009).

28. Wu, Z., Cui, Q. & Yethiraj, A. A new coarse-grained model for water: the importance of electrostatic interactions. *J. Phys. Chem. B* **114**, 10524–10529, https://doi.org/10.1021/jp1019763 (2010).

29. Yesylevskyy, S. O., Schäfer, L. V., Sengupta, D. & Marrink, S. J. Polarizable water model for the coarse-grained MARTINI force field. *PLoS Comput. Biol.* **6**, e1000810, https://doi.org/10.1371/journal.pcbi.1000810 (2010).

30. Van Der Spoel, D. *et al*. GROMACS: Fast, flexible, and free. *J. Comp. Chem.* **26**, 1701–1718, https://doi.org/10.1002/jcc.20291 (2005).

31. Monticelli, L. *et al*. The martini coarse-grained force field: Extension to proteins. *J. Chem. Theory Comput.* **4**, 819–834, https://doi.org/10.1021/ct700324x (2008).

32. de Jong, D. H. *et al*. Improved parameters for the martini coarse-grained protein force field. *J. Chem. Theory Comput.* **9**, 687–697, https://doi.org/10.1021/ct300646g (2013).

33. Marrink, S. J., de Vries, A. H. & Mark, A. E. Coarse grained model for semiquantitative lipid simulations. *J. Phys. Chem. B* **108**, 750–760, https://doi.org/10.1021/jp036508g (2004).

34. Pezeshki, S., Chimerel, C., Bessonov, A. N., Winterhalter, M. & Kleinekathöfer, U. Understanding ion conductance on a molecular level: an all-atom modeling of the bacterial porin OmpF. *Biophys. J.* **97**, 1898–1906, https://doi.org/10.1016/j.bpj.2009.07.018 (2009).

35. Lu, J., Modi, N. & Kleinekathöfer, U. Simulation of ion transport through an N-acetylneuraminic acid-inducible membrane channel: from understanding to engineering. *J. Phys. Chem. B* **117**, 15966–15975, https://doi.org/10.1021/jp408495v (2013).

36. Bhattacharya, S. *et al*. Rectification of the current in α-hemolysin pore depends on the cation type: the alkali series probed by MD simulations and experiments. *J. Phys. Chem. C* **115**, 4255–4264, https://doi.org/10.1021/jp111441p (2011).

37. Vögele, M., Holm, C. & Smiatek, J. Properties of the polarizable MARTINI water model: A comparative study for aqueous electrolyte solutions. *J. Mol. Liq.* **212**, 103–110, https://doi.org/10.1016/j.molliq.2015.08.062 (2015).

38. Roux, B. The membrane potential and its representation by a constant electric field in computer simulations. *Biophys. J.* **95**, 4205–4216, https://doi.org/10.1529/biophysj.108.136499 (2008).

39. Böckmann, R. A., de Groot, B. L., Kakorin, S., Neumann, E. & Grubmüller, H. Kinetics, statistics, and energetics of lipid membrane electroporation studied by molecular dynamics simulations. *Biophys. J.* **95**, 1837–1850, https://doi.org/10.1529/biophysj.108.129437 (2008).

40. Gumbart, J., Khalili-Araghi, F., Sotomayor, M. & Roux, B. Constant electric field simulations of the membrane potential illustrated with simple systems. *Biochim. Biophys. Acta* **1818**, 294–302, https://doi.org/10.1016/j.bbamem.2011.09.030 (2012).

41. Rodriguez, J. R. & García, A. E. Concentration dependence of NaCl ion distributions around DPPC lipid bilayers. *Interdiscip. Sci.* **3**, 272–282, https://doi.org/10.1007/s12539-011-0107-4 (2011).

42. Böckmann, R. A., Hac, A., Heimburg, T. & Grubmüller, H. Effect of sodium chloride on a lipid bilayer. *Biophys. J.* **85**, 1647–1655, https://doi.org/10.1016/S0006-3495(03)74594-9 (2003).

43. Lee, S.-J., Song, Y. & Baker, N. A. Molecular dynamics simulations of asymmetric NaCl and KCl solutions separated by phosphatidylcholine bilayers: Potential drops and structural changes induced by strong Na⁺-lipid interactions and finite size effects. *Biophys. J.* **94**, 3565–3676, https://doi.org/10.1529/biophysj.107.116335 (2008).

44. Song, L. *et al*. Structure of staphylococcal α-hemolysin, a heptameric transmembrane pore. *Science* **274**, 1859–1865, https://doi.org/10.1126/science.274.5294.1859 (1996).

45. Scott, K. A. *et al*. Coarse-grained MD simulations of membrane protein-bilayer self-assembly. *Structure* **16**, 621–630, https://doi.org/10.1016/j.str.2008.01.014 (2008).

46. Guros, N. B., Balijepalli, A. & Klauda, J. B. The role of lipid interactions in simulations of the α-hemolysin ion-channelforming toxin. *Biophys. J.* **115**, 1720–1730, https://doi.org/10.1016/j.bpj.2018.09.009 (2018).

47. Desikan, R., Patra, S. M., Sarthak, K., Maiti, P. K. & Ayappa, K. G. Comparison of coarse-grained (MARTINI) and atomistic molecular dynamics simulations of α and β toxin nanopores in lipid membranes. *J. Chem. Sci.* **129**, 1017–1030, https://doi.org/10.1007/s12039-017-1316-0 (2017).

48. Periole, X., Cavalli, M., Marrink, S.-J. & Ceruso, M. A. Combining an elastic network with a coarse-grained molecular force field: Structure, dynamics, and intermolecular recognition. *J. Chem. Theory Comput.* **5**, 2531–2543, https://doi.org/10.1021/ct9002114 (2009).

49. Tagliazucchi, M. & Szleifer, I. Transport mechanisms in nanopores and nanochannels: can we mimic nature? *Mater. Today* **18**, 131–142, https://doi.org/10.1016/j.mattod.2014.10.020 (2015).

50. van Oeffelen, L. *et al*. Ion current rectification, limiting and overlimiting conductances in nanopores. *PLoS One* **10**, e0124171, https://doi.org/10.1371/journal.pone.0124171 (2015).

51. Liu, Y. & Zhu, F. Collective diffusion model for ion conduction through microscopic channels. *Biophys. J.* **104**, 368–376, https://doi.org/10.1016/j.bpj.2012.11.3826 (2013).

52. Cozmuta, I., O’Keeffe, J. & Stolc, V. Towards an MD simulation of ion currents in the α-hemolysin channel. In *2003 Third IEEE Conference on Nanotechnology, 2003. IEEE-NANO 2003.*, vol. 1, 187–190 vol. 2, https://doi.org/10.1109/NANO.2003.1231747 (2003).

53. Nosé, S. A unified formulation of the constant temperature molecular dynamics methods. *J. Chem. Phys.* **81**, 511–519, https://doi.org/10.1063/1.447334 (1984).

54. Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions. *Phys. Rev. A* **31**, 1695–1697, https://doi.org/10.1103/PhysRevA.31.1695 (1985).

55. Martini force field., http://cgmartini.nl/index.php/example-applications2/lipid-membranes. (2019).

56. Memprotmd., http://memprotmd.bioch.ox.ac.uk. (2019).

## Acknowledgements

This work was performed using HPC resources from GENCI-CINES under the grants 2015 - c2015077139, 2016 - c2016077139 and 2017 - c2017077139. Part of this work was funded by CNRS, Défi InFiniTi for the DYNANO project (2016 and 2017).

## Author contributions

N.B., D.D. and R.R. conducted the simulations and analyzed the results. All authors participated in writing and reviewed the manuscript.

## Competing interests

The authors declare no competing interests.
**Additional information**

**Supplementary information** is available for this paper at https://doi.org/10.1038/s41598-019-51942-y.

**Correspondence** and requests for materials should be addressed to R.R.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_8.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2019