When Peter Agre discovered the first water-conducting channel in human red blood cells, denoted CHIP28 at the time of identification, genomic sequencing revealed that it belonged to the MIP family, which derives its name from the Major Intrinsic Protein, found in lens fiber cell   . Once the water-channel function of CHIP28 was determined, the MIP family was renamed the aquaporin family; CHIP28 was renamed Aquaporin-1, and MIP was renamed AQP0  . Proteins homologous to these water channels are found in all domains of life, which have since been grouped into approximately 30 subfamilies through phylogenetic analysis  . In humans, there are 13 isoforms (AQP0-AQP12), expressed in a tissue-specific manner      . The properties of each aquaporin (selectivity and permeability) correspond to its molecular form.
The general structure for an aquaporin is a homotetramer. Each monomer is made up of 6 alpha-helices that form a barrel-shaped pore, 3 extracellular loops (A, C, and E), 2 intracellular loops (B and D), and cytoplasmic C and N termini   . Two restriction sites regulate the channel. It selects against hydronium ions with highly conserved NPA motifs (Asn-Pro-Ala) located on loops B and E, which form a constriction      . The pore restricts larger solutes through the ar/R constriction site, also known as construction site I, composed of four residues including arginine and aromatic amino acids. The pore diameter at the ar/R constriction site creates an opening 3 Å wide in water-selective aquaporins, allowing water in, but excluding larger solutes     .
Out of the 13 human isoforms of aquaporin channels, 9 are expressed in the human eye  . In the lens fiber, there is a unique isoform that is of special interest. AQP0, which is expressed solely in the mammalian lens fiber, accounts for more than 60 percent by weight of its plasma membrane proteins. Functioning as a junction as well as an aquaporin, AQP0 is crucial for maintenance of the structural integrity of cortical fiber cells, resulting in lens transparency, allowing for light to be focused on the retina.
There are two forms of AQP0 present in a lens. The open configuration functions as a water channel, contributing to the microcirculation of water that nurtures and cleans the lens of waste products     . As the fiber cell matures, AQP0 is modified into the closed configuration. The C-terminal is cleaved through post-transcriptional modification, which causes two AQP0 links to form an intercellular adhesion junction called a thin lens junction    -  . These modifications include phosphorylation at the calmodulin-binding domain  . While the absence of AQP0 from the lens does not prevent interlocking protrusions of young fiber cells, it did cause the loss of function of these protrusions to maintain the structural integrity of the fiber cells, leading to cell separation and cataract formation  . A number of mutations in AQP0 contribute to congenital cataracts in humans and mice by interfering with the protein’s ability to maintain osmotic balance     . Hence, determining the structure of AQP0 has clinical relevance to cataract formation.
The water permeability of AQP0 is 40 times lower compared to AQP1   . This has been attributed to two unique tyrosine residues (Tyr149 and Tyr23). The NPA site in AQP0 is also narrower than in AQP1, resulting from a substitution of Leu151 to Phe141 in loop B  . AQP0 channels are regulated by external pH    , external and internal concentration of calcium ions (Ca2+)    , and phosphorylation. Protonation and deprotonation of His40 has been shown to be caused by the pH regulation of AQP0  . The optimal pH for AQP0 water conduction is at 6.5: an increase in pH from 6.5 to 7.2 decreases water permeability by half    . Calmodulin is a universal regulatory protein that mediates the effects of Ca2+ concentration on the water permeability of AQP0 through allosteric modulation  -  . Calmodulin (CaM) is a small protein (16.7 kDa) with a C-lobe and an N-lobe connected by a flexible linkage     . The binding of four Ca2+ ions to CaM exposes hydrophobic regions that can interact with hydrophobic regions of target proteins    . Recent research on the pseudoatomic model suggests that one CaM can cooperatively bind non-canonically to two open configuration AQP0 monomers and allosterically inhibit the channel by narrowing constriction site II (CSII), a feature absent in other AQPs   . CSII is a located on the cytoplasmic half of the pore and involves Tyr149 extending into the pore and interacting with His66 and Phe75  . The pH and Ca2+ concentration vary inside the layers of the lens, so regulation by these properties may be physiologically significant  . The modifications that affect responsiveness of the fiber cell include phosphorylation at the CaM-binding domain  . An anchoring protein, AKAP2 positions protein kinase A, which phosphorylates AQP0 at Ser235 in the CaM binding site in order to prevent CaM modulation. This may be done in order to preserve fluid circulation in the middle of the lens   . Another study found when either Ser235 or Ser231 was phosphorylated, as it is in posttranslational modification, binding affinity is reduced for dansyl-CaM in a bovine model  .
Previous dynamic simulations suggested that CaM binding to an AQP0 tetramer reduced the dynamics of all AQP0 channels  . Permeability assays carried out on oocytes expressing wild type AQP0 showed that when cytoplasmic Ca2+ increased from 0 to 1.8 mM, water permeability decreases by a factor of 0.4  -  . Another study showed that calcium and protons can regulate AQP0 via single-monomer regulation, as well as in a cooperative fashion. The non-canonical double binding of CaM suggests that the regulation of Ca2+ is through cooperative modulation: the CaM molecule linking the two AQP0 bound monomers to each other  . This leads to the question of whether the dynamic restriction of AQP0 is the only mechanism through which CaM-Ca2+ regulates AQP0’s permeability to water. In this study, we applied various Molecular Dynamics Simulation techniques to investigate the structural and energetic variations among the AQP0 monomers with and without CaM binding.
Molecular visualization programs VMD 1.9.1  and CHARMM36  were used for modeling and visualization. MD Simulations were performed via NAMD2.0  , on local CPU and GPU clusters, the University of Iowa HPC facility, and the University of Texas TACC. Two separate tetrameric models were built with free-AQP0: (PDB ID: 2B6P  ) and CaM-bound AQP0 (PDB ID: 3J41  ). In each model, hydrogen atoms were added to the backbone of an AQP0 tetramer by VMD psfgen. Then the AQP0 tetramer was embedded in 1-palmitoyl-2-oleoyl-sn-glycerol-3-phosphocholine (POPC) bilayer of 150 Å × 150 Å by VMD membrane builder, solvated in TIP3 water molecules using VMD solvate and neutralized with NaCl by VMD autoionize. The final CaM bound AQP0 model had 386,322 atoms, and the open AQP0 simulation model has 284,903 atoms. Both systems were minimized for 2000 steps, and equilibrated for 30 ns before the production runs.
Steered molecular dynamic simulation (SMD) was used to study water conduction energy profile, transport dynamics, as well as determination of the critical residues in the AQP0 pathway. An extracellular water molecule approximately 5 Å away from the opening was aligned with the center z-axis of each channel. In other words, SMD was conducted with four water molecules in a free-AQP0 or a CaM-bound AQP0 model. The water molecules were pulled along the z-axis with a force constant of 1.5 kcal/mol/ Å2 and a constant velocity of 5.0 Å/ns to AQP0 pore. All of the simulations were carried out in a periodic cell of 150 × 150 × 170 Å, at 310 K, and 1 atm. Z-coordinate and force were recorded every 100 fs. The center of mass of a water molecule was defined as the center of mass of the oxygen atom. Instantaneous force was calculated using the equation: F = zi + vt − z, where zi is initial position of the center of mass, v is pulling velocity, t is time and z if position of the center of mass at the point of measurement. Work is the product of the average force between two consecutive data point and the distance between these two data point. The energy profile of each water molecule was calculated based on Jarzynsky’s equality  : exp(−βW) = exp(−βΔG), where β = 1/kBT, kB is the Boltzmann constant and T is temperature. In addition, RMSD was used to keep track of the dynamics of critical residues at constriction site I and constriction site II, especially Tyr149. VMD was used to study water conduction pathway and conformational changes of critical residues.
Umbrella sampling was used to investigate the energy profile of a water molecule in an AQP0 channel. The simulation was carried out for 4 channels simultaneously with a simulation window size of 0.5 Å. First, a water molecule approximately 5 Å directly above each channel was fixed using VMD mdff, and then the systems were minimized for 2000 steps. After that, these water molecules were constrained in z-coordinate, the coordinate along the longitude of the channel, with a force constant of 10 (kcal/mol/Å2). The simulation was conducted at 310 K and 1 atm for 450 ps with data recorded every 500 fs, which resulted in 900 data points. The first 300 data points were used for equilibration and the next 600 data points were used for analysis. The initial position of the water molecule was defined as the fixed position of oxygen during minimization.
HOLE  was used to analyze the channels of AQP0 to confirm the results of the free energy profiles obtained with SMD and Umbrella sampling. The HOLE software creates a visual output and coordinates of the radius of the channel for each monomer. A HOLE profile and visualization was created for each individual monomer, one for the open conformation and another for the closed conformation after over 50 ns of MD simulations. In the end 8 outputs were created. This allowed for the comparison of the channel dynamics when CaM was bound versus when the channel was in an open conformation. Comparisons of the HOLE pathways and SMD trajectories were created by superimposing the channels with their alpha-carbons onto one another using the RMSD tool on VMD for alignment.
RMSD values of the side chains were gathered from the VMD software RMSD Trajectory Tool. Data acquired from equilibration of the two models was used for the analysis. Both of the models were equilibrated and had 200 usable DCD frames over the span of 30 nanoseconds per model.
3. Results and Discussion
3.1. AQP0-Cam System Has Two Distinct Water Conduction Energy Profiles
The various MD simulations applied to CaM bound and unbound AQP0 tetramers revealed three unique energy profiles for water conduction. The first representing the four channels in free-AQP0, the second representing the two channels that were not covered by CaM in CaM-bound AQP0, and the final energy profile representing the two channels that were bound by CaM. For convenience, these three energy profiles are denoted open, closed-open, and closed-closed, respectively. Each of the energy profiles has three prominent peaks.
Figure 1. (Top Left) AQP0 tetramer in the open conformation with no CaM attached to the protein. (Top Right) AQP0 tetramer bound to two CaM. In the closed configuration each CaM is aligned with one of the AQP0 pores. In order to simplify the image, the lipids and water molecules are not shown. (Bottom) Top down perspective of the closed configuration with CaM. In this study, the monomers in dark blue are named as closed-closed, and those in cyan are closed- open.
AQP0 tetramer in the open conformation without CaM. Figure 1―RIGHT shows the two CaM units are bound to the bottom portion of the tetramer. CaM is highlighted in dark blue ribbon below the tetramer, facing outwards towards the intracellular region.
The radius of the open monomer is always at a value above 1 Å (see Figure 2). In contrast with the closed-open and the closed-closed monomers, the channel is very clear and permeable to water. The only slight point of constriction is along Z coordinate 3.53 where the radius dips down to a value of 1.001 angstroms for a short section. This corresponds to a group of residues that constitute the ar/R site (includes His172, Ala181, and Arg187). The degree of constriction along the initial portion of the channel depends on how close His172 is to Ala181 and Arg187. His172 is in a vertical conformation with the channel, allowing water to pass through. Further down the channel, the radius near CSII is at 1.30 angstroms. The residues Tyr149, Phe75, and His66 provide enough room for passage of water molecules to occur.
Next we investigate the changes that occur to the tetramer when two CaM molecules bind to it. The combined use of SMD trajectories and HOLE visuals shows that the path the water molecules take doesn’t change for any of the four monomers when CaM is bound to AQP0 versus when the tetramer is in the open conformation. All of the HOLE
Figure 2. Measured energy profiles for water conduction, and the pore radius calculations (Right) along with the HOLE visualization of the open conformation of an AQP0 monomer with no CaM bound to the system (Left). The black solid line represents the pore radius. The blue-dashed and red-dotted lines represent SMD and umbrella sampling free energy calculations along the pore, respectively.
plots reveal the same pathway. However, the overall conductivity decreases due to the smaller channel radius and less total channel volume. This decrease in conductivity occurs at the second constriction site (CSII), as well as the ar/R site.
Figure 3 shows the major points of constriction along the monomer when it is bound directly to CaM. The radius decreases significantly at several points along the channel. HOLE analysis and radius data give the radius at CSII as only 0.97 Å. At this point in the channel residues Try149, Phe75, and His66 all come together in close proximity to restrict the passage of water. However, the most significant point of constriction in the channel occurs near residues Ala181, Arg187, His172, and Met183. The channel radius decreases to 0.665 Angstroms. This small radius is partly due to His172 being in a horizontal conformation, and the constriction of the ar/R site.
Since only half of the channel is directly bound to CaM, we wanted to see whether or not there are any allosteric changes to the shapes of the two closed-open monomers. To date, there hasn’t been much insight into whether or not the closed-open monomers differ from the closed-closed monomers. Our investigation shows that the closed-open monomers are very restrictive to the passage of water in both the SMD simulations as well as equilibration of the protein along the initial portion of the channel. When comparing between the closed-open and the open conformation, the data makes it clear that there is less permeability in the closed-open, with a higher energy barrier existing along the ar/R site in the closed-open conformation.
Analysis of the closed-open conformation of AQP0 reveals a distinct contrast with
Figure 3. Measured energy profiles for water conduction, and the pore radius calculations (Right) along with the HOLE visualization of the closed-closed conformation of an AQP0 monomer (Middle). The ar/R and CSII residues orientations are shown in the boxes (Left). The black solid line represents the pore radius. The blue-dashed and red-dotted lines represent SMD and umbrella sampling free energy calculations along the pore, respectively.
the open conformation. We predicted the closed-open channel would be similar to the open channel, and it is more permeable at CSII when compared to the closed-closed conformation. Figure 4 shows that the energy barrier around residues His66, Phe75, and Tyr149 is around 2.5 kcal per mole, in comparison to the 5 kcals per mole the closed-closed conformation had at CSII. So it appears that there may indeed be a lesser degree of allosteric modulation that occurs at CSII in the closed-open conformation. Since it is not directly covered by calmodulin, there may be less strain on the protein on the two closed-open conformations, and the residues at CSII may not be brought into as close of proximity as the closed-closed monomers bound directly to CaM.
During SMD trajectories, it has been noted that His172 serves as a sort of “plug” when it is oriented horizontally. The water molecule spends a large portion of time above these residues until the orientation of His172 shifts into a vertical position, allowing for the passage of the water to the rest of the channel. This gating mechanism is visible in Figure 5. Note how the water molecule is blocked initially, and in the second conformation the molecule was allowed passage through the ar/R site. In the closed-closed and the closed-open conformations the radii at this point are 0.665 and 0.95 Angstroms, respectively. The discrepancy between the last two conformations is due to the fact that His172 was vertical in the closed-open monomer. The fluctuation between His172 between a horizontal and vertical orientation may have an effect on
Figure 4. Measured energy profiles for water conduction, and the pore radius calculations (Right) along with the HOLE visualization of the closed-open conformation of an AQP0 monomer (Middle). The ar/R and CSII residues orientations are shown in the boxes (Left). The black solid line represents the pore radius. The blue-dashed and red-dotted lines represent SMD and umbrella sampling free energy calculations along the pore, respectively.
Figure 5. The dynamic motion of the His172 observed to facilitate the water conduction Steered Molecular Dynamics. The closed-closed (Top) and closed-open (Bottom) monomers show distinct difference.
water permeability between the channels and may be the deciding factor in how restrictive the monomer is at the ar/R site.
3.2. Multiple Residues Take Part in Stabilizing CSII
Figure 6 highlights RMSD data collected from a 30 ns equilibration in both the closed and open conformations. Open RMSD values are shown in orange, and closed values are in blue. The values were taken when compared to the first frame, and the atoms on the whole residue were averaged. All of the listed residues, Tyr149-His66-Phe75 are part of CSII. The end result shows that the His66 and Phe75 shifts more on closed configuration to stabilize the CSII site.
However, His172, Arg187, and Ala181, all involved with the ar/R restriction site showed no difference on RMSD analysis for closed and open models.
Figure 6 shows that that the binding of CaM protein stabilizes the AQP0 tetramer allosterically and results in a shift in position in residues that make up the CSII site of the channel. This result confirms other studies claiming that the main form of regulation of AQP0 by CaM is done via restriction of CSII  .
4. Conclusion and Discussion
Our investigation into the AQP0 tetramer was meant to expand previous studies done and to explore the differences between the closed-closed and closed-open conformations of the protein when CaM is bound. From analysis of HOLE plots and the energy
Figure 6. The RMSD distributions of the three critical residues in CSII site, Tyr149-Phe75-His66, during 30 ns equilibration run. Blue line represents the CaM bound AQP0, Orange line is AQP0 tetramer with no CaM. The coordinates are compared to the first frame, and all four monomers results were averaged.
profiles, there appears to be a significant difference between the closed-closed and the closed-open conformation in the channels. Both of these conformations have the same regions of constriction at CSI (ar/R site), but differential constriction at CSII. Energy profiles showed that the closed-closed conformation reached a peak of around 5 kcal/ mol at CSII, while the value for CSII for the closed-open conformation was around 2.5. Both the HOLE plots and the energy profiles show that the different conformations in the AQP0-CaM complex may be differentially modulated, with the most restriction of water occurring in the closed-closed conformations.
Something of potential interest is the large energy peaks at the top of the channel near the ar/R site as opposed to CSII. In all conformations the energy peak is higher at the ar/R site, even when CaM is bound. The ar/R site appears to have the highest energy barrier for all channels and conformations.
Analysis of the RMSD fluctuations of the side chains involved in the ar/R site and the allosterically modulated CSII showed that the binding of CaM didn’t necessarily reduce the internal thermodynamic fluctuations of side chains in the channel. Two simulations produced conflicting results, and we cannot determine whether or not CaM stabilizes the interior of the channel by reducing these fluctuations. Previous literature has suggested that CaM reduces RMSD fluctuations when it is bound. Our findings are inconsistent with this assertion. But possible interpretation may be that binding of CaM causes the side chains involved in CSII to adopt a certain conformation (a closed conformation) during the course of equilibration rather than having lower fluctuations. Further investigation should be done, especially with longer equilibration periods.
More individual frames should be taken over the course of simulations to produce a variety of HOLE visuals for comparison. This should paint a more accurate picture of whether or not there is a clear difference between these two conformations. It would have also been useful if more SMD trajectories were produced. One problem encountered was that in a few cases the water molecule would venture off target, even though the channel was more permeable in the traditional path of water.
In conclusion, the energy profiles show that CSII has a higher energy barrier in closed-closed as opposed to closed-open. Our preliminary findings indicate that they may be significantly different. We also note the importance of His66 and Phe75 in regards to the reduced radius at CSII when CaM is bound, while previous literature focused only on Tyr149. Future investigation should flesh out and fully characterize these differences, but to date no one has determined if these channels have differential permeability. It is assumed that all channels have the same dynamics in the closed AQP0 tetramer.
The authors give our special thanks to Extreme Science and Engineering Discovery Environment (XSEDE) for the computational time awarded to our project. Dean of Coe College, via The Ella Pochobradsky Endowment Award, funded this research. Authors also would like to give special thanks to the Coe College Chemistry Department for their generous support throughout this project.
*These authors contributed equally to this work.