Forced diffusion of water molecules through aquaporin-5 biomembrane; a molecular dynamics study.

Aquaporins (AQPs) are protein channels located across the cell membrane which conduct the water permeation through the cell membrane. Different types of AQPs exist in human organs and play vital roles, as the malfunction of such protein membranes can lead to life-threatening conditions. A specific type of AQP, identified as AQP5, is particularly essential to the generation of saliva, tears and pulmonary secretions. We have adopted Molecular Dynamics (MD) simulation to analyze the water permeation and diffusion in AQP5 structure in a 0.5 microsecond simulation time window. The MD numerical simulation shows the water permeability of the human AQP5 is in the nominal range for other members of human aquaporins family. In addition, we have considered the effect of the osmotic water diffusion and the diffusion occurred by pressure gradient on the protein membrane. The water permeability grows monotonically as the applied pressure on the solvent increases. Furthermore, the forced diffusion increases the minimum radius of Selectivity Filter (SF) region of region AQP5 up to 20% and consequently the permeability coefficients enhance enormously compared to osmotic self-diffusion in AQP5 tetramer. Finally, it is revealed that the MD simulation of human AQP5 provides useful insights into the mechanisms of water regulation through alveolar cells under the different physical conditions; osmotic self-diffusion and forced diffusion condition.

In this study, Molecular Dynamics method is used to analyze the water permeability and diffusion in Aquaporin-5 (AQP5) structure in a 0.5 microsecond simulation time. For the first time, we have considered the effect of the osmotic water diffusion and the diffusion occurred by pressure gradient on the protein membrane. The water permeability grows monotonically as the applied pressure on the solvent increases. Interestingly, the forced diffusion increases the minimum radius of Selectivity Filter (SF) of region AQP5 up to 20% and consequently the permeability coefficients enhance to a great extent in comparison with osmotic self-diffusion in AQP5 tetramer.

Biomembranes are biological structures which act as a selectively permeable barrier against water and other solute permeation. Aquaporins (AQPs) compose a class of integral biomembrane proteins that exist in a wide range of living organisms, e.g. archaea, bacteria, plants and animals. It has been known that, water can be transported more easily across aquaporin cellular membranes rather than across phospholipid bilayer through the diffusion [1,2]. Aquaporins are of significant importance to cell volume control and transcellular water transport. Generally, water transport through AQPs occurs by a passive mechanism, but water regulation and cellular distribution change considerably from one cell and tissue type to another [3]. It should be noted that, the AQPs performance in human body is excessively important, as the malfunction of AQPs can result in several diseases. Genetic defects of AQP genes are notorious for several human diseases, such as nephrogenic diabetes insipidus and neuromyelitis optica [4].

There are different aquaporin types in the human body since there is a diverse range for the water homeostasis regulation in different human cells and organs [5]. Water can be transported through the AQPs of the kidney collecting duct, capillaries in the lung and secretory cells in salivary glands. Also other osmotic water regulation processes such as water homeostasis spinal fluids in brain, maintenance of lens transparency in the eye, secretion of tears, generation of sweat and saliva and etc. are directly affected by AQPs performance. According to their performance, AQPs are divided into two subsets of classical and non-classical AQPs. Classical AQPs which are water permeable, include AQP0, AQP1, AQP2, AQP4, AQP5, AQP6, AQP8 and AqpZ. Nonclassical AQPs comprise AQP3, AQP7, AQP9, AQP10 and GplF. Non-classical APQs are relatively permeated by water and other small solutes, especially glycerol. Therefore, the second subset is typically referred to as aquaglyceroporins [6,7]. Aquaglyceroporins contain two additional peptide spans compared to classical aquaporins, hence they have different behaviors in terms of glycerol permeation. Moreover the possible movement of ions and carbon dioxide is observed through some particular AQPs [8].

In the recent years, AQP transporters have attracted a lot of interest due to their exceptional physical properties [9-11]. AQPs are very hydrophobic, intrinsic membrane proteins which can allow for selective permeation of water through their biochemical structure. These unique properties make AQPs the ideal candidate for various applications. Therefore, it is essential to have a better understanding of dynamics of physical phenomena occurring in AQPs. There is no experimental procedure reported so far that demonstrates the water permeation through AQP membranes at atomic scale. Therefore, computational methods have been widely used to study biomembranes channels. Today with significant growth of computational capacities, numerical methods, such as Molecular Dynamics (MD), are becoming a promising tool to study the biophysical and biochemical processes at molecular length scale [12-15] in the absence of an effective experimental schemes. The MD has been particularly employed to study the biological process, such as gating mechanism of human AQP4 [16]. Also the water movement through AqpZ has been studied by MD simulation [17]. Moreover a recent work demonstrated that MD simulation of human AQP2 gives new insightful information about water permeation mechanisms and urine concentration in the human kidney [18,19].

Some particular kinds of AQPs allow bidirectional water fluxes in response to osmotic water permeation [20]. For instance, the water permeation through AQP0–AQP5 is bidirectional and symmetric. Water permeation in AQPs is a sophisticated process which has to be modeled at molecular length scale. This kind of permeation is first discovered in APQ1 as AQP1-mediated water transport in erythrocytes and the relevant reflection coefficients for urea and glycerol have been measured [20]. AQP1 along with AQP5 majorly participate in osmotic water transport between airspace and capillary compartments [21]. AQP5 has a significant function in lung fluid transport due to submucosal glands structure. Moreover, AQP5 provides the main route for gas transport in alveolar cells and exhibits substantial solute selectivity, which can have an important physiological effect on controlling gas fluxes [22]. Considering the importance of AQP5 performance in human alveolar cells, the main objective of the current study is to analyze water permeation in the AQP5 structure.

Recently, the MD simulation of water permeation through AQP5 under the osmotic diffusion and the associated gating phenomena have been carried out [23]. To have a deeper understanding of osmotic self-diffusive and forced-diffusive water permeability at molecular length scale, here we focus on AQP5 permeability and diffusivity. The MD method simulates the human AQP5 structure embedded in a lipid biological membrane. The biophysical parameter water permeation under osmotic diffusion, is extracted based on numerical calculations and is validated by experimental results related for AQP5. Moreover for the first time, we have studied the effect of external pressure field on the permeation and diffusion of water molecules through the channels of AQP5 crystal structure. We have found that for the forced diffusion condition, the external pressure field has a two-fold effect on the enhancement of the permeability coefficient; first, by directly increasing the velocity of water molecules passing through aquaporin nanopores and second, by enhancement of SF region radius. The present study lays the groundwork for investigating the permeation of other solute; like aqueous gas solution through AQP5 which has its own great physiological importance.

Methods

The System Structure

A high resolution (2 Å) X-Ray crystal structure of the human aquaporin-5 tetramer is studied [24]. The crystal structure is available on the Protein Data Bank (PDB ID: 3D9S). The AQP5 is a membrane protein structure composed of 7390 atoms. The AQP5 is placed into a Palmitoyl-Oleoyl-Phosphatidyl-Ethanolamine (POPE), lipid bilayer using Membrane Builder plugin of VMD 1.9.3 software [25]. The lipid is 120 Å long in both X and Y directions. After alignment of AQP5 structure, the lipid molecules in the range of 0.8 Å of AQP5 protein are deleted to prevent overlapping of protein and lipid molecules. The merged AQP5 and POPE bilayer are surrounded by two 18 Å-thick water layers at the bottom and top of protein-lipid structure. Finally the NaCl ions with 0.15 mol/L concentration are added to the structure to completely model the biological environment. The complete system consisting of 126006 atoms is shown in the Figure 1.

Figure 1

(a) AQP5 in POPE lipid bilayer and water molecules around them from top view. (b) The same system from a side view.

Molecular Dynamics details

To carry out the equilibrium MD simulation of human AQP5, NAMD software package [26] is used. The interaction parameters required for MD calculation, are defined according to CHARMM36 force field [27] for protein, lipid and solvent molecules. The TIP3P model for water molecules is used based on the CHARMM force field. Along with MD algorithm, Langevin dynamics thermostat with a Langevin damping coefficient γ=5 ps−1 is used to maintain the temperature of the system at 300 K. To maintain a constant pressure of 1 bar, the Langevin piston algorithm is employed. Langevin Thermostat is calculated based on the thermostat algorithm applied on thermodynamical ensemble at a constant temperature. Application of the Langevin thermostat and Langevin piston algorithm results in the NPT condition, namely, constant particle number, constant pressure and constant temperature.

Assuming periodic boundary conditions, Particle Mesh Ewald (PME) method are applied to consider finite system size and full long-range electrostatic interactions. For the Van der Waals interactions, a cutoff distance of 12 Å with a switching distance of 10 Å are applied. The time step for MD calculation is 2 fs.

Molecular Dynamics Simulation

To reach an equilibrated system, three different computational phases are performed. First, to equilibrate the POPE lipid layer, all atoms expect lipid tails are assumed to be fixed and the lipid tails are minimized and equilibrated for 0.5 nanosecond. The first phase models the proper disorder of a fluid-like bilayer. At the second phase, minimization and equilibration of the system containing constrained protein is performed for 0.5 nanosecond. The second phase simply guides the system to the nearest local energy minimum in configuration space while the harmonic constrains are applied on the AQP structure. These two phases; minimization and equilibration, cause the lipids to be well packed around the aquaporin. At the third phase the protein is released and the system is further equilibrated. Eventually the system is well equilibrated, and the stability of protein through the MD computation can be detected by monitoring protein RMSD. Finally, 500 nanosecond of NPT molecular dynamics simulation without any restraint is carried out. Also, to model external pressure field in the forced diffusion simulations, a direct force is assumed to be applied to each water molecule [28]. The applied force coefficient is calculated based on the pressure gradient and simulation box dimensions.

Results and Discussion

Aquaporin-5 Channel Radius

The mean pore radius of four AQP5 monomer channels is shown in the Figure 2. The program Hole 2.0 [29] is used to obtain the AQP5 pore size along the channel direction. First, the time average of each monomer pore radius is obtained from the results of MD simulation for a total time of 500 ns. The four monomer pores have similar graphs, so only the mean pore radius is shown in the Figure 2 (b). The indicated blue area in the graph is associated to the Selectivity Filter (SF) region, at the narrowest part of the channel where the water molecules can permeate when SF region is wide, i.e. nearly 1.1 Å pore radius. When the SF radius is less than 1 Å, the water molecules will be restricted in the channel. The red region on the graph indicates the NPA (Asparagine-Proline-Alanine) motif, which provides a high electrostatic barrier that excludes all ion leaks, including the protons located across the bilayer (Fig. 2a). Figure 2 (a) shows one AQP5 monomer with the location of SF region and NPA motif. The five loops constructed of α-helices are represented in different colors.

Figure 2

(a) SF region and NPA motif for an AQP5 monomer, colored by its loop structure. (b) Mean pore radius of four AQP5 monomer channels, calculated by Hole 2.0 program [29]. In the SF region, an arginine residue provides hydrogen bonds which causes water molecules permeation through the selectivity filter. The NPA motif, at the middle of the channel, stabilizes the half helices loops through multiple hydrogen bonds.

Aquaporin-5 Stability Analysis

To analyze the stability of the equilibrium MD simulation, the Root Mean Square Deviation (RMSD) of AQP5 protein tetramer structure versus calculation time is traced. Figure 3 shows the RMSD fluctuation for the AQP5 crystal structure for the first 50 ns of simulation. For each MD time step, RMSD is calculated based on the AQP5 structure configuration. As it can be inferred from Figure 3, at the first step (before 5 ns) considerable fluctuations in the RMSD graph is observed. As the MD simulation progresses, at almost 8 ns, the RMSD fluctuations tend to slow down and the stable equilibrium state for the whole system is achieved which allows for a stable MD analysis.

Figure 3

The RMSD of human AQP5 crystal structure as a function of MD time for free diffusion simulation.

Water Permeability

Figure 4 demonstrates the osmotic water permeability through two monomers of human AQP5 under the equilibration state. This figure shows water transport through AQP5 occurs in a single-file chain, similar to the other aquaporin channels.

Figure 4

The water molecules passing through AQP5 monomers. The AQP5 molecules is represented as a single Surface in VMD.

Moreover, according to the MD results, water molecules cannot pass the middle pore of AQP5. Water molecules can only pass through the four monomers. This observation is in accordance with biochemical structure of the central pore in which a lipid occludes the putative central pore, preventing the passage of water molecules through the center of the tetramer [24]. Figure 5 displays the absence of water molecules in the central pore of AQP5.

Figure 5

Water molecules passing through four AQP5 monomers.

The single channel diffusion permeability constant p’ for each monomer can be calculated from equilibrium MD equilibrium simulation based on the following equation [30];

In which, v is the average volume of a single water molecule and q0 is the number of water molecules that cross the channel per unit time. At the equilibrium MD simulation, q0 can be calculated as the number of water molecules completely crossing the tetramer pores during 500 ns of MD simulation. Furthermore, the osmotic permeability coefficient pf, can be calculated from equilibrium MD simulations using the collective diffusion model [31]. Considering the z-direction displacement of each water molecule located at the channel constriction region, during time step dt, the normalize collective coordinate n(t) is defined as:

where L(t) is the constriction region length, S(t) is the set of the whole water molecules inside the channel constriction region at time t and δt is the time step. Considering the collective diffusion constant defined by Einstein relation;

Also, the osmotic permeability coefficient pf, can be calculated according to the following relation;

It should be noted that 〈n2(t)〉 is the Mean Square Displacement (MSD) of normalized collective coordinate n(t). To calculate the osmotic water permeability coefficient, first collective coordinate n(t) should be evaluated. Equation 3 is employed to calculate n(t) for each AQP5 monomers. The length of pore monomer is important for calculation of n(t). To calculate MSD of the normalized collective coordinate n(t), the averaging is preformed over a 500 time interval for a 500 nanosecond simulation time. Each averaging time interval is 1 ns. Figure 6 shows the MSD of n(t) for each AQP5 monomers.

Figure 6

Mean Square Displacement (MSD) of the normalized collective coordinate for each AQP5 monomers and linear curve fitted on the MSD results.

To calculate single channel diffusion permeability coefficient p, number of water molecules crossing the channel are counted for each AQP5 monomer pore in both positive and negative direction of the channel for over a 500 ns MD simulation time. The total number of water molecules that cross the tetramer channels is 1503. The channel diffusion permeability coefficient p is calculated based on Equation 1. Based on the MD calculation, the value of p is 1.24×10−14 cm3/s. Moreover, the osmotic permeability coefficient p is also calculated based on the Equation 5 which according to the MD results for MSD of n(t) in the human AQP5 (Fig. 5), yields to a p value of 5.11×10−14 cm3/s. An experimental study has reported the value of 5.0±0.4×10−14 cm3/s for p [32]. The result for permeability coefficient calculated over 500 ns, has better accordance with experimental data than the previous MD simulation of water permeation through AQP5 over 200 ns [23]. Therefore, the current MD modeling conducted here for the osmotic permeability coefficient is in the good agreement with existing experimental data for osmotic water permeability through human AQP5. As a result, we can conclude that the ideal single-file model is an appropriate model to describe the water molecules permeation through aquaporins channels.

Pressure Effect

So far, the MD results for osmotic water permeation through human AQP5 have been validated with experimental data, therefore further MD simulation can elucidate the biophysical phenomena related to AQPs performance. To compare the effect of pressure on the water permeation with the water self-diffusion process, an external pressure force is applied on the water molecules. The external pressure gradient is equal to the pressure gradient which exists on sides of cell membrane in vivo condition [28]. Two different cases with two different pressure fields of ΔP/Δz=1 atm/Å and ΔP/Δz=0.5 atm/Å on the water molecules, are assumed. For each case, the MD simulation is carried out for a total time of 20 ns. To stabilize the MD simulation, it is assumed that POPE bilayer membrane is constrained by setting the value of β=1 for the phospholipid membrane molecules.

The RMSD fluctuations of the AQP5 crystal structure for each of two pressure-induced MD simulations, are represented in Figure 7. The same method as explained previously, is applied to calculate RMSD for each MD simulations. As the pressure force is exerted on the water molecules, there is no significant difference between RMSD charts of the two simulations. Also, as it is observed for the equilibrium simulation (Fig. 3), the RMSD chart shows smaller fluctuations as the simulation time elapses. Therefore, after 10 ns the protein structure reaches an equilibrium state.

Figure 7

The RMSD of human AQP5 structure as a function of calculation time.

Figure 8 shows the mean MSD of n(t) for water permeation in human AQP5 tetramer under different physical conditions. As it can be inferred from this figure, applying pressure force on water molecules can considerably increase the water permeation through AQP5 membrane and enhance the biomembrane performance. This is due to the fact that the applied pressure force increases osmotic pressure difference. This result is in consistent with transport phenomenon due to diffusion law which governs the fluid flow from high pressure upstream to low pressure downstream.

Figure 8

Mean MSD of normalized collective coordinate for AQP5 under different water permeation conditions. Linear curves fitted to the MSD results are shown as dashed lines.

In the 500 ns equilibration simulation time and for two different conditions of ΔP/Δz=1 atm/Å and ΔP/Δz=0.5 atm/Å, the total numbers of water molecules that completely cross the AQP5 tetramer channels are 4575 and 3701 respectively. Accordingly, the channel diffusion permeability coefficients for the two conditions are calculated based on Equation 1. Table 1 shows the calculated results for channel diffusion permeability coefficient and the osmotic permeability coefficient for osmotic self-diffusion and force diffusion conditions. Higher pressure gradient results in greater value for permeability coefficients, since a higher pressure gradient is associated to a higher normalized collective coordinate (Fig. 8).

Table 1

MD results for channel diffusion permeability coefficient and osmotic permeability coefficients for three simulation conditions

	Osmotic self-diffusion	Force diffusion ΔP/Δz=0.5 atm/Å	Force diffusion ΔP/Δz=1 atm/Å
pd (cm3/s)	1.24×10−14	2.99×10−14	3.75×10−14
pf (cm3/s)	5.11×10−14	12.60×10−14	15.41×10−14

Forced Diffusion Effect on Pore Radius

As described earlier, Selectivity Filter region is the narrowest part of the aquaporin channels for which its minimum radius can hinder or facilitate the water molecules permeation through the tetramer channels. According to the previous section results, applying pressure force of water molecules considerably increases the permeability coefficients for water molecules passing through AQP5. Apart from the higher momentum which enhances water molecules permeation, the minimum size of SF region can be also affected by applying external pressure field on water molecules. Figure 9 shows the mean minimum radius for AQP5 tetramers calculated from 500 ns MD simulation. The mean minimum radius of SF region is the average value of minimum radius of four AQP5 channels and is calculated for each 5 ns time interval. Obviously, this value for osmotic diffusion simulation is much lower than the value for the forced diffusion simulation.

Figure 9

Minimum radius of AQP5 SF region under different water permeation conditions. The dashed lines show the mean value for each simulation condition.

The results of Figure 9 demonstrate that applying external pressure force on water molecules can affect the minimum radius of SF region up to 20%. This increase in the pore size can be due to the intermolecular interaction between water molecules and AQP5 molecules located at the pore wall. Under forced diffusion condition, the increase in mean minimum radius of SF region along with higher momentum of solvate molecules, result in a greater value for permeability coefficient.

Conclusion

We have studied the water permeability characteristics of crystal structure of the human aquaporin 5 channels within their tetrameric configuration. Since water transport occurs on nanosecond time scales, atomistic molecular dynamics simulations are employed to model the transport phenomena of AQP5 channels. The MD analysis gives the detailed insights on how water conduction occurs within the tetramer. In the current study, the water permeation through AQP5 is modeled under osmotic self-diffusion condition and forced diffusion due to pressure gradient. Two different conditions for pressure gradient are employed in order to study the effect of pressure field on the diffusion rate. The MD simulation for each case is carried out for a 500 nanosecond to reach the equilibration condition.

Initially, the osmotic self-diffusion condition, without any external pressure field, is simulated. For this condition, the osmotic permeability coefficient p for human AQP5, calculated from MD outcomes, is in good agreement with previously reported experimental data. Hence our presumption of ideal single-file fashion model is accurate enough to describe the water molecules permeation through aquaporins monomer channels. Moreover, the MD results show that the water permeability of the human AQP5 is in the range of other human aquaporins family. Besides, there is no water permeation through central pore of AQP5 tetramer which is due to the hydrophobic structure of central pore. Under the condition of forced diffusion due to pressure gradient, applying pressure force on water molecules can considerably increase the water permeation through AQP5 membrane monomers and enhance the biomembrane performance. This is because of the increase in water molecules momentum and size of Selectivity Filter (SF) region of AQP5 nanopores. Based on MD results, a higher pressure gradient results in larger minimum radius of SF region and eventually greater permeability coefficient, as the greater pressure gradient is associated to a higher normalized collective coordinate.

The current study proves that the capability of MD method to model the transport phenomena through aquaporin membrane at the molecular length scale. The numerical results show that the water permeation through AQP5 tetramer and consequently the diffusion rate can be enhanced by imposing adequate pressure gradient. Also, for the first time the permeability coefficients for the forced diffusion condition are precisely calculated using the MD simulation results. Interestingly for the first time, we have found that under the forced diffusion condition, the minimum radius of SF region can be increased up to 20%. Similar MD methods such as what carried out in the current study can lay the groundwork in the understanding of molecular transport of aquaporin biomembranes which would be beneficial to further investigation of other related transport phenomena.