Simulation study of a rectifying bipolar ion channel: Detailed model versus reduced model

We study a rectifying mutant of the OmpF porin ion channel using both all-atom and reduced models. The mutant was created by Miedema et al. [Nano Lett., 2007, 7, 2886] on the basis of the N-P semiconductor diode, in which an N-P junction is formed. The mutant contains a pore region with positive amino acids on the left-hand side and negative amino acids on the right-hand side. Experiments show that this mutant rectifies. Although we do not know the structure of this mutant, we can build an all-atom model for it on the basis of the structure of the wild type channel. Interestingly, molecular dynamics simulations for this all-atom model do not produce rectification. A reduced model that contains only the important degrees of freedom (the positive and negative amino acids and free ions in an implicit solvent), on the other hand, exhibits rectification. Our calculations for the reduced model (using the Nernst-Planck equation coupled to Local Equilibrium Monte Carlo simulations) reveal a rectification mechanism that is different from that seen for semiconductor diodes. The basic reason is that the ions are different in nature from electrons and holes (they do not recombine). We provide explanations for the failure of the all-atom model including the effect of all the other atoms in the system as a noise that inhibits the response of ions (that would be necessary for rectification) to the polarizing external field.


Introduction
Rectification mechanisms in nanopores and ion channels are based on asymmetries in the structure of the pore [1,2]. The asymmetry is either geometrical or electrostatic in nature. In the former, the shape of the pore is asymmetrical as in the case of conical nanopores [3,4].
The latter case, when the charge distribution in the pore is asymmetrical [5], is the subject of this study. This phenomenon is well known in the case of semiconductor diodes [6,7], where the charge asymmetry is achieved by doping different regions of the device differently thus forming an N-P diode, where the majority charge carriers are electrons and holes in the N and P regions, respectively. The N-P junction between these two regions forms a depletion zone for both electrons and holes. An external electric field in forward (ON) and reverse (OFF) bias acts differently on this region by making it even wider in the OFF state and thinner in the ON state. In the ON state, the majority carriers will conduct the current, while in the OFF state, the minority carriers will do the job; hence, the rectification.
In this paper, we consider devices where the charge carriers are ions solvated in a liquid solvent (usually water) that migrate through a pore that is embedded in a membrane. The two major classes of these pores are artificial nanopores and biological ion channels. Nanopores with an N-P charge distribution on their pore walls are called bipolar nanopores [8][9][10][11][12][13][14][15][16][17][18][19]. Nanopores are etched into plastic membranes

Experimental facts for the RREE mutant
In the experimental work of Miedema et al. two filters have been identified inside the pore (see table 1 of reference [47]). In the first filter, the negative amino acids D113 and E117 have been mutated into positive arginines, R113 and R117. In the second filter, the positive arginines, R167 and R168, have been mutated into negative glutamates, E167 and E168. This way, the first filter has been positively doped, while the second filter has been negatively doped, at least, in theory (see figure 1). The point mutations aiming the N-P junction are hard facts, but we do not know whether the protein is folded in the way we want it to fold: crystal structure data are not available for the mutant.
Using 0.1 M NaCl and ±100 mV voltage, the authors found a rectification 0.22 ± 0.02 for the RREE mutant as opposed to the value 1.14±0.03 in the case of the WT channel. Rectification, which is a voltagedependent quantity, is defined as In the case of a 1M NaCl electrolyte, the rectification values are 0.65 ± 0.06 and 0.99 ± 0.01 for the RREE and WT channels, respectively. Rectification, therefore, decreases as concentration increases.
The authors hypothesize in a cartoon (figure 5 (b) in their paper [47]) about the rectification mechanism that is adapted from the case of the semiconductor N-P diodes. The supposed mechanism is that a depletion zone is formed at the junction of the N and P regions that becomes wider and more depleted at the OFF sign of the voltage. It seems to be a widespread assumption that the rectification mechanism is the same in bipolar pores (where ions are the charge carriers) and in semiconductor diodes (where electrons and holes are the charge carriers). In this paper, we show that the mechanism of rectification is different, or, at least, that it can be different in narrow nanopores and ion channels.  [34,35] and its RREE mutant (bottom row) [47] made by the VMD package [52] after changing the indicated amino acids: D113→R113, E117→R117, R167→E167, and R168→E168.
The structure of the OmpF trimer [34,35] was constructed according to the ProteinDataBank database (identifier: 2OMF). The protein/membrane complex was generated with the help of CHARMM-GUI [54], embedding the protein into a DMPC lipid bilayer. We used the VMD program package [52] to mutate the WT channel into the RREE mutant (see figure 1).
We performed all-atom MD simulations with the GROMACS program suite [55,56] using the leap frog integrator with a 2 fs timestep. The system temperature was set with the Nose-Hoover thermostat [57]. Simulations in the N pT ensemble were conducted with a Parrinello-Rahman barostat [58]. We used CHARMM27 force-field based flexible models together with position restraints for the backbone atoms of the protein [59]. The bonds of hydrogen atoms were considered rigid; this allows us to use a slightly larger timestep (larger than that required for an accurate simulation of bond vibration with hydrogen atoms).
In simulations with electric fields we applied a ±200 mV potential (with the ground at the left-hand side).
Periodic boundary conditions were present in all spatial directions.
Most of our simulations were performed in a simulation cell with the size of 105.6 × 105.6 × 114.5 Å 3 in x, y, and z dimensions with z being the transport direction. The solvent phase was constructed of 561 Na + , 528 Cl − , and 29 317 TIP3P water molecules resulting in ≈132 000 atoms including ≈ 15 000 from the protein trimer, and ≈ 28 000 from the DMPC lipid layers.
To check for a system size dependence, we performed two simulations (for 200 and -200 mV) for a larger simulation volume with approximate dimensions 220×220×1 130 Å 3 containing four RREE trimers and ≈ 5 000 000 atoms. The protein and lipid membrane were constructed using four times the smaller simulation volume that was elongated in the direction of the transfer (along axis z) and filled with water and ions.
We followed the simulation procedure of Faraudo et al. [44]. In five preliminary equilibration runs we did not apply an external electric field. We started with an energy minimization run after the construction of the simulation cell. This was followed by a 100 ps NV T run at 100 K and another 100 ps NV T simulation at 296 K. After these steps we turned the barostat on and performed a 1 ns N pT calculation at 296 K and 1 bar with isotropic pressure coupling. The last preliminary equilibration step was to perform a 3 ns N pT simulation at 296 K and 1 bar pressure with semi-isotropic pressure coupling (independent coupling in the direction of transfer).
After we let the system relax, we started the simulations with an applied external field. To achieve a stationary state we did a 10 ns long NV T run at 296 K and with an external electric field corresponding to a 200 mV potential difference across the simulation cell in the z direction. Next, we did the actual production run in which we counted the diffusing particles through the membrane. We have monitored the number of ions that completely crossed the protein by following the individual trajectories of each ion. An ion was considered to cross the channel if it is initially at one side of the membrane, and then ends at the opposite side of the membrane after propagating through the protein channel (some ions enter the channel but instead of crossing it, they return to the bulk where they started).

Results for the all-atom model
The number of counted ion-crossings as a function time (in the final production run) is plotted in figure 2 for the WT channel. We found the channel slightly selective for Cl − at 200 mV. The real channel is known to be slightly cation selective. We also have simulations for KCl, but with much shorter runs and weaker statistics. In this case, we found K + selectivity for 200 mV. No significant rectification was observed.
When the relevant amino acids are mutated (see figure 1), the channel becomes perfectly anion selective, so we plot only the Cl − currents in figure 3. The lack of cation current is probably due to the mutations made in the left-hand side filter; the ring formed by positive amino acids has a very narrow opening that repulses the cations effectively. The negative ring on the other side has a much larger hole in the middle that makes the passage of anions possible. It is more important that we have not found rectification for this model of the RREE mutant. The Cl − current is practically the same for 200 mV and −200 mV within the statistical error of the simulation. These statistical errors can be estimated on the basis of the standard deviations of block averages; we obtained a large number for the error (±50 pA). Even if this large error indicates a weak statistics for the simulations, one thing can be concluded from figure 3 safely: rectification cannot be observed within the applied simulation lengths. From shorter runs for KCl we can draw the same conclusion.

13802-4
If we want to find an explanation for the discrepancy between experiment and simulations, or, at least, we want to get closer to the explanation, we can look at the concentration profiles. Figure 4 shows the concentration profiles, n i (z), which are defined as the average number of ions in a slab divided by the volume of the whole slab (the simulation cell is divided into slabs with a thickness of 2.5 Å in the z direction). An alternative way to plot the results is to show an effective local concentration, c i (z), where the average number of ions is divided by an effective volume. The effective volume is defined as the part of the whole slab, where the ions do not overlap with the body of the protein and the membrane -practically, the region of electrolyte. We will show results for the concentration profiles (in mol/dm 3 ) in the case of the all-atom model, because the effective volume is not a well-defined quantity due to the flexibility of the protein/membrane system.
In figure 4, one of the relevant observations is that the Na + ions are depleted inside the pore (note the logarithmic scale of the concentration axis). This depletion zone acts as a high-resistance segment of the pore that effectively cuts the current of Na + . The other observation is that changing the sign of the voltage has little effect on the concentration profiles of Cl − (the ion that conducts). The effect is that a depletion zone is formed at ≈ −5 Å for −200 mV, while for 200 mV the depletion zone is formed at ≈ 5 Å. Rectification would happen if the depletion zone were deeper at one voltage than at the opposite sign voltage. Here, the depletion zone is just shifted. From the point of view of conductance, the two profiles do not make a difference, therefore, the currents are the same for the two opposite signs of the voltage. Third, the profiles obtained from the simulation for the large system (four trimers) and the small system (one trimer) agree. This justifies the use of the smaller simulation volume and indicates that the results obtained from it can be the basis of analysis.

Reduced model for a bipolar ion channel
The other way of figuring out what is going on in this system is to create a reduced model that takes into account only the "important" degrees of freedom and ignores the noise of the "unimportant" degrees of freedom. The "important" degrees of freedom are those that Miedema et al. [47] manipulated when they created their mutant in order to achieve a rectifying N-P junction in the ion channel. They are the amino acids that form an N-P junction inside the pore as shown in figure 1. To build a reduced model that is appropriate for our purpose, we choose the ion channel model that we used in our previous papers for the L-type calcium channel [50,[60][61][62][63][64][65][66][67][68], the Ryanodine Receptor calcium channel [51,[69][70][71][72], and the neuronal sodium channel [73][74][75]. These reduced models were able to capture the essential features of these channels and reproduce various anomalous selectivity behaviors.

Reduced model
In this model, we work with a reduced representation of the electrolyte, the protein, and the membrane. The ions are charged hard spheres immersed in a dielectric continuum that models the solvent implicitly. The ionic radii are 2 Å for both the cation and the anion (we work with a 1:1 electrolyte), the dielectric constant is = 78.5, the temperature is 298.15 K. The ions electrostatically interact through the screened Coulomb potential if they do not overlap (which is forbidden). The membrane is confined between two hard walls (their distance is 30 Å), with which the ions cannot overlap.
A pore of radius 4 Å penetrates the membrane. The pore has hard walls with which the ions cannot overlap. The central cylindrical portion (of length 20 Å) represents the selectivity filter. ; it is a parameter we can change). In the vestibules the diffusion coefficient is interpolated between these two values in a way described by Boda [51].
The simulation cell is a finite cylinder with hard walls (the 3D cell is obtained by rotating figure 5 around the z-axis). The two cylindrical compartments on the two sides of the membrane represent the two bulk regions between which the ion transport flows. Such a bulk compartment has two parts: one is a transport region that is in non-equilibrium (indicated by a blue line), and the other is an equilibrium bulk region that surrounds it (outside of the blue line). The NP transport equation is solved for the transport region and the boundary conditions are specified on the outer surfaces of the transport regions (two half cylinders).

NP+LEMC method
The ion transport is described by the NP transport equation: where j i (r) is the particle flux density, k B is Boltzmann's constant, D i (r) is the diffusion coefficient profile, c i (r) is the concentration profile, and is the electrochemical potential profile that is the sum of the chemical potential and the interaction with the mean electric potential, Φ(r). In these equations, z i is the ionic valence, e is the elementary charge, µ 0 i is a standard chemical potential, and µ ex i (r) is the excess chemical potential profile. The transport is driven by the gradient of the electrochemical potential, ∇µ i (r).

13802-7
the hard sphere ions studied here, this theory cannot be applied, because it is a mean field approach for point charges. To handle the hard sphere ions, a more developed statistical mechanical theory is needed, for example, the Density Functional Theory of Gillespie et al. [84,85].
Here, we use the LEMC method that is an adaptation of the GCMC method for a non-equilibrium situation [48][49][50][51]. The system is divided into small elementary cells, D k , in which different electrochemical potentials can be assumed [µ i (r k ), where r k is the center of D k ]. Such an elementary cell is assumed to be in local equilibrium that makes it possible to perform particle insertions and deletions with the acceptance criterion of GCMC simulations, but using the particle number in the given cell, N k , its volume, V k , and the electrochemical potential assigned to the cell, µ i (r k ). The energy of the ion insertion/deletion contains the interaction with all the ions in the whole simulation cell and the interaction with the applied field, Φ app (r).
The applied field is computed by solving Laplace's equation for the empty solvation domain (all the charges removed) with the Dirichlet boundary condition that the potential is zero at the half cylinder on the left-hand side and the value of the voltage, U , at the half cylinder on the right-hand side. These surfaces are indicated with a blue line in figure 5. The NP equation is solved inside this surface.
The LEMC simulation provides the concentration profiles as an output, c i (r k ), given an electrochemical potential profile, µ i (r k ). An iteration procedure is used to obtain a self-consistent system in which the flux satisfies the continuity equation, ∇ · j i (r) = 0, namely, the conservation of mass. The heart of the iteration can be summarized as Starting from an electrochemical potential profile in iteration n, the concentration profile for that iteration, c i [n], is obtained from LEMC. The electrochemical potential profile for the next iteration is obtained from writing the integral form of the continuity equation for the elementary cell, D k , as

Results for the reduced model
The electrical current flowing through the pore is obtained from where A is the cross section of the pore. The negative sign makes the current positive for positive voltage.
The current-voltage curves for different values of the filter diffusion constant are shown in figure 6. Rectification is clearly observed; the current is larger in magnitude at positive voltages than at negative voltages [see the rectification curves in the inset; rectification is defined in equation (1)]. Decreasing D filter i decreases the net current, but it has no effect on rectification. If we increase the number of positive/negative structural charges in the filters, rectification improves (results not shown). The fact that rectification is not sensitive to the diffusion constant indicates that rectification is rather determined by another factor in the NP equation: concentration, c i (r).
The effective local concentration profiles are shown in figure 7. The figure illuminates the rectification mechanism observed in the reduced model of a bipolar ion channel. It can be summarized as follows: The rectification mechanism is similar to that observed in semiconductor diodes from the point of view that enhanced depletion in the OFF state produces rectification, but the list of similarities ends here. In the case of semiconductors, the width of the N-P junction is modulated by the voltage. When electrons get into the P zone, they produce a net current even if they are not the majority charge carriers there. The reason is that they recombine with holes arriving from the other direction.
In the case of the ion channel, ions are the charge carriers that cannot recombine. Therefore, the anions, for example, must conduct current in their own depletion zone (the P zone) if we want a net current. The same is true for the cations. The total current is determined by the depletion zone, because that is the largest resistance element of the resistors connected in series if we imagine the slabs as resistors in an equivalent circuit.
The OFF state of the voltage makes its own depletion zone of an ionic species even deeper. The important zone from the point of view of depletion is not the junction zone between the N and P regions, but the N and P regions themselves.
This finding contradicts the usual assumption in the ion channel and nanopore literature where authors assume that the rectification mechanism is the same in semiconductor and ionic devices. We showed here that this is not necessarily true. A deeper discussion of the rectification mechanism observed for the ionic diode follows.

Discussion
In the following, we analyze how the concentration profiles become more depleted in the OFF state. First, we must realize that electrical double layers are formed at the two sides of the membrane. For example, in figure 7 at 100 mV (black curves with full circles) there are more anions at the left-hand side than cations, and vice versa on the right hand side. The important thing is that double layers of the opposite sign are formed in the case of −100 mV.
To understand why these oppositely charged double layers are formed, we need to look at the potential profiles (figure 8). The average electrostatic potential can be computed during simulation by inserting test charges into the system, sampling the potential with them, and averaging. The total electrostatic potential has two components: Φ reduced (r) = Φ app (r) + Φ ion (r). (9) The applied potential created by the electrode charges, Φ app (r), is obtained by solving Laplace's equation  for the ion-free system with the prescribed Dirichlet boundary conditions. The other term is the potential produced by the ions, Φ ion (r), that is related to the ionic charge distribution (including the structural ions) through Poisson's equation.
The slope of the total potential profile is supposed to be small in the bulk solutions because the resistance of the bulk electrolytes is small compared to the ion channel. The potential drop across the membrane region dominates over the drops in the bulk regions (solid curves with full circles). To achieve this, the ions have to arrange into a distribution that imposes an appropriate counterfield (red curves with open squares) against the applied potential. This is the Φ ion (r) term that is zero at the boundaries of the system, as it should be, if we expect it from the total potential to satisfy the boundary conditions. The Φ ion (r) profile is decisively influenced by the double layers shown in figure  This is an important distinction compared to the electron/hole charge carriers. Cations and anions do not recombine, but they are trying to stay close to each other and screen each other's electric field. If the amount of cations in the N zone is already small, it becomes even smaller if the amount of anions (that are eventually responsible for bringing the cations in) decreases.
Therefore, the change in the voltage sign has an indirect effect on the depletion zones of the minority charge carriers in a given zone. The OFF-state voltage creates double layers that deplete the majority charge carriers in the given zone. The further depletion of the minority charge carriers is a consequence of the depletion of the majority charge carriers.
The question arises why the all-atom model does not show the expected behavior. We do not see significant double layers in figure 4 and, what is more important, we do not see a significant effect of the sign change of the voltage. This can be seen even more clearly if we plot the charge profiles (the difference of cation and anion profiles). Figure 9 shows the profiles for both models. To understand the absence of double layers, let us investigate the potential profiles obtained for the all-atom model of the RREE mutant ( figure 10). Now there are more players in the simulation cell, so the potential has more components. In addition to the Φ ion (r) term that was produced solely by the ions in the reduced model, now we have components due to the partial charges in the protein, the membrane, and water: Φ all-atom (r) = Φ app (r) + Φ ion (r) + Φ protein (r) + Φ membrane (r) + Φ water (r).
(10) Figure 10 shows these four terms in the left-hand panels for voltages ±200 mV. The right-hand panels show the total potential with (top) and without (bottom) the applied potential. Our statistics, unfortunately, are quite poor, but a few major conclusions can be drawn nevertheless.
First, the qualitative statement can be laid down that the external field polarizes the system. In this cubes in that part of the system that is attainable for ions (practically, the electrolyte) and averaging over configurations.
case, however, it polarizes not only the ionic distribution, but also everything else in the system that carries partial charges. The external field must exert work to polarize the system. Only a portion of this work is spent on the ions, most of it is "wasted" on the protein, the membrane, and the water molecules. The ionic profiles, therefore, do not respond so sensitively to the polarizing field as in the case of the reduced model. Although the poor statistics prevent us from drawing accurate quantitative conclusions, it seems that water is the component that is chiefly responsible for creating the counterfield to the applied field (bottom-left-hand panel). The ionic profiles also respond to some degree if we look at the potential (topleft-hand panel), but this does not manifest in the change of the ionic distribution that would be sufficient to produce rectification on the basis of the mechanism seen in the reduced model.
The slope of the total potential in the right hand bulk is close to zero (top-right-hand panel) as a result of the counterfield (bottom-right-hand panel) that is added to the applied field (dashed lines in the topright-hand panel). This is clearly seen despite the poor statistics and the small size of the simulation cell.
As a matter of fact, this was the reason that we performed the simulations for the large (four trimer, 5 million atoms) simulation cell. We hoped that in a larger bulk we had more space for the double layers. Potential profiles have not been calculated for the large cell, but for the concentration profiles and the current we obtained the same results as in the small cell.
Summarizing, the presence of all the other atoms and charges in the all-atom system screen the small N-P region inside the pore so effectively that the external field has no observable effect on the ionic profiles, and, thus, its sign change does not produce any observable rectification.
Another possible explanation of the lack of rectification can be deduced from the top panel of figure 9. The charge profile looks like a P-N-P charge distribution rather than an N-P one. It is quite symmetric, albeit rectification requires asymmetry in the charge distribution.
Although these explanations of the failure of the all-atom model make sense, they do not explain why the mutated ion channel of Miedema et al. [47] does rectify. The answer can be some local structural effect that our model cannot capture.
It is strange that (1) Miedema et al. [47] assumed a mechanism of rectification (the N-P junction), (2) created an ion channel based on that assumption with point mutation, (3) showed that the channel really rectifies, but (4) the all atom model of this mutant does not rectify. (5) In the meantime, the reduced model -based on the same assumption Miedema et al. [47] started with -does rectify.
We do not really know where is the flaw in this chain. It can be that the mutant of Miedema et al. [47] folds in a way that has nothing to do with how we imagine its folding. It can be that the all-atom model is not accurate enough due to force field problems. By all accounts, there are several problems with classical force fields. They seem to be more appropriate to study local effects rather than long-range phenomena, including an applied field, screening double layers, and so on. Force fields that handle polarization more realistically are definitely needed. Finally, it can be some kind of problem with the MD methodology, although we think that this is improbable.

Summary
In this work we presented a system, in which a powerful experimental fact (rectification) is studied with all-atom and reduced models. We found the puzzling result that the all-atom model, that is supposed to be more "accurate", cannot reproduce rectification, while the reduced model, that is admittedly simplistic, can. The results show that there are cases when reduced representations can serve our purpose better than detailed representations if our purpose is to understand a given phenomenon.
Rectification is a result of the balance of long-range effects, such as the applied field and the counterfield of the ionic distributions. If we concentrate on these effects in our reduced model, we can better focus on the phenomenon at hand. Building efficient reduced models is far from being trivial. Our earlier works for ion channels [50,51,[60][61][62][63][64][65][66][67][68][69][70][71][72][73][74][75] showed that such models can capture an essential portion of reality that is necessary, and in some cases sufficient, to explain a well-specified phenomenon. All-atom models, however, can guide us in creating these models.
The other main message of this paper is that rectification mechanism in bipolar ionic diodes (biological ion channels or narrow synthetic nanopores) is different from the mechanism in semiconductor N-P diodes. One of the reasons is that ions are different in nature from electrons and holes. Cations and anions do not recombine, so an ion must go through the whole pore all the way, including its own depletion zone. That depletion zone is more depleted at the OFF sign of the voltage than at the ON sign.
The explanation is the effect of the double layers formed at the entrances of the channel as detailed above. These double layers are everywhere. They form at the wall of the nanopore too. If the nanopore is too wide compared to the Debye length, a bulk electrolyte is formed in the center of the pore. In this case, depletion zones do not form and the rectification mechanism described in this paper does not work efficiently. The interesting and efficient pores, therefore, are those whose radius is smaller than the Debye length. Ion channels obviously belong to this category.
The other difference between bipolar ionic and semiconductor diodes, therefore, is that narrow pores are needed in the case of ions as carriers in order to make the formation of depletion zones possible. There is no such a requirement in the case of semiconductors. Furthermore, while the junction region between the N and P regions is important in the case of semiconductors, it is the junction region at the entrances of the pore that has a large impact on the behavior of the system. The double layers extend into the N and P zones and deplete the majority carriers that, in turn, deplete minority carriers further in the OFF state.
Rectification mechanism in long nanopores can be different from that in short nanopores because the resistance of the pore itself dominates over the access resistances at the pore entrances in the case of long pores. This question has been thoroughly discussed by Vlassiouk et al. [14] using both numerical and analytical solutions of PNP. Interestingly, their concentration profiles (figure 2 in reference [14]) do not seem very different from ours (figure 7): in the OFF state, both ions become depleted compared to the ON state. Furthermore, the ions become depleted not only at the junction in the middle (that Vlassiouk et

13802-13
al. call a depletion zone), but also in the entire half zones in the pore (these are the real depletion zones, in our view).
For us, the profiles of Vlassiouk et al. [14] imply a similarity with the mechanism described in the present paper. Although Vlassiouk et al. emphasize that they found "a striking similarity to the corresponding solid-state devices", we suspect that the similarity is limited. It will be fun to sort out these uncertainties in future studies for nanopores. We expect that our NP+LEMC method will provide an additional insight compared to PNP studies due to its improved capabilities to handle ionic correlations in the nanopore and in the ionic double layer.