Anharmonic Vibrational Analysis of Biomolecules and Solvated Molecules Using Hybrid QM/MM Computations.

Quantum mechanics/molecular mechanics (QM/MM) calculations are applied for anharmonic vibrational analyses of biomolecules and solvated molecules. The QM/MM method is implemented into a molecular dynamics (MD) program, GENESIS, by interfacing with external electronic structure programs. Following the geometry optimization and the harmonic normal-mode analysis based on a partial Hessian, the anharmonic potential energy surface (PES) is generated from QM/MM energies and gradients calculated at grid points. The PES is used for vibrational self-consistent field (VSCF) and post-VSCF calculations to compute the vibrational spectrum. The method is first applied to a phosphate ion in solution. With both the ion and neighboring water molecules taken as a QM region, IR spectra of representative hydration structures are calculated by the second-order vibrational quasi-degenerate perturbation theory (VQDPT2) at the level of B3LYP/cc-pVTZ and TIP3P force field. A weight-average of IR spectra over the structures reproduces the experimental spectrum with a mean absolute deviation of 16 cm-1. Then, the method is applied to an enzyme, P450 nitric oxide reductase (P450nor), with the NO molecule bound to a ferric (FeIII) heme. Starting from snapshot structures obtained from MD simulations of P450nor in solution, QM/MM calculations have been carried out at the level of B3LYP-D3/def2-SVP(D). The spin state of FeIII(NO) is likely a closed-shell singlet state based on a ratio of N-O and Fe-NO stretching frequencies (νN-O and νFe-NO) calculated for closed- and open-shell singlet states. The calculated νN-O and νFe-NO overestimate the experimental ones by 120 and 75 cm-1, respectively. The electronic structure and solvation of FeIII(NO) affect the structure around the heme of P450nor leading to an increase in νN-O and νFe-NO.

Introduction

Vibrational spectroscopy is one of the valuable methods to elucidate molecular mechanisms for biological functions.[1] Since molecular vibrational states are highly sensitive to intermolecular interaction, the time-resolved spectrum provides atomically detailed information on biomolecules. Advanced spectroscopic methods, such as resonance Raman,[2,3] infrared (IR) difference,[4] and two-dimensional IR,[5] have further made it feasible to measure the spectrum arising from local structures of biomolecules. These techniques have been used, for example, to reveal conformational changes,[4,6] the structure of recognition sites,[7,8] and intermediate species in enzymatic reactions.[9] However, extensive efforts such as mutation and isotope studies are needed to assign the observed vibrational bands and relate the spectrum with biomolecular structures. Therefore, theoretical calculations that predict the vibrational spectrum with reliable accuracy are of great importance.

Anharmonic vibrational calculations have significantly advanced in recent decades.[10−12] Because the conventional harmonic approximation has limited accuracy, it is necessary to incorporate the anharmonicity of the potential energy surface (PES) for quantitative purposes. The second-order vibrational perturbation theory (VPT2)[12−14] based on a harmonic solution is widely used in conjunction with the quartic force field (QFF). Vibrational self-consistent field (VSCF)[15,16] and post-VSCF methods,[16−23] together with PES generation techniques using electronic structure calculations,[24−30] have been developed as a solver of the vibrational Schrödinger equation for general polyatomic molecules. Yagi et al.[31,32] have extended VSCF based on the quasi-degenerate perturbation theory (VQDPT), and developed a suite of programs, SINDO,[33] to generate the PES and solve the VSCF equations. Hanson-Heine et al.[34] have extended a partial Hessian approach[35,36] for VPT2 calculations for large molecular systems. Bowman and co-workers[37,38] have proposed to use normal modes localized to a water molecule, and performed vibrational configuration interaction (VCI) calculations for water clusters and ice. Methods to transform normal coordinates to local coordinates are also actively studied.[39−47] These methods are applicable to molecules with a single minimum (or a few minima) in the gas phase.

Vibrational calculations for complex systems remain a challenge, because (1) the system visits many minima with different frequencies in finite temperature, and (2) the solvent/protein environment needs to be considered. We have recently developed a weight-averaged method,[48−50] which incorporates the former effect using classical molecular dynamics (MD) simulations and vibrational calculations based on the ab initio PES. In this method, MD simulations are first carried out to sample a large number of structures in the target system. Then, the trajectory is analyzed to extract relevant local structures and their statistical weights, and the vibrational calculations are carried out for each local structure to obtain its spectrum. Finally, the spectra are averaged using the weights of local structures to yield a total spectrum. The method has been used to calculate the IR and Raman spectra of a sphingomyelin bilayer[48,51] and the IR spectrum of Nylon6 with different humidity.[49]

In the present study, we combine anharmonic vibrational calculations with a quantum mechanical/molecular mechanical (QM/MM) method to incorporate the effect of solvent/protein environment on the vibrational spectrum. The QM/MM method, first proposed in seminal papers by Warshel, Levitt, and Karplus,[52,53] is a multiscale model, which treats a local region of chemical interest by the electronic structure theory (QM) and the surrounding environment by molecular mechanics force field (MM).[54−63] Although QM/MM has been extensively used to study enzymatic reactions, it has also been employed to compute the vibrational spectrum of biomolecules. Cui and Karplus[64] have derived an analytical Hessian for QM/MM and carried out normal-mode analyses. Tavan and co-workers[65−68] have combined instantaneous normal-mode analyses (INMA) with QM/MM. However, these methods neglect the effect of anharmonicity. Anharmonic vibrational calculations combined with QM/MM are still scarce; the only report so far, to the best of our knowledge, is the work by Ringer and MacKerell,[69] in which VPT2 calculations have been carried out for nitrile probes in solution using QM/MM calculations. Here, we report VQDPT2 calculations based on a PES derived from QM/MM calculations.

The QM/MM method is implemented in many MD programs. The implementation in CHARMM, which was originally combined with semiempirical electronic structure methods,[55] has been extended to interface with external ab initio programs (CADPACK,[70] Q-Chem,[71] and TURBOMOLE[72]) and DFTB.[73−75] AMBER also includes the QM/MM method.[54,76,77] Recently, a portable interface[78] has been developed and combined with FIREBALL[79] and TeraChem.[80] Other QM/MM-MD programs are also found in the literature: NAMD,[81,82] Dynamo,[83,84] etc. In addition, QM/MM calculations are available in electronic structure programs (Gaussian,[85] Q-Chem,[86] deMon2k,[87] Platypus,[88] etc.) and by interface programs (ChemShell[89] and PUPIL[90]).

Following these works, we have implemented the QM/MM method into GENESIS[91,92] (Generalized-Ensemble Simulation System). A QM/MM module in GENESIS interfaces QM programs (Gaussian, Q-Chem, TeraChem, and DFTB+[93]) with MM calculations based on the CHARMM force field.[94] An efficient method to optimize the geometry, which is necessary prior to vibrational calculations, is also implemented on the basis of a limited memory version of Broyden–Fletcher–Goldfarb–Shanno (L-BFGS-B) developed by Zhu et al.[95−97] and a macro-/microiteration scheme.[98] A partial Hessian of a local region of the system is computed by numerical differentiations to obtain normal coordinates. Then, the anharmonic PES is generated in terms of the normal coordinates by combining GENESIS with SINDO, which is then used for solving the VSCF equations to obtain the spectrum.

We first describe an overview of the QM/MM method in Section , and the vibrational analysis based on the QM/MM potential in Section . Pilot applications are carried out for a phosphate ion in solution and P450 nitric oxide reductase (P450nor). The computational details and the results are given in Sections and 5, respectively. Concluding remarks and outlooks are given in Section .

Overview of the QM/MM Method

QM/MM Potential

The QM/MM method divides a system into a QM region of particular interest and an MM region that provides its surrounding environment. The nuclei and electrons are considered explicitly in the QM region, while the MM region is described in terms of atoms defined in a molecular force field. The electrons are treated quantum mechanically by solving the following Schrödinger equation in atomic unitswhere i, a, and m are indices for electrons, nuclei, and MM atoms, respectively, and Z and q are the charge of nuclei and MM atoms, respectively. r denotes a distance between X and Y. Note that |Ψ⟩ and E are the functions of coordinates of nuclei () and MM atoms ()Then, the potential energy of the system is written aswhere V is the QM energy defined asand VLJQM/MM and VMM are the Lennard-Jones interaction between nuclei and MM atoms and the classical force field for MM atoms, respectively.The derivative of the potential with respect to and gives the forces acting on nuclei and MM atomsThe first term in the right-hand side of eq represents the electrostatic interaction between the MM charge and the electron density of the QM region, which can be rewritten in terms of the electric field, (77)In eqs , 9, and 10, VQM, ∇QM, and are calculated using an external QM program, whereas other terms are calculated by GENESIS.

Link Atom Method

Many QM/MM applications require setting QM and MM boundaries cutting covalent bonds of a molecule, in particular, when dealing with large biomolecules. For this purpose, we have implemented the link atom method.[54,55,57] In this method, (1) the QM atom (nuclei) at the boundary (Q1) is terminated by a hydrogen atom, and (2) MM charges in the vicinity of the boundary are excluded in QM calculations to avoid overpolarization of the electron density (see Figure a). There are two options to exclude the charges: one is only the charge of the MM atom at the boundary (M1), and the other is that of a group of MM atoms that M1 belongs to. Note that the group is defined in the CHARMM force field. The MM potential includes bonded terms between QM and MM atoms at the boundary as shown in Figure b. We follow the implementation in CHARMM, where all possible combinations of QM and MM atoms are included for bond (Vb) and angle (Va) terms, whereas the interaction between M3 – Q1 and M2 – Q2 is included, but that between M1 – Q3 is not for 4-site terms (V14) such as dihedral, improper, and (ϕ, ψ) grid-based energy correction map (CMAP) terms.

Figure 1

(a) Schematic illustration of QM and MM boundary. A hydrogen atom is added to the QM atom at the boundary (Q1), and the MM charges are excluded in the vicinity of the link hydrogen atom (open circles). (b) Bonded terms between QM and MM atoms included in the MM potential.

In the presence of link hydrogen atoms, the QM energy becomes a function of their coordinates, To remove the dependency on the position of link atoms, we introduce constraint to the position of the hydrogen atom.[57,76,99] Let us place the hydrogen atom along a bond between Q1 and M1 with a distance d from Q1where RQ = (RQ, RQ, RQ) and RM = (RM, RM, RM) are the position of Q1 and M1, respectively, and r1 is a distance between Q1 and M1. Then, the gradient of the QM energy is obtained by a chain rulewhere X is either Q or M. From eq , the gradient at Q1 and M1 is obtained aswhere δ is the Kronecker delta.

Further details on the implementation of QM/MM in GENESIS are given in Section 1 of the Supporting Information (SI).

Vibrational Analysis Based on the QM/MM Potential

Geometry Optimization and Normal Coordinates of a Partial System

Let us assume that a local region of the system that yields the spectrum in question is known (denoted VIB region): for example, reactant/product species in enzymatic reactions, hydrogen bonds in a protein, and so on. Then, the QM region is taken to be the VIB region as well as some surrounding molecules to account for the nearest neighbor interaction that affects the vibrational states of the VIB region.

The geometry is first optimized for the QM region and a part of the MM region. An efficient optimizer is implemented using an L-BFGS-B code developed by Zhu et al.[95−97] in terms of Cartesian coordinates. Furthermore, the macro-/microiteration scheme is implemented following Kästner et al.[98] This scheme optimizes the geometry in two loops of cycle. The outer loop (macroiteration) is the usual cycle for all atoms, whereas the inner loop (microiteration) optimizes only the MM atoms keeping the QM atoms fixed. In the microiteration, the electrostatic interaction between QM and MM atoms is approximated using charges of QM atoms, q, asNote that, in the presence of link atoms, the sum over MM atoms excludes MM atoms in the vicinity of the link hydrogen in which the charge is removed. The charges of QM atoms are derived, for example, from the electrostatic potential (ESP) obtained from QM calculations in the macroiteration.[100,101] The scheme drastically reduces the computational cost, because the MM atoms are optimized in microiteration without requiring QM calculations. However, the approximation in eq becomes less accurate as MM atoms move more. For a compensation for such deviations, the gradient of MM atoms is corrected aswhere 0 denotes coordinates of the system in macroiteration. The correction recovers the gradient of the QM energy at = 0 and makes the gradient in macro- and microiterations consistent with each other.

Once the geometry optimization is converged, a partial Hessian matrix[35,36] is computed for VIB atoms in 3N dimension, where N is the number of VIB atoms. The Hessian is obtained from numerical differentiations aswhere x denotes 3N Cartesian components of VIB atoms, δ is a step-size of numerical differentiations, and g(±δ) is the ath component of the gradient at (x ± δ). QM calculations are required at 6N points to construct the Hessian matrix, which is processed in MPI parallel. Then, the Hessian is mass-weighted and diagonalized to obtain normal coordinates (Q) and harmonic frequencies (ω)where m and xeq are the atomic mass and coordinates at the equilibrium structure, respectively, and is a diagonal matrix with M = δm.

Anharmonic Vibrational Calculations

The vibrational Hamiltonian is written in terms of the normal coordinates aswhere f is the number of modes, and V( is the so-called n-mode representation of the potential energy surface (nMR-PES),[102] in which the PES is expanded in terms of mode coupling and truncated at the nth orderNote that the PES is a function of coordinates of all atoms, but that the coordinates of atoms out of the VIB region () are fixed to the equilibrium geometryIn other words, the present method neglects the motion of atoms out of the VIB region.

The nMR-PES can be generated using the techniques matured in the field of vibrational calculations.[24−26,29,30,103−105] Multiresolution PES[29] is a composite of grid potentials[27] and QFF.[28] QFF is written in the form of nMR asThe third- and fourth-order coefficients are generated by numerical differentiations of the gradient (see Section 2 in Supporting Information). The number of grid points to compute 3MR- and 4MR-QFF accumulated over all coupling order is, respectively, given asThe grid potential is based on the harmonic oscillator discrete variable representation (HO–DVR). The number of energy points needed for generating nMR grid potential iswhere M is the number of grid points along each coordinate. Note that the grid potentials of the low-order coupling (n′ < n) are generated simultaneously by setting one of the grid point to be the origin, i.e., by taking M to be an odd number. Although the number of grid points scales steeply with respect to n and f, it may be reduced by selecting the mode of interest such that f < 3N and by varying the size of grid according to the strength of coupling terms.

These grid points are generated by SINDO, and the QM/MM energy and gradients at the points are calculated by GENESIS in conjunction with a QM program. The loop over grid points is parallelized by MPI. The information is used to generate the multiresolution nMR-PES; in the present study, we combine 3MR-QFF and 1MR-grid PES. The PES is subsequently used for solving the vibrational Schrödinger equation by VSCF and post-VSCF methods. Details on VSCF and VQDPT2 methods are described in Section 3 of the Supporting Information.

Computational Details

Phosphate Ion in Solution

H2PO4– was prepared in a box of water (48.6 × 48.6 × 48.6 Å3) by the solvation utility in VMD.[106] The total number of atoms was 11 747. Then, the all-atom MD simulation was performed using spdyn of GENESIS 1.2.2. The force field parameters for H2PO4– and water were, respectively, those in the previous work[107] and TIP3P.[108] Following a 1000-step minimization and an equilibration NPT run for 1.0 ns, a production NVT run was carried out for 10.0 ns. The leapfrog algorithm was used for integration with a time step of 2.0 fs. SHAKE[109] and SETTLE[110] were used to constrain the bond length involving hydrogen atoms. A Langevin thermostat and barostat[111] were used to control the target temperature and pressure at 300 K and 1 atm, respectively. The long-range electrostatic interaction was calculated using the smooth particle mesh Ewald (PME) method.[112,113] The non-bonded interaction was reduced to zero between 10 and 12 Å employing a switching function. The neighbor list was updated every 20 fs. The trajectory was saved every 10 ps, and the resulting 1000 snapshots were used for the analysis.

The snapshots were classified in terms of the number of water molecules near H2PO4– (Nwat). The water molecule was judged to be “near” when the oxygen atom of a water molecule was within 3.5 Å in distance of any of the oxygen atoms of H2PO4–. The number of snapshot structures (Nss) for each Nwat is shown in Table indicating a distribution of Nwat in the range 9–16. Then, one representative snapshot was selected for each class of Nwat by k-means clustering method with respect to a root-mean square deviation (RMSD) of all atoms of H2PO4–. The snapshot was employed as an initial structure of the QM/MM optimization. QM/MM calculations were performed combining a development version of GENESIS/atdyn with Gaussian09.[85] The QM region was set to H2PO4– and the neighboring water molecules, so that the number of QM atoms was 7 + 3Nwat. DFT calculations were performed at the level of B3LYP/cc-pVTZ.[114−116]

Table 1

Number of MD Snapshots Classified in Terms of the Number of Water Molecules around H2PO4–

Nwata	9	10	11	12	13	14	15	16
Nssb	21	70	194	268	264	135	39	9

The number of water molecules, where the oxygen atom is within 3.5 Å of any of the oxygen atoms of H2PO4–.

The number of MD snapshots classified by Nwat.

The geometry was optimized relaxing the atoms within 10 Å of H2PO4– and fixing the others. The optimization was performed using the macro-/microiteration scheme with the Merz–Singh–Kollman charge.[100,101] The convergence threshold was set to 0.36 and 0.54 kcal mol–1 Å–1 for the RMS and maximum gradient, respectively. The optimization converged after 110 cycles on average and 183 cycles in maximum.

The vibrational analysis was performed for the optimized geometry taking H2PO4– as a VIB region. For the calculation of a partial Hessian, 42 points of gradients were required. After the normal-mode analysis, 3MR-QFF and 1MR-grid PES (with 11 grid points to each coordinates) were constructed in 21 normal coordinates of H2PO4–. The generation of 3MR-QFF and 1MR-grid PES required 925 points of gradients and 211 points of energies, respectively. Then, VSCF and VQDPT2 were carried out using the harmonic oscillator basis functions up to n = 10. The quasidegenerate space in VQDPT2 was selected using Ngen = 1 and k = 4. The PES generation and VQDPT2 calculations were carried out using SINDO. The dipole moment surface was constructed from 1MR-grid points and used for calculating the IR intensity. The IR spectrum was constructed by a convolution of vibrational bands constructed by Lorentz functions of 40 cm–1 width (see Section 4 of the Supporting Information). Finally, the IR spectra for each Nwat were weight-averaged using Nss/Ntot as a weight, where Ntot is the total number of snapshots (=1000).

P450 Nitric Oxide Reductase (P450nor)

The structure was taken from a recent X-ray crystal structure[117] (PDBID 5Y5H), which contains P450nor with a NO molecule bound to a ferric iron (FeIII) and crystal water molecules. First, the missing residues of the N terminus were added using GaussView5.0.[118] The hydrogen atoms and protons were added to the system using UCSF Chimera.[119] His96 and Glu361 were protonated based on pKa values at pH 6.5 estimated by PROPKA3.1.[120,121] The CHARMM36 (c36)[94] force field was used for the protein and ions, and TIP3P for water molecules. The parameters for a heme–NO unit were set as follows: (1) By considering the unit as a ferrous heme–NO+, the heme was described by the c36 force field for a ferrous heme, and the atomic charges of N and O atoms were set to +1.272e and −0.272e, respectively. (2) Other parameters for the NO (the N–O bond function and the vdW parameters) were taken from the previous study on myoglobin.[122,123] (3) The parameters for an Fe–NO bond function were set to req = 1.670 Å, = 350.0 kcal mol–1 Å–2; those for Npor–Fe–N (Npor is a nitrogen atom of porphyrin) and Fe–N–O angle functions were set to θeq = 90.0°, k = 100.0 kcal mol–1 rad–2, and θeq = 158.0°, k = 8.0 kcal mol–1 rad–2, respectively, and those for a X–Fe–N–Y dihedral function (X and Y denote an arbitrary atom type) were set to δ = 0.0°, n = 4, k = 0.05 kcal mol–1.

The MD simulation was first carried out fixing heavy atoms of P450nor and crystal water molecules using NAMD 2.12.[81] The energy was minimized for 5000 steps in vacuum to relax the position of hydrogen atoms. The resulting structure is referred to as XtalV. Next, a simulation box of 120 × 113 × 101 Å3 in size was constructed adding 38 284 water molecules, 116 Na+, and 110 Cl–. Then, three simulations were carried out each for 3 ns in NVT, NPT, and NVT conditions. A Langevin thermostat and the Nosé–Hoover Langevin piston methods[124,125] were used to control the temperature at 300 K and the pressure at 1 atm, respectively. The size of the box was relaxed to 115.2 × 108.5 × 97.0 Å3 after the NPT simulation. The last snapshot of these simulations is referred to as XtalW.

Then, the crystal structure (P450nor and water molecules) was gradually relaxed in a series of NVT simulations using harmonic positional restraints. The force constants, initially set to 5.0 and 2.0 kcal mol–1 Å–2 for the backbone and the others, respectively, were reduced to zero in 9 steps for 27.0 ns (see Table S1 in the Supporting Information). Then, the system was further equilibrated for 10.0 ns without restraint, followed by a production run for 70.0 ns. The temperature was controlled at 300 K with Bussi thermostat.[126] The equilibration runs with restraint were carried out with a time step of 2.0 fs using the velocity Verlet integrator, whereas the 10 ns equilibration and the 70 ns production runs were performed using r-RESPA integrator[127] with a time step of 2.5 fs. The long-range electrostatic interaction was calculated using the smooth PME method.[112,113] The non-bonded interaction was reduced to zero between 10 to 12 Å using a switching function. All bonds involving hydrogen atoms were kept rigid using SHAKE[109] and SETTLE.[110] These simulations were performed using spdyn of GENESIS 1.1.6. The snapshot structures of the production run were classified into five clusters using the k-means clustering method by fitting the structure to nitrogen atoms of the heme and measuring the RMSD of Fe, N, and O atoms. The representative structures, i.e., the center of each cluster, were found to be the snapshot at 27.27, 33.59, 35.70, 45.20, and 65.53 ns, which are referred to as NVT1, NVT2, NVT3, NVT4, and NVT5, respectively.

QM/MM calculations were performed starting from the seven snapshots, XtalV, XtalW, and NVT1–5, using a development version of GENESIS/atdyn. A QM subsystem consisted of a heme–NO unit and a side chain of Cys352 (79 atoms). The QM and MM boundary was set between Cα and Cβ atoms of Cys352. The distance between a QM boundary atom and a link H atom was set to 1.0 Å. The charge of group atoms near the MM boundary was excluded in QM calculations. DFT calculations were performed at the level of BP86 or B3LYP-D3 (B3LYP hybrid functional and the D3 version of Grimme’s dispersion with Becke–Johnson damping[128]) using Gaussian09, and ωB97M-V[129] using Q-Chem 4.4. The def2-SVPD[130] basis functions were used for seven atoms (four N atoms of porphyrin, a S atom of Cys352, and NO), and the def2-SVP[131] basis sets for other atoms [denoted def2-SVP(D)]. The ωB97M-V calculation was also carried out using larger basis sets, def2-TZVPPD[130] for Fe, the atoms bound to Fe, and NO, and def2-TZVP[131] for the others [denoted def2-TZVP(PD)]. The spin multiplicity was set to singlet. In addition to a closed-shell singlet (cs-1et) state, an open-shell singlet (os-1et) state was also calculated by a spin-unrestricted Kohn–Sham (KS) method. The initial guess orbitals for the os-1et state were prepared by combining KS orbitals obtained for the heme, the NO ligand, and the side chain of Cys352 by specifying the charge/spin multiplicity as −1e/doublet, 0e/doublet, and −1e/singlet, respectively. An MM subsystem consisted of the remaining part of the protein and water molecules excluding the ions.

The geometry was optimized for the whole protein and water molecules within 20 Å of NO keeping the other atoms fixed. The optimization was performed using the macro-/microiteration scheme based on the Merz–Singh–Kollman ESP charges[100,101] by setting the atomic radii for Fe and H to 1.3 and 1.0 Å, respectively. A convergence threshold was set to 0.30 and 0.45 kcal mol–1 Å–1 for the RMS and maximum gradients, respectively. The optimization converged after 154 cycles on average and 219 cycles in maximum.

The vibrational analysis was performed setting four atoms as a VIB region; Fe, N, and O atoms of NO, and S of Cys352. For the computation of a partial Hessian, 25 points of gradients were required. After the normal-mode analysis, 3MR-QFF and 1MR-grid PES (with 11 grid points to each coordinates) were constructed for 9 normal modes excluding 3 modes related to the motion of the S atom. The generation of 3MR-QFF and 1MR-grid PES required 181 points of gradients and 91 points of energies, respectively. Then, VSCF and VQDPT2 calculations were carried out using the harmonic oscillator basis functions up to n = 10. The quasidegenerate space in VQDPT2 was selected using Ngen = 1 and k = 4. The PES generation and the VQDPT2 calculations were carried out using SINDO.

QM cluster calculations without the MM environment were also carried out for reference. The geometry of the QM subsystem was taken from the X-ray crystal structure. The position of Cβ atom of Cys352 and a dihedral angle, Cβ–S–Fe–NA (NA is one of the nitrogen atoms of a heme group), were frozen during the geometry optimization to avoid rotation of a heme plane around an Fe–S axis. The geometry optimization and the vibrational calculations were performed in the same way as above.

Results

Phosphate Ion in Solution

The previous studies[66,132] have reported a strong solvation effect in the vibrational spectrum of H2PO4– in solution. Comparison of the spectrum in the gas phase and solution has shown that P–O stretching modes (sPO and aPO in Figure ) are red-shifted by 50–170 cm–1, whereas P–O–H bending and P–OH stretching modes (δH and sPO(H)/aPO(H), respectively, in Figure ) are blue-shifted by 120–150 cm–1. The large spectral shift arises from strong hydrogen bonds (HBs) formed between the ion and water molecules, which require an accurate treatment by QM methods. Thus, H2PO4– serves as a test case to assess the accuracy of the proposed method.

Figure 2

Vibrational modes of phosphate ion (H2PO4–). δOPO, OPO bending modes; τOPOH, OPOH torsion; sPO(H) and aPO(H), symmetric and asymmetric P–OH stretching modes of hydroxyl groups; sPO and aPO, symmetric and asymmetric P–O stretching modes; δH, POH bending mode.

From the MD trajectory, the number of water molecules around H2PO4– (Nwat) was found to be distributed in the range 9–16 with Nwat = 12 and 13 being the most abundant, as shown in Table . QM/MM optimized structures for Nwat = 12 and 13, shown in Figure , indicate that the PO– bond accepts three water molecules whereas the hydroxyl group accepts one and donates to another water molecule. Thus, the current QM/MM calculation with Nwat = 9–16 takes into account the HB interaction between H2PO4– and water molecules in the first solvation shell and some in the second solvation shell at the DFT level.

Figure 3

Optimized structures for (a) Nwat = 12 and (b) Nwat = 13. H2PO4– and highlighted water molecules are QM atoms, whereas the surrounding other water molecules are MM atoms. Hydrogen bonds are indicated by blue broken lines, which are drawn when rOO < 3.2 Å and θOHO > 130°.

IR spectra obtained by the harmonic approximation for a QM/MM system with Nwat = 12 (QM/MM-12) and H2PO4– in the gas phase are shown in Figure e and 4f, respectively. The solvation of H2PO4– causes a drastic change in the spectrum as reported in the previous studies.[66,132] The prominent peak shifts, i.e., the red-shift of aPO and sPO and the blue-shift of δH, aPO(H), and sPO(H), are clearly seen in the present result as well. Furthermore, peak positions of τOPOH are found to be dramatically shifted by more than 500 cm–1. This is because the formation of HBs restricts the motion of hydroxyl groups and makes the torsional motion rigid and stiff. Such behavior is similar to bulk water, where the rotational motion of water molecules is restricted by HBs giving rise to libration bands around 400–700 cm–1.

Figure 4

(a) Experimental IR spectrum of H2PO4– in solution.[66] (b) IR spectrum obtained by Born–Oppenheimer (BO) MD based on the B3LYP-D/TZV2P level in the previous work.[132] IR spectra calculated by QM/MM in this work: (c) weight-average of VQDPT2 spectra over Nwat = 9–16, (d) VQDPT2 for Nwat = 12, and (e) harmonic calculation for Nwat = 12. (f) Harmonic spectrum of H2PO4– isolated in the gas phase. The spectrum is constructed using Lorentz functions with a width of 40 cm–1.

To further investigate the strong solvation effect, we have calculated the vibrational frequency of H2PO4– using the polarizable continuum model (PCM)[133] and a QM/MM model which takes only H2PO4– as a QM region (QM/MM-0). The results are listed in Table . Although PCM is capable of predicting the direction of the shift of P–O and P–OH stretching modes, the amount of shift is insufficient. Furthermore, δH and τOPOH are shifted to an opposite direction. The mean absolute deviation (MAD) is 200.5 cm–1 compared to QM/MM-12. In QM/MM-0, the result is significantly improved predicting the shift of δH and τOPOH in the right direction and reducing the MAD to 21.3 cm–1. Nonetheless, QM/MM-0 is semiquantitative yielding a maximum deviation of 52.9 cm–1 for sPO. Therefore, it is essential to take not only H2PO4– but also the water molecules as a QM region to account for the effect of HBs on the vibrational states of H2PO4–.

Table 2

Fundamental Frequencies of H2PO4– Obtained by the Harmonic Approximation and VQDPT2 for an Isolated System (Gas), PCM with Water as Solvent (PCM), a QM/MM System with and without 12 Water Molecules in the QM Region (QM/MM-12 and QM/MM-0, Respectively), and a Weight-Average over 9–16 Water Molecules (QM/MM-wa), together with the Experimental Results

	gas	PCM	QM/MM-0	QM/MM-12		QMMM-wa
mode	harmonic	harmonic	harmonic	harmonic	VQDPT2	VQDPT2	exptla
δOPO	429.8	441.5	505.9	505.6	502.5	517	521
	492.2	481.8	540.0	519.0	513.7	517	521
	497.6	492.8	552.7	529.8	527.3	517	521
τOPOH	203.5	190.1	774.6	810.5	768.0	713
	313.3	275.1	864.7	858.5	829.3	713
sPO(H)	750.1	770.3	841.5	854.1	843.2	848	879
aPO(H)	778.2	789.7	894.9	914.0	899.0	922	944
sPO	1093.9	1090.4	1119.4	1066.5	1048.5	1051	1077
aPO	1326.6	1275.4	1206.6	1176.5	1149.9	1145	1156
δH	1056.4	1020.3	1250.1	1268.7	1225.8	1211	1213
	1071.0	1040.3	1309.1	1324.0	1285.9	1283

ref (66).

Comparison of VQDPT2 and harmonic spectra in Figure d and 4e shows that the anharmonicity makes the frequency of all modes lower, although the overall shape of the spectrum looks similar. Nevertheless, it is notable that the size of anharmonic correction differs by the character of modes. δH and τOPOH make a larger shift than others, whereas δOPO hardly shows any shift. Thus, the scaling of harmonic frequencies by a single factor is not appropriate. It is also interesting to note that Fermi splitting is found for sPO, which is due to the resonance between a fundamental of sPO and an overtone of δOPO. The splitting makes the bandwidth slightly broader in VQDPT2 than in the harmonic spectrum.

The weight-average over Nwat = 9–16 (QM/MM-wa) smears out a peaky structure in a spectrum of QM/MM-12, as shown in Figure c and 4d. It is notable that the frequency of δH and τOPOH is distributed in an extremely wide range, which is a consequence of the motion of hydrogen atoms being sensitive to the hydration structure of H2PO4–. Note that the broadening makes it difficult to identify the peak position of δH and τOPOH in the total spectrum. It is also interesting to see that the relative intensity of sPO and aPO is reversed in Figure c and 4d. Although aPO has larger intensity than sPO for each cluster, the frequency of aPO fluctuates more than that of sPO. Thus, the weight-average yields the band of aPO broader in frequency and smaller in intensity than that of sPO.

Figure a and 4b shows the experimental spectrum[66] and the theoretical one obtained in the previous study,[132] respectively. The previous study carried out Born–Oppenheimer MD (BOMD) at the B3LYP-D level with TZV2P basis sets and obtained the spectrum by a Fourier transform of ensemble-averaged autocorrelation functions of the dipole moment. The two methods reproduce the experimental spectrum with similar quality, which is encouraging given that the computational techniques are totally different. One of the conclusions of the BOMD study was that the use of a nonpolarizable MM force field for water molecules (i.e., TIP3P) is inadequate for predicting the IR spectrum of H2PO4– in solution. The present result is in line with this finding. We observe that treating all water molecules at the MM level incurs a large error in the calculated frequency (see QM/MM-0 in Table ) and that the inclusion of water molecules in the QM region is important for quantitative purposes. Nonetheless, some improvements are worth mentioning. The previous work shows a sharp dip between aPO and δH, and a prominent peak of δH, whereas the present result shows a smooth shoulder in the range 1200–1300 cm–1 that resembles the experimental feature. Also, τOPOH was not identified in the previous work. These shortcomings may be due to a limited sampling time (56 ps) in the BOMD simulation. In the present study, we performed a sampling of the hydration structure for 10 ns and a clustering in terms of Nwat. This procedure, albeit at the MM level, yields essential HB structures between H2PO4– and water molecules that contributes to the spectrum, and leads to a smooth convergence of the spectrum with respect to the number of snapshot as small as 8 (one for each Nwat = 9–16). The result suggests that the weight-average approach[50] combining a structural sampling using MD and vibrational calculations of IR spectra using VQDPT2 and QM/MM is a promising approach to simulate the IR spectrum of solvated molecules.

P450nor is a heme enzyme isolated from the fungus Fusarium oxysporum that reduces a nitric oxide (NO) to a nitrous oxide (N2O) with electrons transferred from a nicotinamide adenine dinucleotide (NADH). The mechanism of the NO reduction has been extensively studied by spectroscopy[117,134,135] and theoretical calculations.[136,137] The initial step of the reaction is a binding of NO to a ferric iron (FeIII) of the heme group, which is characterized by N–O and Fe–NO stretching modes.[138−141] Here, we study these vibrations by VQDPT2 based on QM/MM.

Figure a plots the RMSD from the crystal structure and the radius of gyration (Rg) obtained from the MD trajectory of the production run. Both RMSD and Rg are found to be stable, suggesting that the structure of P450nor is well-maintained during the MD simulation. The RMSD is obtained as 1.74 ± 0.19 Å, and the Rg, which is 21.06 Å at the crystal structure, is slightly increased to 21.38 ± 0.12 Å.

Figure 5

(a) RMSD from the X-ray crystal structure (top) and the radius of gyration (Rg) (bottom) as a function of simulation time in the production run. RMSD is calculated for Cα atoms of residues 4–403. Rg is obtained using backbone heavy atoms of residues 4–403 with mass weighting. (b) QM/MM optimized geometry of P450nor for NVT4 (colored) superimposed with the X-ray crystal structure (gray). The two structures are fit by Cα atoms of residues 4–403. Water molecules within 20 Å of NO (blue dots) are relaxed in QM/MM calculations. (c) Top view of part b without water molecules. (d) Close view of part b around the heme group with NO.

The results of QM/MM calculations are listed in Table . The RMSD for XtalV and XtalW is as small as 0.26 and 0.22 Å, respectively, confirming that the optimized structures are close to the crystal structure. On the other hand, RMSD for NVT is obtained in the range 1.24–1.87 Å, indicating that the optimized structure for NVT corresponds to the one relaxed in solution. The structure of NVT4 is compared with the crystal structure in Figure b,c. Although NVT4 is slightly expanded in accordance with the increase in Rg, and loop regions show some deviations, the overall structure of NVT4 matches well with the crystal structure without undergoing conformational changes or unfolding. Nevertheless, it is notable that the heme group is significantly shifted as shown in Figure d. The distance between Fe atoms of the crystal and NVT4 structures is 1.52 Å. A similar shift is observed for all the NVT snapshots, and thus, the shift is a general trend of P450nor in solution. Figure d also shows that the orientation of NO around the Fe–S axis is different between the crystal and NVT4 structures. In fact, NO is oriented in many different ways in NVT1–5, as evident from the dihedral angles, ϕNA–Fe–N–O and ϕCβ–S–N–O, listed in Table . The result implies that NO rotates around the Fe–S axis in the absence of NADH.

Table 3

RMSD from the Crystal Structure and Geometric Parameters in the Active Site of P450nor [Bond Lengths (r in Å), Bond Angles (θ in deg), and Dihedral Angles (ϕ in deg)], and the N–O and Fe–NO Stretching Vibrations [Harmonic (ω) and VQDPT2 (ν) Frequencies Both in cm–1] Obtained from QM Cluster and QM/MM Calculations Together with Experimental Values

	cs-1eta								os-1eta
	Gasc	XtalV	XtalW	NVT1	NVT2	NVT3	NVT4	NVT5	NVT4	exptlb
Cα RMSDd		0.26	0.22	1.24	1.72	1.87	1.72	1.47	1.72
rFe–NO	1.672	1.680	1.672	1.659	1.665	1.654	1.653	1.661	1.697	1.67
rN–O	1.149	1.146	1.143	1.139	1.140	1.138	1.137	1.137	1.144	1.15
rFe–S	2.282	2.287	2.301	2.305	2.296	2.319	2.310	2.316	2.297	2.33
θFe–N–O	158.6	155.7	158.3	164.7	161.9	167.6	165.0	164.6	163.3	157.9
ϕNA–Fe–N–Oe	48.2	46.3	52.9	46.4	99.5	120.3	130.8	122.6	129.7	60.3
ϕCβ–S–N–O	57.3	55.7	63.3	39.3	93.0	115.0	125.3	123.3	123.8	71.4
ωN–O	1936.5	1946.2	1967.6	1998.1	1988.9	2008.9	2017.4	2011.6	1974.4
νN–O	1894.7	1913.4	1932.1	1966.2	1956.8	1975.4	1984.7	1981.5	1929.3	1853
ΔN–Of	(41.8)	(32.8)	(35.5)	(31.9)	(32.1)	(33.5)	(32.7)	(30.1)	(45.1)
ωFe–NO	603.5	605.1	607.2	633.7	614.9	620.9	628.4	617.8	417.3
νFe–NO	595.3	587.7	585.6	617.4	596.2	600.5	610.9	600.0	327.2	530
ΔFe–NOf	(8.2)	(17.4)	(21.6)	(16.3)	(18.7)	(20.4)	(17.5)	(17.8)	(90.1)

Closed-shell singlet (cs-1et) and open-shell singlet (os-1et) states.

Structural parameters are from X-ray crystal structure[117] and frequencies from resonance Raman spectra in solution.[135]

A QM cluster calculation in the gas phase without the MM region.

The RMSD from the X-ray crystal structure calculated for Cα atoms of residues 4–403.

NA is one of the nitrogen atoms of a heme.

The difference between harmonic and VQDPT2 frequencies (Δ = ω – ν).

The harmonic (ω) and VQDPT2 (ν) frequencies are also listed in Table . The shift in the frequency due to anharmonicity on average over the seven snapshots in the cs-1et state is found to be 32.7 ± 1.5 and 18.5 ± 1.7 cm–1 for the N–O and Fe–NO stretching modes, respectively. The ratio between VQDPT2 and harmonic frequencies for the N–O stretching mode, νN–O/ωN–O, is obtained as 0.984, which is larger than 0.961 used in the previous study to correct the harmonic frequency.[137]

Interestingly, νN–O increases as 1913.4, 1932.1, and 1956.8–1984.7 cm–1 for XtalV, XtalW, and NVT1–5, respectively, as the solvation is incorporated in the model. Such a drastic blue shift may be surprising given that NO is a hydrophobic molecule that does not directly interact with water molecules. The result has inspired us to carry out a QM cluster calculation without the MM environment (denoted Gas in the following), which yield 1894.7 cm–1 for νN–O and a similar geometry to XtalV. Figure shows plots of a bond angle between Fe–N–O (νFe–N–O) and an Fe–S bond length (rFe–S) as a function of νN–O for Gas, XtalV, XtalW, and NVT1–5. It is notable that the increase in νN–O correlates with an expansion of θFe–N–O and an elongation of rFe–S. The behavior is consistent with the so-called thiolate trans effect,[138] where the ligation of thiolate to FeIII(NO) induces an occupation of antibonding orbitals of Fe–NO and N–O that makes Fe–NO and N–O bonds weak and an Fe–N–O angle bent (see Section 5 of the Supporting Information for details). These results suggest that the solvation of P450nor elongates rFe–S and weakens the thiolate trans effect, thereby inducing a strengthening of N–O bond and an increase in νN–O.

Figure 6

Fe–N–O bond angle (θFe–N–O) and the Fe–S bond length (rFe–S) as a function of the NO stretching frequency (νN–O) calculated by QM cluster (gas) and QM/MM (XtalV, XtalW, and NVT).

FeIII(NO) has another spin state, os-1et state, where the radical electrons of NO and Fe remain unpaired, but have an antiparallel spin so that the total system is singlet. The os-1et state is known to be close in energy to the cs-1et state.[136,137,139] Thus, we have repeated the same calculation for the os-1et state starting from NVT4. The results are also listed in Table . It is notable that νN–O and νFe–NO are both decreased in the os-1et state compared to those in the cs-1et state. This is because the radical electron of NO occupies a π* orbital of NO and reduces the d → π* back-donation, and hence, the N–O and Fe–NO bonds are both weakened. The decrease in νN–O, obtained as 55.5 cm–1, is consistent with the previous QM/MM calculation by Riplinger et al.[137]

νN–O and νFe–NO have been measured at 1853 and 530 cm–1, respectively, by the resonance Raman spectroscopy for P450nor in solution.[135] However, the calculated νN–O and νFe–NO for NVT1–5 in the cs-1et state, which are 1973 and 605 cm–1 on average, overestimate the experiment by 120 and 75 cm–1, respectively. In view of such a deviation, we have extended the size of a QM region to include residues and water molecules nearby the heme group, but found that it has little effect on the frequencies (data shown in Table S2 of the Supporting Information). In the os-1et state, νN–O is decreased to 1929.2 cm–1, and thus, the deviation is reduced to 76.2 cm–1. However, νFe–NO is also decreased to 327.2 cm–1 and deviates from the experiment by −202.8 cm–1. Note that the ratio of νN–O/νFe–NO is obtained as 3.25 and 5.90 for the cs-1et and os-1et states, respectively. Judging from the corresponding value for the experiment (3.49), it is plausible that the spin state of FeIII(NO) is the cs-1et state in the experimental condition.

On the basis of these findings, we have carried out QM/MM calculations with the same QM region for NVT4 in the cs-1et state using different DFT functionals and basis sets. The results are listed in Table . The pure functional, BP86, gives νN–O as 1881.3 cm–1 in agreement with the experimental value of 1853 cm–1. However, νN–O/νFe–NO is obtained as 3.09. Furthermore, the variation in the geometry (rFe–S and νFe–N–O) from B3LYP-D3 to BP86, which indicates a weakening of the thiolate trans effect, is inconsistent with the decrease in νN–O. Thus, the agreement of νN–O using BP86 may be fortuitous. We also employed ωB97M-V, which is one of the latest density functionals developed by Mardirossian and Head-Gordon.[129] As shown in Table , both νN–O and νFe–NO are increased compared to B3LYP-D3 resulting in a larger deviation from the experiment. Nonetheless, νN–O/νFe–NO is obtained as 3.27 and 3.20 using the def2-SVP(D) and def2-TZVP(PD) basis sets, and the variations in the geometry and frequency are consistent with each other. Considering that cs-1et and os-1et states lie close in energy,[139] multiconfigurational methods, which can describe the two states on an equal footing, may be needed for the electronic structure calculations. Furthermore, it may be necessary to treat multiple PESs and their coupling terms in the vibrational calculation to achieve quantitative accuracy.

Table 4

Representative Geometric Parameters in the Active Site of P450nor [Bond Lengths (r in Å), Bond Angles (θ in deg)], and the N–O and Fe–NO Stretching Vibrations (ν in cm–1) Obtained from QM/MM Calculationsa Using Different DFT Functionals and Basis Sets Together with Experimental Values

	BP86	B3LYP-D3	ωB97M-V
	def2-SVP(D)c	def2-SVP(D)c	def2-SVP(D)c	def2-TZVP(PD)d	exptlb
rFe–NO	1.655	1.653	1.643	1.641	1.67
rN–O	1.160	1.137	1.121	1.115	1.15
rFe–S	2.337	2.310	2.298	2.300	2.33
θFe–N–O	170.0	165.0	167.5	170.0	158
νN–O	1881.3	1984.7	2097.5	2054.6	1853
νFe–NO	608.0	610.9	641.1	642.0	530

The snapshot is taken from NVT4, and the spin state is set to cs-1et.

Structural parameters are from X-ray crystal structure[117] and frequencies from resonance Raman spectra in solution.[135]

def2-SVPD for the atoms bound to Fe and NO, and def2-SVP for the others.

def2-TZVPPD for Fe, the atoms bound to Fe, and NO, and def2-TZVP for the others.

Summary and Outlook

An anharmonic vibrational method based on a PES from hybrid QM/MM computations is developed for vibrational analyses of biomolecules and solvated molecules. For this purpose, we have implemented into GENESIS an interface for QM/MM calculations and methods to facilitate the normal-mode analysis: the geometry optimization (L-BFGS-B minimizer and macro-/microiteration scheme) and a numerical differentiation scheme to derive a Hessian matrix for a local region. Using the resulting vibrational coordinates, multiresolution nMR-PESs are generated from the QM/MM energy and gradient at grid points, and the PES is used for VSCF and VQDPT2 calculations to obtain the anharmonic vibrational spectrum.

The method is first applied to a phosphate ion in solution. The strong solvation effect is incorporated by treating not only the phosphate ion but also neighboring water molecules as a QM region at the B3LYP/cc-pVTZ level. The calculated spectrum demonstrates a characteristic red-shift of the P–O stretching frequency, and a blue-shift of the P–OH stretching and POH bending frequencies. The calculated band positions agree with experiment with a mean absolute deviation of 16 cm–1. The present calculation also reproduces the shape of the spectrum by taking a weight-average of IR spectra over relevant hydration structures. Second, the method is applied to a NO molecule bound to a ferric (FeIII) heme of P450nor. The spin state of FeIII(NO) is suggested to be a cs-1et state based on a ratio νN–O/νFe–NO calculated for cs- and os-1et states. The absolute values of νN–O and νFe–NO deviate from the experiment by ∼100 cm–1, which may be due to imbalance in the QM and MM interaction or deficiency in DFT to describe the electronic structure of FeIII(NO).

Applications of the proposed vibrational analysis method extend to a wide range. Hydrogen-bonded systems are relevant above all: for example, internal water molecules in a protein, proton transfer systems, amide modes of polypeptides and proteins, etc. To treat such systems, not only the explicit treatment of hydrogen-bonded molecules by QM calculations but also the use of advanced force fields that incorporate polarizability and charge transfer is a promising direction. These subjects will be the scope of future works.