Molecular dynamics study of the archaeal aquaporin AqpM

Background Aquaporins are a large family of transmembrane channel proteins that are present throughout all domains of life and are implicated in human disorders. These channels, allow the passive but selective movement of water and other small neutral solutes across cell membranes. Aquaporins have been classified into two sub-families: i) strict aquaporins that only allow the passage of water and ii) the less selective aquaglyceroporins that transport water and other neutral solutes, such as glycerol, CO2 or urea. Recently, the identification and characterization of a number of archaeal and bacterial aquaporins suggested the existence of a third sub-family; one that is neither a strict aquaporin nor an aquaglyceroporin. The function and phylogeny of this third family is still a matter of debate. Results Twenty nanosecond molecular dynamics (MD) simulation of a fully hydrated tetramer of AqpM embedded in a lipid bilayer permitted predictions to be made of key biophysical parameters including: single channel osmotic permeability constant (pf), single channel diffusive permeability constant (pd), channel radius, potential water occupancy of the channel and water orientation inside the pore. These properties were compared with those of well characterized representatives of the two main aquaporin sub-families. Results show that changes in the amino acid composition of the aromatic/arginine region affect the size and polarity of the selectivity filter (SF) and could help explain the difference in water permeability between aquaporins. In addition, MD simulation results suggest that AqpM combines characteristics of strict aquaporins, such as the narrow SF and channel radius, with those of aquaglyceroporins, such as a more hydrophobic and less polar SF. Conclusions MD simulations of AqpM extend previous evidence that this archaeal aquaporin exhibits hybrid features intermediate between the two known aquaporin sub-families, supporting the idea that it may constitute a member of a novel class of aquaporins.


Background
Aquaporins are a large family of transmembrane channel proteins that are present throughout all domains of life and their malfunction has been implicated in several human disorders [1]. Aquaporin channels allow the passive but selective movement of water and other small neutral solutes such as glycerol, CO 2 , or urea across cell membranes [1][2][3][4].
Structurally, Aquaporins present a homotetrameric organization in which each monomer forms an individual functional pore. The canonical fold of the aquaporin monomer is characterized by a right-handed helical bundle of six transmembrane α-helices (TM1 -TM6) connected by 5 loop regions (loops A to E), in which both amino and carboxyl termini face the cytoplasmic side of the membrane. Loops B and E are formed by a half-membrane spanning helical section (HB and HE respectively) and a non-helical section that contains the highly conserved Asn-Pro-Ala (NPA) motif, considered a signature of this protein family [5][6][7][8][9][10].
The aquaporin channel posses two major constriction zones, the NPA region, located at the centre of the pore and the selectivity filter (SF) containing the aromatic/ arginine region, (ar/R), located~8 Å above the NPA region extending towards the extracellular side of the channel. The ar/R region is formed by a residue from helix TM2 (H2 position), a residue from helix TM5 (H5 position) and two residues from loop E (LE1, LE2 position respectively) [11,12]. This region plays a major role in solute selectivity and permeability [13].
Different permeability to water and other solutes and variations in the amino acid composition of the ar/R region has led to the classification of aquaporins into two main functional sub-families [14,15]. The first are strict aquaporins, that only allow the passage of water molecules, in which the ar/R region contains a Phe residue from TM2 (H2 position), a His residue from TM5 (H5), an Arg residue from the loop E (LE2) and a fourth residue from loop E (LE1) that provides a backbone carbonyl oxygen, usually a Cys, Thr or Ala. In these aquaporins, the His and Arg residues of the SF are believed to provide donor hydrogen bonds for water molecules [14,15]. Examples of strict aquaporins are: AqpZ from Escherichia coli [16,17], and the human Aqp0 [3], Aqp1 [7] and Aqp4 [18]. The second are aquaglyceroporins, less selective aquaporins that can transport water and other solutes, such as glycerol, urea and other uncharged small molecules [19,20]. In these aquaporins, the Phe residue at position H2 is replaced by Trp residue, the His residue at position H5 is replaced by Gly residue and the small residue at position LE1 is replaced by a Phe, giving rise to a wider SF and a so called "hydrophobic corner" formed by Trp (H2) and Phe (LE1) residues [3,[20][21][22][23]. Similar to strict aquaporins, position LE2 is occupied by an Arg residue that is highly conserved in both aquaporin sub-families. Examples of aquaglyceroporins are: GlpF [19][20][21]23] from E. coli, and human aquaporins Aqp7, Aqp3, Aqp9 and Aqp10 [3].
The recent identification of archaeal aquaporins that exhibit a different amino acid composition of their SF, and thus cannot be readily classified in either of the two aforementioned aquaporin sub-families, has led to the suggestion that a third functional sub-family of aquaporins may exist [24,25]. It has been suggested that this family could contain special adaptations for the conduction of solutes specific for the life-style of the organisms or, alternatively, that they could be a primitive non-specialized aquaporin, an ancestor of the more specialized aquaporins [2,24,26]. The SF of these archaeal aquaporins exhibits the highly conserved Arg residue in position LE2, the Phe residue at position H2, and a Ser residue in position LE1 that provides its main chain carbonyl oxygen. However a key difference in the composition of the SF of these aquaporins lies in the presence of a medium size hydrophobic residue such as Ile, Leu or Val instead of the His residue highly conserved in strict aquaporins or the Gly residue present in aquaglyceroporins in position H5 [24,25].
To date, the only representative of this type of aquaporins that has been experimentally studied is AqpM, found in the archaeon Methanothermobacter marburgensis [24,25]. The osmotic water permeability constant, P f , obtained for AqpM, showed that its osmotic permeability was lower than the P f of the strict aquaporins AqpZ from E.coli, and slightly higher than the P f of the aquaglyceroporin from E. coli GlpF [24,25]. (Pf (AqpM) = 57 μms-1 [24]; Pf (AqpZ) = 330 μms-1 [27]; Pf (GlpF) = 49 μms-1 [27]). Moreover, the transient glycerol permeability measured for AqpM resulted in values considerably lower than the glycerol permeability for GlpF [24]. These previous results have led to the suggestion that AqpM could be a water transporter and also with the ability to transport other neutral solutes required by the archaeon such as CO 2 , the only carbon source available in the natural environment of M. marburgensis, or even H 2 S, the terminal electron acceptor in its energy production pathway [24] [25].
In this work we report the results of the first molecular dynamics (MD) simulation study performed on a fully hydrated model of the AqpM tetramer embedded in a lipid bilayer (as shown in Figure 1). Key biophysical parameters were estimated from the MD simulation, including: single channel osmotic permeability constant (p f ), single channel diffusive permeability constant (p d ), channel radius, potential water occupancy of the channel and water orientation inside the pore.
The results of this study extend the information regarding the permeability and selectivity mechanisms of AqpM and provide further insight into the water permeability of aquaporins in general and the particular differences between the sub-families of aquaporins and the implications of these differences for their function.

Results and discussion
Single channel osmotic water permeability constant The single-channel osmotic water permeability constant ( p f ) is a key parameter that allows the characterization of the water transport mechanism of aquaporins [28,29]. An estimation of p f can be obtained from equilibrium MD simulations and can provide information about the water permeation mechanisms of aquaporins at the atomic level [28,30]. Here, we report the calculation and analysis of p f obtained from a 20 ns simulation of the AqpM tetramer. The calculation of p f is based on the collective diffusion model for water permeation through microscopic channels developed by Zhu et al. [30].
Initially, a trajectory for the collective coordinate n(t) was derived (Eq. 5) and the corresponding mean square displacements n (t) 2 were calculated (Eq. 6). Figure  2a displays the trajectory of the collective coordinates n (t) for each monomer of AqpM in the 20 ns simulation and Figure 2b shows n (t) 2 computed by averaging 200 100-ps time windows, each one considered as a time origin (i.e n(t = t`) = 0). The length of the channel is a necessary parameter for computing n(t) and was defined as the average length of the constriction region (CR) of the AqpM channel, L = L(t) = t 19.9¯, considered as the average distance between the atoms G195 (N182/G198):O † and C79(G60/A65):O, over the 20 ns simulation. The single-channel osmotic water permeability constant p f was obtained from n (t) 2 via D n (Eq. 6 and Eq. 7) and yields a value ranging between 5.6x10 -14 cm 3 s -1 and 8.3x10 -14 cm 3 s -1 among the four monomers. Averaging over the four monomers, a value of (6.4±1.4) x10 -14 cm 3 s -1 is obtained. Additional details are provided in Table 1.
Experimental measurements of the osmotic permeability of AqpM have been reported in terms of total membrane osmotic permeability constant (P f ) obtained from AqpM reconstituted into liposomes with values of 57 [24] and 60 μms -1 [26]. In order to compare these values to the single-channel water osmotic permeability (p f ) values obtained in this study, an estimation of p f from P f was made using the method of Borgnia et al.  (b) Mean-square displacement of n <n 2 > (eq. 7) for each AqpM monomer. Each color represents an AqpM monomer (black=M1, red=M2, green=M4, blue=M4). [16] to estimate the p f of AqpZ [24,26]. The estimation consists in dividing P f by the number of channels per unit area (channel density). According to this, an experimentally derived estimate of the p f of AqpM is ≈ 0.7x10 -14 cm 3 s -1 . This suggests that the p f value obtained from the AqpM MD simulation presented here ((6.4±1.4) x10 -14 cm 3 s -1 ) overestimates the osmotic permeability of AqpM by a factor of~9. Interestingly this overestimation of p f from MD simulations has already been observed for AqpZ and other strict aquaporins and most notoriously for GlpF [28,31]. In tables 2 and 3 we present results for water permeability obtained for AqpM, Aqp1, Aqp4 and GlpF from MD simulations studies and experiments, respectively. The results in these tables show an apparent inverse relationship between experimental and MD simulation derived results. GlpF has been shown experimentally to be a slow water transporter. However in MD simulation studies it appears to be a fast water transporter, with p f values equal or greater than those obtained for prototypical and well studied strict aquaporins like AqpZ, Aqp1 or Aqp4 [28,31,32]. The cause of this discrepancy has been discussed and is attributed mainly to assumptions related to the forcefield and water model used for MD simulations and on the difficulty to determine accurately the number of active channels present when estimating p f experimentally [28,31].
An examination of the published MD simulation results obtained for representative aquaporins from both sub-families, including AqpZ [28,31,32] and GlpF [28,[31][32][33] from E.coli, and mammalian Aqp1 [29,31,32] and Aqp4 [31,34], revealed significant discrepancies on the values obtained for the osmotic and diffusive water permeability and the p f /p d ratio for the same aquaporin between different studies (see Table 2). These discrepancies may arise from the use of different force-fields, water models, simulation protocols, simulation times, and insufficient sampling due to the limited time scale of all-atom MD simulations [31]. This issue together with the fact that only one representative from the aquaglyceroporins sub-family, GlpF from E. coli, has been subjected to MD-simulation studies, whereas four representatives from the strict aquaporins have been widely studied (Aqp0, Aqp1, Aqp4, AqpZ) exacerbates the problem of aquaporin sub-family comparisons in terms of water permeability and thus difficult the interpretation of the data obtained for newly described aquaporins such as AqpM in terms of water permeability. In an attempt to overcome these difficulties and to establish a clear point of reference, the present study of AqpM was performed using simulation protocols and analyses methods similar to those previously used in the study of AqpZ and GlpF [28,35], including: i) the use of the same force-field (CHARMM22 for protein and CHARMM27 for lipids), and water model (TIP3P), ii) similar simulation protocol (NPT ensemble, Periodic Boundary Condition, no constraints applied in the production run), iii) similar simulation time (20ns) and temperature (310K) and iv) the same MD engine (NAMD) (see methods). The comparison results in the    following trend: indicating that the single channel osmotic water permeability of AqpM presents a value of p f that is intermediate between those obtained for AqpZ and GlpF. As a whole, this evidence suggests that AqpM exhibits water osmotic permeability that is intermediate between strict aquaporins and aquaglyceroporins and thus could belong to a third sub-family of aquaporins.

Single-channel diffusive water permeability constant
The single channel diffusive water permeability constant p d is another key feature that describes the water permeation mechanism of aquaporins, and is related to the number of water molecules that traverse the channel per unit time (i.e the number of permeation events) [28,36]. In the 20 ns simulation reported here, a total of 69 permeation events were observed (N ± ), with N + = 35 and N -= 34. The accumulated number of permeation events grows linearly with time and the number of permeation events occurring in either direction is as expected, since in an equilibrium MD simulation no net water flow should be present (Fig. 3). Then the unidirectional rate constant (k 0 ) was determined by using: k 0 =N ± /2nmt sim , where N ± is the number of accumulated bidirectional permeation events, t sim is the simulation time discarding the first ns (tsim = 20-1 = 19 ns) and nm represents the number of monomers (nm=4). Using k 0 (k 0 = 0.5ns -1 ) and Eq. 3, a value of p d = 1.4 x10 -14 cm 3 s -1 was obtained for the AqpM tetramer. Additional details and p d estimates for each monomer are provided in Table 1.
As noted for p f , different MD simulation studies of the same aquaporin show significant discrepancies between the values obtained for p d (see table 2). However, despite these discrepancies a trend is apparent, in which GlpF appears with the highest p d values with respect to all other aquaporins including AqpM. The observation that the higher transport rate appears to be in an opposite trend with respect to experimental results can be explained by the wider pore of GlpF that could allow more water-water interchanges and thus more permeation events per unit time, increasing the value of p d (see below). As mentioned above for p f , this fact can be mainly attributed to the nature of MD simulations, and specifically to the water model used, the van der Waals parameters of the forcefield, and non-physical breathing movements of the channel. Following the argument used for the comparison of p f values, we obtained for p d the same trend observed for p f (p d (AqpZ) <p d (AqpM) <p d (GlpF). This result indicates that the single channel diffusive water permeability of AqpM presents values that are intermediate between those obtained for AqpZ and GlpF, using equivalent simulation methods [28]. This observation, together with the results observed for p f , provides the evidence to support the notion that AqpM could be a representative of a third aquaporin sub-family that exhibits properties that are intermediate between those of strict aquaporins and aquaglyceroporins.
According to the continuous time random walk (CTRW) model developed for single-file water transport [37], the p f / p d ratio is related to the number of steps that a water molecule must perform in order to traverse the channel and follows the relation: [37,38] were N is the average number of water molecules inside the channel. Thus the p f /p d ratio can be used as a measure of the "single-fileness" of  * Obtained from P f using the method reported in [16,27]. water movement inside the channel, and for a channel in which the water molecules are perfectly aligned in single file p p < N f d / [28,37,38]. Values of p p < N f d / indicates that the single-fileness of the water molecules inside the channel is interrupted by occasional waterwater interchanges that increase p d , i.e., more water molecules traverse the channel per unit time [28,37,38]. Water-water interchanges do not affect p f because the latter is related to the displacement of individual water molecules inside the channel and is not influenced by the occurrence of water-water interchanges. From the 20 ns simulation of the AqpM tetramer, we obtained 3.0 <p f / p d < 9.1 with an average value over the four monomers of p f / p d ≈ 6.1 and an average channel occupancy equal to N ≈ 7.5 (see below). These results suggest that water molecules move in single file in the AqpM channel with occasional water-water interchange events.
For  28,31]. This indicates that strict-aquaporins, with their narrower pore posses a more idealized single-file water transport mechanism. Interestingly, this relation holds despite the aforementioned discrepancies between the values for p d and p f observed between different MD-simulation studies. Our results shows that, in terms of the p f / p d ratio, AqpM presents a similar level of single-fileness with respect to strict aquaporins (see table 2) but with an average value of p p < N f d / suggesting that AqpM allows more water-water interchanges than strict aquaporins. This could be explained by the slightly wider ar/R region of AqpM due to the replacement of the His residue in position H5 conserved in strict aquaporins by the smaller and apolar Ile residue.

Channel occupancy
Water occupancy histograms from the 20 ns simulation of the AqpM tetramer are shown for each monomer in Figure 4, and the averages are listed in Table.1. Histograms were generated by counting the number of water molecules inside the constriction region (CR) of each AqpM monomer for 10000 snapshots separated by 2 ps taken from the 20 ns MD simulation. The four monomers show a similar behaviour with their occupancy fluctuating closely around the tetramer average water occupancy n N ≈ ≈7.5 (Fig.4). The value of water occupancy reported here for AqpM is in agreement with results obtained from similar MD simulations of other aquaporins (see Table 2) [28,31,35] and for water transport in carbon nanotubes of similar length [37,38] in which a value of N -≈ 7 1 was found.

Channel radius
The channel radius of AqpM was computed using the HOLE program [39]. A total of 10000 configurations separated by 2 ps of each AqpM monomer taken from the MD simulation were considered, yielding an average channel radius of 1.9±0.4 Å for the AqpM constriction region (CR) (-9 Å < z <12 Å). A plot of the radius profile obtained from the 20 ns simulation accompanied by a snapshot of the AqpM constriction region showing the conformation and position of relevant residues are presented in Figure 5a and Figure 5. b, respectively. The channel radius value reported here is in agreement with the channel radius measured from the x-ray structure of AqpM [26]. Comparison of the channel radius profile obtained for AqpM with channel radius measurements with similar methods for other MD-simulated aquaporins, shows that the AqpM channel exhibits an average radius similar to AqpZ ≈ 1.9 Ǻ but is narrower than the average channel radius (≈ 2.5 Ǻ) of GlpF [28,31,32]. As reported for other aquaporins [28,32,35], the minimum radius of the AqpM channel is found at the extracellular side in the selectivity filter (SF) or ar/R region composed of residues R202(R189/R206) † , F62(F43/W48) and I187 (H174/G191) with an average value of 1.4 Å, close to the radius of a single water molecule. This narrowing of the pore at the SF, together with the small fluctuations (RMS) during the MD simulation found on this region (Fig. 5a), suggest that only water and small neutral solutes with similar radius such as CO 2 or H 2 S could potentially pass through the selectivity filter of AqpM [26]. However, it has been shown recently that AqpM is not an H 2 S transporter [25], thus more experimental and theoretical studies of AqpM are required to determine whether this aquaporin can transport other solutes through its channel.

Water orientation
To measure the degree of ordering and orientation of water molecules inside the constriction region of AqpM, two order parameters P (z) = ( ) z 1 cos q and P (z) = ( ) z 2 2 0.5 3cos 1 q − , were measured. In both parameters θ is the angle between the unit vector approximately aligned along the channel axis/membrane normal,n z , and the water dipole vector. Both parameters P 1 (z) and P 2 (z) serve as indications of the average alignment and orientation of the water dipole relative to the channel axis and thus permit the prediction of the water orientation inside aquaporin channels [31][32][33][34][35]. The order parameter P 1 (z) ranges from -1 to 1, and shows the average orientation of the water dipoles with respect to the channel axis ( ) n z . Negative values indicate that the water dipoles are aligned parallel to the channel axis and pointing towards the intracellular side of the channel. Positive values indicate that the water dipoles are aligned parallel to the channel axis and pointing towards the extracellular side of the channel. Values of P 1 (z) near or equal to 0 indicate a perpendicular orientation of the water dipoles with respect to the channel axis. On the other hand P 2 (z), ranges from -0.5 to 1 and positive values indicates an average preferential alignment of the water dipoles parallel to the channel axis ( ) n z , and negative values are an indication of preferential alignment of the water dipole perpendicular tô n z . As can be inferred from an inspection of Figure 6, water molecules appear to be highly ordered within the constriction region (CR) of the AqpM pore and a bipolar orientation of the water dipoles with a dipole inversion at the NPA region can be observed, where both a minimum of P 2 (z) (Fig. 6b) and values near 0 for P 1 (z) (Fig. 6a) are found, indicating that water dipoles are oriented perpendicularly to the channel axis in this region. Water molecules are aligned parallel to the channel axis on both sides of the NPA motifs with their dipoles pointing towards the NPA region. It has been suggested that the bi-polar orientation of water inside the pore of aquaporins is related to the blocking of proton passage through the channel and that could involve a hydrogen donor/hydrogen acceptor pattern that probably induce the dipole inversion around the NPA motifs [40]. The nature of this bi-polar water orientation has been a matter of debate [41][42][43][44] and the current hypothesis is that the bi-polar water orientation itself is not part of the proton blocking mechanism but rather is a signature, an indirect effect of the electrostatic field generated by the two macro-dipoles of the helical part of the re-entrant loops B and E. These macro-dipoles produce a concentration of positive charge in the NPA region of the channel that generates an energetic barrier for protons blocking their passage through the pore [44].
A bi-polar water orientation is clearly observed inside the AqpM constriction region in the course of our MD simulation (Fig.5b and Fig.6) and the hydrogen donor/ hydrogen acceptor pattern appears to be formed by: i) the main chain carbonyl groups exposed to the channel lumen at both sides of the NPA region provided by residues G195(N182/G199), S196(T183/F200), S197(S184/ A201) on the extra-cellular side, and residues H80(H61/ H66), C79(G60/A65) on the intra-cellular side, that could act as hydrogen bond acceptors forming an hydrogen-bond ladder within the pore and ii) the Asn residues at the NPA region N199(N186/N203) and N82 (N63/N68) that act as hydrogen bond donors. An increase in the values of P 1 (z) and negative values of P 2 (z) are observed in the SF of AqpM, indicating that water molecules in this region orient their dipoles perpendicular to the channel axis. This could reflect the interaction of water molecules with the polar hydrogen atoms Arg202(R189/R206):H η and Arg202(R189/R206): H ε that causes a change in the dipole orientation of water molecules in this region. This interaction appears to be enhanced due to the narrowing of the channel in the SF of AqpM. Thus the unique amino acid composition of the SF of AqpM appears to influence the water orientation in this zone of the channel. Calculation of the same order parameters for other MD-simulated aquaporins indicate that, in the case of strict aquaporins like AqpZ and Aqp1, water molecules assume a more perpendicular orientation of their dipoles with respect to the channel axis around the SF due to the narrowing 0.5 3cos 1 q − . θ is the angle between the unit vector, approximately aligned along the channel axis/ membrane normal, n z , and the water dipole vector. Both parameters were measured for each AqpM monomer from 10000 snapshots separated by 2 ps taken from the 20 ns simulation. The black continuous line in each panel represents the average over the 10,000 snapshots and over the four monomers. The error bars represent the standard deviation from the average. of the pore and the strong interaction between water molecules and the Arg (LE2 position) and His (H5 position) residues that comprise the SF of strict aquaporins [28,31,32]. On the other hand, the wider and hydrophobic SF of the aquaglyceroporin GlpF allows more waterwater interactions and thus favours a water dipole orientation parallel to the channel axis [31,32,35].  [47,48]. Fig.7a-c compares the view from the extracellular side of AqpM ( Fig.7a) with that of AqpZ (Fig.7b) and GlpF (Fig.7C); the latter were selected as prototypical representatives of the strict aquaporin and aquaglyceroporin sub-families respectively.

Structure / sequence comparisons
As can be observed in Figure 7(a-d), there are significant differences between the compared aquaporins, regarding the residues that comprise their respective SFs and their spatial arrangement (residues highlighted in blue in Fig. 7d). Of particular interest, is the residue that occupies position H5 (TM5). In AqpZ, b-Aqp1 and h-Aqp4 this position is occupied by a His residue conserved among strict aquaporins [8] (Fig.7d). In GlpF and other aquaglyceroporins position H5 is occupied by a Gly, a feature that allows the formation of the so-called "hydrophobic corner" by Trp48 (replaced by a Phe residue in AqpM and in strict aquaporins) in position H2 and Phe200 (a Ser residue in AqpM, a Tyr residue in AqpZ, a Cys residue in b-Aqp1 and an ala residue in h-Aqp4) in position LE1 that allows the accommodation of a glycerol molecule [19,20,49]. In AqpM, position H5 is occupied by Ile (I187), whereas in the other archaeal aquaporins it is occupied either by Ile, Val or Leu [24,26] generating an SF that is narrower than the SF of aquaglyceroporins but wider and less polar than the SF of strict aquaporins. Thus, the SF of AqpM can be viewed as a slightly wider and more-hydrophobic version of the SF from strict aquaporins. Regarding the hydrophobic residues that line the pore below the SF (residues highlighted in red in red), it is interesting to note that among the compared aquaporins these residues correspond to Val, Ile or Leu, except for AqpZ that has two aromatic residues instead in two of these Horizontal black lines labeled TM1 to TM6 indicate the consensus transmembrane helical regions of the five aquaporins. HB and HE correspond to the half-helices from the re-entrant loops B and E, respectively. NPA motifs are highlighted in green, residues belonging to the selectivity filter (positions H2, H5, LE1 and LE2) are highlighted in blue and hydrophobic residues that line the pore below the SF are highlighted in red. nine positions. These considerations underscore the notion that differences in water permeability between aquaporins are most likely related to the size and polarity of the SF resulting from differences in critically located amino acids. AqpM, with its unique SF, clearly differs from the two other sub-families, exhibiting characteristics that are hybrid between an aquaglyceroporin such as GlpF and a strict aquaporins like AqpZ. Whether this differences permits AqpM to carry out different functions or whether it is just a "molecular fossil" representing an early stage in the evolution of the other aquaporins, as previously suggested [24,26,50], is an issue that remains to be determined. However, the apparent absence of aquaporin genes in a large number of microbial genomes [49], together with the assertion that given the small size and resultant large surface-tovolume ratio of individual microbial cells, lipid membranes are sufficiently permeable to support growth without water channels that enhance water permeability [51,52], have opened the debate about possible alternative roles for microbial aquaporins other than simple water transport. These potential roles include osmotic and turgor sensors [50], involvement in freeze tolerance mechanisms [51] and transport of other neutral solutes such as CO 2 , H 2 S, NH 3 [51].

Conclusions
We have performed the first molecular dynamics simulation of the archaeal aquaporin AqpM, the most studied representative of a new group of aquaporins that are not easily classifiable within the current functional sub-families of aquaporins. From our MD simulation we have provided estimations for key biophysical parameters of AqpM including: single channel osmotic permeability constant (p f ), single channel diffusive permeability constant (p d ), channel radius, potential water occupancy of the channel and water orientation inside the pore. The results obtained from the MD simulations of AqpM were compared with those obtained by similar simulation methods for the well characterized microbial representatives of the two main aquaporin sub-families, AqpZ and GlpF, a strict aquaporin and an aquaglyceroporin respectively. Our results extend the evidence that the size and polarity of the residues that comprise the SF of aquaporins control much of the selectivity and water permeability of aquaporins and supports the notion that AqpM belongs to a third class of aquaporins that exhibits properties that are intermediate between those of the other two classes. However, it is clear that more experimental and theoretical evidence is needed to assess whether this third class of aquaporins differs from the other sub-families in other functional aspect besides water permeability rates i.e., if it exhibits permeability to other uncharged solutes, related to the specific life style of the organisms in which they are found or whether it could represent an earlier evolutionary form, a molecular fossil, which existed before the divergence of the other two subfamilies.

Modelling and simulation
The crystal structure of the AqpM monomer was obtained from the Protein Data Bank [PDB :2F2B] [26]. Hydrogen atoms were added using the psfgen plug-in from VMD v1.86 [53] (considering pH = 7.5 for protonation states). The tetrameric structure of AqpM was generated with VMD v1.86 [53] using the coordinate transformation matrices provided in the PDB file. The simulated periodic cell was constructed using VMD v1.88 [53] and comprised the AqpM tetramer embedded in a pre-equilibrated 1,2-dipalmitoylphosphatidylcholine (DPPC) membrane patch of dimensions 100x100x58 Å. The energy of the crystallographic waters of the monomer was evaluated with the program DOWSER [54], keeping the waters that were located inside the pore. The system was solvated with a water layer of 20Å above and below the membrane. Ions (Na+, Cl-) were randomly placed to neutralize the system, reaching a final concentration of 50mM. The final dimensions of the periodic cell were 100x100x96Å comprising a total of 85,019 atoms. The full system and the structural features of AqpM can be seen in Figure 1.
The system was minimized and subjected to MD for 0.5 ns with fixed protein. The protein was released keeping Cα atoms constrained with a force constant of 5 kcal/mol/Å 2 . The full system was minimized and a slow relaxation procedure was performed in which the constraint applied to Cα atoms of the protein was decreased at a rate of 0.5 kcal/mol/ps until no constraint was applied. Then 21 ns of NPT-MD simulations were performed with the first nanosecond considered equilibration and the last 20 ns used for analysis. The time evolution of Cα -RMSD with respect to the final conformation of the minimization/relaxation procedure is shown for each AqpM monomer in Figure S1(a) of Additional file 1. The fluctuations of Cα atoms during the production run are presented a Cα -RMSF plot for each AqpM monomer in Figure S1(b) of Additional file 1.
The program NAMD [55] with CHARMM27 parameter set [56][57][58] was used for the simulation. Periodic boundary conditions were imposed and the particle mesh Ewald method [59,60] was used for electrostatic forces calculation. Constant temperature (310K) and pressure (1 atm) were maintained by using Langevin dynamics [61].

Single channel water permeability constants
The key quantities that characterize transport properties of a water channel are the single-channel permeability constants: the osmotic permeability constant p f and the diffusive permeability constant p d (measured in cm 3 s -1 ) [28][29][30]37]. In dilute solutions these constants are related to the fluxes j s and j tr due to solute (s) and tracer (tr) concentration differences (ΔC) respectively.
Using equilibrium MD simulations both p f and p d can be calculated. From the total number of complete permeation events (i.e the number of water molecules that traverse the channel lumen per unit time) k 0 [28,29]p d can be computed using: Where v w =V w /N A is the average volume of a water molecule (2.99x10 -23 cm 3 ).
For the estimation of p f from equilibrium MD simulations the collective diffusion model for water permeation through microscopic channels proposed by Zhu et al [30] was used. In this model, a collective coordinate can be defined by accumulating at time t the individual displacement along the channel axis (dz i ) of all water molecules inside the constriction region (CR) of the channel of average length L(t) t relative to t-δt. Thus: By setting n=0 at t=0, n(t) can be determined by the integration of dn: