|Home | About | Journals | Submit | Contact Us | Français|
In order to elucidate how phosphatidylserine (PS6) interacts with AQP5 in a cell membrane, we develop a hybrid steered molecular dynamics (hSMD) method that involves (1) simultaneously steering two centers of mass of two selected segments of the ligand and (2) equilibrating the ligand-protein complex with and without biasing the system. Validating hSMD, we first study vascular endothelial growth factor receptor 1 (VEGFR1) in complex with N-(4-Chlorophenyl)-2-((pyridin-4-ylmethyl)amino)benzamide (8ST), for which the binding energy is known from in vitro experiments. In this study, our computed binding energy well agrees with the experimental value. Knowing the accuracy of this hSMD method, we apply it to the AQP5-lipid-bilayer system to answer an outstanding question relevant to AQP5’s physiological function: Will the PS6, a lipid having a single long hydrocarbon tail that was found in the central pore of the AQP5 tetramer crystal, actually bind to and inhibit AQP5’s central pore under near-physiological conditions, namely, when AQP5 tetramer is embedded in a lipid bilayer? We find, in silico, using the CHARMM 36 force field, that binding PS6 to AQP5 is a factor of 3 million weaker than “binding” it in the lipid bilayer. This suggests that AQP5’s central pore will not be inhibited by PS6 or a similar lipid in a physiological environment.
Accurately computing the free-energy of binding a ligand to a protein is essential in the field of biophysics/biochemistry/quantitative biology because protein-ligand interactions are fundamentally involved in most biological processes such as substrates binding to enzymes, antigens binding to antibodies, antagonists binding to receptors etc. (Olsson et al., 2008, Drake and Klakamp, 2007, General et al., 2012, Leavitt and Freire, 2001) One particular binding problem of our concern is the interaction between the membrane protein AQP5(Krane et al., 2001a, Krane et al., 2001b, Verkman, 2002, Burghardt et al., 2006, Delporte and Steinfeld, 2006, Horsefield et al., 2008, Kang et al., 2008, Musa-Aziz et al., 2009, Watanabe et al., 2009, Törnroth-Horsefield et al., 2010, Zhang et al., 2010, Larsen et al., 2011, Zhang et al., 2011, Geyer et al., 2013b, Qin and Boron, 2013) and phosphatidylserine (PS6), a lipid having a single long hydrocarbon tail that was found (residing in the central pore) in the crystallographic structure of the AQP5 tetramer.(Horsefield et al., 2008)
It is clear from in vitro experiments and in silico simulations that AQP5 facilitates water transport across the cell membrane by lining water molecules in single file throughout its water-conducting pore formed within each monomer protein.(Krane et al., 2001a, Krane et al., 2001b, Verkman, 2002, Burghardt et al., 2006, Delporte and Steinfeld, 2006, Horsefield et al., 2008, Törnroth-Horsefield et al., 2010, Gravelle et al., 2013, Janosi and Ceccarelli, 2013, Zhang and Chen, 2013, Eckhard et al., 2014) It is also clear from in vitro experiments that AQP5 conducts transport of gas molecules as well,(Musa-Aziz et al., 2009, Geyer et al., 2013b, Qin and Boron, 2013) while carbon dioxide transport across the cell membrane is subject to much active debate.(Hub and de Groot, 2006, Wang et al., 2007, Missner et al., 2008a, Missner et al., 2008b, Hub et al., 2010, Wang et al., 2010, Wang and Tajkhorshid, 2010, de Groot and Hub, 2011, Itel et al., 2012, Kaldenhoff, 2012, Geyer et al., 2013a, Geyer et al., 2013c, Endeward et al., 2014, Hulikova and Swietach, 2014, Kaldenhoff et al., 2014) However, how AQP5 facilitates the transport of apolar, hydrophobic gas molecules is currently subject to much discussion.(Wang et al., 2007, Kaldenhoff et al., 2014) The X-ray structure of the AQP5 crystal tells us that AQP5 is tetrameric in conformation and that the four monomers in quasi-four-fold symmetry leave a significant void in the middle along the symmetry axis. This void (the central pore) is lined with hydrophobic residues in quadruplets (one from each of the tetrameric monomers) which interact optimally with the long hydrocarbon tail of PS6. If not plugged up by PS6 or another ligand, the central pore would be an ideal channel for conducting gas molecules that cannot easily pass through the four amphipathic water pores. Therefore, the relevant question is: Does the PS6 (bound in the AQP5 crystal) actually bind to the AQP5 central pore (AQP5c.p.) in a physiological setting, namely, when the protein is embedded in the cell membrane?
The in silico study of (Zhang and Chen, 2013) gave a value of the dissociation constant kDc.p.→aq.=6.1μM for binding PS6 in the AQP5c.p., indicating strong PS6-AQP5 interaction. This question remains: Is this binding strong enough in comparison with “binding” PS6 in a lipid bilayer (LBL)? In this paper, we develop a new method capable of resolving both binding PS6 in AQP5c.p. and binding PS6 in LBL. We answer the afore-stated question by computing the binding energies of the PS6-in-AQP5c.p. complex and of the PS6-in-LBL system. We conduct steered molecular dynamics (SMD)(Isralewitz et al., 1997, Gullingsrud et al., 1999, Jensen et al., 2002, Tajkhorshid et al., 2003, Li and Makarov, 2004, Park and Schulten, 2004, Jensen et al., 2007, Wang and Zhang, 2007, Minh and McCammon, 2008, Chen, 2011, Trinh et al., 2011, Fukunishi et al., 2012, Moradi and Tajkhorshid, 2013, Nicolini et al., 2013, Trinh et al., 2013, Velez-Vega and Gilson, 2013, Yu et al., 2013) simulations to steer/pull the long-tailed PS6 from its bound state in the central pore of AQP5 to the dissociated state in the bulk region and to pull PS6 from its equilibrium residence in the lipid bilayer to the dissociated state in the bulk. The simulation is not the widely used, regular SMD of pulling one center of mass of one selection of the ligand atoms but a hybrid SMD (hSMD) of pulling two centers of mass of two selected segments of the ligand. The hSMD also involves necessarily equilibration processes with and without partially biasing the pulling centers. To validate the hSMD approach, we apply it to computing the absolute binding energy of the vascular endothelial growth factor receptor 1 (VEGFR1) in complex with N-(4-Chlorophenyl)-2-((pyridin-4-ylmethyl)amino)benzamide (8ST).(Furet et al., 2003, Tresaugues, 2013) The computed binding energy of the VEGFR1-8ST system well agrees with the in vitro experimental data.(Furet et al., 2003) Knowing its accuracy, we use the hSMD method to study the PS6-in-AQP5-c.p. complex and the PS6-in-LBL system. The computed binding affinity of PS6-in-lipid-bilayer is about 3 million times the binding affinity of PS6-in-AQP5-central-pore, suggesting that PS6 does not inhibit AQP5 in a physiological environment. AQP5 can facilitate transport of apolar gas through its central pore practically unhindered even with PS6 or another single-long-tailed lipid present in the aqueous environments sandwiching the cell membrane.
We use all-atom model systems with explicit solvent (TIP3P water). All the atomistic interactions (bonded or non-bonded) are represented by the CHARMM 36 force field.(Brooks et al., 2009, Vanommeslaeghe et al., 2010) The equilibrium molecular dynamics (MD) runs and the nonequilibrium SMD runs (Langevin dynamics in both cases) were implemented with NAMD. The temperature was kept at 293.15K and the pressure at 1 bar (using Langevin piston). The Particle mesh Ewald was used for computing the Coulomb interactions in an exact manner with periodic conditions applied in all three dimensions. The time step of 1.0 fs was used for short-range interactions while 4.0 fs for long-range interactions. The Langevin damping was chosen as 5.0/ps.
With theoretical details given in the first part of the supplemental information (SI), the gist of our hSMD is the following: To efficiently steer a long ligand out of its binding site, we choose two (mutually exclusive) segments of the molecule. The two centers of mass of the two chosen segments are used as pulling centers. We divide the dissociation path from the binding site to the bulk region into multiple sections. In each section, we pull the pulling centers forward and backward to sample a number of forward and reverse pulling paths. Work done to the system is computed along every pulling path. The work curves are used to extract the free energy or reversible work or Potential of mean force (PMF)(Kirkwood, 1935, Chandler, 1978, Pratt et al., 1994, Roux, 1995, Allen et al., 2006) as a function of the ligand displacement from its binding site. The PMF difference between the bound state and the dissociated state gives a large part of the absolute binding energy. Computing the ligand fluctuations at the binding site gives another part of the binding energy. Computing the rotation and fluctuations of the ligand in the bulk region gives the final part of the binding energy. The three parts put together give the absolute binding energy that should be compared the experimental data.
We applied this hSMD method to the VEGFR1-8ST complex whose binding affinity was experimentally measured. The agreement between the computed result (details shown below in this subsection) and the experimental value validates our method. The all-atom model system of the VEGFR1-8ST complex is illustrated in Fig. 1, which was built as follows: We took the atomistic coordinates of the X-ray crystal structure of the vascular endothelial growth factor receptor 1 in complex with N-(4-Chlorophenyl)-2-((pyridin-4-ylmethyl)amino)benzamide (PDB: 3HNG(Tresaugues, 2013)) and set up the all-atom interactions using the CHARMM force field.(Brooks et al., 2009, Vanommeslaeghe et al., 2010) We placed this complex in the center of a TIP3P(Jorgensen et al., 1983) water box, neutralized the system with eight Cl− ions, and ionized it with 85 Na+ and 85 Cl− ions to a concentration of 150 mM. A long (30 ns) equilibrium MD run was conducted to fully equilibrate the system. During the MD and SMD runs, we fixed some alpha carbons on the helices that are farther than 10Å away from the ligand to fully respect the crystal structure of the protein. We chose two segments of the ligand as shown in Fig. 1(d) for steering/pulling.
The PMF along the pulling path from the bound state to the dissociated state is shown in Fig. 2(a), which was computed from the work curves (shown in supplementary Fig. S1) using SI Eq.(2). Note that the two steering centers (the center of mass of Seg. 1 and the same of Seg. 2) are subject to identical displacement along the z-axis Δz = vdΔt with a pulling speed of vd=2.5Å/ns. The PMF in Fig. 2(a)
is a function of Δz with Δz= 0 at the bound state and Δt representing a time interval. It needs to be pointed out that the path, along which two centers of mass of two selections of atoms are pulled, should not be mistaken as the reaction path from the dissociated state to the bound state. The high barrier value of about 40 kcal/mol does not imply the prediction of a vanishingly small reaction (association) rate but indicates the actual association path is very different from the straight lines of our SMD pulling. However, free energy is a function of state. The free-energy difference between any two states (in particular, the bound and the dissociated states here) is independent of the paths connecting the two states. The PMF curve in Fig. 2(a) is only and correctly used for computing the binding energy, namely, the free-energy difference between the bound and the dissociated states.
In the dissociated state, Δz =25Å, a second set of SMD runs were conducted: Fixing the center of mass of Seg. 2 and steering the center of mass of Seg. 1 from and to Seg. 2. The work done to the system along the pulling paths are shown supplemental Fig. S2, from which we used SI Eq. (2) to extract the PMF, W(r), as a function of the distance
between the Seg. 1 and Seg. 2 centers. The PMF is shown in Fig. 2(b). Note that the same reference point is used for both PMFs:
And W(r0) = 0 with
Using these two PMFs and conducting equilibrium MD for the bound state deviations (shown in Supplemental Information) involved in SI Eq.(3), we obtained the dissociation constant of kD ≈ 215nM in comparison with the in vitro experimental data of 180 nM.(Furet et al., 2003, Tresaugues, 2013) In terms of absolute binding energy, we obtained the free-energy of binding in the amount of −9.2±0.6 kcal/mol that well agrees with the in vitro value of −9.3 kcal/mol. All this indicates the validity of our hSMD approach of pulling two centers of mass of two segments of a ligand.
The all-atom model system of the AQP6-PS6 complex embedded in a Phosphatidylethanolamine (POPE) bilayer is illustrated in Fig. 3. Briefly, we took the crystal coordinates of AQP5 (PDB: 3D9S) to build the AQP5-tetramer in complex with PS6 (residing in the central pore). We embedded the AQP5-PS6 complex in the center of a POPE-bilayer. We also embedded a PS6 in the bilayer but away from the protein. The AQP5-PS6-POPE complex was solvated with TIP3P waters, neutralized with 22 Cl− ions, and ionized with Na+ and Cl− ions to the concentration of 150 mM. The system was equilibrated for 300 ns while the alpha carbons on the transmembrane helices were fixed to their crystal coordinates in full respect for the crystallographic structure. The illustrations in Fig. 3 are from the fully equilibrated system.
The binding affinities were computed by conducting two sets of SMD runs pulling PS6 out of the AQP5 central pore (Fig. 3(c)) and pulling PS6 out of the POPE-bilayer (Fig. 3(d)). Following the schemes described in the Method section and using two steering centers (the centers of mass of Seg. 1 and Seg. 2 shown in Fig. 3(b)), we obtained the PMFs along the dissociation path from the AQP5 central pore and along the dissociation path from inside the POPE-bilayer. The SMD run details are identical to that for 8ST described in the previous subsection. The work curves along the pulling paths are shown in the supplemental figures Figs. S4 and S5. The PMFs shown in Fig. 4 clearly tell us that the binding affinity of PS6 in the central pore of AQP5 is much weaker (by a factor of 1 over 3 million) than the “binding” affinity of PS6 residing inside the POPE-bilayer. Therefore, we suggest that, under physiological conditions when AQP5 resides in the cell membrane, PS6, if present, will reside alongside the membrane lipids instead of binding to the AQP5 central pore. The central pore of the AQP5 will remain uninhibited by PS6 in a physiological environment.
In this case, it is unnecessary to compute the absolute binding energies because the equilibrium deviations of PS6 residing inside the POPE bilayer are far greater than its deviations in the AQP5 central pore. In fact, PS6 is free to diffuse in the 2 dimensions of the POPE-bilayer just like other lipids in the bilayer. These effects further reduce the PS6-AQP5-central-pore affinity relative to the PS6-POPE-bilayer affinity and strengthen our conclusion that the AQP5 central pore will not be inhibited by PS6 in presence of the cell membrane.
It needs to be noted that Fig. 4 shows an apparently very high barrier (~19 kcal/mol) on the pathway for PS6 to dissociate from AQP5c.p. into the bulk region and then move to take residence in the lipid bilayer. The PMF, however, is only one of the three factors contributing to the Gibbs free energy along the ligand dissociation path. When the other two factors (fluctuations at the binding site and rotation and fluctuations in the bulk) are fully accounted for, the barrier will be approximately 7.2 kcal/mol.(Zhang and Chen, 2013) Therefore, the time it takes PS6 to dissociate from AQP5c.p. would be on the order of submicroseconds.
In terms of computational approach, we have developed a hybrid steered molecular dynamics approach for computing absolute binding energy from the PMF along a dissociation path implemented by simultaneously steering/pulling two center(s) of mass of two segments of a ligand. The segments chosen to be steered need to be the most stable parts of the ligand in the bound state because the deviations of these centers were approximated as Gaussian. This approximation is expected to be valid for a wide range of ligand-protein complexes because there should always be some part of a ligand that does not deviate greatly from its bound-state coordinates but is tightly bound to the protein except the cases of very weak binding.
In terms of biological functions, we have ascertained the following conclusion based on the well accepted CHARMM force field: PS6 will practically not bind to the central pore of AQP5 in a physiological environment but practically certainly integrate itself into the cell membrane. It is worth noting that PS6, unlike other lipids, has only one long hydrophobic tail (hydrocarbon chain) interacting with the hydrophobic side-chains that line the central pore of AQP5. These non-specific interactions are responsible for binding PS6 to AQP5 in absence of a cell membrane. It is thus speculated that other single-tailed lipids would bind there too when AQP5 is not embedded in the cell membrane. For systems with AQP5 embedded in a lipid bilayer, PS6 or another lipid would reside in the bilayer, leaving the AQP5 central pore open for whatever functions it may have. This conclusion is in line the in vitro findings that AQP5 facilitates transport of CO2 across the membrane.
The author gratefully acknowledges support from the NIH (GM 084834) and the Texas Advanced Computing Center.
DECLARATION OF INTEREST
The author reports no declarations of interest.