Systematic QM Region Construction in QM/MM Calculations Based on Uncertainty Quantification
While QM/MM studies of enzymatic reactions are widely used in computational chemistry, the results of such studies are subject to numerous sources of uncertainty, and the effect of different choices by the simulation scientist that are required when setting up QM/MM calculations is often unclear. In particular, the selection of the QM region is crucial for obtaining accurate and reliable results. Simply including amino acids by their distance to the active site is mostly not sufficient as necessary residues are missing or unimportant residues are included without evidence. Here, we take a first step towards quantifying uncertainties in QM/MM calculations by assessing the sensitivity of QM/MM reaction energies with respect to variations of the MM point charges. We show that such a point charge variation analysis (PCVA) can be employed to judge the accuracy of QM/MM reaction energies obtained with a selected QM region, and devise a protocol to systematically construct QM regions that minimize this uncertainty. We apply such a PCVA to the example of catechol O-methyltransferase, and demonstrate that it provides a simple and reliable approach for the construction of the QM region.Our PCVA-based scheme is computationally efficient and requires only calculations for a system with a minimal QM region. Our work highlights the promise of applying methods of uncertainty quantification in computational chemistry.
Introduction<!>Methodology<!>Sensitivity analysis for QM/MM energies<!>Computational Details<!>Modification of the input files concerning point charge variation and analysis of the results<!>Systematic QM Region Construction<!>Size<!>PCVA<!>Conclusion
<p>Since the introduction of the QM/MM approach by Warshel and Levitt in 1976 1 it has evolved to a broadly used tool in computational chemistry. Dividing a large biomolecular system into two subsystems, one smaller subsystem treated quantum-mechanically (QM) and one larger subsystem treated with molecular mechanics (MM), allows one to investigate the mechanisms and energetics of enzymatic reactions in an efficient way. [2][3][4] Consequently, the application of the QM/MM approach to large biomolecules such as enzymes has become a common practice in the last two decades. [5][6][7][8][9][10][11] However, in practice setting up QM/MM calculations is far from trivial and requires many manual choices by the simulation scientist. 4,8 Besides the selection of a suitable QM method and an MM force field, the most important decision is the choice of the QM region.</p><p>The convergence of QM/MM results with the choice and particularly the size of the QM region has been investigated in several studies, [12][13][14][15][16] which underline that it generally has a large effect on the quality of the final results, and that often rather large QM regions are required for reaching converged results. To alleviate this problem, schemes for the systematic construction of QM regions have been proposed, usually with the aim of obtaining medium-sized QM regions that provide reliable QM/MM reaction energies. Examples of such schemes include free energy perturbation analysis, 17 charge deletion analysis, 18 charge shift analysis (CSA), 15 and Fukui shift analysis (FSA). 19 For their recently developed selfparametrizing system-focused atomistic models (SFAM), Brunken and Reiher proposed an automatic scheme for the construction of hybrid QM/SFAM models, including a systematic determination of the QM region based on the energy gradient. 20 For a given QM region, there are different possible choices of the embedding method, 21 of a suitable coupling scheme, 22 and for the treatment of the boundary region 23 if covalent bonds cross the border between the subsystems. For the commonly used electrostatic embedding scheme, 4 the choice of the MM point charges will influence the resulting QM and QM/MM energies, and further parameters need to be chosen in advanced polarizable embedding 24,25 or flexible embedding schemes. 26,27 Different coupling schemes are available such as IMOMM, 28 ONIOM, 29 or the AddRemove 30 model. The most common approach for treating covalent bonds across the QM-MM boundary is saturating the QM region with capping atoms. 4 Here, the position and the type of these atoms is crucial for the calculations. An alternative are frozen orbitals, 31 e.g., in the Localized SCF (LSCF) 32 or the Generalized Hybrid Orbital (GHO) 33 approach. The uncertainties introduced by all these different choices and the corresponding parameters are interconnected and will again depend on the choice of the QM region.</p><p>Consequently, there is a need for rigorous uncertainty quantification for QM/MM calculations, i.e., to systematically assess the sensitivity of the QM/MM energy with respect to these technical choices and empirical parameters and to ultimately provide rigorous error bounds on the results of QM/MM calculations (compared to a full QM treatment). Mathematical and computational tools for quantifying uncertainties in computer simulations have been developed intensively in the past decades (for textbooks, see, e.g. Refs. 34,35) and are employed in many areas of simulation science, 36,37 but their application is just starting to emerge in computational chemistry. 38 Recently, we have applied such tools for analyzing the sensitivity of calculated spectra with respect to distortions of the molecular structure. 39,40 Here, we aim at taking a first step towards uncertainty quantification for QM/MM methods by analyzing the sensitivity of QM/MM reaction energies with respect to variations of the MM point charges. While there are other relevant empirical parameters entering in QM/MM calculations, most importantly those related to the treatment of the QM-MM boundary (e.g., to the placement of the link atoms), we expect the MM charges to be a key factor influencing the final QM/MM results.</p><p>The ability to quantify the sensitivity of QM/MM calculations with respect to parameters of the MM environment provides a natural starting point for guiding the systematic choice of the QM region. With increasing size of the QM region and approaching a full QM calculation, one can expect this sensitivity to decrease. Therefore, choosing the QM region such that the uncertainty is reduced implies that the QM/MM calculation approach those of a full QM calculation. Here, we exploit this idea by proposing a simple and efficient scheme for the systematic construction of the QM region that is guided by uncertainty quantification. This work is organized as follows. In Section 2 we recall the necessary theoretical background of QM/MM approaches (Sect. 2.1, introduce our point charge variation analysis (PCVA) for analyzing the sensitivity of QM/MM energies (Sect. 2.2, and give the computational details (Section 2.3). In Section 3, we introduce the model system used in this work and discuss the convergence of ligand charges and reaction energies for QM regions of increasing size. The sensitivity of these quantities with respect to global point charge variations is analyzed in Section 4. This is followed in by the evaluation of the energy sensitivity for single amino acids in Section 5, which are used to devise a scheme for the systematic construction of QM regions based on a PCVA. The QM regions obtained with this PCVAbased scheme are assessed for QM regions of increasing size and for atom-economical QM regions in Sections 6 and 7, respectively. Finally, conclusions and an outlook can be found in Section 8.</p><!><p>2.1 QM/MM energy partitioning QM/MM is based on the partitioning of the full target system into a QM region (A) including the interesting part of the system, such as the active center of an enzyme, and an MM region (B) containing all other atoms, i.e., the active center's environment. The total energy can be expressed as the sum of the QM energy of subsystem A, E QM (A), the MM energy of subsystem B, E MM (B), and an interaction energy between the two subsystems, E int (A, B),</p><p>The exact definition of these three energy contributions varies between different implementations of QM/MM schemes. 4,8 In the following, we focus on the general principles as far as they are relevant for the current work. Details on the QM/MM implementation employed here are given in Section 2.3.</p><p>The energy of the QM region, E QM (A), is obtained from a quantum-chemical calculation of the corresponding subsystem A. To allow for covalent bonds to cross the boundary between the QM and MM regions, the QM subsystem is usually saturated using capping atoms. The energy of the MM region, E MM (B), is calculated using a classical force field and contains the usual bonding and non-bonding force-field energy contributions. Finally, the interaction energy E int (A, B) contains both electrostatic interactions [E int,el (A)] and non-electrostatic interactions [E int,ne (A, B)] between the two subsystems.</p><p>In the electrostatic embedding scheme, which is commonly used when applying QM/MM to enzymatic reactions, the electron density of subsystem A is polarized by the MM point charges of subsystem B. The MM point charges are included in the QM Hamiltonian 4 to compute the interaction between the electron density ρ A (r) of subsystem A and the electrostatic potential V B (r) derived from the point charges of subsystem B,</p><p>where N B is the number of atoms in subsystem B and q I,B represents the MM point charge and R i,B the position of the I-th atom. This electrostatic interaction energy is usually included in the QM energy of subsystem A, i.e.,</p><p>Altogether, the QM/MM energy can be expressed as</p><p>and the non-electrostatic interactions between the two subsystems.</p><p>When investigating enzymatic reactions with QM/MM calculations, the main quantity of interest (QoI) is generally the reaction energy,</p><!><p>The reaction energy calculated within a QM/MM model is subject to numerous sources of uncertainty (see Introduction). One important element of uncertainty quantification 34,35 is the analysis of the sensitivity of the simulation results with respect to its input parameters. 41 Here, we consider the QM/MM reaction energy ∆E reaction QM/MM as our quantity of interest (QoI) and analyze how sensitively it depends on parameters of the QM/MM model.</p><p>We consider the MM point charges q MM as one of the most important sources of uncertainty and want to systematically analyze the effect of variations in the MM point charges on our QoI, i.e., the reaction energy ∆E reaction QM/MM (q MM ). To this end, we follow our earlier work on the sensitivity of calculated spectra with respect to distortions of the molecular structure 39 and consider a collective variation of the MM point charges, i.e.,</p><p>where q MM is a vector of size N B containing all MM point charges, q 0 MM is the vector of the undistorted MM point charges as provided by the employed force field, and ∆q MM is a collective variation of these point charges, which depends on a parameter ∆q that controls the size of the variation. We chose the collective variations of the MM point charges such that I ∆q MM,I = 0, i.e., the sum of the MM point charges is preserved.</p><p>In the following, we will consider two types of collective point-charge variations. First, we change the charges of all protein MM atoms simultaneously by an equal magnitude ∆q, while changing all solvent MM charges equally such that the total charge is preserved, i.e,</p><p>where N protein B and N solvent B are the numbers of protein and solvent atoms in subsystem B, respectively. Second, we consider variations of the MM charges of the i-th amino acid,</p><p>where N aa,i is the number of atoms in the i-th amino acid. The first will provide an estimate of the overall sensitivity of the QM/MM reaction energy to variations of the MM point charges, while the second will allow us to assess the effect of the individual single amino acids.</p><p>For these collective point-charge variations, we perform a local sensitivity analysis 41 and consider the derivative of ∆E with respect to a the parameter ∆q,</p><p>that is, the derivative is taken in the direction of the collective point-charge variation.</p><p>Of the components of the QM/MM energy [see Eq. ( 4)], the non-electrostatic interaction energy does not depend on the MM point charges. Therefore, the sensitivity of the QM/MM energy of the reactants or the products is given by,</p><p>and the sensitivity of the reaction energy can be calculated as</p><p>where the superscripts R and P designate the reactants and products, respectively. If the protein environment is similar for the reactants and the products, which is usually the case for enzymatic reactions, the second term can be expected to be small and could possibly neglected, i.e,</p><p>The simplest way of obtaining the derivatives necessary for the calculation of the sensitivity δ∆E reaction QM/MM is their numerical evaluation, either using a symmetric two-point finite difference formula,</p><p>or using a forward two-point finite difference formula,</p><p>An analytical evaluation of these derivatives is also possible, but would require modifications of the quantum-chemical software packages used for subsystem A.</p><!><p>Molecular dynamics calculations were performed using GROMACS 2019.3 42,43 with the AMBER99SB-ILDN 44 force field. The already equilibrated initial structure provided by Ku-lik et al. in the Supporting Information of Ref. 15 was solvated in TIP3P 45 water molecules in a cubic simulation box with 1 nm distance between the borders and the enzyme. The system was neutralized by adding six sodium cations. The positions of the solvent molecules and ions were minimized using the force field, while the enzyme structure was held fixed.</p><p>Subsequently, a droplet was extracted including COMT with substrates, sodium ions and all water molecules within 33 Å from the COMT center of mass for the following QM/MM calculations.</p><p>All QM/MM calculations were performed using the Amsterdam Modeling Suite (AMS Version 2020.203). 46 The Amsterdam Density Functional (ADF) engine 47 was used for the QM part applying density functional theory (DFT) with the PBE exchange-correlation functional 48 employing a DZ and a TZP Slater-type orbital basis set 49 for all geometry optimizations and single point calculations, respectively. For the MM region the ForceField engine of AMS was used with the AMBER95 force field, 50 which was extended by parameters for SAM and catecholate using Antechamber 51,52 and acpype. 53,54 Electrostatic embedding as implemented in AMS 55 was applied for the interaction be- All QM/MM geometry optimizations were performed using the FIRE minimization algorithm 56 with all solvent molecules fixed to their initial coordinates. All charges evaluated for charge convergence tests are calculated from the Voronoi deformation density (VDD) 57 of the reactant structure only.</p><!><p>were achieved using Python. Plots were generated with Matplotlib 58,59 and structures were visualized using Vmd. 60 A data set containing PDB files of the reactant and product starting structures, a modified AMBER95 force field file, AMS fragment files for the ligands and ions, and the AMS input files for all geometry optimizations and single point calculations is available at Ref. 61.</p><p>3 Catechol O-methyltransferase as model system for QM/MM calculations As model system for investigating the sensitivity of QM/MM calculations on point-charge variations and for exploring automatic QM region selection schemes, we chose the enzyme catechol O-methyltransferase (COMT), 62 which plays a crucial role in the regulation of neurotransmitters in the human body. In a previous study Kulik et al. 15 investigated the convergence of catalytic properties with increasing QM region size and established it as a test case for benchmarking QM/MM approaches, in particular of schemes for the systematic determination of the QM region. 14,19,64 We use the structural model of COMT introduced by Kulik et al., starting from the initial equilibrated MM structure taken from the Supporting Information of Ref. 15 based on the protein crystal structure (PDB: 3BWM). 63 The COMT active site includes the neutral Sadenosyl methionine (SAM), the catecholate anion (CAT), and the catalytically active Mg 2+</p><p>(see Fig. 1).  All QM calculations were performed using the GGA exchange-correlation functional PBE. In contrast to Ref. 15, we did not encounter a closing of the HUMO-LUMO gap with increasing size of the QM region. Instead, the HOMO-LUMO gap remained constant at about 1 eV for the larger QM regions (see Supporting Information, Fig. S1). Note that the results presented in Ref. 15 that will be discussed in the following, have been obtained with the range-separated hybrid functional ωPBEh, which also avoids a spurious closing of the HOMO-LUMO gap.</p><p>To test the overall QM region size convergence we first evaluated the distances between SAM methyl and the catecholate oxygen atom in the reactant structures (see Supporting Information, Fig. S2). The SAM to catechol distance in the initial MM-equilibrated structure before QM/MM optimization is 3.11 Å which is similar to the most likely distance in the underlying distance distribution. 64 With increasing QM region size the distance significantly decreases for regions 2 and 3 from about 3.05 Å to 2.7 Å and 2.8 Å, respectively. With a distance of 2.97 Å for region 4 and 5 the SAM-CAT distance starts converging for region 6 and larger to between 2.8 and 2.9 Å. Given the differences in the QM treatment, this behavior is in reasonable agreement with the results of Ref. 15, particularly for the larger QM regions. The absence of only one of these three residues causes the charge of +0.55.</p><p>Overall, it can be stated that small QM regions are not sufficient for reproducing the ligand charges found for large QM regions because important residues coordinating the ligands might be missing. When using a distance-based construction of the QM region, rather larger QM regions need to be reached before all relevant residues are included in the QM region.</p><p>Finally, the QM/MM reaction energy for QM regions of increasing size is shown in Fig. 2C.</p><p>For small QM regions, the reaction energy shows large oscillations, while starting from QM Overall, our results confirm the slow convergence of the reaction energy with increasing size of the QM region found in earlier studies and emphasize the need for systematic protocols for the construction and selection of the QM region. We note that Jindal and Warshel 14 found that in COMT, the activation barrier shows a much smaller sensitivity with respect to the choice of the QM region than the reaction energy. Therefore, we will not consider activation barriers and focus on the reaction energies in the present work. In addition to the VDD charges obtained with the undistorted MM charges for the different QM regions considered above, Fig. 2A includes the VDD charges obtained with collective point-charge variations ∆q tot MM with ∆q between +0.005 and −0.005. The corresponding sensitivities (defined in analogy to Eq. ( 9) for the VDD charges and evaluated numerically using a symmetric two-point finite difference formula [Eq. ( 13)] with ∆q = 0.005) are plotted in Fig. 2B.</p><p>It can be seen that the variation of all MM point charges slightly affects the VDD charges.</p><p>Significant deviations from the unvaried curve are visible for ∆q = ±0.005. Overall, the variations do not affect the charge convergence behavior, but lead to different ligand-dependent observations concerning the sensitivity. For CAT and Mg 2+ , the VDD charge sensitivity decreases with increasing QM region size (with the exception of region 1 for CAT) and converges for region 4 and larger. For the largest QM regions 8 and 9, the sensitivity of the Mg 2+ charge is reduced to almost zero.</p><p>For SAM, in contrast to CAT and Mg 2+ , the VDD charge sensitivity is initially increasing until a QM region size of about 300 atoms is reached (region 4), before the sensitivity starts to slightly decrease when further enlarging the QM region. This different behavior arises because SAM is a much larger molecule than CAT and thus offers numerous possibilities for charge redistribution and more contact sites to adjacent residues. For the small QM regions, the possibilities of charge redistribution to adjacent QM residues are reduced, which results in a smaller VDD charge sensitivity in these cases. In region 4, five residues (GLY65, TYR67, TYR70, SER71 and ILE90) which are part of the SAM coordination sphere are added to the QM region. This enables a wide variability for the SAM charge to be redistributed, resulting in an increase in sensitivity for this QM model. In larger QM regions, the SAM charge sensitivity is then gradually decreasing because only single residues being part of the SAM coordination sphere are added, such as MET39 in region 5 or TRP142 in region 6.</p><p>The plot of the QM/MM reaction energies for the different QM regions in Fig. 2C also includes the reaction energies obtained for collective point-charge variations ∆q tot MM , whereas Fig. 2D shows the corresponding sensitivities δ∆E reaction QM/MM . The point charge variation leads to significant changes in the reaction energy especially for a variation of ∆q = ±0.005, for which changes of up to 8 kcal/mol could be observed for small QM regions. The convergence behavior of the energy itself is not affected by the point charge variation. The sensitivity starts decreasing strongly with region 4 while the reaction energy converges towards the one found for the largest QM region 9. An exception is found for QM region 8, for which the reaction energy does not follow this trend. This can be attributed to the inclusion of four rather critical charged residues in region 8, namely the negatively charged GLU5, GLU63, and ASP168 as well as the positively charged LYS45. Remarkably, for this outlier the sensitivity is also strongly increased.</p><p>Overall, these first point charge variation tests show that small changes in the MM point charges have an impact on the QM region. While VDD charges and the reaction energy slowly converge for larger QM regions, the corresponding sensitivities to global point-charge variations decrease, whereas outliers are accompanied by an increased sensitivity. Generally, sensitivities are smaller for reaction energies closer to our best estimate (i.e., the reaction energy obtained for the largest QM region). This indicates that the sensitivity to global point-charge variations might indeed be useful as an indicator for the reliability of QM/MM calculations, and that systematically reducing this sensitivity could be a promising strategy for the systematic construction of the QM region.</p><!><p>Motivated by the results of a global point charge variation analysis presented in the previous section, we set out to develop a protocol of the systematic construction of the QM region that aims at minimizing the sensitivity of the QM/MM reaction energy. To this end, we consider the sensitivity of the QM/MM reaction energy with respect to variations of the point charges in single amino acids, i.e., we perform a single amino acid PCVA. The resulting protocol for systematic QM region construction is summarized in Fig. 3 and will be described in the following.</p><p>As starting point, we consider the minimal QM/MM model with QM region 1, i.e., for our COMT test case only the ligands and the catalytically active magnesium ion are included in the QM region. For single-point energy calculations of the geometry-optimized reactant and product structures, we calculate the sensitivity of the QM/MM energy with respect to variations of a single amino acid ∆q aa,i MM,I [cf. Eq. ( 8)]. Here, we use ∆q = −0.5, i.e., the total charge of the considered amino acid is decreased by 0.5 while an equal charge of opposite sign is distributed over all other MM atoms.</p><p>To reduce the computational effort for the evaluations of this sensitivity for each amino acid, we assessed several possible simplifications compared to the global PCVA in the previous section. The different levels of simplification are referred to as PCVA-A to PCVA-E, and the sensitivities obtained in these different approximations are shown in Fig. S3. First, instead of the full QM/MM reaction energy (PCVA-A) we use only the QM contribution (PCVA-B), i.e., the approximation of Eq. ( 12) is employed. This is justified as the effect on the QM energy is the main focus of the analysis and possible undesired MM-only effects will be excluded. Furthermore, no significant differences in the calculated sensitivities are observed comparing the full and QM energy approach. Second, instead of the reaction energy we could use the sensitivity of the product (PCVA-C) or reactant (PCVA-D) energy instead.</p><p>While this does change the values of the sensitivities, it leads to overall similar trends (cf. Fig. S3). Using the product or reactant energy sensitivities could also be advantageous for avoiding error cancelation that might be present when considering the sensitivity of the reaction energy only. Consequently, we choose to use only the reactant energies, which reduces the number of calculations to be performed by half. The impact for this approximation on the selection of the QM region will be discussed below. Finally, instead of using a symmetric two-point formula for the numerical differentiation, it turns out to be sufficient to use a forward finite-difference formula (PCVA-E), which again reduces the number of QM calculation by another factor of two.</p><p>Altogether, we employ the PCVA-E approximation and calculate the sensitivity with respect to point-charge variations for the i-th amino acid as,</p><p>where we employ a point charge variation of ∆q = −0.5 (PCVA-E). For our COMT test case, this requires 214 QM calculations for the minimal QM region with different sets of MM point charges.</p><p>The resulting sensitivities δ i E QM/MM for all single amino acid point-charge variations are shown in Fig. 4A, in which the amino acids are sorted according to their center of mass (COM) distance from the substrates in the active site of the reactant structure (QM region 1).</p><p>Naively, it could be expected that residues closer to the active site show higher sensitivities to point charge variations than distant ones. However, there are also several high-sensitivity amino acids at medium distances to the active site and low-sensitivity residues very close to the substrates. This observation confirms that an exclusively distance-based approach to include residues into the QM region is not able to detect all important amino acids and furthermore includes residues which are probably not necessary to obtain consistent QM regions. To find a compromise between including amino acids that show a high sensitivity and those that are close to the active site, we define an empirical QM region indicator Θ i for each amino acid by dividing the sensitivity by the COM distance between the amino acid and the active site, i.e.,</p><p>This definition ensures that distant residues with high sensitivities are considered, but also residues close to the active site which may show medium sensitivities are not overlooked.</p><p>The resulting indicator is plotted in Fig. 4B. A comparison of the indicators for the schemes PCVA-A to PCVA-D is shown in the Supporting Information (Fig. S4), and an additional Table S3 lists the 16 amino acid residues with the highest QM region indicators for the schemes PCVA-A to PCVA-E (see also Section 7). When using the QM-only instead of the full QM/MM reaction energy (i.e., comparing PCVA-A and PCVA-B), the only change among the highest-ranked amino acids is that ILE88 is replaced by CYS68. Similarly, only two changes among the 16 amino acids with the highest indicators are found when using the reactant instead of the product energies (i.e., comparing PCVA-C and PCVA-D). Here, LEU64 and TRP142 are replaced by VAL41 and CYS172. All these amino acids show very similar sensitivities and indicators in the different schemes. Larger differences are found going form the reaction energies (PCVA-B) to using the reactant energies (PCVA-D), where six differences appear among the top-16 amino acids. However, these changes mostly occur for amino acids with very close values of the QM region indicator. Furthermore, the use of the reactant or product energies avoids error cancelation that might be present when using the reaction energy. Finally, only two differences (VAL41 and TYR67 instead of LEU64 and TYR146) are found among the 16 highest-ranked amino acids when going to a simplified numerical differentiation formula (i.e., from PCVA-D to PCVA-E). Overall, this confirms that the PCVA-E approximation seems to be a reasonable choice.</p><p>By including the amino acids with the highest QM region indicators Θ i (see Fig. 5), we are now able to systematically construct QM regions that should reduce the sensitivities of the QM/MM reaction energy to point charge variations.</p><p>6 Assessment of PCVA-Based QM Regions of Increasing</p><!><p>As a first test, we assess the PCVA-based construction of QM regions with increasing size.</p><p>Ideally, by using PCVA for a systematic construction of QM regions, the convergence towards the results obtained with very large QM regions should be accelerated compared to an exclusively distance-based construction of the QM region. To this end, we construct QM regions consisting of just as many residues as in the distance-based approach (see Section 4)</p><p>and label these regions as 2', 3' etc. (e.g., QM regions 3 and 3' both contain seven amino acid residues). Fig. 5 and Tab. S4 in the Supporting Information show which residues are included for each QM region. Again, geometry optimizations of the reactant and product structures were performed for each of these QM regions.</p><p>Fig. 6 shows the convergence of the VDD charges and the QM/MM reaction energy as well as the corresponding sensitivities to global point charge variations (cf. Fig. 2 for the same plot PCVA-constructed QM regions deliver overall lower sensitivities regarding VDD ligand charges (see Fig. 6B), which marks an expected behavior because high-sensitivity residues are included in the QM region by PCVA. The convergence of the charge sensitivities is similar to the analysis for distance-based QM regions (cf. Fig. 2) with a fast convergence for Mg 2+ , a constant behavior for CAT, and jumps for SAM with a slightly decreasing trend.</p><p>The reaction energy (see Fig. 6C) decreases for the smallest QM regions, starts to stabilize for QM regions 4', and oscillates around our best estimate of about −11 kcal/mol for larger QM regions. However, these oscillations are smaller than for the distance-based construction of the QM regions. Note that compared to Fig. 2, we updated our best estimate to correspond to QM region 9' in Fig. 6C.</p><p>In contrast to the distance-based case, a slightly increasing trend is observed for the reaction energy sensitivity (Fig. 6D). For QM region 8', at which the reaction energy is further from our best estimate, the sensitivity also shows a marked increase. Except for the smallest QM regions, the global PCVA sensitivity can thus be used to judge the accuracy of the calculated QM/MM reaction energy. The overall larger sensitivity for the PCVAbased QM regions compared to the distance-based construction could be an effect of a larger contact area between the MM and QM regions because in the PCVA case, MM residues close to the active site with low sensitivities remain in the MM part and can affect adjacent QM residues.</p><p>7 Assessment of Atom-Economical QM Region Based on</p><!><p>We now consider the determination of an atom-economical QM region consisting of 16 amino acids for COMT based on PCVA. This approach follows the work of Kulik et al., who constructed such QM regions using their CSA and FSA approach, 15,19 and thus allows for a direct comparison to these approaches. To this end, we included the 16 amino acids with the highest QM region indicators Θ i in the QM region (see Fig. 5), and again performed QM/MM geometry optimizations of the reactant and product structures. Similarly, we performed QM/MM geometry optimizations using the QM regions obtained in Refs. 15,19   which also each include 16 amino acids.</p><p>Table 1 compares the amino acid residues included in such an atom-economical QM region for a distance-based QM region construction, the CSA and FSA approaches, and for our PCVA-based QM-region construction.</p><p>Compared to the CSA and FSA approaches applied by Kulik et al. 15,19 our PCVA approach detects the majority of the residues that should be part of the 16-residue QM region, as it can be seen in Table 1 (e.g., VAL41, GLU89, or ASP140). Moreover, also more distant residues are included that would not be considered in a 16-residue QM region in the distance-based case (e.g., MET39, LYS45, ALA66, or ASP168). This indicates the ability of PCVA to even detect high-impact residues located relatively distant from the active site COM. Nevertheless, there are several residues included by CSA and/or FSA that are not under the 16 highest-ranked residues when applying the PCVA approach. However, most of these residues are assigned a PCVA rank close to 16 such as ASN40, GLU63, or ILE90.</p><p>An extreme case with a very low PCVA rank (178) compared to CSA (10) and FSA (10-11) is SER71. This residue is located very close to the active site and thus potentially important for ligand binding, but it shows a very low electrostatic effect on the substrates and is thus not detected by PCVA. Another case is SER118, which is ranked about 30 places Furthermore, even between CSA and FSA major differences in the evaluation of residues can be observed (e.g. for GLU63, GLY65, TYR70, or ALA72) leading to differently composed QM regions.</p><p>In Fig. 6, we compare the VDD charges of SAM, CAT, and the catalytically-active Mg 2+ cation as well as the QM/MM reaction energy obtained for the atom-economical QM regions constructed using CSA, FSA, and PCVA to those obtained for PCVA-constructed QM regions of increasing size discussed in Section 6. Note that the size of the atom-economical QM regions of 16 amino acids is in between those of QM regions 4' (13 amino acids) and 5'</p><p>(19 amino acids).</p><p>Regarding VDD charges (see Fig. 6A), PCVA performs as well as CSA and FSA. The change in the Mg 2+ charge to about 0.3 with the distance-based region 8 (34 residues) and larger is achieved already in the QM regions containing 16 residues. For SAM, all three atom-economical QM regions also provide the converged charge. In the SAM case, none of the methods is able to detect the charge change to 0.25 for region 7 and larger. This is a fundamental limitation of QM regions of a given size and can probably only be rectified by using larger QM regions.</p><p>Concerning the reaction energy (Fig. 6B</p><!><p>As a first step towards systematically quantifying the uncertainties of QM/MM calculations, we have presented a point charge variation analysis for assessing the sensitivity of QM/MM reaction energies to changes of the MM point charges. To this end, we considered the derivative of the QM/MM reaction energy with respect to selected distortions of the MM charges.</p><p>This derivative can be calculated numerically by performing QM/MM calculations with varied MM point charges, and different efficient approximations can be employed. Generally, the most simple approximation (PCVA-E), in which a forward finite difference is used and only the reactant structure is considered, turns out to be sufficient for a qualitative assessment of uncertainties.</p><p>A global PCVA, in which all protein point charges are varied simultaneously, can be used as a simple indicator of the accuracy of the resulting QM/MM reaction energy as well as other properties of the active site, such as ligand charges. For the considered test cases we found that when comparing different QM regions, the sensitivity in a global PCVA is smaller for those QM regions that yield results closer to the best estimate obtained for the largest QM regions, i.e., the sensitivity generally decreases when approaching the limit of a full QM calculation. However, this only holds once the QM regions have reached a reasonable size and the correlation between the sensitivity in a global PCVA and the deviation from the best estimate is not always clear.</p><p>Nevertheless, our results demonstrate that the analysis of the sensitivity with respect to the MM point charges is a good starting point for the investigation of uncertainties in QM/MM calculations. We are planning to extent this work in the future by considering not only selected distortions of the MM point charges, but performing a full sensitivity analysis that allows one to identify the collective point charge variations that have the largest influence on QM/MM reaction energies, following our earlier work on structural distortions in theoretical spectroscopy. 40 While the reliability of considering only one global distortion of the MM point charges is limited, such a more comprehensive assessment of the sensitivity can be expected to overcome the limitations of our present approach.</p><p>In addition to a global PCVA, we have considered a single amino acid PCVA, in which the MM charges of each amino acid residue are varied. This makes it possible to assess the contribution of each amino acid to the uncertainty in the QM/MM reaction energy, and can be used to guide the systematic construction of the QM region. By including amino acids with a high sensitivity to point charge variations in the QM region, the overall sensitivity of the QM/MM reaction energy can be reduced. Here, we devised a PCVA-based scheme for the systematic construction of the QM region.</p><p>For the considered test case, our scheme leads to a faster and more reliable convergence with the size of the QM region compared to distance-based QM region construction. Comparing to atom-economical QM regions of the same size provided by the other common approaches, in particular CSA and FSA, our PCVA-based approach performed well and yields similar QM regions.</p><p>The huge advantage of PCVA is its much lower computational cost compared to the CSA or FSA approach (see Supporting Information, Tab. S5). Our PCVA-based approach requires only a geometry optimization of the target system with a minimal QM region including substrates, which is followed by single-point calculations for the point-charge variation of each amino acid. In contrast, CSA is a very expensive approach based on large QM regions with up to 1000 atoms. For these large systems, geometry optimizations have to be performed</p>