Molecular Dynamics Simulations of Mite Aquaporin DerfAQP1 from the Dust Mite Dermatophagoides farinae (Acariformes: Pyroglyphidae)

Abstract

Aquaporins are a large family of transmembrane channel proteins that facilitate the passive but highly selective transport of water and other small solutes across biological membranes. House dust mite (Dermatophagoides farinae) is the major source of household immunogens, and we have recently reported six cDNA sequence encoding aquaporins from this mite species. To better understand the structure and role of mite aquaporin, we constructed a tertiary structure for DerfAQP1 by homology modeling from the X-ray structure of malaria aquaporin PfAQP (Protein Data Bank code No. 3C02) and conducted molecular dynamics simulation. The simulation arranged seven water molecules in a single file through the pores of the DerfAQP1. Further, two conserved Asn-Pro-Ala motifs were located on Asn203 and Asn77; residues Arg206, Trp57, Met190, Gly200, and Asp207 constituted an extracellular vestibule of the pore; and residues His75, Val80, Ile65, and Ile182 constituted the cytoplasmic portions. The overall free energy profile for water transport through DerfAQP1 revealed an energy barrier of ~2.5 kcal/mol. These results contribute to the understanding of mite physiology and pathology.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

House dust mites are a worldwide health problem and are a major risk factor for asthma, rhinitis, dermatitis, and other allergic disorders. More than 15%–20% of the populations in industrialized countries are affected by house dust mite sensitization [1]. The two most important elicitors of mite allergies are house dust mites Dermatophagoides pteronyssinus and D. farinae, which are both members of the Pyroglyphidae family. The prevalence of these two major species of mites depends on the geographical region, with D. pteronyssinus as the most prevalent species in Europe, North and South America, Africa, and Australia. Meanwhile, D. farinae is often misleadingly described as the American house dust mite, although it is dominant only in the northeastern parts of North America [2]. Chemicals have been tested in the laboratory for their ability as acaricides to kill house dust mites, but most have poor control or 100% mortality at very high concentrations, thus posing a threat to the environment and humans.

Aquaporins (AQPs) are a family of small, integral membrane proteins that play critical roles in facilitating water transport across cellular plasma membranes in response to osmotic gradients. AQP performs fundamental tasks in the physiology of virtually every living organism, including higher mammals, other vertebrates, invertebrates, plants, eubacteria, archaebacteria, and other microbes. Proteins are involved in diverse biological processes throughout the natural world, making them attractive targets for the development of novel drug therapies or inhibitors of membrane water permeability, such as mercury, gold, silver, copper, phloretin, tetraethyl ammonium salts, and acetazolamide compounds [3]. Thus, studies of mite aquaporins allow both precise biophysical measurements of their functions and exploration of high-affinity inhibitors for control of these insects.

AQPs are highly conserved proteins across diverse species. The completed Human Genome Project acknowledges 13 AQP members (AQP0–12), which are divided into three subgroups based on their structure, function, cellular location, and permeability: (1) orthodox AQPs (AQPs 0, 1, 2, 4, 5, 6, and 8), which have exclusive water permeability; (2) aquaglyceroporins (AQPs 3, 7, 9, and 10), which are permeable to water, glycerol, and possibly other small polar molecules; and (3) superaquaporins (AQPs 11 and 12) [4]. More classes of AQPs can be found in older lineages of mammals (Metatheria and Prototheria), including AQPs 13 and 14 [5]. The construction of an insect phylogenetic tree by using clustalX2.1, MEGA7.0, Jalview, and Mesquite software classified 45 insect AQPs into four major paralogs, with two conserved Asp–Pro–Ala (NPA) motifs located in the central and C-terminal domains [6].

The AQP primary sequences contain six membrane-spanning segments (TM1–6) with five connecting loops (A–E). Each AQP monomer consists of an extracellular and cytoplasmic vestibule connected by a central amphipathic pore containing hydrophilic and hydrophobic surfaces. The hydrophilic surface comprises alpha-carbonyl groups from the polypeptide backbone, including NPA motifs, and the hydrophobic residues of the pore prevent the permeation of larger molecules. The extracellular vestibule contains an ar/R constriction region or selectivity filter (SF), which comprises aromatic and arginine residues and is also the narrowest region of the pore [7].

Three-dimensional (3D) structural research provides further insights into the atomic-level mechanisms of water and solute conduction and of proton and/or ion exclusion. However, the AQP structure is classically difficult to characterize because purification and crystallization are impeded in a membrane environment. Further, no current experimental methods can monitor permeation through AQPs on a molecular level with sufficient spatial and temporal resolution. Molecular dynamics simulations (MDS) can go beyond static structures, thus offering an ideal method to understand the working mechanism of water and glycerol transport through AQPs at atomic resolution. On the basis of our reports on human AQP4 [8], PfAQP from the malarial parasite Plasmodium falciparum [9], and Aqy1 from Pichia pastoris [10], the setup for these simulations typically involves a box containing the protein embedded in a simple membrane surrounded with water.

Elucidating the structure of AQPs provides a unique opportunity to understand the mechanisms of these biomolecular systems. By using RNA sequencing and RT-PCR, we identified six AQPs from the common house dust mite Dermatophagoides farinae, namely, DerfAQP1–4, DerfAQP5.01, and DerfAQP5.02 (GenBank Nos. KY231248, KY231249, KY231250, KY231251, KY231252, and KY231253) [11]. To better understand its structure and role in mite physiology and pathology, we constructed the tertiary structure and conducted MDS for DerfAQP1.

2. Materials and Methods

2.1. System Preparation

Templates of the DerfAQP1 (GenBank Accession No. ATL75744.1) protein were analyzed for homology modeling in the Protein Data Bank (PDB) [12] by using Protein BLAST. The appropriate templates were selected based on high score, low $e$ value, and maximum sequence identity. The tertiary structure of DerfAQP1 was constructed with Modeller 9.0 (https://salilab.org/modeller). PROCHECK [13] was employed to determine the stereochemical quality of the DerfAQP1 protein structure by analyzing the residue-by-residue geometry and overall structural geometry. VERIFY 3D [14] was used to detect the compatibility of an atomic model (3D) of DerfAQP1 protein with its own amino acid sequence (1D). ERRAT [15] was used to analyze the statistics of nonbonded interactions between different atomic types. QMEAN [16] was used to derive both global (the entire structure) and local (per residue) absolute quality estimates based on a single model. The GMQE [17] score was expressed as a number between 0 and 1, which reflects the expected accuracy of a model built with that alignment and template, and the coverage of the target wherein higher numbers indicate higher reliability. The interactive visualization and analysis of molecular structures was completed using UCSF Chimera 1.13.1 [18].

We used the automatic protein structure file (PSF) builder plugin of the VMD software to generate the PSF and PDB files for the monomeric DerfAQP1 and crystallographic water molecules and ions. The protonation states of all residues were selected to conform to the physiological pH range of 7.4 according to the pKa values of their side chains using PSF builder. To solvate the protein, we used the SOLVATE program [19] to fill any empty space inside the pores and to surround the protein with water. By using the VMD software and our previous scripts [20], the all-atom model system for the DerfAQP1 monomer was formed approximately along the $z$-axis [21] at the center of an $80\times 80\times 100$ ų water box. We then added 52 Na⁺ and 55 Cl^- ions to neutralize and salinize the system to 150 mM, with a total of 59,914 atoms.

2.2. DerfAQP1 Simulations

Simulations were performed using the Chemistry at Harvard Macromolecular Mechanics 36 (CHARMM36) force field [21, 22], TIP3P water model [23], and molecular dynamics program NAMD version 2.9 software [24] with periodic boundary conditions at a constant temperature of 298 K and constant pressure of 1 atm (NPT ensemble) to conform to standard laboratory conditions. The vDW interactions were calculated with a cutoff distance of 1.0 nm, switching distance of 0.9 nm, and pair-list distance of 1.2 nm. The Langevin dynamics and piston methods [25] were used to maintain constant temperature and pressure. Full electrostatics were employed with the particle mesh Ewald method [26]. Time steps of 2.0 and 1.0 fs were used for long-range electrostatic and for short-range nonbonded and bonded interactions, respectively. The energy of the bulk water molecules of the system was minimized, while protein and ions were fixed. Next, we fixed alpha-carbons on the central parts of the helices $\left(-10\text{Å}<z<10\text{Å}\right)$ during the equilibrium molecular dynamics run to preserve the backbone structure [20].

2.3. Computation of Water Density in Each Interval through DerfAQP1

To compute the density of water $\rho \left(z_i\right)$ at every interval $\left[z_i,z_i+\Delta z\right]$ along the channel, we first obtained a time average of the amount of water molecules within that interval using the following equation: $$\begin{array}{}\text{(1)}& N\left(z_i\right)=\frac{{\sum }_{j=0}^{t_f/\Delta t}N\left(z_i,t_j\right)}{t_f/\Delta t},\end{array}$$where $t_f=20$ ns is the duration of the simulation and $\Delta t=1$ ps is the time interval between each frame of the trajectory. We then divided by the length of the interval using the following equation: $$\begin{array}{}\text{(2)}& \rho \left(z_i\right)=\frac{N\left(z_i\right)}{\Delta z}.\end{array}$$

Thus, we obtained the number of water molecules per unit Angstrom.

2.4. Free Energy Computation

As stated by Wambo et al. [27], considering that the simulation is performed in the isobaric−isothermal ensemble, the density function of water along the channel can be considered proportional to the probability distribution function of the ensemble: $$\begin{array}{}\text{(3)}& \rho \left(z_i\right)\propto e^{\frac{-\Delta G\left(z_i\right)}{k_{B}T}},\end{array}$$where $\Delta G\left(z_i\right)$ is the Gibbs free energy at point $z_i$. If $\Delta G\left(z_A\right)\equiv 0$ at a reference point $z_A$, the free energy profile along the channel may then be expressed as follows: $$\begin{array}{}\text{(4)}& ∆G\left(z_i\right)=k_{B}T\ln \left[\frac{\rho \left(z_A\right)}{\rho \left(z_i\right)}\right].\end{array}$$

3. Results

3.1. Preparation and Equilibrium of the MDS System

After searching with Protein BLAST, 6QZI, 6F7H, and 3C02 were selected as templates for the DerfAQP1 homology modeling with sequence identities of 35.59%, 35.63%, and 34.89%, respectively, which meet the requirements of homology modeling. The DerfAQP1 protein’s SAVES evaluation results show that its PROCHECK, ERRAT value, and VERIFY 3D results are all ideal, indicating that the results of the DerfAQP1 homology modeling are acceptable (Table 1). The GMQE value of the DerfAQP1 protein is between 0 and 1 and is close to 1, indicating a high reliability of the three protein structures (Table. 1). The QMEAN $Z$-scores showed a good agreement between the DerfAQP1 protein structure and the template protein structure (Table 1). Accordingly, we determined the reliability of the tertiary structure of DerfAQP1 obtained by homologous modeling. By using the automatic PSF builder and Add Solvation Box plugins in the VMD software, the monomer DerfAQP1 was placed in an $80\times 80\times 100$ ų box, and the TIP3P water model and 52 Na⁺ and 55 Cl^- ions were added throughout the system. The system consisted of a total of 59,914 atoms (Figure 1).

Table 1

Evaluation of the tertiary structure of DerfAQP1.

Protein	Structural assessment method	PROCHECK (%)	ERRAT value	VERIFY 3D (%)	GMQE score	QMEAN $Z$-score
DerfAQP1	SAVES	88.9 core	93.25	86.93
		9.6 allowed, 1.2 generously allowed, 0.2 disallowed
	GMQE				0.73
	QMEAN					−3.9

[figure omitted; refer to PDF]

According to the changes of the root mean square deviation (RMSD) of the whole protein, protein without hydrogen atoms, and protein side chains (Figure 2), DerfAQP1 monomers exhibited significant fluctuation at the beginning of simulations but reached a conformational steady state after 30 ns with a RMSD of 3.8–4.2 Å. We saved the files for free energy computation.

[figure omitted; refer to PDF]

3.2. Conduction Channel, SF, and NPA Motifs

Comparison of the radii of DerfAQP1 at the equilibrated state to the crystal structures of bAQP1 (PDB code No. 1J4N), hAQP4 (PDB code No. 3GD8), and PfAQP (PDB code No. 3C02) revealed differences between the DerfAQP1 and PfAQP radii, although the radius of DerAQP1 was similar to that of bAQP1. Through the channel, the narrowest site was located on the SF. The two most conserved NPA motifs formed a canonical “fireman’s grip-like” structure (Figure 3).

[figure omitted; refer to PDF]

3.3. Water Arrangement in the Channel

In the last frame, seven water molecules formed a single file in the channel of the DerfAQP1 monomer (Figure 4(a)). Before the water molecules entered the channel, they stayed in the extracellular vestibule, which contained Arg206, Trp57, Met190, Gly200, and Asp207 (Figure 4(b)). The middle section of the DerfAQP1 channel contained two conserved NPA motifs similar to that in bAQP1 and hAQP4 (Figure 4(c)), resulting in a narrower radius than the channel of PfAQP containing Asn–Leu–Ala and Asn–Pro–Ser [30]. In the channel near the cytoplasm, the water molecules interacted with residues His75, Val80, Ile65, and Ile182 (Figure 4(d)).

[figures omitted; refer to PDF]

3.4. Free Energy Computation

Considering that the single file in the channel of the DerfAQP1 monomer is located between –10 Å and 10 Å $\left(–10\text{Å}<z<10\text{Å}\right)$, the free energy profile showed an energy barrier of ~2.5 kcal/mol with respect to the SF structure (Figure 5) [31]. As we described previously [20], the error associated to our free energy calculation originates from the thermal fluctuation energy $\text{kBT}=0.6$ kcal/mol inherent to molecular dynamics simulations. Since we are subtracting two measurements of energy at different points of the channel, this becomes ±1.2 kcal/mol.

[figure omitted; refer to PDF]

4. Discussion

The current study investigated the working mechanisms of house dust mite AQP to provide information about the future development of AQP inhibitors. MDS showed that seven water molecules were arranged in a single file through the pores of DerfAQP1. As the water molecules trickled through the pores, their polarization changed around the NPA motifs in Asn203 and Asn77. Residues Arg206, Trp57, Met190, Gly200, and Asp207 constitute the extracellular vestibule of the pore, while residues His75, Val80, Ile65, and Ile182 constitute the cytoplasmic portions. Therefore, these key amino acid residues should be the targets for the development of new inhibitors, which can be designed and tested in silico against the docking energies of AQP4 [32].

In 2004, Law and Sansom performed MDS (totaling almost 40 ns) by using the X-ray structures of Aqp1, a homology model of Aqp1 based on the GlpF structure, and the intermediate resolution structures of Aqp1 derived from the electron microscopy [33]. These results suggest that the homology models based on bacterial homologues could be used to derive functionally relevant information on the structural dynamics of mammalian transport proteins. Here, we performed MDS of DerfAQP1 with homology modeling constructed from the X-ray structure of malaria PfAQP (PDB code No. 3C02) [30]. Subsequently, we placed the monomer DerfAQP1 in an $80\times 80\times 100$ ų box and added the TIP3P water model and 52 Na⁺ and 55 Cl^- ions throughout the system. The RMSD change in the last few nanoseconds of the MDS is insufficient to demonstrate our success.

All identified AQP structures have a common homotetrameric arrangement, and each AQP monomer folds into a right-handed alpha-barrel architecture. Our monomer AQP model in simulation achieved structural equilibrium after 30 ns according to the RMSD. Moreover, the free energy profile showed an energy barrier of ~2.5 kcal/mol, and this condition was observed in the SF structure, which was the narrowest part of the pore. Computational free energy profiles deliver valuable structure–function relationships in regard to selectivity. For example, the free energy profiles for the computational mutants of two channels, Aqp1 and GlpF, where the bipolar field was eliminated by artificially discharging backbone atoms [34], suggested that mutating the energy barrier (i.e., SF structure) may effectively destroy mite AQPs.

5. Conclusions

We conducted MDS of the mite AQP DerfAQP1 and computed its free-energy profile for the first time. We are conducting similar simulations for other mite AQPs. The analysis of water molecules through the water channel can further elucidate the AQP structure and functional implications. The structural comparison of mite AQPs, such as mutating the key amino acid residues, can also be helpful in the development of inhibitors for mite control.

Authors’ Contributions

LW conducted homology modeling. LY conducted the molecular dynamics simulations and free energy computation. LW, YZ, MW, FT, and NW conducted the structural analysis. YC devised the experimental plan, interpreted the results, and wrote the manuscript. All authors have commented on the drafts of the manuscript and have read and agreed to the submitted version. Li-lei Wang and Li-li Yu contributed equally to this work.

Acknowledgments

The authors acknowledge the helpful discussions with Prof. Liaoyuan Chen of the Department of Physics at the University of Texas at San Antonio. This study was supported by the National Natural Sciences Foundation of China (Grant No. NSFC31272369), the 333 Project of Jiangsu Province in 2017 (Grant No. BRA2017216), the Medical Innovation Team of Jiangsu Province (Grant No. CXTDB 2017016), the Major Program of Wuxi Health and Family Planning Commission (Grant No. Z201701), the Primary Research and Development Plan of Jiangsu Province (Grant No. BE2018627), the Jiangsu Provincial Medical Youth Talent (Grant No. QNRC2016804), and by the Youth Scientific Research Project of Jiangsu Provincial Population and Family Planning (Grant No. Q201310).

References

[1] J. P. Zock, J. Heinrich, D. Jarvis, G. Verlato, D. Norback, E. Plana, J. Sunyer, S. Chinn, M. Olivieri, A. Soon, S. Villani, M. Ponzio, A. Dahlman-Hoglund, C. Svanes, C. Luczynska, "Distribution and determinants of house dust mite allergens in Europe: the European Community Respiratory Health Survey II," The Journal of Allergy and Clinical Immunology, vol. 118 no. 3, pp. 682-690, DOI: 10.1016/j.jaci.2006.04.060, 2006.

[2] W. R. Thomas, T. K. Heinrich, W. A. Smith, B. J. Hales, "Pyroglyphid house dust mite allergens," Protein and Peptide Letters, vol. 14 no. 10, pp. 943-953, DOI: 10.2174/092986607782541169, 2007.

[3] A. S. Verkman, M. O. Anderson, M. C. Papadopoulos, "Aquaporins: important but elusive drug targets," Nature Reviews Drug Discovery, vol. 13 no. 4, pp. 259-277, DOI: 10.1038/nrd4226, 2014.

[4] K. Ishibashi, S. Hara, S. Kondo, "Aquaporin water channels in mammals," Clinical and Experimental Nephrology, vol. 13 no. 2, pp. 107-117, DOI: 10.1007/s10157-008-0118-6, 2009.

[5] T. Laloux, B. Junqueira, L. C. Maistriaux, J. Ahmed, A. Jurkiewicz, F. Chaumont, "Plant and mammal aquaporins: same but different," International Journal of Molecular Sciences, vol. 19 no. 2,DOI: 10.3390/ijms19020521, 2018.

[6] W. Xia, P. Zhao, Z. Yi, Y. Cui, "Phylogenetic and specific sequence analysis of four paralogs in insect aquaporins," Molecular Medicine Reports, vol. 16 no. 4, pp. 4903-4908, DOI: 10.3892/mmr.2017.7148, 2017.

[7] M. Kourghi, J. V. Pei, M. L. De Ieso, S. Nourmohammadi, P. H. Chow, A. J. Yool, "Fundamental structural and functional properties of aquaporin ion channels found across the kingdoms of life," Clinical and Experimental Pharmacology and Physiology, vol. 45 no. 4, pp. 401-409, DOI: 10.1111/1440-1681.12900, 2018.

[8] Y. Cui, D. A. Bastien, "Water transport in human aquaporin-4: molecular dynamics (MD) simulations," Biochemical and Biophysical Research Communications, vol. 412 no. 4, pp. 654-659, DOI: 10.1016/j.bbrc.2011.08.019, 2011.

[9] Y. Cui, D. A. Bastien, "Molecular dynamics simulation and bioinformatics study on yeast aquaporin Aqy1 from Pichia pastoris," International Journal of Biological Sciences, vol. 8 no. 7, pp. 1026-1035, DOI: 10.7150/ijbs.4703, 2012.

[10] Y. Cui, D. A. Bastien, "Molecular dynamics simulations of PfAQP from the malarial parasite Plasmodium falciparum," Molecular Medicine Reports, vol. 5 no. 5, pp. 1197-1201, DOI: 10.3892/mmr.2012.822, 2012.

[11] J. Peng, Y. Zhou, H. Jia, L. Li, J. Qian, F. Han, H. Yin, Y. Cui, "Transcriptomics-based identification of aquaporin diversity in the house dust mite Dermatophagoides farinae (Acariformes: Pyroglyphidae)," Journal of Insect Science, vol. 18 no. 5,DOI: 10.1093/jisesa/iey099, 2018.

[12] S. K. Burley, H. M. Berman, G. J. Kleywegt, J. L. Markley, H. Nakamura, S. Velankar, "Protein Data Bank (PDB): the single global macromolecular structure archive," Methods in Molecular Biology, vol. 1607, pp. 627-641, DOI: 10.1007/978-1-4939-7000-1_26, 2017.

[13] R. A. Laskowski, J. A. Rullmannn, M. W. MacArthur, R. Kaptein, J. M. Thornton, "AQUA and PROCHECK-NMR: programs for checking the quality of protein structures solved by NMR," Journal of Biomolecular NMR, vol. 8 no. 4, pp. 477-486, DOI: 10.1007/BF00228148, 1996.

[14] J. U. Bowie, R. Lüthy, D. Eisenberg, "A method to identify protein sequences that fold into a known three-dimensional structure," Science, vol. 253 no. 5016, pp. 164-170, DOI: 10.1126/science.1853201, 1991.

[15] C. Colovos, T. O. Yeates, "Verification of protein structures: patterns of nonbonded atomic interactions," Protein Science, vol. 2 no. 9, pp. 1511-1519, DOI: 10.1002/pro.5560020916, 1993.

[16] P. Benkert, S. C. Tosatto, D. Schomburg, "QMEAN: a comprehensive scoring function for model quality assessment," Proteins, vol. 71 no. 1, pp. 261-277, DOI: 10.1002/prot.21715, 2008.

[17] M. Biasini, S. Bienert, A. Waterhouse, K. Arnold, G. Studer, T. Schmidt, F. Kiefer, T. G. Cassarino, M. Bertoni, L. Bordoli, T. Schwede, "SWISS-MODEL: modelling protein tertiary and quaternary structure using evolutionary information," Nucleic Acids Research, vol. 42 no. W1, pp. W252-W258, DOI: 10.1093/nar/gku340, 2014.

[18] E. F. Pettersen, T. D. Goddard, C. C. Huang, G. S. Couch, D. M. Greenblatt, E. C. Meng, T. E. Ferrin, "UCSF Chimera--a visualization system for exploratory research and analysis," Journal of Computational Chemistry, vol. 25 no. 13, pp. 1605-1612, DOI: 10.1002/jcc.20084, 2004.

[19] H. Grubmüller, V. Groll, SOLVATE, 1996.

[20] L. Yu, O. D. Villarreal, L. L. Chen, L. Y. Chen, "1,3-Propanediol binds inside the water-conducting pore of aquaporin 4: does this efficacious inhibitor have sufficient potency?," Journal of Systems and Integrative Neuroscience, vol. 2 no. 1, pp. 91-98, DOI: 10.15761/JSIN.1000117, 2016.

[21] K. Vanommeslaeghe, E. Hatcher, C. Acharya, S. Kundu, S. Zhong, J. Shim, E. Darian, O. Guvench, P. Lopes, I. Vorobyov, A. D. Mackerell, "CHARMM general force field: a force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields," Journal of Computational Chemistry, vol. 31 no. 4, pp. 671-690, DOI: 10.1002/jcc.21367, 2009.

[22] B. R. Brooks, C. L. Brooks, A. D. Mackerell, L. Nilsson, R. J. Petrella, B. Roux, Y. Won, G. Archontis, C. Bartels, S. Boresch, A. Caflisch, L. Caves, Q. Cui, A. R. Dinner, M. Feig, S. Fischer, J. Gao, M. Hodoscek, W. Im, K. Kuczera, T. Lazaridis, J. Ma, V. Ovchinnikov, E. Paci, R. W. Pastor, C. B. Post, J. Z. Pu, M. Schaefer, B. Tidor, R. M. Venable, H. L. Woodcock, X. Wu, W. Yang, D. M. York, M. Karplus, "CHARMM: the biomolecular simulation program," Journal of Computational Chemistry, vol. 30 no. 10, pp. 1545-1614, DOI: 10.1002/jcc.21287, 2009.

[23] W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, M. L. Klein, "Comparison of simple potential functions for simulating liquid water," The Journal of Chemical Physics, vol. 79 no. 2, pp. 926-935, DOI: 10.1063/1.445869, 1983.

[24] L. Kale, R. Skeel, M. Bhandarkar, R. Brunner, A. Gursoy, N. Krawetz, J. Phillips, A. Shinozaki, K. Varadarajan, K. Schulten, "NAMD2: greater scalability for parallel molecular dynamics," Journal of Computational Physics, vol. 151 no. 1, pp. 283-312, DOI: 10.1006/jcph.1999.6201, 1999.

[25] P. Langevin, "On the theory of Brownian motion," Comptes Rendus de l'Académie des Sciences, vol. 146, pp. 530-533, 1908.

[26] U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, L. G. Pedersen, "A smooth particle mesh Ewald method," The Journal of Chemical Physics, vol. 103 no. 19, pp. 8577-8593, DOI: 10.1063/1.470117, 1995.

[27] T. O. Wambo, R. A. Rodriguez, L. Y. Chen, "Computing osmotic permeabilities of aquaporins AQP4, AQP5, and GlpF from near-equilibrium simulations," Biochimica et Biophysica Acta - Biomembranes, vol. 1859 no. 8, pp. 1310-1316, DOI: 10.1016/j.bbamem.2017.04.022, 2017.

[28] W. Humphrey, A. Dalke, K. Schulten, K. VMD, "VMD: Visual molecular dynamics," Journal of Molecular Graphics, vol. 14 no. 1, pp. 33-38, DOI: 10.1016/0263-7855(96)00018-5, 1996.

[29] O. S. Smart, J. M. Goodfellow, B. A. Wallace, "The pore dimensions of gramicidin-A," Biophysical Journal, vol. 65 no. 6, pp. 2455-2460, DOI: 10.1016/S0006-3495(93)81293-1, 1993.

[30] Z. E. R. Newby, J. O'Connell, Y. Robles-Colmenares, S. Khademi, L. J. Miercke, R. M. Stroud, "Crystal structure of the aquaglyceroporin PfAQP from the malarial parasite Plasmodium falciparum," Nature Structural & Molecular Biology, vol. 15 no. 6, pp. 619-625, DOI: 10.1038/nsmb.1431, 2008.

[31] B. Roux, T. Allen, S. Berneche, W. Im, "Theoretical and computational models of biological ion channels," Quarterly Reviews of Biophysics, vol. 37 no. 1, pp. 15-103, DOI: 10.1017/S0033583504003968, 2004.

[32] V. J. Huber, M. Tsujita, T. Nakada, "Identification of aquaporin 4 inhibitors using in vitro and in silico methods," Bioorganic & Medicinal Chemistry, vol. 17 no. 1, pp. 411-417, DOI: 10.1016/j.bmc.2007.12.040, 2009.

[33] R. J. Law, M. S. P. Sansom, "Homology modelling and molecular dynamics simulations: comparative studies of human aquaporin-1," European Biophysics Journal, vol. 33 no. 6, pp. 477-489, DOI: 10.1007/s00249-004-0398-z, 2004.

[34] H. Chen, B. Ilan, Y. Wu, F. Zhu, K. Schulten, G. A. Voth, "Charge delocalization in proton channels, I: the aquaporin channels and proton blockage," Biophysical Journal, vol. 92 no. 1, pp. 46-60, DOI: 10.1529/biophysj.106.091934, 2007.