Article pubs.acs.org/JPCA
Toward Molecular Mechanism of Xenon Anesthesia: A Link to Studies of Xenon Complexes with Small Aromatic Molecules Natalya N. Andrijchenko,† Alexander Yu. Ermilov,† Leonid Khriachtchev,‡ Markku Ras̈ an̈ en,‡ and Alexander V. Nemukhin*,†,§ †
Department of Chemistry, M. V. Lomonosov Moscow State University, Leninskie Gory 1/3, Moscow 119991, Russian Federation Department of Chemistry, University of Helsinki, P.O. Box 55, Helsinki FI-00014, Finland § N. M. Emanuel Institute of Biochemical Physics, Russian Academy of Sciences, Kosygina 4, Moscow 119334, Russian Federation ‡
ABSTRACT: The present study illustrates the steps toward understanding molecular mechanism of xenon anesthesia by focusing on a link to the structures and spectra of intermolecular complexes of xenon with small aromatic molecules. A primary cause of xenon anesthesia is attributed to inhibition of N-methyl-D-aspartate (NMDA) receptors by an unknown mechanism. Following the results of quantum mechanics/molecular mechanics (QM/MM) and molecular dynamics (MD) calculations we report plausible xenon action sites in the ligand binding domain of the NMDA receptor, which are due to interaction of xenon atoms with aromatic amino-acid residues. We rely in these calculations on computational protocols adjusted in combined experimental and theoretical studies of intermolecular complexes of xenon with phenol. Successful reproduction of vibrational shifts in molecular species upon complexation with xenon measured in low-temperature matrices allowed us to select a proper functional form in density functional theory (DFT) approach for use in QM subsystems, as well as to calibrate force field parameters for MD simulations. The results of molecular modeling show that xenon atoms can compete with agonists for a place in the corresponding protein cavity, thus indicating their active role in anesthetic action.
■
INTRODUCTION The role of relatively weak noncovalent interactions between atoms and molecules in many physical, chemical, and biological phenomena cannot be underestimated. Numerous experimental and theoretical studies of intermolecular complexes are related to distinctively different directions, among which we emphasize those devoted to interactions of rare gas xenon with small organic molecules. The results of these investigations find important biomedical applications such as using xenon as a marker in protein cavities1 or understanding mechanism of xenon anesthesia.2 A link between the molecular mechanism of xenon anesthesia and subtle details of intermolecular complexes of xenon may be outlined as follows. A primary cause of xenon action is attributed to inhibition of N-methyl-D-aspartate (NMDA) receptors, multidomain membrane proteins, three-dimensional atomic structures of which have become available only recently.3 If xenon atoms are capable of occupying places in a ligand binding cavity and competing with the natural agonist, glycine, then the mechanism of action is getting clearer. In turn, the walls of such protein cavities contain amino acid residues with aromatic rings, tyrosine (Tyr), phenylalanine (Phe), and tryptophan (Trp), an interaction with which is a fingerprint of the resolved xenon−protein structures deposited in the Protein Data Bank. A typical motif of these complexes is a pattern with a xenon atom residing aside the plane of an aromatic ring; the distances between xenon and carbon atoms vary between 3.5 and 5 Å. For brevity, we call these complexes the π-type © XXXX American Chemical Society
structures; they represent a part of noncovalently bound systems with aromatic rings.4 From the theoretical side, straightforward simulations of xenon atoms surrounded by organic molecules using molecular dynamics (MD) calculations5−7 face a well-known problem of reliability of force field parameters8,9 required to estimate the energy of intermolecular complexes involving xenon. Quantum-based models seem to be more promising in this respect; however, the necessity to include large fractions of protein cavity walls to the model systems encourages one to use the quantum mechanics/ molecular mechanics (QM/MM) approaches.10 Apparently, xenon and the nearest molecular groups must be assigned to the quantum subsystem, and a proper accuracy level11,12 at a reasonable computational cost must be provided to describe the corresponding complexes. Experimentally, complexes of rare gas atoms with small organic molecules can be studied in the gas phase13 and in the low-temperature matrices.14 In the latter case, subtle changes in vibrational bands of a molecule due to its complexation with xenon atoms can be accurately measured. In this computational work, we use the data obtained in matrix experiments with xenon and aromatic molecules to select a proper functional form in the density functional theory Special Issue: Markku Räsänen Festschrift Received: August 31, 2014 Revised: September 30, 2014
A
dx.doi.org/10.1021/jp508800k | J. Phys. Chem. A XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry A
Article
(DFT) approach for quantum-based characterization of π-type xenon complexes, as well as to calibrate force field parameters for MD simulations. The QM/MM and MD modeling tools allow us to determine the trapping sites of xenon in the ligand binding domain of N-methyl- D -aspartic acid (NMDA) receptors.
in spite of the presence of water molecules competing for the spots near amino acid residues forming the walls of the cavity. The analysis of conformations along MD trajectories shows at least three trapping sites occupied by xenon inside the cavity. To characterize these sites, we performed QM/MM calculations with the DFT(BHHLYP)-D2/6-31G* approach to estimate the energies and forces in the QM subsystems. The NWChem program suite was employed in these calculations.23 Specification of the QM parts for each trapping site is given in the Results. The starting coordinates for the QM/MM minimization were generated from appropriate frames along MD trajectories. This strategy did not guarantee that all sound equilibrium geometry configurations could be located; however, our goal was to find at least a few stable structures with a xenon atom trapped inside the cavity when the glycine molecule was replaced.
■
MODELS AND METHODS To model interactions of the rare-gas element xenon with small biologically related organic molecules, it is recommended to apply approaches capable to provide a well balanced description of noncovalent interactions, such as the second-order Møller− Plesset perturbation theory, MP2.11,15 However, attempts to use MP2 in QM/MM simulations of systems modeling xenon in protein clefts would lead to prohibitively expensive calculations. Therefore, we applied a two-step strategy to select a modeling tool. First, we showed that the MP2/aug-cc-pVTZ approach allowed us to reproduce vibrational spectra of the phenol···Xen (n = 2−4) complexes examined by infrared absorption spectroscopy in neon matrices.14 Second, we used the results of equilibrium geometry configurations and binding energies of the complexes of xenon with small aromatic molecules obtained in the MP2 calculations as reference data to select an inexpensive DFT-based method. It was demonstrated16 that DFT with a dispersion correction of type D217 with the functional BHHLYP and a modest basis set 6-31G* was able to approach the MP2 data. In particular, the DFT-D2 computed binding energy for the Xe−phenol π-type complex was 2.7 kcal mol−1, which is similar to 3.1 kcal mol−1 attainable in the MP2/cc-pVTZ calculations. The deviations of the distances between Xe and carbon atoms obtained in the DFTD2 calculations from the reference values did not exceed 0.05 Å. The DFT(BHHLYP)-D2 method also reproduced the differences between the two structures (π-type and H-bonded) of the Xe−phenol complex, as well as the trend toward an increase in binding energy of Xe in the series of aromatic amino acids Phe, Tyr, Trp.16 This theoretical approach was used for to calculate energies and forces in the QM subsystems in the QM/ MM calculations. To construct molecular models of a NMDA ligand binding domain, we used the coordinates of heavy atoms of the protein structure PDBid: 2A5T18 corresponding to the crystal of the ligand-binding core complex containing glycine. After hydrogen atoms and sodium and chloride ions with a concentration of up to 0.15 M were added, via the VMD software,19 energy minimization and a few MD runs were carried out. Water molecules recognized in the crystal structure were included in the model system, and the TIP3P water molecules20 formed a solvation shell for the protein. The MD simulations were performed with the NAMD program21 using the CHARMM force field parameters.22 During MD simulations of about 20 ns length, the systems were kept at a constant temperature of 298 K and pressure of 1 atm (NPT ensemble); the periodical boundary conditions were imposed. At the next step, the glycine moiety was replaced by a xenon atom, and energy minimization and MD runs (up to 50 ns) were repeated. In these calculations the CHARMMM force field parameters were corrected by improving the Lennard-Jones interaction potential (ε = 1.7 kcal mol−1) for the pairs containing xenon in such a manner that the equilibrium geometry configurations of the Xe−phenol complex obtained in the MP2/aug-cc-pVTZ calculations could be reproduced in MM simulations. Importantly, the xenon atom did not escape from the cavity
■
RESULTS MD Simulations. Ten classical molecular dynamics trajectories extended up to 50 ns were executed from the initial coordinates generated by replacement of glycine by a xenon atom in the model system constructed by motifs of the crystal structure PDBid: 2A5T.18 Eight of ten MD runs showed that Xe migrated between three local areas inside the protein cavity called below trapping sites I, II, and III. In two cases, a firm assignment of a xenon atom to a particular spot could not be accomplished during the simulation period. Figure 1 shows
Figure 1. Trapping sites of a xenon atom migrating inside the cavity of the ligand-binding domain along a representative MD trajectory.
the location of trapping sites I, II, III in the cleft; Figure 2 illustrates migration of a xenon atom, initially placed near the Phe92 residue (similarly to the glycine molecule in the crystal structure PDBid: 2A5T), along a representative trajectory. An important conclusion that can be drawn from these MD calculations is that xenon does not escape from the cavity despite the latter being filled by water molecules (not shown in Figure 1), which also occurs for the binding sites with amino acid residues forming walls of the cleft. In qualitative agreement with the results of previous classical MD simulations5,6 our data evidence that xenon can compete with glycine in the ligand binding domain of the NMDA receptor. To characterize possible trapping sites of xenon atoms inside the cleft, we applied a higher level theoretical approach, namely, the QM/ MM theory. B
dx.doi.org/10.1021/jp508800k | J. Phys. Chem. A XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry A
Article
arrangement is typical for complexes of xenon with phenol studied in our previous papers.14,16 According to the MD simulations (Figure 2), a xenon atom can migrate to another site leaving the initial position near Phe92. In the calculations for trapping site II, the QM subsystem included Xe, side chains of Gln13, Phe92, Pro124, Ser180, Trp223, Asp224, Val227, Phe250, and one water molecule. The equilibrium geometry configuration corresponding to trapping site II is illustrated in Figure 4. In this case, xenon atom is still bound to Phe92 but tends to find an additional binding partner, this time, another residue of aromatic nature, Trp223.
Figure 2. Root mean square deviations (RMSD) values along a MD trajectory computed from the initial position of Xe, which corresponds to replacement of glycine in the crystal structure PBDid: 2A5T by a xenon atom.
QM/MM Simulations. According to the MD results, a xenon atom migrates from one local environment to another; therefore, it is impossible to select molecular groups constituting a quantum subsystem for QM/MM calculations in the same manner for three possible trapping sites. In every selection of the QM subsystem, we followed a goal to completely surround a xenon atom by environmental molecular groups keeping an amount of quantum atoms to about 100. To generate initial coordinates for the QM/MM minimization, we used appropriate frames of MD trajectories assigned to areas I, II, or III. For trapping site I, the QM subsystem included Xe, side chains of Phe92, Leu125, Thr126, Ser180, Trp223, Phe250, and one water molecule. The equilibrium geometry configuration is illustrated in Figure 3. Apparently, the xenon atom forms a πtype complex with the aromatic ring of Phe92. This
Figure 4. Structure of trapping site II.
For trapping site III shown in Figure 5, the QM subsystem included Xe, side chains of Gln13, Phe16, Phe92, Pro124, Gln144, Asp224, Ala226, Val227, Phe250, and three water molecule. In this site, the side chain of Phe16 is the primary partner to form the π-type complex with a xenon atom.
Figure 3. Structure of trapping site I. Here and in all figures, the atoms of the QM subsystem are shown in balls and sticks; carbon atoms are green, oxygen red, nitrogen blue; and the distances are given in Å. C
dx.doi.org/10.1021/jp508800k | J. Phys. Chem. A XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry A
Article
the NMDA protein cavity can be identified, in which xenon is located. The motif of the π-type complex between the xenon atom and the aromatic molecule is typical for these trapping sites. Correspondingly, xenon atoms can compete with agonist species (glycine) for a place in the ligand binding domain of a protein receptor, thus indicating its role in anesthetic action. The computational method was calibrated by using the experimental data on the xenon−phenol complexes in lowtemperature matrices, which demonstrates the benefits of matrix-isolation spectroscopy for biology-oriented research.
■
AUTHOR INFORMATION
Corresponding Author
*A. V. Nemukhin. Phone: +7-495-939-10-96. E-mail: [email protected]; [email protected]. Notes
Figure 5. Structure of trapping site III.
The authors declare no competing financial interest.
■
■
CONCLUDING REMARKS The results of MD and QM/MM simulations show, first, that a xenon atom can be trapped in the ligand binding domain replacing glycine and, second, that the motif of the π-type complex is typical for these trapping sites. Analysis of biochemical aspects of xenon anesthesia is far beyond the scope of the present work. We rely on the knowledge that, most likely, xenon serves as an antagonist of the NMDA receptor inhibiting the receptor at its glycine binding site. Analysis of the available crystal structure of the ligand binding domain of the NMDA receptor shows that xenon can be bound to the side chains of aromatic amino acid residues, consistent with the molecular structures of xenon− protein complexes deposited to the protein data bank. Therefore, detailed knowledge on the structure and properties of intermolecular complexes of xenon atoms with the corresponding parent aromatic molecules like phenol, cresol, and toluene is essential when the interaction of xenon with protein receptors is modeled. We have shown that the subtle changes in vibrational frequencies of xenon complexes measured in low-temperature matrices are well reproduced by the MP2/aug-cc-pVTZ approach. In turn, the equilibrium geometry configurations and binding energies of this theoretical level are well reproduced using the much less expensive DFT-D2 approach17 with the functional BHHLYP and the basis set 6-31G*. Therefore, we can rely on the QM/MM calculations describing xenon binding sites in the NMDA receptor using the DFT(BHHLYP)-D2/6-31G* approximation in relatively large quantum subsystems. The results of these QM/MM calculations demonstrate that xenon can reside in several trapping sites showing patterns of π-type complexes. The force field parameters for the MD simulations were also calibrated by the structures of the π-type complexes of xenon with phenol. Specifically, the parameter ε of the Lennard-Jones interaction potential was adjusted to reproduce the MP2/augcc-pVTZ data. The large-scale MD calculations of the ligandbinding domain of the NMDA receptor will be reported elsewhere because their description requires a specific presentation typical for biochemistry oriented journals, e.g., ref 24. In this work, we used relatively short MD trajectories prior the expensive QM/MM minimization procedures to account for suitable arrangement of protein conformations. To conclude, the results of the present MD and QM/MM calculations show that several plausible xenon binding sites in
ACKNOWLEDGMENTS We thank Bella Grigorenko for discussion and advice. This study was supported by the Russian Science Foundation (project 14-13-00124). We acknowledge the use of supercomputer resources of the M. V. Lomonosov Moscow State University25 and of the Joint Supercomputer Center of the Russian Academy of Sciences.
■
REFERENCES
(1) Oros, A. M.; Shah, N. J. Hyperpolarized Xenon in NMR and MRI. Phys. Med. Biol. 2004, 49, R105−R153. (2) Derwall, M.; Coburn, M.; Rex, S.; Hein, M.; Rossaint, R.; Fries, M. Xenon: Recent Developments and Future Perspectives. Minerva Anestesiol. 2009, 75, 37−45. (3) Sobolevsky, A. I.; Rosconi, M. P.; Eric Gouaux, E. X-Ray Structure, Symmetry and Mechanism of an AMPA-Subtype Glutamate Receptor. Nature 2009, 462, 745−758. (4) Li, S.; Xu, Y.; Shen, Q.; Liu, X.; Lu, J.; Chen, Y.; Lu, T.; Luo, C.; Luo, X.; Zheng, M.; et al. Non-Covalent Interactions with Aromatic Rings: Current Understanding and Implications for Rational Drug Design. Curr. Pharm. Des. 2013, 19, 6522−6533. (5) Dickinson, R.; Peterson, B. K.; Banks, P.; Simillis, C.; Martin, J. C.; Valenzuela, C. A.; Maze, M.; Franks, N. P. Competitive Inhibition at the Glycine Site of the N-Methyl-D-Aspartate Receptor by the Anesthetics Xenon and Isoflurane: Evidence from Molecular Modeling and Electrophysiology. Anesthesiology 2007, 107, 756−767. (6) Liu, L. T.; Xu, Y.; Tang, P. Mechanistic Insights into Xenon Inhibition of NMDA Receptors from MD Simulations. J. Phys. Chem. B 2010, 114, 9010−9016. (7) Weinrich, M.; Worcester, D. L. Xenon and other Volatile Anesthetics Change Domain Structure in Model Lipid Raft Membranes. J. Phys. Chem. B 2013, 117, 16141−16147. (8) Paton, R. S.; Goodman, J. M. Hydrogen Bonding and PiStacking: How Reliable are Force Fields? A Critical Evaluation of Force Field Descriptions of Nonbonded Interactions. J. Chem. Inf. Model. 2009, 49, 944−955. (9) Zgarbová, M.; Otyepka, M.; Sponer, J.; Hobza, P.; Jurecka, P. Large-Scale Compensation of Errors in Pairwise-Additive Empirical Force Fields: Comparison of AMBER Intermolecular Terms with Rigorous DFT-SAPT Calculations. Phys. Chem. Chem. Phys. 2010, 12, 10476−10493. (10) Warshel, A.; Levitt, M. Theoretical Studies of Enzimic Reactions: Dielectric Electrostatic and Steric Stabilization of the Carbonium Ion in the Reaction of Lysozyme. J. Mol. Biol. 1976, 103, 227−249. (11) Hobza, P. Calculations on Noncovalent Interactions and Databases of Benchmark Interaction Energies. Acc. Chem. Res. 2012, 45, 663−672. D
dx.doi.org/10.1021/jp508800k | J. Phys. Chem. A XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry A
Article
(12) Sedlak, R.; Janowski, T.; Pitoňaḱ , M.; Rezác,̌ J.; Pulay, P.; Hobza, P. The Accuracy of Quantum Chemical Methods for Large Noncovalent Complexes. J. Chem. Theory Comput. 2013, 9, 3364− 3374. (13) Dessent, C. E. H.; Müller-Dethlefs, K. Hydrogen-Bonding and van der Waals Complexes Studied by ZEKE and REMPI Spectroscopy. Chem. Rev. 2000, 100, 3999−4022. (14) Cao, Q.; Andrijchenko, N.; Ahola, A. E.; Domanskaya, A.; Rasanen, M.; Ermilov, A.; Nemukhin, A.; Khriachtchev, L. Interaction of Phenol with Xenon and Nitrogen: Spectroscopic and Computational Characterization. J. Chem. Phys. 2012, 137, 134305. (15) Riley, K. E.; Hobza, P. Assessment of the MP2 Method, along with Several Basis Sets, for the Computation of Interaction Energies of Biologically Relevant Hydrogen Bonded and Dispersion Bound Complexes. J. Phys. Chem. A 2007, 111, 8257−8263. (16) Andrijchenko, N. N.; Ermilov, A.Yu. Using the DFT-D Method to Describe Dispersion Interactions in Systems of Weakly-Bonded Xe−Aromatic Molecules. Russ. J. Phys. Chem. A 2012, 87, 1342−1348. (17) Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787−1799. (18) Furukawa, H.; Singh, S. K.; Mancusso, R.; Gouaux, E. Subunit Arrangement and Function in NMDA Receptors. Nature 2005, 438, 185−192. (19) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graph. 1996, 14, 33−38. (20) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (21) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kalé, L.; Schulten, K. Scalable Molecular Dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781− 1802. (22) Best, R. B.; Zhu, X.; Shim, J.; Lopes, P. E. M.; Mittal, J.; Feig, M.; MacKerell, A. D. Optimization of the Additive CHARMM AllAtom Protein Force Field Targeting Improved Sampling of the Backbone ϕ, ψ and Side-Chain χ1 and χ2 Dihedral Angles. J. Chem. Theory Comput. 2012, 8, 3257−3273. (23) Valiev, M.; Bylaska, E. J.; Govind, N.; Kowalski, K.; Straatsma, T. P.; Van Dam, H. J. J.; Wang, D.; Nieplocha, J.; Apra, E.; Windus, T.; et al. NWChem: A Comprehensive and Scalable Open-Source Solution for Large Scale Molecular Simulations. Comput. Phys. Commun. 2010, 181, 1477−1489. (24) Vemparala, S.; Domene, C.; Klein, M. L. Computational Studies on the Interactions of Inhalational Anesthetics with Proteins. Acc. Chem. Res. 2010, 43, 103−110. (25) Voevodin, Vl. V.; Zhumatiy, S. A.; Sobolev, S. I.; Antonov, A. S.; Bryzgalov, P. A.; Nikitenko, D. A.; Stefanov, K. S.; Voevodin, Vad. V. Practice of “Lomonosov” Supercomputer. Open Systems J. - Moscow 2012, 36−39.
E
dx.doi.org/10.1021/jp508800k | J. Phys. Chem. A XXXX, XXX, XXX−XXX
Toward Molecular Mechanism of Xenon Anesthesia: A Link to Studies of Xenon Complexes with Small Aromatic Molecules Natalya N. Andrijchenko,† Alexander Yu. Ermilov,† Leonid Khriachtchev,‡ Markku Ras̈ an̈ en,‡ and Alexander V. Nemukhin*,†,§ †
Department of Chemistry, M. V. Lomonosov Moscow State University, Leninskie Gory 1/3, Moscow 119991, Russian Federation Department of Chemistry, University of Helsinki, P.O. Box 55, Helsinki FI-00014, Finland § N. M. Emanuel Institute of Biochemical Physics, Russian Academy of Sciences, Kosygina 4, Moscow 119334, Russian Federation ‡
ABSTRACT: The present study illustrates the steps toward understanding molecular mechanism of xenon anesthesia by focusing on a link to the structures and spectra of intermolecular complexes of xenon with small aromatic molecules. A primary cause of xenon anesthesia is attributed to inhibition of N-methyl-D-aspartate (NMDA) receptors by an unknown mechanism. Following the results of quantum mechanics/molecular mechanics (QM/MM) and molecular dynamics (MD) calculations we report plausible xenon action sites in the ligand binding domain of the NMDA receptor, which are due to interaction of xenon atoms with aromatic amino-acid residues. We rely in these calculations on computational protocols adjusted in combined experimental and theoretical studies of intermolecular complexes of xenon with phenol. Successful reproduction of vibrational shifts in molecular species upon complexation with xenon measured in low-temperature matrices allowed us to select a proper functional form in density functional theory (DFT) approach for use in QM subsystems, as well as to calibrate force field parameters for MD simulations. The results of molecular modeling show that xenon atoms can compete with agonists for a place in the corresponding protein cavity, thus indicating their active role in anesthetic action.
■
INTRODUCTION The role of relatively weak noncovalent interactions between atoms and molecules in many physical, chemical, and biological phenomena cannot be underestimated. Numerous experimental and theoretical studies of intermolecular complexes are related to distinctively different directions, among which we emphasize those devoted to interactions of rare gas xenon with small organic molecules. The results of these investigations find important biomedical applications such as using xenon as a marker in protein cavities1 or understanding mechanism of xenon anesthesia.2 A link between the molecular mechanism of xenon anesthesia and subtle details of intermolecular complexes of xenon may be outlined as follows. A primary cause of xenon action is attributed to inhibition of N-methyl-D-aspartate (NMDA) receptors, multidomain membrane proteins, three-dimensional atomic structures of which have become available only recently.3 If xenon atoms are capable of occupying places in a ligand binding cavity and competing with the natural agonist, glycine, then the mechanism of action is getting clearer. In turn, the walls of such protein cavities contain amino acid residues with aromatic rings, tyrosine (Tyr), phenylalanine (Phe), and tryptophan (Trp), an interaction with which is a fingerprint of the resolved xenon−protein structures deposited in the Protein Data Bank. A typical motif of these complexes is a pattern with a xenon atom residing aside the plane of an aromatic ring; the distances between xenon and carbon atoms vary between 3.5 and 5 Å. For brevity, we call these complexes the π-type © XXXX American Chemical Society
structures; they represent a part of noncovalently bound systems with aromatic rings.4 From the theoretical side, straightforward simulations of xenon atoms surrounded by organic molecules using molecular dynamics (MD) calculations5−7 face a well-known problem of reliability of force field parameters8,9 required to estimate the energy of intermolecular complexes involving xenon. Quantum-based models seem to be more promising in this respect; however, the necessity to include large fractions of protein cavity walls to the model systems encourages one to use the quantum mechanics/ molecular mechanics (QM/MM) approaches.10 Apparently, xenon and the nearest molecular groups must be assigned to the quantum subsystem, and a proper accuracy level11,12 at a reasonable computational cost must be provided to describe the corresponding complexes. Experimentally, complexes of rare gas atoms with small organic molecules can be studied in the gas phase13 and in the low-temperature matrices.14 In the latter case, subtle changes in vibrational bands of a molecule due to its complexation with xenon atoms can be accurately measured. In this computational work, we use the data obtained in matrix experiments with xenon and aromatic molecules to select a proper functional form in the density functional theory Special Issue: Markku Räsänen Festschrift Received: August 31, 2014 Revised: September 30, 2014
A
dx.doi.org/10.1021/jp508800k | J. Phys. Chem. A XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry A
Article
(DFT) approach for quantum-based characterization of π-type xenon complexes, as well as to calibrate force field parameters for MD simulations. The QM/MM and MD modeling tools allow us to determine the trapping sites of xenon in the ligand binding domain of N-methyl- D -aspartic acid (NMDA) receptors.
in spite of the presence of water molecules competing for the spots near amino acid residues forming the walls of the cavity. The analysis of conformations along MD trajectories shows at least three trapping sites occupied by xenon inside the cavity. To characterize these sites, we performed QM/MM calculations with the DFT(BHHLYP)-D2/6-31G* approach to estimate the energies and forces in the QM subsystems. The NWChem program suite was employed in these calculations.23 Specification of the QM parts for each trapping site is given in the Results. The starting coordinates for the QM/MM minimization were generated from appropriate frames along MD trajectories. This strategy did not guarantee that all sound equilibrium geometry configurations could be located; however, our goal was to find at least a few stable structures with a xenon atom trapped inside the cavity when the glycine molecule was replaced.
■
MODELS AND METHODS To model interactions of the rare-gas element xenon with small biologically related organic molecules, it is recommended to apply approaches capable to provide a well balanced description of noncovalent interactions, such as the second-order Møller− Plesset perturbation theory, MP2.11,15 However, attempts to use MP2 in QM/MM simulations of systems modeling xenon in protein clefts would lead to prohibitively expensive calculations. Therefore, we applied a two-step strategy to select a modeling tool. First, we showed that the MP2/aug-cc-pVTZ approach allowed us to reproduce vibrational spectra of the phenol···Xen (n = 2−4) complexes examined by infrared absorption spectroscopy in neon matrices.14 Second, we used the results of equilibrium geometry configurations and binding energies of the complexes of xenon with small aromatic molecules obtained in the MP2 calculations as reference data to select an inexpensive DFT-based method. It was demonstrated16 that DFT with a dispersion correction of type D217 with the functional BHHLYP and a modest basis set 6-31G* was able to approach the MP2 data. In particular, the DFT-D2 computed binding energy for the Xe−phenol π-type complex was 2.7 kcal mol−1, which is similar to 3.1 kcal mol−1 attainable in the MP2/cc-pVTZ calculations. The deviations of the distances between Xe and carbon atoms obtained in the DFTD2 calculations from the reference values did not exceed 0.05 Å. The DFT(BHHLYP)-D2 method also reproduced the differences between the two structures (π-type and H-bonded) of the Xe−phenol complex, as well as the trend toward an increase in binding energy of Xe in the series of aromatic amino acids Phe, Tyr, Trp.16 This theoretical approach was used for to calculate energies and forces in the QM subsystems in the QM/ MM calculations. To construct molecular models of a NMDA ligand binding domain, we used the coordinates of heavy atoms of the protein structure PDBid: 2A5T18 corresponding to the crystal of the ligand-binding core complex containing glycine. After hydrogen atoms and sodium and chloride ions with a concentration of up to 0.15 M were added, via the VMD software,19 energy minimization and a few MD runs were carried out. Water molecules recognized in the crystal structure were included in the model system, and the TIP3P water molecules20 formed a solvation shell for the protein. The MD simulations were performed with the NAMD program21 using the CHARMM force field parameters.22 During MD simulations of about 20 ns length, the systems were kept at a constant temperature of 298 K and pressure of 1 atm (NPT ensemble); the periodical boundary conditions were imposed. At the next step, the glycine moiety was replaced by a xenon atom, and energy minimization and MD runs (up to 50 ns) were repeated. In these calculations the CHARMMM force field parameters were corrected by improving the Lennard-Jones interaction potential (ε = 1.7 kcal mol−1) for the pairs containing xenon in such a manner that the equilibrium geometry configurations of the Xe−phenol complex obtained in the MP2/aug-cc-pVTZ calculations could be reproduced in MM simulations. Importantly, the xenon atom did not escape from the cavity
■
RESULTS MD Simulations. Ten classical molecular dynamics trajectories extended up to 50 ns were executed from the initial coordinates generated by replacement of glycine by a xenon atom in the model system constructed by motifs of the crystal structure PDBid: 2A5T.18 Eight of ten MD runs showed that Xe migrated between three local areas inside the protein cavity called below trapping sites I, II, and III. In two cases, a firm assignment of a xenon atom to a particular spot could not be accomplished during the simulation period. Figure 1 shows
Figure 1. Trapping sites of a xenon atom migrating inside the cavity of the ligand-binding domain along a representative MD trajectory.
the location of trapping sites I, II, III in the cleft; Figure 2 illustrates migration of a xenon atom, initially placed near the Phe92 residue (similarly to the glycine molecule in the crystal structure PDBid: 2A5T), along a representative trajectory. An important conclusion that can be drawn from these MD calculations is that xenon does not escape from the cavity despite the latter being filled by water molecules (not shown in Figure 1), which also occurs for the binding sites with amino acid residues forming walls of the cleft. In qualitative agreement with the results of previous classical MD simulations5,6 our data evidence that xenon can compete with glycine in the ligand binding domain of the NMDA receptor. To characterize possible trapping sites of xenon atoms inside the cleft, we applied a higher level theoretical approach, namely, the QM/ MM theory. B
dx.doi.org/10.1021/jp508800k | J. Phys. Chem. A XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry A
Article
arrangement is typical for complexes of xenon with phenol studied in our previous papers.14,16 According to the MD simulations (Figure 2), a xenon atom can migrate to another site leaving the initial position near Phe92. In the calculations for trapping site II, the QM subsystem included Xe, side chains of Gln13, Phe92, Pro124, Ser180, Trp223, Asp224, Val227, Phe250, and one water molecule. The equilibrium geometry configuration corresponding to trapping site II is illustrated in Figure 4. In this case, xenon atom is still bound to Phe92 but tends to find an additional binding partner, this time, another residue of aromatic nature, Trp223.
Figure 2. Root mean square deviations (RMSD) values along a MD trajectory computed from the initial position of Xe, which corresponds to replacement of glycine in the crystal structure PBDid: 2A5T by a xenon atom.
QM/MM Simulations. According to the MD results, a xenon atom migrates from one local environment to another; therefore, it is impossible to select molecular groups constituting a quantum subsystem for QM/MM calculations in the same manner for three possible trapping sites. In every selection of the QM subsystem, we followed a goal to completely surround a xenon atom by environmental molecular groups keeping an amount of quantum atoms to about 100. To generate initial coordinates for the QM/MM minimization, we used appropriate frames of MD trajectories assigned to areas I, II, or III. For trapping site I, the QM subsystem included Xe, side chains of Phe92, Leu125, Thr126, Ser180, Trp223, Phe250, and one water molecule. The equilibrium geometry configuration is illustrated in Figure 3. Apparently, the xenon atom forms a πtype complex with the aromatic ring of Phe92. This
Figure 4. Structure of trapping site II.
For trapping site III shown in Figure 5, the QM subsystem included Xe, side chains of Gln13, Phe16, Phe92, Pro124, Gln144, Asp224, Ala226, Val227, Phe250, and three water molecule. In this site, the side chain of Phe16 is the primary partner to form the π-type complex with a xenon atom.
Figure 3. Structure of trapping site I. Here and in all figures, the atoms of the QM subsystem are shown in balls and sticks; carbon atoms are green, oxygen red, nitrogen blue; and the distances are given in Å. C
dx.doi.org/10.1021/jp508800k | J. Phys. Chem. A XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry A
Article
the NMDA protein cavity can be identified, in which xenon is located. The motif of the π-type complex between the xenon atom and the aromatic molecule is typical for these trapping sites. Correspondingly, xenon atoms can compete with agonist species (glycine) for a place in the ligand binding domain of a protein receptor, thus indicating its role in anesthetic action. The computational method was calibrated by using the experimental data on the xenon−phenol complexes in lowtemperature matrices, which demonstrates the benefits of matrix-isolation spectroscopy for biology-oriented research.
■
AUTHOR INFORMATION
Corresponding Author
*A. V. Nemukhin. Phone: +7-495-939-10-96. E-mail: [email protected]; [email protected]. Notes
Figure 5. Structure of trapping site III.
The authors declare no competing financial interest.
■
■
CONCLUDING REMARKS The results of MD and QM/MM simulations show, first, that a xenon atom can be trapped in the ligand binding domain replacing glycine and, second, that the motif of the π-type complex is typical for these trapping sites. Analysis of biochemical aspects of xenon anesthesia is far beyond the scope of the present work. We rely on the knowledge that, most likely, xenon serves as an antagonist of the NMDA receptor inhibiting the receptor at its glycine binding site. Analysis of the available crystal structure of the ligand binding domain of the NMDA receptor shows that xenon can be bound to the side chains of aromatic amino acid residues, consistent with the molecular structures of xenon− protein complexes deposited to the protein data bank. Therefore, detailed knowledge on the structure and properties of intermolecular complexes of xenon atoms with the corresponding parent aromatic molecules like phenol, cresol, and toluene is essential when the interaction of xenon with protein receptors is modeled. We have shown that the subtle changes in vibrational frequencies of xenon complexes measured in low-temperature matrices are well reproduced by the MP2/aug-cc-pVTZ approach. In turn, the equilibrium geometry configurations and binding energies of this theoretical level are well reproduced using the much less expensive DFT-D2 approach17 with the functional BHHLYP and the basis set 6-31G*. Therefore, we can rely on the QM/MM calculations describing xenon binding sites in the NMDA receptor using the DFT(BHHLYP)-D2/6-31G* approximation in relatively large quantum subsystems. The results of these QM/MM calculations demonstrate that xenon can reside in several trapping sites showing patterns of π-type complexes. The force field parameters for the MD simulations were also calibrated by the structures of the π-type complexes of xenon with phenol. Specifically, the parameter ε of the Lennard-Jones interaction potential was adjusted to reproduce the MP2/augcc-pVTZ data. The large-scale MD calculations of the ligandbinding domain of the NMDA receptor will be reported elsewhere because their description requires a specific presentation typical for biochemistry oriented journals, e.g., ref 24. In this work, we used relatively short MD trajectories prior the expensive QM/MM minimization procedures to account for suitable arrangement of protein conformations. To conclude, the results of the present MD and QM/MM calculations show that several plausible xenon binding sites in
ACKNOWLEDGMENTS We thank Bella Grigorenko for discussion and advice. This study was supported by the Russian Science Foundation (project 14-13-00124). We acknowledge the use of supercomputer resources of the M. V. Lomonosov Moscow State University25 and of the Joint Supercomputer Center of the Russian Academy of Sciences.
■
REFERENCES
(1) Oros, A. M.; Shah, N. J. Hyperpolarized Xenon in NMR and MRI. Phys. Med. Biol. 2004, 49, R105−R153. (2) Derwall, M.; Coburn, M.; Rex, S.; Hein, M.; Rossaint, R.; Fries, M. Xenon: Recent Developments and Future Perspectives. Minerva Anestesiol. 2009, 75, 37−45. (3) Sobolevsky, A. I.; Rosconi, M. P.; Eric Gouaux, E. X-Ray Structure, Symmetry and Mechanism of an AMPA-Subtype Glutamate Receptor. Nature 2009, 462, 745−758. (4) Li, S.; Xu, Y.; Shen, Q.; Liu, X.; Lu, J.; Chen, Y.; Lu, T.; Luo, C.; Luo, X.; Zheng, M.; et al. Non-Covalent Interactions with Aromatic Rings: Current Understanding and Implications for Rational Drug Design. Curr. Pharm. Des. 2013, 19, 6522−6533. (5) Dickinson, R.; Peterson, B. K.; Banks, P.; Simillis, C.; Martin, J. C.; Valenzuela, C. A.; Maze, M.; Franks, N. P. Competitive Inhibition at the Glycine Site of the N-Methyl-D-Aspartate Receptor by the Anesthetics Xenon and Isoflurane: Evidence from Molecular Modeling and Electrophysiology. Anesthesiology 2007, 107, 756−767. (6) Liu, L. T.; Xu, Y.; Tang, P. Mechanistic Insights into Xenon Inhibition of NMDA Receptors from MD Simulations. J. Phys. Chem. B 2010, 114, 9010−9016. (7) Weinrich, M.; Worcester, D. L. Xenon and other Volatile Anesthetics Change Domain Structure in Model Lipid Raft Membranes. J. Phys. Chem. B 2013, 117, 16141−16147. (8) Paton, R. S.; Goodman, J. M. Hydrogen Bonding and PiStacking: How Reliable are Force Fields? A Critical Evaluation of Force Field Descriptions of Nonbonded Interactions. J. Chem. Inf. Model. 2009, 49, 944−955. (9) Zgarbová, M.; Otyepka, M.; Sponer, J.; Hobza, P.; Jurecka, P. Large-Scale Compensation of Errors in Pairwise-Additive Empirical Force Fields: Comparison of AMBER Intermolecular Terms with Rigorous DFT-SAPT Calculations. Phys. Chem. Chem. Phys. 2010, 12, 10476−10493. (10) Warshel, A.; Levitt, M. Theoretical Studies of Enzimic Reactions: Dielectric Electrostatic and Steric Stabilization of the Carbonium Ion in the Reaction of Lysozyme. J. Mol. Biol. 1976, 103, 227−249. (11) Hobza, P. Calculations on Noncovalent Interactions and Databases of Benchmark Interaction Energies. Acc. Chem. Res. 2012, 45, 663−672. D
dx.doi.org/10.1021/jp508800k | J. Phys. Chem. A XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry A
Article
(12) Sedlak, R.; Janowski, T.; Pitoňaḱ , M.; Rezác,̌ J.; Pulay, P.; Hobza, P. The Accuracy of Quantum Chemical Methods for Large Noncovalent Complexes. J. Chem. Theory Comput. 2013, 9, 3364− 3374. (13) Dessent, C. E. H.; Müller-Dethlefs, K. Hydrogen-Bonding and van der Waals Complexes Studied by ZEKE and REMPI Spectroscopy. Chem. Rev. 2000, 100, 3999−4022. (14) Cao, Q.; Andrijchenko, N.; Ahola, A. E.; Domanskaya, A.; Rasanen, M.; Ermilov, A.; Nemukhin, A.; Khriachtchev, L. Interaction of Phenol with Xenon and Nitrogen: Spectroscopic and Computational Characterization. J. Chem. Phys. 2012, 137, 134305. (15) Riley, K. E.; Hobza, P. Assessment of the MP2 Method, along with Several Basis Sets, for the Computation of Interaction Energies of Biologically Relevant Hydrogen Bonded and Dispersion Bound Complexes. J. Phys. Chem. A 2007, 111, 8257−8263. (16) Andrijchenko, N. N.; Ermilov, A.Yu. Using the DFT-D Method to Describe Dispersion Interactions in Systems of Weakly-Bonded Xe−Aromatic Molecules. Russ. J. Phys. Chem. A 2012, 87, 1342−1348. (17) Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787−1799. (18) Furukawa, H.; Singh, S. K.; Mancusso, R.; Gouaux, E. Subunit Arrangement and Function in NMDA Receptors. Nature 2005, 438, 185−192. (19) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graph. 1996, 14, 33−38. (20) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (21) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kalé, L.; Schulten, K. Scalable Molecular Dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781− 1802. (22) Best, R. B.; Zhu, X.; Shim, J.; Lopes, P. E. M.; Mittal, J.; Feig, M.; MacKerell, A. D. Optimization of the Additive CHARMM AllAtom Protein Force Field Targeting Improved Sampling of the Backbone ϕ, ψ and Side-Chain χ1 and χ2 Dihedral Angles. J. Chem. Theory Comput. 2012, 8, 3257−3273. (23) Valiev, M.; Bylaska, E. J.; Govind, N.; Kowalski, K.; Straatsma, T. P.; Van Dam, H. J. J.; Wang, D.; Nieplocha, J.; Apra, E.; Windus, T.; et al. NWChem: A Comprehensive and Scalable Open-Source Solution for Large Scale Molecular Simulations. Comput. Phys. Commun. 2010, 181, 1477−1489. (24) Vemparala, S.; Domene, C.; Klein, M. L. Computational Studies on the Interactions of Inhalational Anesthetics with Proteins. Acc. Chem. Res. 2010, 43, 103−110. (25) Voevodin, Vl. V.; Zhumatiy, S. A.; Sobolev, S. I.; Antonov, A. S.; Bryzgalov, P. A.; Nikitenko, D. A.; Stefanov, K. S.; Voevodin, Vad. V. Practice of “Lomonosov” Supercomputer. Open Systems J. - Moscow 2012, 36−39.
E
dx.doi.org/10.1021/jp508800k | J. Phys. Chem. A XXXX, XXX, XXX−XXX
