Does CO₂ Permeate through Aquaporin-1?

INTRODUCTION

Aquaporins are passive integral membrane channels facilitating efficient, yet selective permeation of water across biological membranes (1,2). Also, glycerol and other small linear alcohols are known to permeate via related aquaglyceroporins like GlpF from Escherichia Coli (3). So far, thirteen different aquaporins and aquaglyceroporins with different specificities have been identified in humans. They are expressed in tissues as diverse as the kidney, lung, red blood cells, brain, and the eye lens. Malfunction of aquaporins has been shown to be associated with a number of pathological conditions (4–8). The elucidation of the structure of different aquaporins (9–13) and aquaglyceroporins (14,15) formed the basis for the determination of the molecular mechanism of water permeation and proton exclusion by molecular dynamics simulations (16–21).

It has been a long-standing and controversially discussed question whether also gas transport across biological membranes is facilitated by aquaporins, or whether gases like carbon dioxide can freely permeate lipid bilayer membranes. The role of aquaporin-1 (AQP1) for CO₂ permeation in mammalian lungs and red blood cells, for example, remains unclear (22–29). Part of the controversy is apparently due to the fact that i), CO₂ permeation is usually measured indirectly, via acidification through the action of carbonic anhydrase (CA); ii), that, depending on the pH and on the concentration of CA, the CO₂ permeation is limited by unstirred water layers on both sides of the membrane (30); iii), that in erythrocytes this acidification is additionally influenced by the action of the -Cl⁻ transporter; and that iv), inhibitors like 4,4′-diisothiocyanato-stilbene-2,2′-disulfonic acid (DIDS) and mercurial compounds may not be specific for the -Cl⁻ transporter and AQP1, respectively (26–29). Moreover, there is a large variability in the CO₂ permeability of different membranes (26,29–31).

Another process for which aquaporin-mediated CO₂ permeation has been suggested to play a physiological role is photosynthesis. In a recent study it was shown that the leaf growth of tobacco plants was dependent on the level of NtAQP1 expression, an aquaporin homologous to human AQP1 (32).

Here, we study the barrier for CO₂ permeation through AQP1 as well as through a pure palmitoyloleoylphosphatidylethanolamine (POPE) bilayer, using extensive atomistic molecular dynamics simulations. We address the question whether CO₂ is likely to permeate through AQP1 and close homologs like the plant aquaporin NtAQP1 by comparing the free energy for CO₂ permeation through AQP1 to the corresponding profiles for a model lipid bilayer membrane. Moreover, we present profiles for CO₂ interactions along possible permeation pathways and estimate permeation coefficients.

METHODS

An equilibrated simulation box of a tetramer of human AQP1 embedded in a solvated lipid bilayer was chosen as starting configuration of the simulations. The starting structure of human AQP1 was modeled based on the x-ray structure of bovine AQP1 (PDB code 1J4N (10)) by mutating differing residues using the WHAT IF modeling software (33). The two sequences are 91% sequence identical, making a structural overlay straightforward. The periodic simulation box contained the AQP1 tetramer, 271 POPE lipids and 19769 TIP4P (34) water molecules. The OPLS all-atom force field (35,36) was employed for the protein, lipid parameters were taken from Berger et al. (37). Four chloride ions were added to neutralize the simulation cell. The system was equilibrated for 2 ns before the production simulations. A typical simulation box is shown in Fig. 1 a.

All simulations were carried out using the GROMACS simulation software (38,39). Electrostatic interactions were calculated with the particle-mesh Ewald method (40,41). Short-range repulsive and attractive dispersion interactions were described together by a Lennard-Jones potential, which was cut off at 1.0 nm. The Settle (42) algorithm was used to constrain bond lengths and angles of water molecules, and LINCS (43) was used to constrain all other bond lengths, allowing a time step of 2 fs. The simulation temperature was kept constant by weakly (τ = 0.1 ps) coupling the protein, lipids, and solvent separately to a temperature bath (44) of 300 K. Likewise, the pressure was kept constant by weakly coupling the system to a pressure bath of 1 bar. The xy (membrane plane) direction was pressure coupled with a coupling constant τ of 1 ps, whereas the box length in z-direction was kept fixed to avoid artifacts in the umbrella sampling simulations (see below).

The CO₂ simulations were set up as follows: the AQP1 channel was divided into 0.5-Å wide equidistant sections with the center of each section representing an umbrella center. For each section the equilibrated trajectory was searched for snapshots with water molecules as close as possible to the chosen position of the umbrella potential. This water molecule was subsequently replaced by CO₂. Umbrella sampling calculations were carried out by applying a harmonic restraint force along the pore coordinate with force constants between 100 and 6400 kJ/mol/nm². The umbrella sampling simulations were performed with CO₂ with the restraint acting on the carbon atom, and—for comparison and to validate the procedure—with water by restraining the oxygen atom. No restraint along the lateral directions was applied.

We observed that a possible flexibility of Arg-195 has strong impact on the barrier in the aromatic/arginine (ar/R) constriction region (see Results). Since these flexibilities are sampled slowly in the simulations we performed additional umbrella samplings in the ar/R region with starting positions taken from an independent equilibrium run. This way, eight independent channels contributed to our profile. Additionally, the protonation state of His-180 opposing the Arg-195 might influence the barrier height. In six out of eight channels simulated the δ-nitrogen was protonated, in two channels the ε²-nitrogen. In accordance with a recent finding for water permeation, no significant influence of the protonation state of His-180 on the CO₂ barrier was found (45).

To reduce CPU time, and hence improve sampling, CO₂ molecules were placed in each pore and, to further improve statistics, at three or four positions along the channel coordinate within one simulation run. A minimal distance of 25 Å between CO₂ molecules along the pore was imposed to ensure that interactions between CO₂ molecules (if any) are negligibly small. The umbrella sampling histograms from all monomeric channels were combined to compute a free energy profile using the weighted histogram analysis method (WHAM) (46). To account for the periodicity of the system, we implemented a cyclic version of the WHAM procedure. In total, 3912 histograms from 400 ps simulations were obtained (taking only the last 300 ps for analysis and using the first 100 ps for equilibration), extracted from a total of 334 ns of simulation time of the aquaporin system.

Additionally, CO₂ molecules were placed along the central channel surrounding the tetrameric axis. Umbrella positions were chosen as before, either by replacing a water molecule, or by adding it to the structure in case of an empty position. For the central pore along the tetrameric axis, 888 histograms were obtained.

Umbrella simulations to obtain a free energy profile for water along the AQP1 pore were performed the same way, except that water molecules close to the desired minimum of the artificial potential were not replaced but were directly restrained. Like for CO₂, the individual umbrella simulations were carried out for 400 ps. In total, 3092 histograms were obtained. Force constants were chosen between 100 and 1600 kJ/mol/nm².

Umbrella simulations for the lipid bilayer were performed by inserting CO₂ molecules at various positions into random snapshots of a 20-ns equilibrium run of a bilayer patch of 128 POPE and 4777 TIP4P water molecules. Up to 12 CO₂ molecules were inserted into one simulation keeping a minimum distance of 25 Å in the lateral direction and 30 Å perpendicular to the bilayer. A total of 3200 histograms were extracted from 260 simulations of 500 ps each, using the last 300 ps for analysis and the first 200 ps for equilibration. Force constants between 100 and 800 kJ/mol/nm² were used.

To account for fluctuations of each of the channels within the tetramer, the umbrella positions were corrected with respect to the center of the corresponding monomer. As a robust measure for the monomer position we chose the center of mass of the backbone atoms of the transmembrane helices.

Partial charges for CO₂ were obtained from electrostatic fitting using the CHELPG procedure to wavefunctions obtained at the RHF/6-31G* level, resulting in a charge of 0.9378e on the carbon atom and −0.4689e on the oxygen atoms. The Lennard-Jones parameters σ and ε for the CO₂ carbon were 0.375 nm and 0.439 kJ/mol, respectively, and for the CO₂ oxygen 0.296 nm and 0.877 kJ/mol, respectively. The CO₂ parameters were validated by comparison to the oil/water partition coefficient of CO₂. The experimentally determined value of near unity (47) compares favorably to the small free energy difference for CO₂ between the bulk water phase and the lipid tail region of the bilayer (cf. Fig. 2, shaded curve).

Profiles for enthalpic interactions of CO₂ (Fig. 3) were calculated as the sum of short-range Coulomb interactions (≤1 nm) and Lennard-Jones interactions. Note that these do not correspond to the total enthalpy for a given CO₂ position since they only include interactions directly involving CO₂. CO₂-induced alterations in the surroundings (like protein-water interactions) are not included in this component of the enthalpy. The profiles do, however, give more insight into the origin of the free energy barriers.

Permeability coefficients P_f for the POPE membrane were estimated by counting barrier crossings in simulations where CO₂ molecules were allowed to diffuse freely through a POPE bilayer (see below). After equilibration, we observed 3.4 barrier crossings per barrier, direction, and nanosecond. The permeability coefficient for one barrier of ΔG^POPE = 4 kJ/mol was calculated by where Φ = 3.4/ns denotes the transition rate, S = 30 nm² our membrane surface and c_s = 0.52 M the CO₂ concentration, giving . Since half of the CO₂ molecules crossing the first barrier will return and not permeate through the bilayer this refers to an estimated permeability coefficient for POPE of .

Using the permeability coefficient for a single POPE barrier we estimated the coefficients for the water pore and the central cavity via exp(− (ΔG′ − ΔG^POPE)/kT). Here, ΔG′ denotes the barrier of the water pore or the central cavity, T the temperature, and k the Boltzmann constant. In terms of rate theory, this estimate assumes the same prefactor for CO₂ permeation across a single POPE barrier as through a channel. Additionally, it rests on the assumption that the permeation rate is limited solely by the main barrier along the reactive coordinate. Along the central cavity two barriers of the same height are present, giving rise to an additional factor of 1/2. Unit permeability coefficients were calculated by p_f = P_fS_mono since the free energy barriers were normalized to an area of one AQP monomer. We estimated the area of one monomer as S_mono = 10.3 nm². Note that the area per monomer in vivo can be expected to be lower than the area per monomer in AQP1 crystals (11.5 nm²) (48).

Statistical errors of the free energy profiles were estimated using bootstrap analysis. To account properly for limited long-time sampling, we performed bootstrap analysis by considering complete histograms as independent data points: From the N original histograms we randomly selected N histograms (allowing duplication) which were subsequently used for WHAM analysis. This procedure was performed 50 times for each profile allowing the calculation of standard deviations. The standard deviations at the main barriers were 1kJ/mol for the monomeric channel, 2 kJ/mol for the central cavity, and 1 kJ/mol for the POPE membrane.

The barrier at the ar/R region in the monomeric channel is influenced strongly by the flexibility of Arg-195, which is sampled slowly (see above). Therefore, we calculated the barrier for each of the eight simulated channels separately. Considering these eight barriers as independent data points gave an additional uncertainty of 3 kJ/mol.

RESULTS AND DISCUSSION

Fig. 2 shows the free energy profile obtained as a potential of mean force from a set of umbrella sampling simulations for three potential pathways for CO₂ permeation through AQP1 embedded in a lipid bilayer: the pathway through an aquaporin water pore (black, solid line), the central channel (black, dashed line), and through the POPE membrane (shaded line). In the monomeric water channel CO₂ molecules encounter a major barrier of ∼23 ± 4 kJ/mol in the ar/R constriction region (also termed selectivity filter, see Fig. 1 c) and a lower one of ∼9 ± 1 kJ/mol around the region with the two conserved NPA motifs. The corresponding enthalpic interactions of CO₂ to water and to the protein, respectively, are shown in Fig. 3 b. As can be seen, there is no significant loss of interactions involving CO₂ in the ar/R region, rendering the free energy barrier either entropic in nature, or, more likely, due to indirect enthalpic effects like an unfavorable configuration of surrounding water molecules. This is illustrated in Fig. 5, which shows a typical pathway of CO₂ through the monomeric AQP1 channel, as derived from the umbrella sampling simulations. As can be seen, water molecules form strong hydrogen bonds to Arg-195 in the ar/R region, which are broken upon CO₂ passage. Fig. 3 a shows 115 pathways for CO₂ molecules that were not restrained by an umbrella potential. The trajectories show that the CO₂ molecules behave as expected from the potential of mean force, not spontaneously passing the ar/R constriction region within the simulated time of 500 ps.

The barrier in the ar/R is strongly influenced by the flexibility of Arg-195. In the crystal structure, its guanidinium group points straight into the pore creating the narrowest part of the channel. Two water molecules forming H-bonds to the ε- and η¹-nitrogens are displaced by a permeating CO₂ molecule (compare Fig. 5, second and third picture from bottom). In this configuration, the barrier is likely to be higher than 23 kJ/mol. However, during simulation we observed some flexibility of Arg-195 probably caused by flexibility in Loop C. Arg-195 is stabilized sterically by the neighboring Asn-127 and by a frequent H-bond to the carbonyl oxygen of Gly-125 (see Fig. 1). In such cases, no breakage of H-bonds between water and Arg-195 is necessary lowering the barrier for CO₂ passage significantly. The effect of the flexibility of Arg-195 on the CO₂ barrier is reflected by the reported uncertainty of 4 kJ/mol.

As described in the methods section we roughly estimated the unit permeability coefficient for CO₂ permeation through AQP1 to be p_f = 1 × 10⁻¹⁵ cm³/s. The error of 4 kJ/mol in the main barrier corresponds to factor of 5 in the permeabilities, i.e., a range of 2 × 10⁻¹⁶ to 5 × 10⁻¹⁵. The unit permeability coefficient for water permeation through a two-dimensional AQP1 crystal is (48), thus, the resistance of AQP1 to CO₂ is between one and two orders of magnitude higher in comparison to water. This result can also be expected by comparison of the corresponding free energy barriers of 12 and 23 kJ/mol, respectively.

The single channel permeability estimated from our results can be compared with the results from Nakhoul et al. (22), Prasad et al. (24), and Yang et al. (25) from experiments with oocytes, erythrocytes, and proteoliposomes containing AQP1. As discussed by Yang et al. (25), the results reported by Nakhoul et al. (22) refer to a single channel CO₂ permeability of 1−2 × 10⁻¹⁴ cm³/s, whereas the results reported by Prasad et al. (24) give ∼5 × 10⁻¹⁵ cm³/s. Yang et al. (25) did not observe AQP1-mediated CO₂ permeation and reported an upper limit of 3 × 10⁻¹⁶ cm³/s. Therefore, our estimate of ∼1 × 10⁻¹⁵ cm³/s most closely matches the result by Prasad et al. and the upper limit by Yang et al.

A second possible pathway through AQP1 would be the central pore along the fourfold axis of the tetramer. Due to its hydrophobicity the central pore is neither filled by water nor does it conducts water (see Fig. 1 c). Two major barriers for CO₂ permeation can be identified: first, a barrier of ∼12 ± 2 kJ/mol at z = 7.5 Å is located near the extracellular entrance to the central cavity, surrounded by the 4 Val-50 of the 4 monomers. In this region CO₂ molecules lose favorable interactions to neighboring water molecules (compare Fig. 3 c, green and red lines). A second barrier of the same height is found ∼7 Å above the first one. Here, the CO₂ is surrounded by the 4 Asp-48, and hence, it is likely that CO₂ competes with water for hydrogen bonds to Asp-48. The barrier translates into a permeability coefficient of P_f = 0.1 cm/s for the central cavity, referring to the area of an AQP1 tetramer, or a unit permeability coefficient of p_f = 4 × 10⁻¹⁴ cm³/s.

The free energy barrier for CO₂ permeation through the central cavity is significantly smaller than for the monomeric channel. When taking into account that each tetramer contributes four monomeric channels and one central pore, these barriers would imply that ∼10% of AQP1-permeated CO₂ would be expected to permeate via the monomeric channels. Note, however, that in these calculations we have assumed an empty central channel. Should an ion or organic molecule be bound in the central cavity under physiological conditions (like in the recently resolved structure of an aquaporin-Z tetramer (49)), this would be expected to drastically increase the barrier for CO₂ permeation through the central pore. Experiments with mercurial AQP1 blockers and with the C189S mutant suggest that AQP1-mediated CO₂ permeation takes place through the monomeric pore (23,24), which would imply that, indeed, the central channel is blocked under the experimental conditions.

The tobacco aquaporin NtAQP1 was recently reported to facilitate physiologically relevant CO₂ permeation (32). The pore region of NtAQP1 is highly similar to human AQP1, with the residues surrounding the major barriers in the NPA and ar/R region fully conserved among the two proteins (32). Our results suggest that—if the central channel is blocked—aquaporin-mediated CO₂ permeation is expected to play a significant role in photosynthesis only if the plant cell membranes have a low intrinsic CO₂ permeability, with an activation barrier well above 20 kJ/mol. However, since the structure of the central pore of NtAQP1 is unknown, we cannot exclude that permeation via the central pore lowers the effective barrier.

CONCLUSIONS

Taken together, these results suggest that CO₂ permeation through AQP1 can be expected to play a physiological role only in membranes with a low intrinsic CO₂ permeability, such as membranes of cells exposed to an inhospitable environment (26,31). Membranes with similar physicochemical characteristics to POPE are highly permeable to CO₂, rendering a physiological role for AQP1-mediated CO₂ permeation in such membranes unlikely. The role of AQP1 mediated CO₂ permeation in red blood cells, mammalian lungs and other tissues therefore depends on the intrinsic CO₂ permeability of the surrounding membrane. Only for membranes with an activation barrier well above 20 kJ/mol or an empty central channel can AQP1-mediated CO₂ permeation be expected to play a major role, if a high enough expression level is provided. An interesting test of this hypothesis would be an evaluation of the CO₂ permeability of vesicles with different lipid composition and different levels of embedded AQP1.