Non-equilibrium molecular dynamics study of human aquaporin-2 in the static external electric fields

Abstract In this paper, non-equilibrium MD simulations (NEMD) of human aquaporin-2 (AQP2) in the presence of an external static electric field have been performed along + z and − z directions of the pore axis. The impacts of the electric field direction on the gating mechanism corresponding to the selectivity filter (SF) region of AQP2 have been studied. Besides, the effects of applied external electric field on the PMF profile of water molecules translocation, water permeability, and molecules dipole orientation are investigated. Our results showed that when the external electric field is implemented along the + z direction of the channels, the selectivity filter region is kept in the wide conformation for the majority of the time. Therefore, a remarkable increase in the overall water permeability can be seen compared to the case without any external electric field. This is in contrast to the effects of − z-directed electric field on the conformations of the selectivity filter, which induces mostly narrow conformations in this constriction region. A substantial higher energy barrier emerged in the middle of the AQP2’s pores under the effect of –z-directed electric field in comparison with the zero and + z-directed electric field strengths, which is mainly ascribed to the deviation from bipolar dipole orientation within the AQP2’s pores. Communicated by Ramaswamy H. Sarma


Introduction
Aquaporins (AQPs) exist in the cell membranes of some tissues, acting as semi-permeable conducting pores. The presence of AQPs is crucial for the body's water homeostasis and controlling the volume of the cells (King et al., 2004). water passes through cell membranes either by diffusion within the lipid bilayer or under an osmotic pressure difference through AQPs. The transport of water molecules through AQP pores is recognized to be the fastest among biological processes (Hub et al., 2009). Besides, translocation of other small molecules and gases such as ammonia, carbon dioxide, and nitric oxide within AQP channels have been observed (Payne & Forbush, 1994;Wang & Tajkhorshid, 2010). Malfunctionality of AQPs may lead to various diseases, including diabetes insipidus, skin disorder, obesity, and cancer (Kreida & T€ ornroth-Horsefield, 2015;Verkman, 2012). Human AQP2 located mainly in the principal cells of the collecting ducts in the kidney, is an important member of the AQPs family and it carries out the body's water balance tuning (Frick et al., 2014). Based on the reports of experimental studies, AQP's monomers in which four of them constitute a tetramer, are mainly composed of six membrane alpha-helices with a central water translocation pore and two cytoplasmic loops (Murata et al., 2000). Each monomer has two conserved Asn-Pro-Ala (NPA) motifs, located in the middle of the pore, acting as a constriction to the translocation of water molecules and other solutes (Frick et al., 2014).
Moreover, there is another important conserved constriction inside the pore known as the aromatic/Arg (ar/R) selectivity filter made of arginine and an aromatic residue.
Detection and analysis of water transport through biological transporters such as aquaporins on a molecular level have become possible in recent years. Molecular dynamics (MD) simulations are known as the main computational tool used for the study of biological processes such as water conduction within AQPs (Binesh & Kamali, 2015).
After the first MD simulation by our group on AQP2 in which the water transport properties within AQP2 channels were evaluated (Binesh & Kamali, 2015), another paper clarified some other issues including the energetics mechanism (Padhi & Priyakumar, 2017).
MD simulations have been helpful in the investigation of the gating mechanism in aquaporins (Alberga et al., 2014;Garate et al., 2011;Janosi & Ceccarelli, 2013;Yamamoto et al., 2014). Janosi and Ceccarelli (Janosi & Ceccarelli, 2013) showed that the channels of AQP5 can be open or closed because of different conformations of histidine at the CE. Similar behavior was reported for AQP4 (Alberga et al., 2014). Recently, the gating impacts of selectivity filter (SF) as wide/narrow conformations on the mechanism and dynamics of water passage through the wild-type AQP2 and two NDI causing mutants AQP2-V168M and AQP2-G64R were investigated .
The development of an effective method of inducing water flow at the nanoscale is of critical importance to the advancement of several nanotechnologies such as lab-on-achip devices, drug release, delivery systems, and desalination. Recently, it has been shown that the application of static and dynamic electric fields can affect the rate of permeation in lipid membranes and aquaporins (Garate et al., 2009;. Considering the increasing focus on mechanisms underpinning the biological effects of electric and electromagnetic (e/m) fields (English & Waldron, 2015), as well as focused attention on applying intense electric fields to disparate applications, e.g. nanotechnology for biosensing (Mulero et al., 2010), there is a great interest in investigating the effect of external electric fields on molecular/water transport in confined geometries (such as nanopores or aquaporins) . The ability to model the real shape of biomolecules makes molecular dynamics a powerful tool in numerical investigations of this field of study. However, because of the complexities of modeling the effect of external forces such as electric forces in the dynamics of biomolecules in aquaporins, only a few researchers applied non-equilibrium MD for simulations regarding the effect of electric fields on water transport through biological membranes.
In 2011, Garate et al. (2011) performed non-equilibrium MD (NEMD) simulations of h-AQP4 embedded in a solvated lipid bilayer and considered the effects of continuously applied static and alternating electric fields on the watertransport process and key features such as single-channel water permeabilities. Reale et al. (2013) carried out a similar simulation of embedded-bilayer h-AQP4, albeit in the absence and presence of nanosecond-scale static and alternating electric-field impulses. Their observations showed that water self-diffusion increases within aquaporin channels during, and immediately following field impulses. Marracino (2016) established water-permeation effects under fields applied perpendicularly to the pores. Their results illustrated decreased levels of water osmotic permeability within aquaporin channels during orthogonally-oriented field impulses. Bernardi et al. (2019) employed non-equilibrium MD simulations in external electric fields to observe and control Cl À conductivity across electro-pores formed in transmembrane proteins. It was found that the intensity of the electric field is required to be modulated carefully to maintains the electro-pores of the transmembrane in a quasi-stable state and control the permeability flux without disrupting the cell-wall structural stability.
In this work, non-equilibrium MD simulations (NEMD) of human AQP2 in the presence of an external electric field have been carried out in both directions of the pores' axis. The effects of the electric field direction on the gating mechanism associated with the selectivity filter of AQP2 have been investigated. Furthermore, the variations of the PMF profile of water molecules transport, water permeability, and molecules dipole orientation in response to an external electric field are explored.

Simulation details
In this study, a high resolution (2.75 Å) X-ray crystal structure of the human AQP2 tetramer was utilized (Frick et al., 2014). This initial structure was downloaded from the Protein Data Bank (PDB ID: 4NEF). Then, the tetramer structure of AQP2 was placed into a POPC (palmitoyl-oleoyl-phosphatidylcholine) lipid bilayers using the Membrane Builder plugin of VMD software v1.9.3 (Humphrey et al., 1996). To avoid any overlapping, the aquaporin and POPC membrane were aligned and the lipid segments and water molecules in the range of 0.8 Å of the protein were removed. Hydration of the system, including the protein and membrane, was done by adding water molecules layers of 18 Å thickness on the two sides of the system with the aid of the Solvate plugin of VMD software. Furthermore, the system was neutralized by adding Na þ and Clions with the 0.15 mol/L ionic concentration using VMD's auto-ionize plugin. Figure 1 depicts the entire system, utilized in our MD simulations. Figure 1 depicts a graphical depiction of the simulated system indicating the þ z and À z directions.

Molecular dynamics setup
All the Molecular dynamics simulations are carried out using the NAMD 2.13 simulation package (Phillips et al., 2005). The CHARMM36 force field was utilized for proteins (MacKerell et al., 1998) and phospholipids (Feller & MacKerell, 2000). Moreover, water molecules were modeled with the TIP3P model (Jorgensen et al., 1983). All of our simulations were carried out in an isothermal isobaric ensemble (1 atm, 310 K) by employing the Langevin piston method and the Langevin dynamics thermostat with a Langevin damping coefficient c ¼ (Feller et al., 1995). For minimizing the finite-size effects, periodic boundary conditions were applied in all three directions. Particle Mesh Ewald (PME) method (Essmann et al., 1995) was implemented to obtain full long-range electrostatic interactions and the SHAKE algorithm was used to constrain all hydrogen atoms (Miyamoto & Kollman, 1992). Furthermore, to calculate the van der Waals interactions, a cut-off distance of 12 Å with a switching distance of 10 Å was used.

Molecular dynamics simulation
Prior to production simulations, MD simulation of the system was performed in three steps. In the first step, the whole system except lipid tails was fixed to perform the minimization and equilibration of the lipid tails for 0.5 ns. In order to minimize the system in configuration space, minimization and equilibration of the system were carried out in the second step for 10 ns by applying harmonic constraints only to the protein atoms. Finally, by removing all of the restraints, MD simulations of the systems were performed at constant temperature and pressure of 310 K and 1 atm, respectively.
After 130 ns of equilibrium MD simulation of AQP2 under zero-field, employing the model used in our previous paper , non-equilibrium MD simulations have been performed under the effects of an external electric field along þ z and À z directions of the pore axis. The intensity of the implemented electric field has been selected as 0.07 V/ nm, similar to the value used in (Garate et al., 2011). Such a field intensity is in the order of those applied in the experimental studies, and the bilayer has been found to be stable under such electric field intensities.
As has been reported in the previous studies concerning electroporation within lipid-bilayer, the field strength required for electroporation induction is at least 0.2 V/nm (Bernardi et al., 2019;B€ ockmann et al., 2008;Marracino et al., 2018) which is about three times greater than the value implemented in our study.
In another study by Tieleman (English & Waldron, 2015;Tieleman, 2004), the field intensity threshold for pore formation was reported to be even higher, as strong as 0.5 V/nm.
Applying an electric field of this magnitude causes the lipid bilayer to become unstable and even influences the stability of the aquaporin. As a result, it causes permanent changes in the membrane permeability and may even result in cell death English & Waldron, 2015). For this reason, we chose a strength of 0.07 V/nm to make sure that the structural stability of the membrane is not altered and that the membrane permeability has not been affected by the electroporation.
The electric field (E) is applied to each partial-charge site i, with charge q i by:

Water permeability analysis
To compute water permeability constants, the Continuoustime random walk (CTRW) model was applied. This model, presented by Berezhkovskii and Hummer (Berezhkovskii & Hummer, 2002), is a known model used for studying water molecules passing within biological pores. It is assumed that nano-pore is filled with water molecules in a single-file fashion and these molecules pass concurrently through the channel. In order to compute the osmotic permeability coefficient (p f ) from equilibrium molecular dynamics simulations, the collective diffusion model proposed by Zhu et al. (Zhu et al., 2004) was used. The normalized collective coordinate n(t) can be calculated based on the displacement of each water molecule through the channel in the z-direction (dz i ) during the time step d t : where indicates the size of the time step. S(t) and L(t) represent the set of passing water molecules and the channel lumen length, i.e. the constriction region containing the water molecules in single-file fashion transport, respectively. The collective diffusion constant can be obtained by the following equation (Binesh & Kamali, 2015): where hn 2 (t)i presents the mean square displacement (MSD) of n(t). Then, the osmotic permeability (p f ) can be obtained by: Here, represents the average volume of a single water molecule.

Calculation of the potential of mean force (PMF)
The potential of mean force (PMF) was computed within the channel by dividing the channel into zones of 2 Å along the z-axis and then calculating the average water density in each zone (Sasseville et al., 2011): where K, and T represent the Boltzmann constant, the bulk density, and the absolute temperature, respectively.

Results and discussion
3.1. Effects of electric field on the RMSD profiles of the pores in the human AQP2 In this section, the RMSD profiles of AQP2 under the effect of an electric field in þ z and À z directions are plotted in Figure 2. Each of the profiles indicates averages in the two independent MD simulations. As shown, it takes about 150 ns for the RMSD profile of AQP2's atoms to reach adequate stability and converge to about 3.7 Å, when the transmembrane is under the effect of -z-directed electric field. But, when the þ z-directed electric field is implemented, the system reaches sooner the stability after about 70 ns and its RMSD converges to the value of 2.2 Å. Moreover, remarkably higher oscillations can be seen for the case in which the À zdirected electric field is applied along with the aquaporin pores compared to the case with the þ z-directed electric field. It means that applying an electric field in the þ z direction perturbs the overall structural stability less than À zdirected electric field. To make sure of the structure's stability, for both cases, the last 350 ns of MD simulation was selected for subsequent data analysis.

Effects of electric field on the conformations of SF region as narrow/wide in AQP2
In this section, the effects of the electric field on the indicator distances defined between the residues constructing the pore at the SF region have been investigated. In order to examine the gating mechanism of AQP2 in the SF region under the effect of an external electric field, the local constriction indicator (D) has been used, which indicates the distance between the NE2 nitrogen atom of His172 and NE nitrogen atom of Arg187 at the SF region. Figure 3 indicates the variation of this indicator within AQP2 under the effect of the À z-directed electric field during the 350 ns of simulation. As can be observed, the transitions between wide and narrow states in AQP2's SF region have been significantly reduced compared to the zero-field case studied in our previous work . This can be due to the interaction between electric force and the side chain of the residue His172 and the reduction of its fluctuations.
Similar behavior was observed when the variations of the SF region were illustrated for independent simulations with 350 ns duration for the negative and positive directed electric fields, as shown in Figures S1 and S2 (Supplementary material), respectively. For convenience, the variations of the SF region associated with AQP2 without an external electric field are given in Figure S3 (Supplementary material).
Moreover, it is evident from Figure 3 that when the À zdirected electric field is applied, the constriction region of SF remains in its narrow confirmation in most of the monomers. This is in contrast to the opposite effect for fields along the þ z-direction, as shown in Figure 4. When the electric field is applied along þ z direction within AQP2, the SF region remains wide for the majority of the time. It is noteworthy that in the case that no external electric field is applied through the aquaporin, there is a protein intrinsic electric field pointing toward the þ z-axis of the channels (English & Garate, 2016), which keeps the SF region open for the majority of the time.
Therefore, exerting the electric field in the same direction of this intrinsic field induces wide conformations for longer times. Figure 5 depicts the average radius profiles for the cases with negative and positive directed external electric fields plotted by HOLE 2 (Smart et al., 1993). Besides, the standard deviation for the profiles during MD runs has been computed and can be observed as shadows on this graph. The distribution of the radius profile indicates that the pore at the SF region for the negative directed electric field is limited to approximately 0.5 Å making the transport of water molecules very difficult. On the other hand, a positive directed electric field induces wide conformation in the SF region. Our calculated smallest radii for the negative and positive directed electric fields are in very good agreement with the values corresponding to narrow and wide states reported in our previous work .

Effects of electric field on the PMF profiles of the pores in AQP2
Variations of the average free energy profile of water molecules transport through all the four monomers of AQP2 under the effect of an electric field for zero-field and statically applied field conditions (along both þ z and -z-directions) have been demonstrated in Figure 6. The profiles corresponding to the cases with applied electric fields are averages over the monomers in the two independent MD simulations. The convergence of the PMF profiles for all the cases was checked, employing the block averaging method (Grossfield & Zuckerman, 2009). For instance, the convergence plots of the zero-field case are given in Figure  S4 (Supplementary material).
A higher energy barrier in the middle of the pores of AQP2 under the effect of À z-directed electric field can be seen, which can be due to the deviation from bipolar dipole orientation inside the aquaporin's pores discussed in detail in the following sections This energy barrier peak impedes the transport of water molecules within the NPA region leading to lower water permeability (see Table 1). Furthermore, the energy barriers within the SF region of the channel for þ z and À z directions are consistent with the conformations of the SF region, which was shown in Figures 3 and 4, such that mostly narrow states of the SF region have resulted in higher energy barrier in this zone for the À zdirected electric field.

Water transport properties of AQP2 with and without the effect of an external electric field
In this section, the effect of the electric field on the water permeability within the monomers of AQP2 has been investigated.  The single-channel osmotic permeability (P f) of AQP2 under the effect of an electric field with an intensity of 0.07 V/nm for wide/narrow conformations and also for the total simulation time can be seen in Tables 1 and 2. Since applying þ z directed external electric field has induced the channels to be wide most of the time (as shown in Figure 4), a substantial increase in the overall water permeability () can be observed compared to the case that no external electric field has been applied () .
However, when the electric field is applied along the À zdirection of the pores' axis, as can be seen in Table 1, the overall water permeability () decreases compared to the zero electric field configuration. The reason is that the conformation of the SF region has become mostly narrow in monomers A, C, and D, as discussed in the previous section. Interestingly, when the permeability is compared in two cases for an identical conformation of the SF region, it is found that the permeability is similarly reduced. The reason behind this can be a deviation from the well-known bipolar dipole orientation of water molecules within aquaporins discussed in the following section in detail. This is also consistent with Figure 5, where the higher energy barrier due to the perturbation of bipolar dipole orientation of water molecules can be observed near the NPA region. Figure 7(a) and (b) illustrate the cumulative permeation events in both directions of the monomers during 350 ns of the MD trajectory time under the effect of E=-0.07 V/nm and E¼ þ0.07V/nm, respectively. As can be seen, when the positive directed electric field is applied, the water molecule permeations within all the monomers increase with time almost linearly, indicating the wide state in these pores. For the negative electric field, however, the permeation of water molecules decreases substantially. It is noteworthy that in monomer 2 under the effect of the negative electric field, in the early times when the SF region is in the open state, the slope of the graph is higher and then decreases. These observations are compatible with this monomer's wide and narrow conformations, as shown in Figure 3.

Effects of electric field on the water dipole orientation
It is reported that the conformational change of the pore is not the only factor that affects the water dynamics within the channel. Dipole orientation of water molecules can affect the water permeability inside the constricted channels remarkably (Su & Guo, 2011;Wan et al., 2009). A well-defined dipole orientation of water molecules is observed within the lumen of the AQP's pores (de Groot & Grubm€ uller, 2001;Jensen et al., 2003). As observed in previous molecular dynamics simulations (Hashido et al., 2005;Tajkhorshid et al., 2002) (see Figure 8), while water molecules transport within the NPA region, a reverse in the direction of water dipoles occurs. The hydrogen bonding within the NPA motifs and the helix macrodipoles caused by the B and E helices (Tajkhorshid et al., 2002) determine the water dipole orientation.
Inside the constriction region, the water dipoles align with the helix macrodipoles (de Groot & Grubm€ uller, 2001) induced by the B and E helices. The center of the rotation is located within the NPA region, where the helices B and E meet.
This effect is essential for selectivity. Since only the molecules with a large dipole moment have the ability of rapid inversion within this region (de Groot & Grubm€ uller, 2001).
The effect of the electric field on the distribution of water dipole orientation within AQP2 is shown in Figure 9. The profiles corresponding to the cases with electric fields are averages over the monomers in the two independent MD simulations. As can be observed, applying an electric field in the À z-direction causes the distribution of water dipole Figure 5. Comparison of the average pore radius profiles (R) along the pore axis (z) for the negative and positive directed electric field in wild-type AQP2. Each profile has been averaged over two independent MD simulations. The shadows represent the standard deviation. Figure 6. PMF profiles of water molecules permeation as a function of position along pore axis (z) through wt-AQP2 with and without the effect of the external electric field. Each of the profiles represents averages over two independent MD simulations. orientation to shift to lower values. It can be due to the impact of À z directed electric field on the water molecules making them align with the field. This leads to deviation from the well-defined dipole orientation within the constriction region to some extent. In particular, water dipole orientation on both sides of the NPA region has become negative. This effect on the water osmotic permeability could be seen in Tables 1 and 2. As mentioned in the previous section, when the À z-directed electric field has applied, the water permeability diminishes in the similar conformation of the SF region compared to the zero electric field configuration, which can be due to the deviation of water dipole orientation from its optimal configuration when no external electric field is applied across the aquaporin. On the other hand, when the channels are affected by þ z directed external electric field, no significant change in the water dipole orientation can be detected. As mentioned before, without an external electric field, the protein intrinsic electric field still exists along the þ z-axis direction of the pores.
Therefore, implementing the electric field in the same direction of this intrinsic field does not alter the dipole orientation of the water molecules remarkably. In particular, the sign of the dipole orientation of the water molecules differs between two sides of the NPA region indicating the bipolar    . Dipole orientation of water molecules within the AQP2 pore under the effect of electric fields E ¼ 0 and 0.07 V/nm. h represents the angle between a water dipole inside the pore and axis z. Each of the profiles shows averages over two independent MD simulations.
characteristic of this dipole orientation distribution through the aquaporin's channels.

Conclusion
In the present work, non-equilibrium MD simulations of human AQP2 in the presence of an external static electric field have been carried out for the first time along þ z and À z directions of the pore axis. The intensity of the applied electric field, as 0.07 V/nm, is in the order of those applied in the experimental studies, and the bilayer has been reported to be stable under such electric field strength. The results revealed that when the external electric field is implemented along the þ z direction of the channels, which is the same direction for the intrinsic electric field within the protein's pores, the SF region remains in the wide conformation for the majority of the time. Consequently, a substantial increase in the overall water permeability can be observed compared to the case with zero external electric fields. This is in contrast to the effects of À z-directed electric field on the conformations of the selectivity filter, which induces mostly narrow states in this constriction region. A higher energy barrier has been observed in the middle of the AQP2's pores under the effect of -z-directed electric field compared to zero and þ z-directed electric field strengths, which mainly corresponds to the deviation from bipolar dipole orientation inside the aquaporin's pores. According to the observations, applying an electric field in the À z-direction, which is opposite to the direction of the intrinsic electric field within the protein, causes the distribution of water dipole orientation to shift to lower values. It can be due to the impact of -z-directed electric field on the water molecules making them align with the field and leads to perturbation of the well-defined dipole orientation from its optimal bipolar distribution within the constriction region to some extent, i.e. water dipole orientation on both side of NPA region becomes negative.
Disclosure statement