v2i3-rohani

Research Paper

Vibrational Frequency Calculations of Phylloquinone in the A1 Binding Site:

Layer Selection in ONIOM Methods

Leyla Rohani and Gary Hastings*

Department of Physics and Astronomy, Georgia State University, Atlanta, Georgia 30302, USA.

*Corresponding author, Gary Hastings, Email: ghastings@gsu.edu

Received 20 April 2018, revised 31 May 2018, accepted 1 June 2018

Publication Date (Web): June 1, 2018

© Frontiers in Science, Technology, Engineering and Mathematics

Abstract

Multi-layer ONIOM methods have been used to model pigment-protein interactions for the phylloquinone (PhQ) molecule that occupies the A1 binding site in photosystem I (PSI).  In the molecular models considered, all amino acids within 10 Å of both carbonyl oxygen atoms of the PhQ were included.  Two classes of multi-layer methods were considered. In the first, PhQ and several nearby amino acid residues were considered at a single quantum mechanical (QM) level of theory.  The rest of the atoms were considered at a molecular mechanics (MM) level.  Second, PhQ and several nearby amino acid residues were modeled using two separate QM levels of theory.  The remainder of the atoms again being considered at the MM level.  Vibrational frequency calculations were undertaken for the optimized QM layers, for PhQ in both the neutral and anion states, allowing the construction of so-called [Anion – Neutral] infrared difference spectra (DS).  The calculations suggest that the three-layer method is superior to the two-layer method. Calculated spectra obtained using the three-layer method are considerably simpler to assess since only vibrational modes associated with the quinone need be considered.  In addition, in the three-layer method, vibrational mode potential energy distributions can be calculated allowing a direct assessment that individual molecular groups of the quinone molecule make to the normal mode in question.

Keywords

Photosystem I, Phylloquinone, A1, density functional theory, ONIOM, QM:MM methods, FTIR difference spectra.

Introduction

A PhQ molecule (2-methyl-3-phytyl-1,4-naphthoquinone) occupies the A1 binding site in PSI. PhQ in PSI is one of the most reducing quinones in biology (Srinivasan et al 2009). This unique characteristic is a result of interactions of PhQ with its surrounding protein environment. The crystal structure suggests that the C1=O group of PhQ in the A1binding site is free from hydrogen bonding (H-bonding), whereas the C4=O group is H-bonded to the backbone NH group of a leucine residue (Fig. 1). In this study, we have used two and three-layer ONIOM (Vreven et al 2006) methods to geometry optimize molecular models associated with PhQ in the A1 binding site.  Bond lengths and angles associated with atoms involved in H-bonding are compared in the two methods. Calculated [Anion – Neutral] FTIR DS are also compared. 

Previously, we have shown that the calculated vibrational properties of pigments in the gas-phase or in solution do not adequately model the vibrational properties of the same pigments in a protein binding site (Lamichhane et al 2011; Lamichhane et al 2013; Makita et al 2017). Any realistic calculation needs to involve the protein environment at some level.  On the other hand, despite the development of computational chemistry methods using multi-core “supercomputers”, it is still challenging and time-consuming (expensive) to calculate properties of a large protein molecular system using high-level quantum mechanical (QM) methods.  One solution to this problem is to consider the important molecular features of a system at a QM level of theory, while treating more peripheral atoms and molecular groups at a lower, less-expensive level of theory, usually using a molecular mechanics (MM) method such as UFF (Rappe et al 1992). 

 is the energy of a two-layer ONIOM model,  is the energy of the entire system at the MM level, 

 is the calculated energy for the intermediate layer at MM level,  is the calculated energy for the model at QM-high level, and  is the calculated energy for the model system at QM-low level.

Fig. 2 shows schematically the layer selection in the ONIOM-2 and -3 methods. The PhQ molecule and several selected residues were treated at a QM level in the ONIOM-2 method.  Extending upon the ONIOM-2 method, the PhQ molecule by itself is treated at a higher QM level, with selected nearby residues treated at a lower QM level. After optimization the QM layer of ONIOM-2 and the highest QM layer of ONIOM-3 were extracted and used in harmonic normal mode vibrational frequency calculations. Calculations were undertaken for the PhQ in both the neutral and anion states allowing [Anion – Neutral] infrared DS to be calculated and compared to corresponding FTIR DS.

Figure 2. Partitioning of the system and layer definition in ONIOM-2 and -3 methods (Chung et al 2015).

Materials and Methods

Molecular Model Construction

   Molecular models were constructed using the x-ray structure of cyanobacterial PSI at 2.5 Å resolution (PDB entry 1JB0) (Jordan et al 2001). All atoms within 10 Å of both oxygen atoms of PhQ in the A1A binding site were selected. PhQ is the only cofactor present in the constructed model. The phytyl tail of PhQ was truncated to a 5-carbon CH2C2H(CH3)2 unit. Hydrogen atoms were added to the molecular model using GaussView 6 (Gaussian 2016). From this molecular structure, PhQ and parts of its immediate protein environment were selected to be treated at a QM level. The residues were chosen so as to be similar to that used in ONIOM calculations of the magnetic spectroscopic parameters of PhQ in PSI (Lin et al 2011). LeuA722, TrpA697, PheA689, SerA692 and MetA688, in addition to part of AlaA721 and SerA723 are considered at the QM level in ONIOM-2 calculations.  The protonation state of amino acid side chains (neutral pH and room temperature) were checked and adjusted where necessary.

Computational Details

      For geometry optimization, heavy atoms of the protein side-chains and backbone are constrained.  The heavy atoms of the quinone, as well as all hydrogen atoms, are unconstrained. Electronic embedding was employed (Thompson et al 2014). The link atom approach was used to link the QM and MM layers (Chung et al 2015). For ONIOM-2, following geometry optimization, the QM layer was extracted, and missing hydrogen atoms were added. This extracted layer with added hydrogens was re-optimized at the same QM level with all heavy atoms constrained.  The final optimized structure was then used in vibrational frequency calculations at the same QM level as the optimization.  In ONIOM-2, QM level calculations were undertaken using density functional theory (DFT) employing the B3LYP functional and the 6-31G(d) basis set.

      For ONIOM-3 calculations, the quinone, which is treated at the highest QM level, is extracted and vibrational frequency calculations are undertaken at the same QM level.  For the semiquinone (anion) state, the charge on the high and intermediate QM layers is set to -1.  In ONIOM-3 calculations the quinone was treated using DFT employing the B3LYP functional and the 6-31G+(d) basis set. The lower QM layer was treated using the B3LYP functional and the 6-31G(d) basis set.  All other atoms were treated at the MM level using UFF (Rappe et al 1992).

In ONIOM-3 calculations the protein backbone and sidechain atoms are constrained.  However, the calculated spectra based on ONIOM-3 method do not depend much on whether the heavy atoms of the amino acid side-chains are unconstrained or constrained (Makita et al 2017). This is a useful observation indicating one possible approach to decrease computational cost.

      All calculations described here were undertaken using Gaussian16 software (Gaussian 2016). The atomic displacements calculated for each normal mode can be animated using GaussView 6.  The most prominent contribution of molecular groups to these normal modes can be assessed visually.

Results and Discussion

Calculated Geometric Parameters

In the x-ray crystal structure, the distance between the carbonyl oxygen atom of PhQ (attached to C4-Fig. 1) and the backbone nitrogen atom of LeuA722 is 2.694 _.  Given this distance such H-bond would be considered as moderately strong, and contain mostly contributions from electrostatic interactions (Minch et al 1999). In the following we will use the nomenclature ‘A-H…B’ to describe the three atoms involved in H-bonding, with A/B being hydrogen donor/acceptor, respectively.  ‘A-H…B’ in our case is N-H…O4.

Figure 3. ONIOM-2 (A) and -3 (B) molecular models used in calculations. Hydrogen atoms are not shown. The two/three layers are obvious in (A)/(B).

For H-bonding the bonding angles should also be considered (Minch et al 1999). Bond lengths and angles calculated using ONIOM-2 and -3 methods are listed in Table 1.

The distance between O4 (of PhQ) and N (of LeuA722) is increased in both neutral optimized structures compared to the distance found in the crystal structure, although less so in ONIOM-3 calculations.  In ONIOM-3 optimized structures the N…O4 distance decreases 0.012/0.109 Å for PhQ_/PhQ, respectively.  For ‘A-H…B’ type H-bonds, bond angles between 130° and 180° usually indicate moderately strong H-bonds, while angles greater than 175° are classified as strong H-bonds (Minch et al 1999). Both ONIOM methods, for both neutral and reduced states would therefore indicate strong H-bonding.

Table 1. Bond lengths (in Å) and angles for the PhQ N-H…O4 H-bond calculated using ONIOM-2 and -3 methods. Neutral/anion state are in black/red text, respectively.

Calculated Spectra

     Vibrational frequency calculations were used to construct [Anion – Neutral] FTIR DS. The calculated FTIR DS obtained using ONIOM-2 and -3 methods are compared in Fig. 4.  Because atoms of the protein were constrained in the geometry optimization, several negative frequencies were calculated. In general, when negative frequencies are calculated the optimized geometry is at a saddle point on the potential energy surface. However, if the calculated negative frequency is low, the geometry is converged to a local minimum.Calculated negative frequencies and their intensities from both ONIOM models are listed in Table 2.  The number of calculated negative frequencies (and their intensity) is dramatically decreased in ONIOM-3 calculations. The highest calculated negative frequency for PhQ_ in the ONIOM-3 calculation is at -180 cm-1, and has very low intensity. It is mainly due to a twisting mode of a methyl group and does not involve the PhQ C=O or C=C groups (Table 2). 

     It is unlikely that employing methods to produce an optimized structure that is free of negative frequencies will alter the spectra in Fig. 4, and it is therefore reasonable to ignore these calculated negative frequencies. Previously we have undertaken calculations in which different sets of atoms were constrained or unconstrained during geometry optimization, and have found that differing numbers of negative frequencies are calculated.  However, the calculated infrared spectra in the spectral regions of interest (1300-1800 cm-1) were unaltered, reinforcing the idea that it is reasonable to ignore calculated negative frequencies, particularly if the frequencies are low and one is interested in high frequency, high intensity, normal modes.

Table 2. Calculated negative frequencies in both ONIOM methods. neutral/anion state are in black/red text, respectively.

Figure 4. Calculated FTIR DS obtained using ONIOM-2 and ONIOM-3 methods. Calculated frequencies are scaled by 0.966/0.977 for the neutral/reduced state in ONIOM-3, respectively. Frequencies are scaled by 0.958/0.964 for the neutral/reduced state in ONIOM-2, respectively.  These scale factors fix the position of the C1=O mode of the neutral/anion state at 1654/1495 cm-1, respectively. 

Figure 5.  Visual assessment of the selected vibrational modes of PhQ using GaussView 6. The red arrow shows the C=O group's vibrations of PhQ. Only the quinone molecule is included at the highest QM level in ONIOM-3. Vibrational frequencies of quinone are coupled to many vibrational groups of the protein (presented by the thin lines) that were treated at the QM level. Neutral/anion state are in black/red text, respectively (See Table 3).

In both ONIOM-2 and -3 calculations, for both neutral and reduced PhQ, the C4=O mode intensity is considerably less than the C1=O mode intensity. 

In ONIOM-2/3 calculations for neutral PhQ, the C1=O mode is 17/25 cm-1 higher than the C4=O mode, respectively. For PhQ_ the C=O mode separation is 63-65 cm-1.  Why it is so much higher for the anion state compared to the neutral state is not entirely clear but one suggestion seems to be that the H-bonding comes on strong in the anion state.

In ONIOM-3 calculations only the quinone molecule is included at the highest QM level.  This QM level involves the use of diffuse functions (the “+” in the 6-31G(+)d basis set), which are known to be important for accurate calculation of vibrational modes of anion states (Bandaranayake et al 2006). In ONIOM-2 calculations, because many atoms in the vicinity of PhQ are included at the QM level, the calculation of vibrational frequencies becomes more expensive, and thus a lower level basis not involving diffuse functions is used.

Fig. 5 indicates that vibrational frequency analysis is difficult using the ONIOM-2 method as the calculated vibrational frequencies of quinone are coupled to many vibrational groups of the protein that were treated at the QM level.  We do note, however, that this complication can be minimized when ONIOM-2 calculations are used to construct isotope edited FTIR DS.  This was demonstrated previously where ONIOM-2 methods were used to calculate [18O – 16O] isotope edited FTIR double difference spectra for PSI with PhQ incorporated (Hastings et al 2015).

Another advantage of the ONIOM-3 method that also stems from the fact that the quinone atoms are treated separately at a higher QM level, is that it is possible to calculate potential energy distributions (Jamroz et al 2013) to quantitatively assess the contribution that different molecular groups of the quinone make to the normal modes. Normal mode analysis bases on calculated potential energy distributions for PhQ in PSI were presented previously (Makita et al 2017).

Conclusions

For neutral PhQ in the A1 binding site in PSI, ONIOM-3 calculations predict a C1=O mode near 1654 cm-1, with a lower intensity C4=O mode downshifted 25 cm-1, due to H-bonding.  In ONIOM-2 calculations the downshift is only 17 cm-1.  Thus the asymmetry in H-bonding is higher in ONIOM-3 calculations.  The mode compositions also differ in the two types of calculations, however.

For PhQ_ in PSI, both types of calculation predict a 63-65 cm-1 separation between the C1=O and C4=O groups.  Thus both calculations predict a very strong H-bond to the C4=O group.

Despite similarities in the results from both types of calculations the ONIOM-3 approach offers advantages allowing easier implementation and assessment of data.  The ONIOM-3 approach also allows for a quantitative assessment of calculated normal modes, rather than just a visual assessment.

Acknowledgements

This work was supported by grant number DE-SC0017937 from the Department of Energy to GH.  We acknowledge the use of Georgia State's research computing resources that are supported by Georgia State's Research Solutions.

References

Bandaranayake, K., Sivakumar, V., Wang, R. & Hastings, G. Modeling The A1 Binding Site In Photosystem I. Density Functional Theory For The Calculation Of “Anion – Neutral” FTIR Difference Spectra of Phylloquinone. Vib. Spectrosc. 42, 78-87, (2006).

Chung, L. W. et al. The ONIOM Method and Its Applications. Chemical Reviews 115, 5678-5796, (2015).

Gaussian 16 Rev. B.01 (Wallingford, CT, 2016).

Hastings, G. Vibrational spectroscopy of photosystem I. Biochimica et Biophysica Acta (BBA) - Bioenergetics 1847, 55-68, (2015).

Jamroz, M. H. Vibrational energy distribution analysis (VEDA): Scopes and limitations. Spectrochimica Acta Part a-Molecular and Biomolecular Spectroscopy 114, 220-230, (2013).

Jordan, P. et al. Three-dimensional structure of cyanobacterial photosystem I at 2.5 angstrom resolution. Nature 411, 909-917, (2001).

Lamichhane, H. P. & Hastings, G. Calculated vibrational properties of pigments in protein binding sites. Proceedings of the National Academy of Sciences of the United States of America 108, 10526-10531, (2011).

Lamichhane, H. P. & Hastings, G. (2013) Calculated Vibrational Properties of Ubisemiquinones. Computational Biology Journal, 11.

Lin, T. J. & O'Malley, P. J. Binding Site Influence on the Electronic Structure and Electron Paramagnetic Resonance Properties of the Phyllosemiquinone Free Radical of Photosystem I. Journal of Physical Chemistry B 115, 9311-9319, (2011).

Makita, H., Rohani, L., Zhao, N. & Hastings, G. Quinones in the A1 binding site in photosystem I studied using time-resolved FTIR difference spectroscopy. Biochimica et Biophysica Acta (BBA) - Bioenergetics 1858, 804-813, (2017).

Minch, M. J. An Introduction to Hydrogen Bonding (Jeffrey, George A.). Journal of Chemical Education 76, 759, (1999).

Rappe, A. K., Casewit, C. J., Colwell, K. S., Goddard, W. A. & Skiff, W. M. Uff, a Full Periodic-Table Force-Field for Molecular Mechanics and Molecular-Dynamics Simulations. J Am Chem Soc 114, 10024-10035, (1992).

Srinivasan, N. & Golbeck, J. H. Protein-cofactor interactions in bioenergetic complexes: The role of the A(1A) and A(1B) phylloquinones in Photosystem I. Biochimica et Biophysica Acta (BBA) - Bioenergetics 1787, 1057-1088, (2009).

Thompson, L. M. et al. Analytical Harmonic Vibrational Frequencies for the Green Fluorescent Protein Computed with ONIOM: Chromophore Mode Character and Its Response to Environment. J Chem Theory Comput 10, 751-766, (2014).

Vreven, T. et al. Combining quantum mechanics methods with molecular mechanics methods in ONIOM. J Chem Theory Comput 2, 815-826, (2006).

Citation:

Leyla Rohani and Gary Hastings (2018) Vibrational Frequency Calculations of Phylloquinone in the A1 binding site: Layer Selection in ONIOM Methods, Frontiers in Science, Technology, Engineering and Mathematics, Volume 2, Issue 3, 130-137