Recent Applications of Hybrid Ab Initio Quantum Mechanics – Molecular Mechanics Simulations to Biological Macromolecules

To understand the functional mechanisms of biological macromolecular systems, investigations of both the three-dimensional geometric and electronic structures are important. Accordingly, quantum mechanical (QM) calculations are essential to obtain the properties relevant to the electronic structures of the macromolecular systems. However, because of the high computational costs, only a few hundred atoms can be included in the usual ab initio QM calculations. Therefore, even though biological macromolecular systems include a huge number of atoms, isolated QM model calculations have been performed to date, in which only the reaction centres are considered. However, this may lead to serious misunderstandings of the properties relevant to the electronic structures of the macromolecular systems. To overcome these problems, the ab initio QM calculation is combined with a molecular mechanics (MM) calculation method; this is referred to as a QM/MM calculation. In this strategy, the QM calculation is adopted for the active site (QM regions), and for the remainder of the system, the MM calculation is adopted (MM regions). This QM/MM methodology was originally developed by Warshel and Levitt (Warshel&Levitt, 1976), and since their pioneering work, great progress has been made in the development of QM/MM algorithms and their applications to biological systems (Bruice&Kahn, 2000; Field, 2002; Field et al., 1990; Gao et al., 2006; Lin&Truhlar, 2007; Mulholland, 2005; Ryde, 2003; Senn&Thiel, 2007, 2009). In this chapter, we discuss some useful QM/MM schemes and their recent applications to analyze the functional mechanisms of biological macromolecular systems. Here, we will utilize three molecular systems as examples of applications of QM/MM calculations to biological macromolecules. We will also introduce the theoretical background to investigate the functional mechanisms in proteins, such as electron transfer.


Introduction
To understand the functional mechanisms of biological macromolecular systems, investigations of both the three-dimensional geometric and electronic structures are important. Accordingly, quantum mechanical (QM) calculations are essential to obtain the properties relevant to the electronic structures of the macromolecular systems. However, because of the high computational costs, only a few hundred atoms can be included in the usual ab initio QM calculations. Therefore, even though biological macromolecular systems include a huge number of atoms, isolated QM model calculations have been performed to date, in which only the reaction centres are considered. However, this may lead to serious misunderstandings of the properties relevant to the electronic structures of the macromolecular systems.
To overcome these problems, the ab initio QM calculation is combined with a molecular mechanics (MM) calculation method; this is referred to as a QM/MM calculation. In this strategy, the QM calculation is adopted for the active site (QM regions), and for the remainder of the system, the MM calculation is adopted (MM regions). This QM/MM methodology was originally developed by Warshel and Levitt (Warshel&Levitt, 1976), and since their pioneering work, great progress has been made in the development of QM/MM algorithms and their applications to biological systems (Bruice&Kahn, 2000;Field, 2002;Field et al., 1990;Gao et al., 2006;Lin&Truhlar, 2007;Mulholland, 2005;Ryde, 2003;Senn&Thiel, 2007. In this chapter, we discuss some useful QM/MM schemes and their recent applications to analyze the functional mechanisms of biological macromolecular systems. Here, we will utilize three molecular systems as examples of applications of QM/MM calculations to biological macromolecules. We will also introduce the theoretical background to investigate the functional mechanisms in proteins, such as electron transfer.

Energy schemes of QM/MM calculations
In the QM/MM calculation, the entire system is separated into the QM and MM regions, as shown in figure 1. Since the QM and MM atoms interact with each other, the total energy can be written as follows: Here, E (ES), E (QM), E (MM), and E (QM/MM) represent the energy of the entire system, the energy of the QM region, the energy of the MM region, and the interaction energy between the QM and MM regions, respectively. In this section, we discuss how to merge the energies that were evaluated by each of the QM and MM calculation programs. Fig. 1. Division of the entire system (ES) into the QM (pink) and MM (light blue) regions.

Additive scheme
The additive scheme can integrate the QM and MM interactions, and this is referred to as a hybrid QM/MM calculation. The energy of the entire system in the additive scheme, E additive (ES), is given by the following equation: The superscript of E represents the evaluation method Here, the subscript of E represents the type of interaction. For example, elec represents an electronic interaction, bond represents a bond interaction, and vdW represents a van der Waals interaction. One can calculate the electronic interaction in the QM Hamiltonian with an MM partial charge as follows: The first term on the right hand side of eq. (4) represents the coulomb interactions between the electron density at the r i and the effective partial charges (q  ) of MM atoms at their position R  , and the second term represents the coulomb interactions between the effective partial charges (q  ) of MM atoms at their positions R  and the QM nuclear charges (Q  ) at their positions R  (figure 2). Therefore, the effects of the partial charges of MM atoms are considered as an external electrostatic field (i.e., one-electron integral term), and this can polarize the electronic structure of the QM region. The consideration of the effects of the environments (MM region) in the QM calculation is a merit of the additive scheme, and is impossible in the subtractive scheme (described in the next subsection). This advantage is clearer in some polarizable systems, such as a transition metal binding protein contained in the system.

Subtractive scheme
In the subtractive scheme, the energy of the entire system is calculated by three steps. First, the energy of the entire system is calculated by the MM calculation. Second, the energy of the QM region is calculated by the QM calculation. Finally, to correct the double counted energy, the energy of the QM region is calculated by the MM calculation. From these steps, the energy of the entire system, E sub (ES), is given in the following form: The superscript of E represents the evaluation method (QM or MM) of the system, and the region on which the calculation is performed is described in parentheses. An advantage of the subtractive scheme is the simplicity of the energy calculation (figure 2). Table 1 shows the advantages and disadvantages of the two calculation schemes. The overpolarization problem is discussed in the next section.

Treatment of a covalent bond located on the boundary between the QM and MM regions
Note that a problem occurs in the case where a covalent bond crosses the boundary between the QM and MM regions. In most biological systems, the active site contains several amino acid residues, which are connected to other amino acid residues of the protein. Therefore, termination of the QM region should be required (for example, the side chain of an amino acid residue is assigned to the QM region, and terminates at the C  ). Since the QM calculation system should be self-consistent with the QM Hamiltonian, several methods have been developed to control the boundary region between the QM and MM regions. Among these methods, the link atom method is generally used (Lin&Truhlar, 2005). In this method, the boundary atoms between the QM and MM regions are capped by a hydrogen atom, which is called a link atom, along the bond direction between the QM and MM atoms at a distance of ~1 Å. These new atoms are considered in the QM calculations. Although an advantage of this method is its simplicity, overpolarization problems occur, since this method introduces artificial atoms that are located too close to the boundary MM atoms (Lin&Truhlar, 2005;Senn&Thiel, 2009). This problem is more serious when one uses a large basis set that can describe an electronic structure more flexibly, such as plane waves (Senn&Thiel, 2009). Therefore, to overcome the overpolarization issue, several schemes have been developed. In one method (Field, 2002;Field et al., 1990), the charges of the MM atoms were not computed in the one-electron integrals. In another (Reuter et al., 2000), the charges of the MM atoms were not included in the one-electron integrals. In yet another strategy (Amara&Field, 2003;Eichinger et al., 1999), the Gaussian functions at the positions of the MM atoms that are located close to the QM region were introduced.
To resolve the disadvantages of the link atom method, the Generalized Hybrid Orbitals (GHO) method was developed. GHO is a major method to overcome the disadvantages of the link atom method, and originated from the Local Self-Consistent Field (LSCF) method (Amara et al., 2000;Assfeld&Rivail, 1996;Ferre et al., 2002;Gao et al., 1998;Monard et al., 1996;Pu et al., 2004;Thery et al., 1994). In this scheme, the "auxiliary" and "active" orbitals are located on the MM boundary atoms, respectively. By incorporating the active orbital, which is directed towards the QM boundary atoms, within the SCF calculation, the valency of the QM subsystem is satisfied, while the auxiliary orbitals are fixed during the SCF calculation.

Software for QM/MM calculations
In this section, we introduce several QM/MM calculation software packages. Various QM calculation packages, such as Gaussian (Frisch et al., 2003), include the MM module in the QM/MM calculation. In addition, several MM calculation packages, such as AMBER (Case et al., 2005) and CHARMM (Woodcock et al., 2007), can calculate QM/MM schemes by using the included QM modules. The merit of these calculation program packages is simple, and only one calculation job is needed to perform the QM/MM calculation. However, such calculation packages have some restrictions. For examples, in the Gaussian, it is impossible to perform high-level statistical calculation schemes, such as umbrella sampling. On the other hand, AMBER can only perform semi-empirical QM calculations, but not the first principle QM calculations.
To avoid such problems, several interface programs were developed, including ChemShell (Sherwood et al., 2003), QoMMM (Harvey, 2004), and PUPIL (Torras et al., 2008). With these types of methodologies, one can perform QM/MM calculations effectively and couple them to statistical sampling schemes. For example, extended conformational sampling schemes are also available by using molecular dynamics simulations implemented in MM packages, and algorithms to find optimal reaction paths and stationary points on the potential energy surface (PES) are also available, by using the program modules implemented in the QM packages.
Recently, in our laboratory, a new QM/MM interface program was developed by connecting the AMBER (Case et al., 2005) and GAMESS (Schmidt et al., 1993) packages (Kang et al., 2009). Thereby, one can couple various levels of QM calculations, such as Hartree-Fock, density functional theory (DFT), and more advanced schemes, with conformational sampling and free energy estimation techniques, such as a standard molecular dynamics (MD) calculation, replica exchange MD, and potential mean force (PMF).
The roles of the interface program are i) to control the sequence of steps for the QM/MM optimization and the MD simulation, and ii) to exchange information between the two programs. At the first stage of a calculation cycle, a single-point calculation is performed in GAMESS concerning the coordinate set of a system. This yields two types of forces: i) the forces on the QM atoms from the QM-QM and QM-MM interactions, and ii) the forces on the MM atoms derived from the QM-MM interactions. To calculate the second force, the contribution from the QM atoms to the MM atoms is calculated, by integrating the interaction between the partial charges of the MM atoms and the electron density at each grid point of the quantum region, as shown in equation (6).
Here, q j is the partial charge of an MM atom, and N is the total number of electron density points. Here, dq i designates a small volume on the electron density grid, dq i =  i dxdydz. With respect to the coordinate set, AMBER calculates (1) forces from van der Waals (vdW) interactions between all atoms of the system, and (2) forces on MM atoms from MM-MM interactions. The interface program combines the calculated forces. This step requires the mapping of the relationship between the QM and MM atoms, since the atom labeling in one program package may differ from that in the other one. In order to reduce the overhead cost of this process (in particular, the required computational time would increase when treating solvated biological macromolecules; i.e., the number of atoms included in the calculation system could exceed 100,000 atoms), we employ a combination of a UNIX shell and C programming. For combining the QM and MM calculations, we should mention another use of (advanced) QM calculations (such as Coupled Cluster); i.e., the development of effective potentials for MM calculations. A recently developed scheme can be applied to a wide range of interactions. For example, the description of cation- interactions, such as the complex of Na + and an aromatic ring, was difficult to achieve using conventional MM potentials. The new scheme mentioned above enabled us to perform accurate MD simulations involving a cation- interaction that is present in the catalytic site of an enzyme (Hagiwara et al., 2009a;Hagiwara et al., 2011).

Azurin
Biological systems possess various types of metalloenzymes, which are involved in biological functions such as electron transfer, metal storage, dioxygen binding, substrate turnover, and protein structure configuration (Holm et al., 1996). An important application of QM/MM methods for metalloenzymes is to investigate the electronic structures of the active sites in blue copper proteins (Paraskevopoulos et al., 2006;Ryde&Olsson, 2001;Ryde et al., 2000;Sinnecker&Neese, 2006;Solomon, 2006), which function in electron transfer (ET) through the bound Cu ion(s) (Gray et al., 2000). X-ray crystallographic analyses have been reported for such enzymes, including azurin, plastocyanin, and stellacyanin. The common features of the active site are one-cysteine coordination and two-histidine coordination to the copper centres, forming an approximate trigonal plane, and their coordination environments are completed with weak axial ligand(s) (Gray et al., 2000;Paraskevopoulos et al., 2006). Azurin includes a five-coordinated copper, and possesses methionine (Met) and the backbone carbonyl oxygen of glycine (Gly) as weak axial ligands (figure 3). Stellacyanin and plastocyanin have a four-coordinated copper site coordinated with glutamine (Gln) and Met as weak axial ligands, respectively. These weak coordinative bonds are presumed to contribute critically to the redox potential related to electron transfer, thereby regulating the biological functions. Thus, blue copper proteins have been used as standard metalloenzymes to examine the accuracy of calculations (Ryde&Olsson, 2001;Ryde et al., 2000;Solomon, 2006;Solomon et al., 2004).

Evaluation of QM/MM calculations: Spin polarization
The above-mentioned QM/MM schemes, the additive (termed simulation A here) and subtractive (simulation B) schemes, were adopted by exploiting our interface program that connects QM and MM calculation engines implemented in parallel supercomputers (see section 2.4) (Hagiwara et al., 2009d;Kang et al., 2009). An advantage of this strategy is that the interface program enables us to compare the results of QM/MM calculations using distinct schemes; thus, the long-range effects on the copper active site were evaluated in this theoretical analysis. Here, an ab inito DFT calculation was utilized as the QM part in the following QM/MM calculations. The QM/MM calculations were evaluated by a comparison to the experimental data, as follows. The singly occupied molecular orbital (SOMO) of the copper binding site of blue copper proteins is known as an antibonding feature, as a combination of a Cu 3d-orbital and a cysteine sulphur 3p-orbital (Solomon, 2006;Solomon et al., 2004;Sugiyama et al., 2005). In simulations A and B, the 128 and 127 -spin orbitals are occupied, and the SOMOs of simulations A and B were found in the 126 -spin orbitals, as shown in figure 4a and b, respectively. The Cu 3d-orbital and the cysteine sulphur 3p-orbital contribute to the SOMO in both simulations A and B, while some differences in the orbital patterns are seen at the Gly45 and His117 residues (figure 4). The electron density of the S atom in Cys112 is significantly delocalized onto the copper dorbital (this was confirmed by the calculation of the Mulliken charge of the S atom) (Kang et al., 2009), resulting in a large spin density on the S atom. In simulation A, the spin polarization is significantly improved (49.9%), as compared to simulation B (56.8%), since the former calculation result (simulation A) exhibits better agreement with the spectroscopic experiment; i.e., 45% S 3p-orbital character was experimentally observed in the SOMO of azurin (George et al., 2003;Solomon et al., 2004). This result indicates that the explicit inclusion of an electrostatic interaction by the QM Hamiltonian is essential for the accurate evaluation of the copper site.

Evaluation of the QM/MM calculation: Three-dimensional geometrical structure
Next, with respect to the three-dimensional geometrical structure, the atomic distances relevant to the Cu coordination are shown in Table 2. It was pointed out that the intense Xray beam reduced the blue copper site during data collection, and thus ambiguities still remain in the experimental data. Accordingly, computational investigations are crucial to elucidate the outstanding electronic and geometrical features of the weak axial Cu-S(Met121) and Cu-O(Gly45) ligands and the relatively strong Cu-S(Cys112) ligand observed in azurin. For the Cu-S(Met121) coordination, a comparison of simulations A and B shows almost the same Cu-S(Met121) distances of ~3.5 Å, indicating that in the QM/MM energy schemes, the weak Cu-S(Met121) coordination is somewhat insensitive to the polarization effect. Note that for the Cu-O(Gly45) distance (2.75-3.16 Å in the crystal structures), Hasnain and coworkers obtained the remarkably short distances of 2.55 and 2.49 Å, using the ONIOM module in the Gaussian package as a QM/MM program, in the presence and absence of the electrostatic energy term in the QM Hamiltonian, respectively (Paraskevopoulos et al., 2006). Thus, the results of the above-mentioned calculations showed much better agreement with the experimental data ( a. Experimental data of the crystal structure (Nar et al., 1991) b. EE (ME) optimized geometries (Paraskevopoulos et al., 2006) Table 2. Comparison of the three-dimensional geometrical structures of the active site of azurin.
In addition, significantly different Cu-O(Gly45) distances, 3.0 and 2.8 Å, were found between the calculated structures for simulations A and B, respectively, indicating that this weak axial Cu-O(Gly45) coordination is sensitive to the treatment of the electrostatic interaction in the QM Hamiltonian. This is a distinctive feature between the Cu-S(Met121) and Cu-O(Gly45) bonds. Note here that the Cu-O(Gly45) bond is fairly well polarized, whereas the Cu-S(Met121) bond is just slightly polarized (this was confirmed by the calculation of the Mulliken charges) (Kang et al., 2009). Next, with respect to the electronic structures of the active site of azurin, the molecular orbitals (MOs) that predominantly include the electrons involved in S(Met121)/O(Gly45) atoms were compared between simulations A and B. Here, note that O(Gly45) is included in a peptide group of the protein backbone, and that the electrons of the atom are delocalized onto the peptide group. Therefore, the contributions of the delocalized electrons should be summed up for the peptide groups of the protein backbone. As a result of the analysis, for the Cu-S(Met121) coordination, the highest occupied molecular orbitals (HOMOs) of the electrons, which are the 128 MOs in both models, were revealed to include the electrons of the S(Met121) atoms predominantly in both models, and they are equivalent to each other. The energy difference between the two HOMOs is as small as 1.8 kcal/mol, which can be considered to be marginal (

Effect of the long-range interaction
The results indicated that the treatment of the electrostatic interactions of the metal site, the surrounding protein moiety, and the solvent in the QM Hamiltonian is important for an accurate description of the electronic structure and the geometry of the blue copper site in azurin. Thus, it is suggested that the electronic structures of active sites are actually modulated by the surrounding regions through long-range electrostatic interactions, resulting in contributions to the biological functions of metalloenzymes, such as electron transfer.
In addition, for the other issues, such as the evaluation of solvent effects on the electronic structures of biological macromolecular systems, the QM/MM calculation schemes are also useful. For example, with respect to a DNA-protein complex, the recognition mode between the DNA and PU.1 (a transcriptional factor protein) influences the solvent accessibility (accessible surface area) of the DNA bases, through masking by the protein. The QM/MM calculations revealed that the differences in the solvent accessibilities of the DNA bases are closely relevant to the MO energy levels of the bases, which may contribute to the regulation of the biological functions of the molecular systems, such as the reactivity of each base (Hagiwara et al., 2010b).

Cu A site of cytochrome c oxidase
Cytochrome c oxidase (CcO), which is the terminal enzyme of the electron transport system, reduces an oxygen molecule to water, thereby generating the proton concentration gradient between the matrix and the intermembrane space of mitochondria or the periplasmic space of bacteria (proton pumping). The Cu A si t e o f C c O ( f i g u r e 6 ) r e c e i v e s e l e c t r o n s f r o m cytochrome c, thereby providing the electrons with heme a. This triggers various chemical reactions occurring from the reduced state in CcO, which are cooperatively induced, followed by oxygen binding (Kaila et al., 2010;Namslauer&Brzezinski, 2004;Papa et al., 2004). Thus, the Cu A site acts as the "gate" of the subsequent cooperative reactions in CcO. The Cu A site possesses two Cu ions and six coordinated ligands (figure 6). These two Cu ions were found to form a covalent bond characterized by a mixed-valence state (DeBeer www.intechopen.com George et al., 2001;Gamelin et al., 1998;Kroneck et al., 1988Kroneck et al., , 1989, which produces a lower reorganization energy and a higher electron transfer rate than those of the blue Cu proteins, such as plastocyanin and azurin (DeBeer George et al., 2001). In bovine CcO, the distances of Cu-S (Met207) and Cu-O (the Glu198 backbone) are slightly longer than the other covalent bonds, and thus the S and O atoms are referred to as the axial-coordinated atoms.

Controversial experimental data involving discrepancies
Each of the Cu-binding sites of plastocyanin and azurin also possesses an axial Met ligand, and the functional roles of the Cu-binding sites in the blue Cu proteins were considered to be similar to those of the Cu A site in CcO (Covello&Gray, 1990;DeBeer George et al., 2001).
In plastocyanin and azurin, this Met residue is important to control their reduction potentials (Berry et al., 2003;Garner et al., 2006). By analogy to the blue Cu proteins, Met207 of the bovine CcO may also be involved in regulating the electron transfer rate. In fact, this Met is highly conserved in various CcO proteins (Paumann et al., 2004). However, several distinct experimental results were reported. For example, in the engineered azurin, in which the Cu A site was incorporated to mimic that of CcO, the roles of the Met residue were analyzed by measuring the redox potentials of several mutants (Hwang et al., 2005;Robinson et al., 1999). This analysis revealed only slight changes in the redox potential of the mutants. In contrast, significant changes in the redox potential were observed in the Met mutants of the bacterial CcO enzymes. In fact, for the CcO from Rhodobacter sphaeroides, the replacement of the axial Met with Leu remarkably increased the redox potential, by 118 mV (Wang et al., 2002;Zhen et al., 2002). With respect to the Thermus thermophilus ba 3 oxidase, systematic mutagenesis studies of the Cu A site revealed various ranges of redox potentials in the mutants (Ledesma et al., 2007). In this manner, the experimental results relevant to the Cu A site were apparently controversial and also involved some discrepancies. In addition, the intense X-ray beam used for the structural analysis reduces the metal binding site during the data collection, thereby delocalizing the electron density of the transition metal binding site. Thus, the X-ray structure of the Cu A site may also include some ambiguities, which is a similar case to that of azurin, as discussed below. These led to some difficulties in the detailed and accurate elucidation of the functional roles of the axial ligands in the Cu A site.

Electronic structure of Cu A site
To clarify the functional roles of the amino acid residues in the Cu A site, intensive theoretical a n a l y s e s h a v e b e e n p e r f o r m e d , a n d s o m e w e r e c o u p l e d t o e x p e r i m e n t s , s u c h a s spectroscopic measurements. In a recent QM/MM analysis of the Cu A site, experimental artifacts were found in its atomic coordinates. The conformation involving the artifacts was improved in the theoretical analysis, and then explicit solvent water molecules were arranged in the calculation model. The accurate descriptions of the electronic structure of the active site were intended to elucidate the functional roles of the axial ligands, where some experimental discrepancies were found, as mentioned above. As a consequence, an interesting feature was revealed in the electronic structure of the Cu A site of CcO, which may be closely relevant to the roles of the axial residues . This finding is briefly discussed in this section.
In the Cu A structure of bovine CcO, a steric clash exists between C  of Met207 and C  of Cys196, which may be derived from the above-mentioned ambiguity of the crystallographic analysis. Accordingly, for the structural improvement, extended sampling techniques coupled to QM/MM geometry optimization were adopted for the Cu A site of bovine CcO. Thereby, the steric clash was completely removed without inducing any structural transitions, as compared with the crystallographic electron density distribution map of the Cu A site. Thus, such computational methodologies are indispensable to minimize the experimental artefacts induced by the X-ray beam in transition metal binding proteins.
The following section includes some descriptions of the electronic structure obtained by the QM/MM calculations, and comparisons with experimental data: the SOMO of the oxidized state of the Cu A site contains Cu and S (Cys) atomic orbitals (figure 7). In this hybrid QM/MM calculation of bovine CcO, the two Cu and S atoms of the Cu A site exhibited 47 % and 45 % spin densities, respectively . On the other hand, to minimize the difficulties in the experimental analysis of the Cu A site of CcO, Hay et al. incorporated the two Cu ions and the coordinated amino acid residues into azurin, to mimic the original Cu A structure and the function of CcO (Hay et al., 1996). With respect to this engineered azurin, DeBeer George et al. performed Cu L-edge and S K-edge X-ray absorption spectroscopic (XAS) experiments. This analysis revealed that the d-orbitals of the Cu atom (44 %) and the p-orbitals of the S atom (46 %) predominantly occupy the SOMO of the Cu A site (DeBeer George et al., 2001). In this manner, the results of the hybrid QM/MM calculation of CcO are quite consistent with the experimental data. In previous studies, similar comparisons were performed with respect to the calculated and experimental data, using the engineered azurin for both analyses (Xie et al., 2008). The above-mentioned QM/MM calculation results are better than the previous ones, in comparison with the experimental data. Next, to evaluate the effects of the axial Met residue, a simple model was also calculated. This model was composed of only the atoms involved in the Cu A site where the Met residue is excluded. The results were compared with the above-mentioned QM/MM calculations of the CuA site involving CcO in the calculation model. This analysis revealed that the axial Met ligand increases the energy of an MO that predominantly contains the Cu d zx orbital, whereas the other MO energy levels were influenced to a much lower extent (please note that the local coordinate system, such as d zx , depends on the calculation programs used in the analysis, and thus is not substantial) . This elevation of the energy level of the Cu d zx orbital can be understood by the hybridization of the atomic orbitals between the Cu ion and the S atom of the axial ligand, as shown in figure 8. Thus, from the view of the electronic structure of the Cu A site, the role of the axial Met ligand may be the selective modulation of the energy level of a Cu d orbital. www.intechopen.com

Inner-sphere reorganization energy of Cu A site
Furthermore, the inner-sphere reorganization energy, which is one of the factors to determine the electron transfer rate in the Marcus theory, was calculated within the framework of the QM/MM calculation, and compared with the experimental data . The Marcus theory gives the relationship between the electron transfer rate and the reorganization energy, as follows (Marcus&Sutin, 1985;Moser et al., 1992): Here, |H DA | is the electronic coupling between the initial and final states, k B is the Boltzmann constant, and ħ is the reduced Planck constant, which is defined as the Plank constant divided by 2. The G 0 represents the standard Gibbs free energy change, and the  represents the reorganization energy, which is the sum of the inner-sphere and outersphere reorganization energies. In this report, the inner-sphere reorganization energy is evaluated by using the following equation, as mentioned previously (Gorelsky et al., 2006;Kelterer et al., 2001) : Here, E represents the energy of the QM region, involving the electrostatic interaction energy between the MM1 and QM regions. The subscript, chg, refers to the charge state of the system. The other subscript, g, refers to the geometry of the redox site, which is obtained by the geometric optimization of the hybrid QM/MM calculations. Thus, each E value should be obtained by QM/MM geometry optimization. The following two calculation models were used to evaluate the axial Met residue in terms of the inner-sphere reorganization energy: A) QM-Met and B) MM-Met. Model A is the same as the QM/MM model that was described above. In model B, the axial Met residue was removed from the QM region (i.e., the Met residue was incorporated into the MM region), and the atomic partial charges of the Met reside were set to zero. The computational (model A) and experimental values of the inner-sphere reorganization energy were consistent (i.e., 228 and 250 meV, respectively). Next, the comparison of the two computational models revealed that the inner-sphere reorganization energy of model A is larger than that of model B, by as little as 44 meV. This indicates that the axial Met residue does not significantly perturb the electron transfer rate of the Cu A site of bovine CcO. In fact, the mutant in which the Met residue is replaced with Leu exhibited a lower electron transfer rate than that of the wild type, by ten-fold (Wang et al., 2002), which seems to be quantitatively consistent with the above-mentioned calculation result. Therefore, it is concluded that the axial Met residue may be the fine-modulator of the electron transfer rate of the Cu A site.

Use of techniques in information science
In this manner, the hybrid QM/MM calculations indicated that the axial Met ligand functions as a fine-modulator of the electronic structure (see 4.3) and the electron transfer rate (see 4.4) of the Cu A site of bovine CcO. Still, one cannot solve the above-mentioned experimental discrepancy in the differences of the redox potentials of the Cu A site, depending on the species. The X-ray structure of the axial Met residue of bovine CcO is shown in figure 9. The Met backbone is located in a short turn structure between an -helix and a -sheet that form a rigid conformation, and thereby, the residue side chain is also fixed together with the Cu-coordination of the S atom. To solve the experimental discrepancy, the structural changes of the Cu A sites were examined with respect to the mutants of oxidases from various species, by combining amino acid multiple sequence alignment and computer modeling techniques. The amino acid sequence alignment of the CcOs revealed that the Rhodobacter and bovine oxidases (large changes were observed in their redox potential values) possess inserted modular structures close to the Cu A sites, whereas the Thermus thermophilus ba 3 oxidase and the engineered azurin (small changes in the redox potentials) do not. In fact, in the crystal structures of the Rhodobacter and bovine oxidases, the axial Met residues contact the inserted modular structural elements. In contrast, in the crystal structures of Thermus thermophilus ba 3 oxidase and the engineered azurin, the Met residues can apparently move more flexibly, and so the replacement of the Met residue may cause only small perturbations of the properties relevant to the Cu A sites. Computer-assisted structural modeling of the various mutants confirmed these findings .
Thus, the use of information science techniques (bioinformatics) is also crucial to perform thereotical analyses employing ab inito QM/MM electronic structure calculations of www.intechopen.com biological macromolecular systems. The development of methodologies in both fields is required, and the methodologies should be unified to analyze biological systems.

Double-stratified architecture of the Cu A site
On the other hand, with respect to the Cu-coordinated Cys/His residues, their replacement influenced the geometric and electronic structures of the Cu A site much more drastically Zhen et al., 2002;Zickermann et al., 1997), resulting in structural disruptions and significant increases in the redox potential . For instance, the replacement of each cysteine residue (the mixed valence state) with serine in the engineered azurin (C112S and C116S) induced the different electronic states (i.e., the C112S mutation induced two distinct Cu centres, each of which is type II, while the C116S mutation induced one Cu centre, which is type I or II depending on pH) . Thus, these Cys/His residues play a critical role in establishing the geometric and electronic structures of the Cu A site. In this manner, the functional roles of the amino acid residues of the Cu A site have been revealed. The two Cys/His residues build the fundamental geometric and electronic structures of the Cu A site; together, these are considered as the "platform" of the Cu A structure, since the replacements of those residues exert critical effects on the Cu A system (figure 10). In contrast, the axial Met ligand functions in the fine modulation of the electronic structure of the "platform"; thereby, this residue may regulate the electron transfer rates from Cyt c to heme a, thus inducing the various subsequent cooperative reactions occurring in CcO (figure 10) . This means that the Cu A site forms the double-stratified architecture, comprising the platform (i.e., Cu 2 , Cys S atoms, His N atoms) and its fine-modulators (Met S and O atoms). This is not trivial. In fact, the weakly coordinated axial Met ligand in the active site of native (wild type) azurin seems to act as a member of the "platform" residues (see 4.1). In addition, the critical roles of the axial Met residues in the Cu A sites had been proposed for the oxidases and the engineered azurin, as mentioned in section 4.2; this suggests that the Met ligand is a member of the "platform" residues. However, it has been concluded that the Met ligand is not a "platform" residue, but a "fine-modulator", on the basis of the analysis using ab initio QM/MM electronic structure calculation and computer modeling techniques coupled to multiple sequence alignment. This type of hierarchical structure relevant to the function may be widely utilized in biological systems. For example, some protein structures seem to be established by a similar architecture; i.e., fundamental folds are common in various proteins, and the addition of specific inserted modular elements into such fundamental folds leads to new biological functions of the proteins in molecular evolution. Such multi-stratified architecture and molecular evolution can be found in various biological systems, as observed in the Cu A sites of various oxidases.

Aminoacyl-tRNA synthetases
Aminoacyl-tRNA synthetases (aaRSs) are responsible for the correct attachment of their cognate amino acid to the 3'-end of the specific tRNA (aminoacylation). Thereby, these enzymes are critical for the conversion of the nucleotide sequences of mRNAs into the amino acid sequences of proteins, in the translational process of every known organism. Translational fidelity is determined by the accuracy with which specific amino acids are attached to their appropriate tRNAs. This reaction proceeds as follows: first, an amino acid is activated to an aminoacyl adenylate by ATP transfer, with the generation of pyrophosphate; second, the amino acid moiety of the aminoacyl adenylate is transferred to the 3'-end of the specific tRNA. The aaRSs are divided into classes I and II, according to their primary and tertiary structures, and are further subdivided into subclasses, a, b and c, within each class (Sankaranarayanan&Moras, 2001). The class I aaRSs are characterized by the nucleotide binding (Rossmann) fold on which the active site is located, whereas the active sites of the class II aaRSs comprise a characteristic seven-stranded  sheet with flanking  helices (Cusack et al., 1990;Eriani et al., 1990). As mentioned above, the fidelity of translation is assured by the strict discrimination of cognate from non-cognate amino acids. However, for the leucine, isoleucine, valine, threonine, alanine and phenylalanine systems, which each share structural similarity to some other systems, their cognate enzymes, i.e., the leucyl-(LeuRS), isoleucyl-(IleRS), valyl-(ValRS), threonyl-(ThrRS), alanyl-(AlaRS) and phenylalanyl-(PheRS) tRNA synthetases, have difficulties in the strict discrimination of their specific amino acids, thus producing mis-activated amino acids or mis-aminoacylated tRNAs. For example, when isoleucine is replaced with valine in the IleRS system, the error rate is approximately 1 in 5 (Pauling, 1957). This is caused by the similarity in the chemical structures of the amino acids, and results in misactivated amino acids or misaminoacylated tRNAs (Fukai et al., 2000;Fukunaga et al., 2004;Fukunaga&Yokoyama, 2006;Lincecum et al., 2003;Mursinna et al., 2001;Nureki et al., 1998;Silvian et al., 1999;Zhai&Martinis, 2005).
To avoid such incorrect products, these enzymes catalyse editing reactions. Two types of editing, pre-transfer editing and post-transfer editing, correct mis-activated amino acids and mis-aminoacylated tRNAs, respectively. A misactivated amino acid is hydrolysed to the amino acid and AMP by the pre-transfer editing pathway, and a mis-aminoacylated tRNA is hydrolysed to the amino acid and tRNA by the post-transfer editing pathway (Fukai et al., 2000;Fukunaga et al., 2004;Fukunaga&Yokoyama, 2005;Lincecum et al., 2003;Mursinna et al., 2001;Nordin&Schimmel, 2002;Nureki et al., 1998;Sankaranarayanan&Moras, 2001;Silvian et al., 1999;Zhai et al., 2007). In this manner, the overall error occurring in the abovementioned isoleucine system is thus reduced to only 1 in 40,000 (Freist et al., 1985).

QM/MM molecular dynamics simulation
Ab initio QM/MM molecular dynamics (MD) simulation is a state-of-the-art theoretical methodology in current computer simulation techniques for huge molecular systems. To investigate the mechanism of the post-transfer editing reaction by the class Ia Thermus thermophilus leucyl-tRNA synthetase (LeuRS) complexed with valyl-tRNA Leu , hybrid QM/MM MD simulations were performed coupled to ab initio DFT calculations as the QM simulation (Hagiwara et al., 2010a). Several possible reaction pathways were explored by combining umbrella sampling techniques with hybrid QM/MM MD simulations to obtain the PESs of the reaction pathways, and the calculated activation barrier led us to identify the most preferable reaction pathway. The result obtained was consistent with the biochemical experimental data, and thus a novel enzymatic reaction mechanism was revealed by the theoretical strategy, as discussed below.
In contrast, reaction pathway exploration by experimental techniques is very difficult. Even if such methodologies are available to identify the enzymatic mechanisms, theoretical approaches are also indispensable to elucidate the electronic structure changes of the active site of the enzymatic systems. This is a recent successful case. figure 11 depicts the entire calculation system of the LeuRS•valyl-tRNA Leu complex, including solvent water molecules (in total ~165,000 atoms), as well as the active site of the editing reaction (Hagiwara et al., 2009b;Hagiwara et al., 2009c). To obtain the PES relevant to the enzymatic reaction, an adiabatic mapping approach was adopted. This method is based upon conformational/configurational sampling by QM/MM MD simulation and QM/MM geometry optimization. The role of the QM/MM MD simulation is the enhancement of such sampling. www.intechopen.com The reaction space was divided into windows, each defined by a particular set of values of the reaction coordinate variables, and then, for each window, a QM/MM MD simulation was performed, followed by a geometry optimization. Both the MD simulations and geometry optimizations were performed with functions that constrained the reaction coordinate variables, to keep them close to their window reference values. The final PESs were then reconstructed from the energies of the geometrically optimized structures, but without the constraint energy terms.

Editing mechanism: A novel enzymatic reaction by a hybrid ribozyme/protein catalyst
This computational analysis revealed that the water molecule that acts as the nucleophile in the editing reaction is activated by a 3'-hydroxyl (3'-HO) group at the 3'-end of tRNA Leu , and that the O2' atom of the leaving group of the substrate is capped by one of the water's hydrogen atoms. This was shown to be consistent with experimental data. Thus, it is the tRNA itself, and not the protein, that drives the editing reaction. Unexpectedly and surprisingly, editing was found to be a ribozymal self-cleavage reaction. The protein does, however, have an important stabilizing effect on certain high-energy intermediates along the reaction path, and thereby plays a critical role in promoting the reaction and providing specificity in ligand recognition. These results led us to call the LeuRS•valyl-tRNA Leu complex a "hybrid ribozyme/protein catalyst", which is more efficient than the ribozyme alone (Hagiwara et al., 2010a). Furthermore, examinations of previously reported structural and biochemical data revealed that this RNA-driven catalytic mechanism seen in the LeuRS•valyl-tRNA Leu complex appears to be quite widespread, as it occurs in the deacylation and peptidyl-transferase reactions of the ribosome as well as the editing reactions in the other class (class II) of aaRSs. The existence of hybrid ribozyme/protein catalysts, exemplified by the LeuRS editing mechanism, also provides insight into the types of transitional forms that could have been important in the RNA world hypothesis for the origin of life (Cusack et al., 1990;Eriani et al., 1990). In this manner, by employing hybrid QM/MM MD simulations, we found a new biological catalyst. This system is distinct from conventional RNA•protein complexes such as the group I intron: the proteins (aaRSs) directly contribute to the reaction through hydrogenbonding with RNA moieties, thus stabilizing the high-energy intermediates. Recently, a similar theoretical analysis of the editing mechanism was conducted by employing our fully-solvated structural model (figure 11), and the author reported almost the same reaction mechanisms as described above (Boero, J. Phys. Chem., 115, 12276, 2011). In this report, an additional possible initial reaction pathway was proposed; however, this process is not feasible, as shown in the previous study (see Simulations 2 and 4 in Hagiwara, et al., 2010a). The difference (artefact) could be due to the insufficient QM regions lacking the crucial atoms, the functional used in the calculations (see Tuckerman, et al., 2006), etc.

Dynamic electronic structure rearrangements in hybrid ribozyme/protein catalysis
The dynamic changes that occur in the electronic structure of the catalytic site were investigated in the hybrid ribozyme/protein catalysis by the LeuRS•valyl-tRNA Leu complex. As a result of this theoretical analysis, dramatic, functionally important rearrangements of the MOs were observed (Kang et al.). We will briefly describe these dynamic rearrangements of the electronic structure in the biological macromolecular systems. The left and right panels represent the initial and transition states, respectively, in the editing reaction. In the initial state of the editing reaction, the orbital that would be reactive for the bond formation is not HOMO, but HOMO-14. In fact, the energy level of HOMO-14 is elevated as the enzymatic reaction proceeds (i.e., nucleophilic attack occurs). In the transition state, this MO becomes HOMO, which leads to the formation of a new covalent bond between the water and the carbonyl group.
Two changes are particularly noteworthy. One concerns the MO that contains a contribution from the nucleophile (a water molecule). It initially has a much lower energy than the HOMO, but it is activated as the reaction proceeds until it becomes the reactive HOMO ( figure 12). The other involves the reactive LUMO with anti-bonding character, which emerges in the bond rupture that leads to the products. We refer to these processes as the dynamic induction of the reactive HOMO (DIRH) and LUMO (DIRL), respectively. Interestingly, the induction of the reactive HOMO is enhanced by the formation of a low-barrier hydrogen bond (LBHB), which represents a novel role for LBHBs in enzymatic systems.
Previous reports of computational studies of enzymatic reactions have discussed the electronic structures of the catalytic sites in terms of frontier orbital theory (FOT). However, the above-mentioned DIRH/DIRL are unusual, in that we have highlighted the dynamic rearrangements of the MOs that are relevant to catalysis along the reaction path. In small compound systems, the reactive HOMO and LUMO are already present in the reactant structures, since the majority of the atoms are directly relevant to the chemical reaction, and so dramatic rearrangements in the MOs are not required. By contrast, in the complex environment of a biological macromolecule, the reactive HOMO and LUMO are much more likely to be hidden, and so these systems must implement the DIRH and DIRL mechanisms that expose them. Therefore, DIRH and DIRL are probably widespread in the reactions of complex and flexible molecular species, such as enzymes, and thus their analysis should be an essential part of the investigation of many reaction mechanisms. In this manner, state-of-the-art methodologies in computational science are essential to facilitate the elucidation of these complex electronic processes.

Conclusion
In this chapter, we introduced two QM/MM schemes, the additive and subtractive schemes, and discussed their advantages and disadvantages. We then showed three types of examples of their applications to biological macromolecular systems. Here, the ab inito DFT calculation was utilized as the QM part. In the first application, we compared the calculation results using the two schemes, to evaluate the effects of the long-range electrostatic interactions to the QM Hamiltonian (i.e., the effects of polarizability). As a consequence, the additive scheme was found to be more reliable to reproduce the experimental data. Therefore, it was concluded that the environmental structures of the protein and solvent water molecules crucially affect the geometric and electronic structures of the active site of a protein. Thereby, biological macromolecular systems may regulate the properties of their active sites. In the second application, a computational technique that extensively explored the conformational space (involving hybrid QM/MM calculations) was important to remove the incorrect structures (e.g. steric clashes) included in the experimental structure. Utilizing the corrected atomic coordinates of the Cu A site of CcO, the functional roles of the axial Met ligand were investigated, where some properties related to the Met residue showed experimental discrepancy. A comparison of the electronic structure obtained using the hybrid QM/MM calculation of CcO with that of a simple model system lacking the axial Met ligand revealed that the effect of this Met residue is not significant, but "selective". Thus, the axial Met ligand was theoretically found to be a fine modulator of the electronic structure of the Cu A site, although the experimental data were controversial concerning the role of this residue. This fine modulation may contribute to a more effective oxygen reduction process, by optimizing the electron transfer reaction to the environment. In this manner, hybrid QM/MM approaches are crucial to resolve the experimental discrepancy related to the electronic properties of biological macromolecules.
In the final application, we showed an example of a hybrid QM/MM MD simulation that was used to elucidate the editing mechanism of the LeuRS•valyl-tRNA Leu complex. The results of this analysis revealed a novel biological catalyst, a hybrid ribozyme/protein enzyme. In this catalysis, the dynamic rearrangement of the electronic structures was discovered. Thus, hybrid QM/MM MD simulation is currently one of the most powerful theoretical methodologies to elucidate the functional mechanisms of enzymatic systems.