MM vs Pure MM Molecular Dynamics for Evaluating Water

Apr 16, 2019 - Hybrid QM/MM vs Pure MM Molecular Dynamics for Evaluating Water Distribution within p21N-ras and the Resulting GTP Electronic Density...
0 downloads 0 Views 4MB Size
Subscriber access provided by UNIV OF LOUISIANA

B: Biophysics; Physical Chemistry of Biological Systems and Biomolecules

Hybrid QM/MM vs Pure MM Molecular Dynamics for Evaluating Water Distribution within p21 and the Resulting GTP Electronic Density N-ras

Ruth Helena Tichauer, Gilles Favre, Stéphanie Cabantous, and Marie Brut J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.9b02660 • Publication Date (Web): 16 Apr 2019 Downloaded from http://pubs.acs.org on April 19, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Hybrid QM/MM vs Pure MM Molecular Dynamics for Evaluating Water Distribution within p21N −ras and the Resulting GTP Electronic Density Ruth H. Tichauer,∗,† Gilles Favre,‡ Stéphanie Cabantous,‡ and Marie Brut∗,† †LAAS-CNRS, Université de Toulouse, CNRS, UPS, Toulouse, France ‡Cancer Research Center of Toulouse, INSERM U1037, 31037 Toulouse, France, Université de Toulouse, Toulouse, France E-mail: [email protected]; [email protected] Phone: +33 (0)5 61 33 63 04. Fax: +33 (0)5 61 33 62 08

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract p21ras proteins activity, regulated by GTP hydrolysis, constitutes an active field of research for the development of cancer targeted therapies that would concern ∼ 30% of human tumors to which specific mutations have been associated. Indeed, the catalyzing mechanisms provided by the protein environment during GTP hydrolysis, and how they are impaired by specific mutations, remain to be fully elucidated. In this article, we present results from Molecular Mechanics (MM) Molecular Dynamics (MD) simulations and Density Functional Theory (DFT) calculations carried out for wild-type p21N −ras and six Gln 61 mutants. In a first part, we present water distribution within the active site of the wild-type protein according to MM MD. Significant differences are observed when comparing the results to the previous distribution assessed through Quantum Mechanics/Molecular Mechanics (QM/MM) MD. Such method-dependent results highlight the importance of accounting for the electrostatic coupling between the protein complex and the solvent molecules to identify hydration sites. In a second part, we present results from DFT calculations performed to determine the electronic distribution of the GTP ligand, considering the wild-type active site arrangement according to both classical and hybrid approaches. Only in the QM/MM based configuration, the ligand electronic density is similar to that of a GDP-like state observed experimentally. For this reason, in a last set of calculations carried out for p21N −ras Gln 61 mutants, the active site structural conformations obtained through hybrid MD are only considered. Through the analysis of GTP electronic density, we conclude that the wild-type active site arrangement according to QM/MM MD is closer to a catalytically efficient conformation of the protein than the arrangement according to MM MD. Hence, water distribution according to the hybrid approach must correspond to the optimal placement of solvent in the active site. Within all the studied Gln 61 substituted proteins, p21ras major catalyzing effect, which consists in stabilizing a more GDP-like state, is lost.

2

ACS Paragon Plus Environment

Page 2 of 35

Page 3 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

1

Introduction

p21ras (referred to as Ras in this article) is a small GTPase protein known for its critical role in signal transduction. 1 Through GTP hydrolysis, 2,3 it switches between an active GTP bound state, which enables interactions with downstream effectors, 4 and an inactive GDP bound state. Furthermore, Ras activity is downregulated by the binding of GTPase Activating Proteins (GAPs) 5–9 that increase the hydrolysis rate by up to 5 orders of magnitude through steric 10–14 and electrostatic 12–17 effects. Ras mutations at sites 12, 13 and 61 are found in approximatively 30% of human tumors. 18 All of them induce a decrease in hydrolysis rate, 19,20 leaving Ras constitutively active. The elucidation of Ras mechanisms for catalyzing GTP hydrolysis remains consequently a milestone in the development of cancer targeted therapies. Diverse reaction pathways have been proposed by experimental 21–30 and theoretical 12–14,31–40 studies carried out in the last 40 years. In particular, distinct roles have been attributed to Gln 61 residue as its substitution results in a defective intrinsic 41 and GAP enhanced 5 hydrolysis, suggesting that it plays a major role in the reaction. Among these assumptions, Gln 61 residue could i) participate in proton transfers, 21,34,35,38–40 ii) stabilize the transition state, 11,22–24,31 iii) precisely position relevant water molecules 13,26,27,30,32,33,42 or iv) have an indirect steric effect in stabilizing the preorganized catalytic configuration of Ras active site. 10,12,14,41,43 In our previous study 43 of p21N −ras (referred to as NRas), we used combined quantum mechanics / molecular mechanics (QM/MM) molecular dynamics (MD) to study the GTP bound conformation of wild-type (WT) NRas and six oncogenic mutant proteins reported in the Catalogue Of Somatic Mutations In Cancer (COSMIC) 44 (i.e. Q61E, Q61P, Q61H, Q61L, Q61K and Q61R, in increasing order of prevalence). We described the structural rearrangements undergone by the protein active site upon Gln 61 substitutions and the resulting water distribution. A mapping of water presence probability density allowed to identify favorable and unfavorable hydration sites and to determine the number of water 3

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

molecules present in each complex active site. In the present study, first, differences in structural elements and water distribution are evaluated when employing a pure molecular mechanics (MM) description. The comparison between MM and QM/MM results is done for WT and Q61R NRas, since they were the most dissimilar in the QM/MM study. This comparison shows that accounting for the electrostatic coupling between the protein complex and the solvent molecules is crucial for determining water distribution within the active site, a prerequisite for hydrolysis. Subsequently, we estimate GTP electronic density through DFT/MM calculations in which this ligand, the Mg2+ ion and water molecules present in the active site are treated at the quantum level. Beforehand, the protein active site conformation and water positioning according to both MM and QM/MM results are considered for WT NRas. These calculations reveal that GTP electronic distribution induced by the wild-type QM/MM structure is the only one that can be correlated to experimental observations 25,26,28,29 where Ras stabilizes a GDP-like state through electrostatic effects. On this basis, for Gln 61 mutant proteins, only the active site conformations extracted from QM/MM dynamics are considered. We notably show that GTP electronic density undergoes drastic modifications upon Gln 61 substitutions. This suggests that i) the QM/MM conformation of the WT active site together with related water distribution is closer to a preorganized catalytic state, which is necessary for hydrolysis to occur, than the conformation and water distribution resulting from MM MD and ii) the GDP-like induced distribution of GTP electronic density is impaired by Gln 61 substitutions. As a consequence, their reduced hydrolysis rate could be explained, at least in part, by the loss of this major catalyzing effect.

4

ACS Paragon Plus Environment

Page 4 of 35

Page 5 of 35

2

Methods

2.1

MM Molecular Dynamics

The initial structures of GTP-bound WT and Q61R NRas (residues 1-166), in complex with GAP arginine binding loop and the cofactors, were obtained from PDB entries 3CON 45 (GDP-bound NRas), 1QRA 27 (GTP-bound HRas) and 1WQ1 11 (HRas-GDP-AlF3 -GAP complex) as previously described. 43 The nucleotide binding pocket 46 is depicted in figure 1(a). From NRas active site, where hydrolysis takes place, constituted of the P-loop (residues 10-17), the Switch I (residues 30-40) and the Switch II (residues 57-67), only residues that were treated quantum mechanically in our previous work 43 are depicted in figure 1(b). Motions undergone by residues belonging to the highly flexible Switch I and Switch II regions during MM MD are described in the Results section and confronted to previous reports. 15,43,47 ``` ```

Gly 13 Gly 12 Gln 61 Pβ Pα



Tyr 32 c c c

(a)

c c c

Switch I

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

P-loop

Switch II

Mg2+

(b)

Gly 60

Thr 35

Figure 1: View of (a) the G domain, a zoom is made on NRas regions that interact with the phosphate moiety of the ligand, and (b) the positioning of NRas residues 12-13, 32, 35, 60-61, GTP, M g 2+ ion and water molecules. All MM MD simulations were carried out using AMBER package, 48 using the ff99bsc0 force field. 49 WT and Q61R NRas proteins were studied in explicit solvent, using the TIP3P water model 50 under periodic boundary conditions, with neutralizing Na+ ions. To get rid of possible vacuum bubbles created during the solvation process, 10000 cycles of energy 5

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

minimization were first performed with a 500 kcal·mol−1 · Å−2 restraint applied to the solute. Subsequently, the entire system was gently heated to 300K using the Langevin thermostat with a collision frequency of 2 ps, at constant volume, with weak restraints of 5 kcal · mol−1 · Å−2 applied to the entire protein complex. Finally, the production set-up was implemented in the NPT ensemble with 1 atm target pressure, using the Berendsen barostat with a pressure relaxation time of 1 ps. Long-range electrostatic interactions were treated using the particlemesh Ewald method with a 12.0 Å cutoff. Covalent bonds involving hydrogen atoms were constrained using the shake algorithm 51 so that a 2 fs step could be used. The obtained trajectories are 10 ns long, that is 10 times longer than those of the previous QM/MM study, 43 enabling a more extended sampling of conformational space. System properties (energy, density, temperature, pressure, volume, backbone root mean square deviation (RMSD)) were monitored to assess their quality. The corresponding plots can be found in figures S1 and S2 of the supporting information (SI).

2.2

DFT Calculations

The starting structures of WT and Gln 61 mutant NRas were extracted from the previous QM/MM trajectories. 43 They correspond to the most representative conformation of each NRas-GTP complex, identified by a clustering algorithm 52 implemented in the AMBER package. These protein-ligand complexes together with the magnesium ion, its two coordinated water molecules and any water molecule present within 5Å of GTP Pγ atom were subsequently energy minimized in implicit solvent, using the generalized Born model 53 until the convergence criterion was reached (root mean square gradient < 10−6 kcal · mol−1 ·Å−1 ). Following this procedure, water molecules occupy positions that correspond to high water density regions identified through the 2D mapping performed in the QM/MM study (see figure 2). DFT calculations of NRas reaction center were carried out using the Orca package 54,55 interfaced to Amber, 56 so that the protein environnement together with the solvent effects 6

ACS Paragon Plus Environment

Page 6 of 35

Page 7 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Water molecules coordinated to the Mg2+ ion appear as gray sticks while transient water molecules are depicted in a colored ball-andstick representation.

(a)WTQM/M M

(b)WTM M

(c)Q61E

(d)Q61P

(e)Q61H

(f )Q61L

(g)Q61K

(h)Q61R

Figure 2: Active site conformation of (a) WTQM/M M , (b) WTM M , (c) Q61E, (d) Q61P, (e) Q61H, (f) Q61L, (g) Q61K and (h) Q61R NRas which properties are further investigated through DFT/MM calculations. could be taken into account in a DFT/MM MD scheme. 57 The quantum portion consisted in GTP, water molecules placed within 5Å of its Pγ atom, Mg2+ ion and its two coordinated water molecules. They were described at the DFT level using the Generalized Gradient Approximation (GGA) functional PBE and the Def2-SVP 58–60 basis set in single point energy calculations performed during short DFT/MM MD simulations. Indeed, the primary purpose is to get insight into the electronic distribution of GTP, bound to WT and mutant NRas, when the active site is in the most representative QM/MM conformation.

7

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

3

Results

For comparison purposes, the trajectories obtained from MM MD were analyzed as those obtained from the previous QM/MM simulations, that is to say i) the intrinsic mobility of individual NRas residues was evaluated by plotting their side chain RMSD against the simulated time, ii) possible interaction partners of GAP Arg 789 residue were assessed by computing its native/non-native contacts and iii) water presence in the active site was evaluated by computing the Radial Distribution Function (RDF) and the 2D RDF of water molecules around GTP Pγ and Pα atoms. This enabled the identification of NRas properties that require a quantum treatment to be accurately described. DFT/MM calculations were carried out to assess the impact on GTP electronic density of the protein active site arrangement together with water placement in it. First, this influence was evaluated for WT NRas active site conformation and water distribution determined from MM and QM/MM MD, revealing that the QM/MM based arrangement is the only one that induces an electronic distribution of GTP characteristic of a precatalytic conformation that is necessary for hydrolysis to occur. In a second place, the effect of Gln 61 substitutions was evaluated considering the active site structure from the hybrid approach only.

3.1 3.1.1

MM vs QM/MM Molecular Dynamics Active Site Stability

Details concerning the structural properties of WT and Q61R NRas active site during MM MD can be found in the supporting information (SI). Overall, it appears that, during the simulated time, an MM description does not enable to observe certain structural changes induced by Q61R mutation observed when employing a QM/MM approach, namely i) modifications in Gly 60 residue interaction network, ii) Thr 35 residue destabilization and iii) GAP Arg 789 residue loss of contact with GTP Pγ atom. Such structural differences between WT and Q61R NRas could be expected at least because

8

ACS Paragon Plus Environment

Page 8 of 35

Page 9 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

of the additional positive charge introduced by the mutation.

3.1.2

Water Occupancy within the Active Site

As summarized above, classical MD does not enable certain structural changes of NRas active site to occur upon Gln 61 substitution. Along this line, we investigated if a classical treatment allows water molecules to explore the same hydration sites as identified during the QM/MM study. Details can be found in the SI.

RDF Concerning RDF of water molecules around GTP Pγ and Pα atoms, similar conclusions are reached when employing a pure MM description as when adopting a more demanding QM/MM approach. Indeed, for both WT and Q61R NRas, a peak of water density arises close to GTP Pγ atom (see figure S7 in the SI). Yet, in a 5Å radius from this atom, as concluded in the previous QM/MM study, there are less water molecules within WT than within Q61R NRas (1.27 vs 2.32 water molecules, respectively, according to table S1 in the SI). Moreover, a more apparent difference in water presence between both proteins is observed in the RDF around GTP Pα atom as a high amplitude peak arises near this atom for Q61R only (see figure S8 in the SI), alike the corresponding curves computed from QM/MM dynamics. The integration of this RDF curves shows that, as reported in the previous QM/MM work, there are much less water molecules in a 7Å radius from GTP Pα atom for WT than for Q61R NRas (1.10 vs 6.64 water molecules, respectively, as indicated in table S1 in the SI) . Nevertheless, a first method-dependent difference in the RDF curves is the more structured water densities observed during MM dynamics than what was observed when using a QM/MM approach. Indeed, as depicted in figures S7 and S8 in the SI, a second peak arises in both WT and Q61R NRas Pγ RDF plots (and a third as well for Q61R only), such that this second peak is also clearly defined when plotting the curve for a MM simulated time equal to that of the QM/MM study (i.e. after 1 ns, see figure S10 in the SI).

9

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

A second method-dependent difference in the RDF curves arises when comparing MM and QM/MM RDF curves of Q61R NRas. It appears that when resorting to a pure MM description rather than a hybrid approach, the first peak observed in the vicinity of both GTP Pα and Pγ atoms is shifted in the same direction, by nearly the same amount. Since the positioning of these peaks is systematically downshifted according to MM dynamics, water molecules can get closer to GTP phosphorus atoms according to this level of theory than predicted by a QM/MM approach. A thorough investigation shows that this difference is firstly due to P − O bond lengths within the α and γ phosphate groups that are longer according to the semi-empirical Hamiltonian used during QM/MM dynamics (PM3) (see table S2 in the SI). Clearly, during MM MD, their length is dictated by the GTP parameters 61 (as their values are almost identical within WT and Q61R NRas) while they elongate under the influence of the surrounding environment when using a quantum description. It is important to mention that the longer bonds are in accordance with a Raman spectroscopy study 26 that showed that upon binding to Ras, the P − O bonds of GTP phosphate groups are loosened by interactions with the protein and the Mg2+ ion. To assess if the treatment of the solvent at either the MM or QM level also contributes to the previously discussed shift between MM and QM/MM RDF plots, this quantity was evaluated around GTP γ oxygen atoms within WT and Q61R NRas (see figure S9 in the SI). Because both approaches (MM and QM/MM) yield very dissimilar water distributions within NRas active site (see 2D RDF plots in figure S11 in the SI), it is difficult to assess if the differences observed in the RDF curves around GTP Oγ atoms are purely methoddependent. MM and QM/MM calculations of GTP in solution might shed light on this matter but this is beyond the scope of the present study. We can only conclude that, the method-dependent positioning of the first peak in Q61R RDF curves around Pα and Pγ atoms is, at least in part, due to the different P − O bond lengths during MM and QM/MM dynamics. It remains to be elucidated if treating the solvent at different levels of theory i) contributes to the shift, ii) reduces it or iii) has no influence on it i.e. if the solute-solvent

10

ACS Paragon Plus Environment

Page 10 of 35

Page 11 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

average minimum distance when treating water molecules at the QM level is the same than the one dictated by the Lennard-Jones terms when employing MM parameters only.

2D RDF Concerning the precise mapping of water molecules in NRas active site, a pure MM description of solute-solvent interactions leads to noticeable differences in water distribution between WT and Q61R i.e. the number and positioning of water density peaks (see figure S11 and related description in the SI). Nevertheless, what appeared to be characteristic of each NRas form according to QM/MM results, is not revealed by a classical approach. In particular, for WT NRas, residues 12 and 60 do not participate in holding transient water molecules in a precise region of the active site that, during the QM/MM study, lead to the formation of the highest water density peak in this region. For Q61R, water molecules delocalization is lost. Furthermore, the 2D RDF plots clearly show how pure MM dynamics do not allow structural rearrangements of the active site to take place upon Gln 61 substitution. Indeed, the projections of the Mg2+ ion and nitrogen backbone atoms of residues 12, 13, 35, 59, 60 and 61 are similar for WT and Q61R NRas. During QM/MM dynamics, these projections were remarkably different both in average positioning and standard deviation. The most noticeable difference was observed for Thr 35 nitrogen backbone atom which, within Q61R NRas, was displaced beyond the threshold distance, so its projection could not be observed. Structural reorganizations of the active site could be expected in the case of Q61R, at least because of the positive charge introduced by such mutation. In order to further characterize the impact of Gln 61 substitutions on the protein active site structure together with water positioning in it, we proceed to examine GTP electronic density. Indeed, the differences observed in the 2D RDF plots of WT and Gln 61 mutant NRas during QM/MM dynamics, gave rise to the hypothesis that water distribution has a direct impact on the electronic density of the bound ligand. Hence, in the next section, first GTP electronic density considering WT NRas active site structure and water distribution

11

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

inferred from both MM and QM/MM MD is evaluated. Subsequently, the impact of Gln 61 mutations, considering the active site conformations and water positioning obtained from QM/MM calculations only, is assessed.

3.2

GTP Electronic Density

Regardless of the numerous reaction schemes that have been proposed for GTP hydrolysis within Ras proteins, it has been shown 25,62 that its binding to the enzyme results in the elongation of the bond to be broken (between Pγ and O3β atoms). In this manner, the associated energy barrier is lowered by several kcal/mol compared to the same reaction in solution. Among other attributed catalyzing effects, several theoretical and experimental studies have shown a strong binding between the protein and GTP β phosphate group, 26 leading to the stabilization of a GDP-like state. Furthermore, in a theoretical study of the enzymatic reaction catalyzed by GAP, Rudack et al. 17 have proposed that the role of the Mg2+ ion is to temporarily store electrons from GTP γ phosphate group. This process would render it more positive while the so-called GAP arginine finger would render the β phosphate group more negative . To determine GTP electronic distribution when bound to WT and Gln 61 mutant NRas, DFT/MM MD calculations of the protein-ligand complexes were performed. The individual contributions of each atom to the overall electron density were then assessed through a Löwdin population analysis 63,64 that was carried out for the portion treated at the DFT level of theory. The reduced atomic charges in this manner attributed1 are only considered for comparing the different environments and, thus, the relative changes they induce on electron distribution. The Kohn-Sham Highest Occupied Molecular Orbital (KS HOMO) for each protein complex was also investigated to assess the distribution of the highest energy 1

Charges held by individual atoms (or groups of atoms) reported all along this section, refer to the reduced Löwdin charges, i.e. the amount of electron density ρ each atom A holds, noted ρA . This quantity needs to be distinguished from the partial charge of the considered atom (or group of atoms) which is usually defined as δqA = ZA − ρA where ZA is the atomic number. In this manner, when ρA < 0 (respectively ρA > 0), δqA is larger (respectively smaller) than the charge held by isolated free atom A

12

ACS Paragon Plus Environment

Page 12 of 35

Page 13 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

electrons which are the most likely to react. Results from both analyses are depicted in figures 3 and 4 for WT and mutant NRas, respectively. 3.2.1

Wild-type NRas Wild-type structure from MM vs QM/MM H H O H H O H H O H H O O2’ O2’ H2’ C2’ H3’ O3’ O3’ C3’ C4 N3 H22 H21 N2 C2 H1 N1 O6 C6 C5 N7 H8 C8 N9 H1’ C1’ O4’ H4’ C4’ 5’’ H5’ C5’ O5’ O2a O1a Pa O3a O2b O1b Pb O3b O3g O2g Pg O1g Mg

WT with 2 H2O (QM/MM) WT with 1 H2O (MM)

(a)

-0.8

(b)

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

Lowdin Reduced Atomic Charges

Figure 3: Löwdin population analysis of GTP, M g 2+ and water molecules within WT NRas active site. The conformation of this site together with water molecules placement is that inferred from QM/MM (black rectangles) and MM (red rectangles) MD. In the insets are shown the corresponding KS HOMOs i.e. according to (a) QM/MM and (b) MM MD. The isovalue selected for representing the delimiting surface of the KS HOMO is |0.02|e− a−3 0 . The positive lobes are depicted in blue and their negative counterparts in red.

Löwdin Population Analysis WT active site structural organization along with water molecules inside it, inferred from MM (1 water) and QM/MM (2 waters) dynamics, culminate in significant differences of the reduced charges2 assigned to GTP atoms and Mg2+ ion (see figure 3). The latter holds a bigger portion of the overall electron density ρ within the active site conformation deduced from the QM/MM study (referred to as WTQM/M M ) than within the conformation obtained from the present MM results (referred to as WTM M ). Indeed, its Löwdin reduced charge ρM g2+ is more negative within WTQM/M M than within WTM M : -0.33 a.u. within WTQM/M M vs -0.18 a.u. within WTM M . 2

determined from a Löwdin population analysis

13

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

GTP γ phosphate group is less negatively charged (i.e. more positive) within WTQM/M M than within WTM M . Its Pγ atom holds a more positive charge (0.86 a.u. vs 0.81 a.u respectively) and two among the three Oγ atoms are less negatively charged. In particular, O3γ atom, the one coordinated to Lys 16 residue, holds 0.24 a.u. less negative charge within WTQM/M M than within WTM M . The bridging oxygen atom between GTP γ and β phosphate groups has slightly the same reduced charge in both environments (∼ 0.002 a.u. charge difference). In the same manner but to a lesser extent, GTP β phosphate group is less negatively charged within WTQM/M M than within WTM M . Its Pβ atom holds a slightly less positive charge (0.82 a.u. vs 0.83 a.u respectively) and both Oβ atoms are less negatively charged by 0.08 and 0.02 a.u. for O1β and O2β , respectively. The bridging oxygen atom between the β and α phosphate groups holds a slightly smaller amount of negative charge within WTQM/M M than within WTM M (-0.41 a.u. vs -0.43 a.u respectively). Finally, the phosphorous atom belonging to the α phosphate group has about the same reduced charge in both environments (∼ 0.003 a.u. of charge difference). O1α and O2α oxygen atoms are respectively more and less negatively charged within WTQM/M M than within WTM M . The bridging oxygen atom between the α phosphate group and the guanosine moiety, O5 , holds a slightly larger amount of negative charge within WTQM/M M than within WTM M (-0.31 a.u. vs -0.29 a.u respectively). Overall, the β phosphate group holds a more negative charge than the γ phosphate group within WTQM/M M (-0.57 a.u. vs -0.50 a.u. respectively) while the opposite is observed within WTM M NRas (-0.66 a.u. vs -0.81 for GTP β and γ phosphate groups, respectively). GTP electronic distribution within WTQM/M M is hence in agreement with numerous experimental and theoretical investigations 25,26,28,29,62 concluding that Ras catalyzes GTP hydrolysis by stabilizing a more GDP-like state. Rudack et al. 17 showed that this is achieved by strong

14

ACS Paragon Plus Environment

Page 14 of 35

Page 15 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

interactions that lead to charge shifts that render the γ phosphate group more positive while the β phosphate group is rendered more negative, compared to GTP in solution. We presume that GTP electronic distribution within WTM M does not correspond to that of a GDP-like state due to the absence of the water molecule that is held by nitrogen backbone atoms of residues 12 and 60 within WTQM/M M . We propose that this water molecule has a structuring role. Its absence from this precise location in WTM M active site disables the finely tuned interactions between NRas and its ligand. As a consequence, the charge shifts that would be necessary to accommodate GTP electronic density to a more GDP-like state cannot take place. In such a scenario, the presence of two water molecules in the active site appears to be crucial. This is in accordance with a recent study of HRas 42 that supports a two water molecules scheme although with an albeit different role and positioning for the second water molecule than what we propose here.

HOMO The previous Löwdin population analysis aimed at evaluating individual atoms electronic contribution to the overall electron density ρ of the DFT system they form. In this part, we study the KS HOMO to assess how the highest energy electrons are distributed as they are the most likely to react. For WTQM/M M NRas, electronic density from the KS HOMO appears prominently on GTP α oxygen atoms, followed by β oxygen atoms and finally γ oxygen atoms. According to the chosen isovalue (|0.02|e− /a30 ) for representing the delimiting surface of the KS HOMO, the highest energy electrons are distributed following the same trend observed when attributing partial charges to GTP phosphate groups by a Löwdin population analysis, i.e.:

αP O− electronic density > βP O− electronic density > γP O2− electronic density 4

4

4

The agreement between these two analyses suggests that the KS HOMO makes large contributions when assigning reduced charges to GTP phosphate groups following the Löwdin scheme, although electronic distribution on γ oxygen atoms does not follow the same trend according to both approaches. Indeed, O1γ atom holds no electronic density according to 15

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

the KS HOMO analysis while it is the most negatively charged γ oxygen atom according to the Löwdin scheme. For WTM M NRas, using the same isovalue, electronic density from the KS HOMO appears mainly concentrated on the presumed lytic water molecule placed in the vicinity of Thr 35 residue. GTP γ oxygen atoms are the only ones that also share this electronic density. O3γ atom appears as the most negatively charged among the three according to both approaches, while O1γ and O2γ are the least negative according to the KS HOMO and the Löwdin population analysis, respectively. The γ phosphate group holds a larger amount of negative charge than the β phosphate group according to both approaches, and than the α phosphate group according to the KS HOMO analysis only. The highest energy electron distribution within WTM M suggests that, when attributing Löwdin partial charges, other lower energy occupied MO contribute more significantly than the depicted KS HOMO.

3.2.2

Gln 61 Mutant Proteins

Löwdin Population Analysis According to a Löwdin scheme for attributing reduced atomic charges, NRas active site structural rearrangements together with the corresponding modifications on water distribution, induced by Gln 61 substitutions, impact the reduced charges held by GTP atoms (see figure 4). Some common features arise for the six studied mutant proteins, namely: i) the Mg2+ ion is more negatively charged within WTQM/M M than within any of the mutant proteins, in particular within Q61P, the reduced Löwdin charge attributed to the cation is positive, ii) GTP Pγ atom is more positively charged within WTQM/M M than within any of the mutant proteins and iii) O3γ and O1β atoms hold less negative charge within WTQM/M M than within any of the mutant proteins (even if for Q61R the charge difference on O1β is small, about ∼ 0.005 a.u.). It is also observed that O1γ and O2γ atoms hold the same negative Löwdin reduced charge within Q61E and Q61P mutants as within WTQM/M M NRas. They are less negatively

16

ACS Paragon Plus Environment

Page 16 of 35

Page 17 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(a) WT vs Q61E H H O H H O H H O H H O O2’ O2’ H2’ C2’ H3’ O3’ O3’ C3’ C4 N3 H22 H21 N2 C2 H1 N1 O6 C6 C5 N7 H8 C8 N9 H1’ C1’ O4’ H4’ C4’ 5’’ H5’ C5’ O5’ O2a O1a Pa O3a O2b O1b Pb O3b O3g O2g Pg O1g Mg

(b) WT vs Q61P H H O H H O H H O H H O O2’ O2’ H2’ C2’ H3’ O3’ O3’ C3’ C4 N3 H22 H21 N2 C2 H1 N1 O6 C6 C5 N7 H8 C8 N9 H1’ C1’ O4’ H4’ C4’ 5’’ H5’ C5’ O5’ O2a O1a Pa O3a O2b O1b Pb O3b O3g O2g Pg O1g Mg

Wild-type Q61E

(i)

-0.8

-0.6

(ii)

-0.4

-0.2

0

0.2

0.4

0.6

0.8

Wild-type Q61P

(i)

-0.8

-0.6

(ii)

-0.4

Lowdin Reduced Atomic Charges

-0.2

(c) WT vs Q61H H H O H H O H H O H H O O2’ O2’ H2’ C2’ H3’ O3’ O3’ C3’ C4 N3 H22 H21 N2 C2 H1 N1 O6 C6 C5 N7 H8 C8 N9 H1’ C1’ O4’ H4’ C4’ 5’’ H5’ C5’ O5’ O2a O1a Pa O3a O2b O1b Pb O3b O3g O2g Pg O1g Mg

(i)

-0.6

(ii)

-0.4

-0.2

0

0.2

0.4

0.6

0.8

-0.8

-0.6

(ii)

-0.4

-0.2

H H O H H O H H O H H O O2’ O2’ H2’ C2’ H3’ O3’ O3’ C3’ C4 N3 H22 H21 N2 C2 H1 N1 O6 C6 C5 N7 H8 C8 N9 H1’ C1’ O4’ H4’ C4’ 5’’ H5’ C5’ O5’ O2a O1a Pa O3a O2b O1b Pb O3b O3g O2g Pg O1g Mg

(ii)

-0.2

0

0

0.2

0.4

0.6

0.8

(f) WT vs Q61R

(i)

-0.4

0.8

Lowdin Reduced Atomic Charges

Wild-type Q61K

-0.6

0.6

(i)

(e) WT vs Q61K

-0.8

0.4

Wild-type Q61L

Lowdin Reduced Atomic Charges

H H O H H O H H O H H O O2’ O2’ H2’ C2’ H3’ O3’ O3’ C3’ C4 N3 H22 H21 N2 C2 H1 N1 O6 C6 C5 N7 H8 C8 N9 H1’ C1’ O4’ H4’ C4’ 5’’ H5’ C5’ O5’ O2a O1a Pa O3a O2b O1b Pb O3b O3g O2g Pg O1g Mg

0.2

(d) WT vs Q61L H H O H H O H H O H H O O2’ O2’ H2’ C2’ H3’ O3’ O3’ C3’ C4 N3 H22 H21 N2 C2 H1 N1 O6 C6 C5 N7 H8 C8 N9 H1’ C1’ O4’ H4’ C4’ 5’’ H5’ C5’ O5’ O2a O1a Pa O3a O2b O1b Pb O3b O3g O2g Pg O1g Mg

Wild-type Q61H

-0.8

0

Lowdin Reduced Atomic Charges

0.2

0.4

0.6

0.8

Wild-type Q61R

(i)

-0.8

-0.6

Lowdin Reduced Atomic Charges

(ii)

-0.4

-0.2

0

0.2

0.4

0.6

0.8

Lowdin Reduced Atomic Charges

Figure 4: Löwdin population analysis of GTP, M g 2+ and water molecules within WT NRas (black rectangles) vs (red rectangles) (a) Q61E, (b) Q61P, (c) Q61H, (d) Q61L, (e) Q61K and (f) Q61R, along with the corresponding KS HOMOs. The positive lobes of the KS HOMOs are depicted in blue and the negative ones in red. charged within Q61H and Q61L. Finally, within Q61K and Q61R, O1γ is more negatively charged while O2γ less negatively charged than within WTQM/M M . 17

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 35

The bridging oxygen atom between the γ and β phosphate groups holds a slightly larger negative charge (∼ 0.02) within Q61P, Q61H, Q61L and Q61K mutants than within WTQM/M M NRas. This same atom is slightly less negatively charged (∼ 0.01 a.u.) within Q61E and Q61R than within WTQM/M M . The Löwdin reduced charge held by Pβ , O2β and Oαβ atoms is practically the same within WTQM/M M and the studied Gln 61 mutant proteins (differences of ∼ ±0.01 a.u.). Finally, the α phosphorus atom also holds approximately the same positive reduced charge within WTQM/M M and the six studied Gln 61 mutant proteins. O1α and O2α atoms are less negatively charged within Q61E, Q61P, Q61H, Q61L and Q61R than within WTQM/M M NRas. Within Q61K, O1α is also less negatively charged than within the WTQM/M M while O2α holds 0.05 a.u. more negative charge than within WTQM/M M . The overall reduced Löwdin charges attributed to GTP α, β and γ phosphate groups and to the Mg2+ ion, when bound to WT and mutant NRas, are summarized in table 1. From this table, it can be seen that WTQM/M M NRas active site arrangement, along with water molecules inside it, is the only structural conformation that leads to a bound ligand with a more negatively charged β phosphate group than its γ phosphate group is. These relative charge distribution is in accordance with numerous experimental and theoretical investigations of WT Ras proteins, 25,26,28,29,62 as stated in the previous section. Table 1: Löwdin reduced charges (in a.u.) for WT NRas and six Gln 61 mutants. WTQM/M M

WTM M

Q61E

Q61P

Q61H

Q61L

Q61K

Q61R

γ P O42−

-0.50

-0.81

-0.76

-0.77

-0.65

-0.74

-0.86

-0.77

β P O4−

-0.57

-0.66

-0.64

-0.59

-0.61

-0.70

-0.63

-0.57

α P O4−

-0.78

-0.82

-0.75

-0.74

-0.76

-0.69

-0.82

-0.70

M g 2+

-0.33

-0.18

-0.02

0.12

-0.04

-0.13

-0.04

-0.29

Furthermore, for five of the six studied Gln 61 mutant proteins (i.e. Q61E, Q61P, Q61L, Q61K and Q61R), the γ phosphate group not only is more negatively charged than the β phosphate group, but also holds more negative charge than the α phosphate group. Hence, 18

ACS Paragon Plus Environment

Page 19 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

for these five NRas mutants, GTP phosphate groups reduced charge follows the trend:

γP O2− electronic density > αP O− electronic density > βP O− electronic density 4

4

4

while it is: αP O4− > βP O4− > γP O42− for WTQM/M NRas, as showed in the previous section (and αP O4− > γP O42− > βP O4− for Q61H and WTM M NRas). These results suggest that Gln 61 substitutions hinder the catalyzing effect of WT NRas that consists in modifying GTP electronic density, further stabilizing a GDP-like transition state than the GTP reactant. The loss of this effect could be correlated to the alteration of the active site optimal electrostatic field necessary for an efficient catalysis of the reaction, as recently demonstrated by Novelli et al. 42 HOMO As for WT NRas, the KS HOMO of the DFT region3 , within the environment provided by Gln 61 substituted proteins studied here, is depicted for each mutant in inset (ii) of the corresponding panel of figure 4. As it can be seen, not only the distribution of the highest energy electrons does not systematically follow that inferred from the Löwdin population analysis, notably concerning reduced charges on GTP phosphate groups, but also the reactivity of the ligand appears to be peculiar to specific Gln 61 mutated proteins. For Q61E and Q61L, using the same cutoff (i.e. the same isovalue) for the KS HOMO delimiting surface as for WT NRas (|0.02|e− /a30 ), the highest energy electrons are mainly localized on, respectively, one and two of the additional water molecules identified in our previous work. For Q61E, a part of the highest energy electron density is also located on the three γ oxygen atoms such that the O3γ atom holds a larger share than the two others. For Q61L, O1γ atom is the only one to also hold a portion of electron density, while significantly smaller than the amount located on the two water molecules. Because within the selected cutoff for delimiting the KS HOMO surface, only GTP γ oxygen atoms appear to hold some electron density, the γ phosphate group appears to be more negatively charged than the α consisting in the GTP ligand, a Mg2+ ion with two coordinated water molecules and other water molecules present in the active site 3

19

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

and β phosphate groups within these two particular mutants. This is in accordance with the Löwdin population analysis. The Q61E GTP γ phosphate group appears to be more negatively charged than the same chemical group within Q61L, again in accordance with the previous reduced charge attribution (see table 1). Transient water molecules within the two mutant proteins appear equally reactive towards an electrophile. For Q61P, electron density is prominent on GTP γ oxygen atoms. A significantly smaller amount is located on (in decreasing order): i) one of the water molecules coordinated to the Mg2+ ion, ii) one of the transient water molecules from the active site and iii) the γ − β bridging oxygen atom. As within Q61E and Q61L, according to the KS HOMO, GTP γ phosphate group within Q61P is more negatively charged than the two other phosphate groups, as determined by the Löwdin population analysis. The intake of electron density by one of the Mg2+ coordinated water molecules could explain why a positive Löwdin reduced charge is attributed to this ion within Q61P. For Q61H and Q61K, the highest energy electron density is located on the guanosine moiety of the ligand. The KS HOMO depicted within the selected cutoff does not allow to assess the relative charges held by GTP phosphate groups. When attributing reduced Löwdin charges to these groups, lower energy occupied MO must contribute more importantly than the KS HOMO. Finally, within Q61R, the highest energy electron distribution according to the KS HOMO is the most similar to that observed within WTQM/M M . It appears prominently on GTP α oxygen atoms, followed by Oαβ and O1β atoms that display a similar amount, succeeded by Oγβ and O1γ atoms after which, O2γ atom holds a smaller share. Lastly, one of the water molecules coordinated to the Mg2+ ion holds the smallest part, in the scope of the selected isovalue. In addition to the Mg2+ coordinated water molecule that appears to take a non negligible amount of KS HOMO electron density, two important differences that arise within Q61R compared to WTQM/M M NRas should be pointed out: i) O2β atom holds no high-energy electron density contrary to what is observed within the WTQM/M M and ii)

20

ACS Paragon Plus Environment

Page 20 of 35

Page 21 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

O1γ atom has the biggest share among the γ oxygen atoms while it has none within the WTQM/M M .

4

Conclusions

The comparison between results obtained from MM and QM/MM MD for studying WT and Q61R NRas, leads to the conclusion that the employed MM force field together with the solvent model are suitable to qualitatively determine the RDF of water molecules around particular GTP atoms. This level of theory predicts the same relative amounts of solvent in the vicinity of atoms of interest as when using a hybrid approach. It remains to be determined if the minimum average distance between solute and solvent atoms is method-dependent, as is the shift observed between peaks in the MM and QM/MM RDF curves. Nevertheless, within the simulated time, a higher level of theory is needed for enabling structural rearrangements of the active site to take place upon Gln 61 substitutions. Moreover, it is also required for accurately mapping water molecules in the protein active site. The electrostatic coupling between atoms forming the latter and solvent molecules appears to be necessary for the establishment of the finely tuned interactions that lead to a precise positioning of reactant and structuring water molecules. Indeed, the examination of GTP electronic density, through a Löwdin population analysis, supports water placement in the WT protein as observed when employing a hybrid QM/MM approach: WTQM/M M is the only active site conformation that leads to GTP β phosphate group being more negatively charged than the γ phosphate. This is in good agreement with experimental studies 25,26,28,29 that showed that GTP binding to Ras proteins induces a GDPlike distribution of its electronic density. We presume that the water molecule identified in the QM/MM study, held by nitrogen backbone atoms of residues 12 and 60, has a structuring role. The conformation ensured by this water molecule, enables the establishment of proteinligand interactions that result in the stabilization of a GDP-like transition state.

21

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Moreover, Rudack et al 17 suggested that this process is partially achieved by the Mg2+ ion that temporarily stores electrons from GTP γ phosphate group, rendering it more positive. This mechanism is more importantly observed within WTQM/M M . Indeed, the Löwdin population analysis shows that this conformation leads to the cation having the largest reduced negative charge, hence storing the most electrons, as well as the γ phosphorus atom being the most positively charged, hence conceding the most electrons. Gln 61 substitutions provoke the loss of this catalyzing feature, i.e. that of stabilizing a GDP-like state. Indeed, within the six studied mutant proteins, the γ phosphate group appears to hold more negative charge than the β phosphate. Nevertheless, to our knowledge, the inability of Gln 61 mutated proteins to stabilize a product-like state remains to be experimentally verified. Such investigation could be carried out employing the same techniques that led to the detection of these charge transfers that are now considered to constitute one of Ras major catalyzing effects. Finally, when analyzing GTP electronic density by depicting the KS HOMO, the highest energy electrons within WTQM/M M appear to be distributed in the same order as determined by the Löwdin population analysis i.e. αP O4− > βP O4− > γP O42− . For all Gln 61 mutant proteins, except Q61R, it appears that other lower energy electrons contribute to the reduced charge held by GTP phosphate groups, thereby shedding light on the most reactive regions of the ligand within the different environments. The impact of GTP modified electronic distribution on NRas GTPase activity should be evaluated. This is in progress within our group by simulating the reaction within the active site of NRas forms studied here.

Acknowledgement The authors thank CALMIP Supercomputer Center for CPU resources and support on project P1237.

22

ACS Paragon Plus Environment

Page 22 of 35

Page 23 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Supporting Information Available The file supplied as Supporting Information (SI) contains: - Figures representing general system properties monitored during MM MD of WT and Q61R NRas - A detailed description of results obtained from these MM MD simulations The Supporting Information is available free of charge on the ACS Publications website at DOI:

References (1) Wittinghofer, A.; Herrmann, C. Ras-Effector Interactions, the Problem of Specificity. FEBS Letters 1995, 369, 52–56. (2) Bourne, H.; Sanders, D.; McCormick, F. The GTPase Superfamily: a Conserved Switch for Diverse Cell Functions. Nature. 1990, 348, 125–132. (3) Carvalho, A.; Szeler, K.; Vavitsas, K.; Aqvist, J.; Kamerlin, S. Modeling the Mechanisms of Biological GTP Hydrolysis. Arch. Biochem. Biophys. 2015, 582, 80–90. (4) Lu, S.; Jang, H.; Muratcioglu, S.; Gursoy, A.; Keskin, O.; Nussinov, R.; Zhang, J. Ras Conformational Ensembles, Allostery, and Signaling. Chem. Rev. 2016, 116, 6607–6665. (5) Trahey, M.; McCormick, F. A Cytoplasmic Protein Stimulates Normal N-ras p21 GTPase, but Does Not Affect Oncogenic Mutants. Science. 1987, 238, 542–545. (6) McCormick, F. ras GTPase Activating Protein: Signal Transmitter and Signal Terminator. Cell. 1989, 56, 5–8. (7) McCormick, F. The World According to GAP. Oncogene. 1990, 5, 1281–1283.

23

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(8) Cherfils, J.; Zeghouf, M. Regulation of Small GTPases by GEFs, GAPs, and GDIs. Physiol. Rev. 2013, 93, 269–309. (9) Mishra, A.; Lambright, D. Small GTPases and their GAPS. Biopolymers 2016, 105, 431–448. (10) Kraulis, P.; Domaille, P.; Campbell-Burk, S.; Van Aken, T.; Laue, E. Solution Structure and Dynamics of ras p21-GDP Determined by Heteronuclear 3-Dimensional and 4Dimensional NMR-Spectroscopy. Biochemistry. 1994, 33, 3515–3531. (11) Scheffzek, K.; Ahmadian, M.; Kabsch, W.; Wiesmüller, L.; Lautwein, A.; Schmitz, F.; Wittinghofer, A. The Ras-RasGAP Complex: Structural Basis for GTPase Activation and Its Loss in Oncogenic Ras Mutants. Science. 1997, 277, 333–338. (12) Glennon, T.; Villà, J.; Warshel, A. How Does GAP Catalyze the GTPase Reaction of Ras? A Computer Simulation Study. Biochemistry. 2000, 39, 9641–9651. (13) Resat, H.; Straatsma, T.; Dixon, D.; Miller, J. The Arginine Finger of RasGAP Helps Gln-61 Align the Nucleophilic Water in GAP-Stimulated Hydrolysis of GTP. Proc. Natl. Acad. Sci. USA. 2001, 98, 6033–6038. (14) Shurki, A.; Warshel, A. Why Does the Ras Switch Break by Oncogenic Mutations? Proteins. 2004, 55, 1–10. (15) Geyer, M.; Schweins, T.; Herrmann, C.; Prisner, T.; Wittinghofer, A.; Kalbitzer, H. Conformational Transitions in p21ras and in Its Complexes With the Effector Protein Raf-RBD and the GTPase Activating Protein GAP. Biochemistry. 1996, 35, 10308– 10320. (16) Ahmadian, M.; Stege, P.; Scheffzek, K.; Wittinghofer, A. Confirmation of the ArginineFinger Hypothesis for the GAP-Stimulated GTP-Hydrolysis Reaction of Ras. Nat. Struct. Biol. 1997, 4, 686–689. 24

ACS Paragon Plus Environment

Page 24 of 35

Page 25 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(17) Rudack, T.; Xia, F.; Schlitter, J.; Kötting, C.; Gerwert, K. The Role of Magnesium for Geometry and Charge in GTP Hydrolysis, Revealed by Quantum Mechanics/Molecular Mechanics Simulations. Biophys. J. 2012, 103, 293–302. (18) Bos, J. ras Oncogenes in Human Cancer: a Review. Cancer Res. 1989, 49, 4682–4689. (19) McGrath, J.; Capon, D.; Goeddel, D.; Levinson, A. Comparative Biochemical Properties of Normal and Activated Human ras p21 Protein. Nature. 1984, 310, 644–649. (20) Sweet, R.; Yokoyama, S.; Kamata, T.; Feramisco, J.; Rosenberg, M.; Gross, M. The Product of ras is a GTPase and the T24 Oncogenic Mutant is Deficient in This Activity. Nature. 1984, 311, 273–275. (21) Pai, E.; Krengel, U.; Petsko, G.; Goody, R.; Kabsch, W.; Wittinghofer, A. Refined Crystal Structure of the Triphosphate Conformation of H-ras p21 at 1.35 A Resolution:Implications for the Mechanism of GTP Hydrolysis. EMBO J 1990, 9, 2351–2359. (22) Privé, G.; Milburn, M.; Tong, L.; deVos, A.; Yamaizumi, Z.; Nishimura, S.; Kim, S.-H. X-ray Crystal Structures of Transforming p21 ras Mutants Suggest a Transition-State Stabilization Mechanism for GTP Hydrolysis. Biochemistry. 1992, 89, 3649–3653. (23) Schweins, T.; Geyer, M.; Scheffzek, K.; Warshel, A.; Kalbitzer, H.; Wittinghofer, A. Substrate-Assisted Catalysis as a Mechanism for GTP Hydrolysis of p21ras and Other GTP-Binding Proteins. Nat. Struct. Biol. 1995, 2, 36–44. (24) Schweins, T.; Warshel, A. Mechanistic Analysis of the Observed Linear Free Energy Relationships in p21ras and Related Systems. Biochemistry. 1996, 35, 14232–14243. (25) Cepus, V.; Scheidig, A.; Goody, R.; Gerwert, K. Time-Resolved FTIR Studies of the GTPase Reaction of H-Ras p21 Reveal a Key Role for the β - Phosphate. Biochemistry. 1998, 37, 10263–10271.

25

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(26) Wang, J.; Xiao, D.; Deng, H.; Webb, M. R.; Callender, R. Raman Difference Studies of GDP and GTP Binding to c-Harvey ras. Biochemistry. 1998, 37, 11106–11116. (27) Scheidig, A.; Burmester, C.; Goody, R. The Pre-Hydrolysis State of p21(ras) in Complex With GTP: New Insights into the Role of Water Molecules in the GTP Hydrolysis Reaction of Ras-Like Proteins. Structure. 1999, 7, 1311–1324. (28) Du, X.; Frei, H.; Kim, S.-H. The Mechanism of GTP Hydrolysis by Ras Probed by Fourier Transform Infrared Spectroscopy. J. Biol. Chem. 2000, 275, 8492–8500. (29) Allin, C.; Ahmadian, M.; Wittinghofer, A.; Gerwert, K. Monitoring the GAP Catalyzed H-Ras GTPase Reaction at Atomic Resolution in Real Time. Proc. Natl. Acad. Sci. USA. 2001, 98, 7754–7759. (30) Buhrman, G.; Holzapfel, G.; Fetics, S.; Mattos, C. Allosteric Modulation of Ras Positions Q61 for a Direct Role in Catalysis. Proc. Natl. Acad. Sci. USA. 2010, 107, 4931–4936. (31) Schweins, T.; Langen, R.; Warshel, A. Why Have Mutagenesis Studies Not Located the General Base in ras p21. Nat. Struct. Biol. 1994, 1, 476–484. (32) Soares, T.; Miller, J.; Straatsma, T. Revisiting the Structural Flexibility of the Complex p21ras-GTP: the Catalytic Conformation of the Molecular Switch II. Protein Struct. Funct. Genet. 2001, 45, 297–312. (33) Topol, I.; Cachau, R.; Nemukhin, A.; Grigorenko, B.; Burt, S. Quantum Chemical Modeling of the GTP Hydrolysis by the RAS-GAP Protein Complex. Biochim. Biophys. Acta. 2004, 1700, 125–136. (34) Grigorenko, B.; Nemukhin, A.; Topol, I.; Cachau, R.; Burt, S. QM/MM Modeling the Ras-GAP Catalyzed Hydrolysis of Guanosine Triphosphate. Proteins. 2005, 60, 495–503. 26

ACS Paragon Plus Environment

Page 26 of 35

Page 27 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(35) Grigorenko, B.; Nemukhin, A.; Shadrina, M.; Topol, I.; Burt, S. Mechanisms of Guanosine Triphosphate Hydrolysis by Ras and Ras-GAP Proteins as Rationalized by AbInitio QM/MM Simulations. Proteins. 2007, 66, 456–466. (36) Prasad, B.; Plotnikov, N.; Warshel, A. Addressing Open Questions about Phosphate Hydrolysis Pathways by Careful Free Energy Mapping. J. Phys. Chem. B. 2013, 117, 153–163. (37) Prakash, P.; Gorfe, A. Overview of Simulation Studies on the Enzymatic Activity and Conformational Dynamics of the GTPase Ras. Mol. Simul. 2014, 40, 839–847. (38) Khrenova, M.; Grigorenko, B.; Kolomeisky, A.; Nemukhin, A. Hydrolysis of Guanosine Triphosphate (GTP) by the Ras-GAP Protein Complex: Reaction Mechanism and Kinetic Scheme. J. Phys. Chem. B. 2015, 119, 12838–12845. (39) Khrenova, M.; Grigorenko, B.; Nemukhin, A. Theoretical Vibrational Spectroscopy of Intermediates and the Reaction Mechanism of the Guanosine Triphosphate Hydrolysis by the Protein Complex Ras-GAP. Spectrochim. Acta, Part A 2016, 166, 68–72. (40) Grigorenko, B.; Khrenova, M.; Nemukhin, A. Amide-Imide Tautomerization in the Glutamine Side Chain in Enzymatic and Photochemical Reactions in Proteins. Phys. Chem. Chem. Phys. 2018, 20, 23827–23836. (41) Der, C.; Finkel, T.; Cooper, G. Biological and Biochemical Properties of Human ras Genes Mutated at Codon 61. Cell. 1986, 44, 167–176. (42) Novelli, E.; First, J.; Webb, L. Quantitative Measurement of Intrinsic GTP Hydrolysis for Carcinogenic Glutamine 61 Mutants in H-Ras. Biochemistry. 2018, 57, 6356–6366. (43) Tichauer, R.; Favre, G.; Cabantous, S.; Landa, G.; Hemeryck, A.; Brut, M. Water Distribution Within Wild-Type NRas Protein and Q61 Mutants During Unrestrained QM/MM Dynamics. Biophys. J. 2018, 115, 1417–1430. 27

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(44) Forbes, S.; Bhamra, G.; Bamford, S.; Dawson, E.; Kok, C.; Clements, J.; Menzies, A.; Teague, J.; Futreal, P.; Stratton, M. The Catalogue of Somatic Mutations in Cancer (COSMIC). Curr. Protoc. Hum. Genet. 2008, (45) Nedyalkova, L.; Tong, Y.; Tempel, W.; Shen, L.; Loppnau, P.; Arrowsmith, C.; Edwards, A.; Bountra, C.; Weigelt, J.; Bochkarev, A. et al. Brookhaven Protein Data Bank entry 3CON. (46) Bourne, H.; Sanders, D.; McCormick, F. The GTPase Superfamily: Conserved Structure and Molecular Mechanism. Nature. 1991, 349, 117–127. (47) Spoerner, M.; Herrmann, C.; Vetter, I.; Kalbitzer, H.; Wittinghofer, A. Dynamic Properties of the Ras Switch I Region and its Importance for Binding to Effectors. Proc. Natl. Acad. Sci. USA. 2001, 98, 4944–4949. (48) Case, D.; Cerutti, D.; Cheatham, T.; Darden, T.; Duke, R.; Giese, T.; Gohlke, H.; Goetz, A.; Greene, D.; Homeyer, N. et al. AMBER 2017. University of California: San Francisco, 2017. (49) Pérez, A.; Marchán, I.; Svozil, D.; Sponer, J.; Cheatham, T.; Laughton, C.; Orozco, M. Refinement of the AMBER Force Field for Nucleic Acids: Improving the Description of alpha/gamma Conformers. Biophys. J. 2007, 92, 3817–3829. (50) Price, D.; Brooks, C. A Modified TIP3P Water Potential for Simulation With Ewald Summation. J. Chem. Phys. 2004, 121, 10096–10103. (51) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. Numerical Integration of the Cartesian Equations of Motion of a System With Constraints: Molecular Dynamics of N-Alkanes. J. Chem. Phys. 1977, 23, 327–341. (52) Roe, D.; Cheatham, T. E. PTRAJ and CPPTRAJ: Software for Processing and Analysis of Molecular Dynamics Trajectory Data. J. Chem. Theory Comput. 2013, 9, 3084–3095. 28

ACS Paragon Plus Environment

Page 28 of 35

Page 29 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(53) Hawkins, G.; Cramer, C.; Truhlar, D. Parametrized Models of Aqueous Free Energies of Solvation Based on Pairwise Descreening of Solute Atomic Charges from a Dielectric Medium. J. Phys. Chem. 1996, 100, 19824–19839. (54) Neese, F. The ORCA Program System. Wiley Interdiscip Rev Comput Mol Sci 2012, 2, 73–78. (55) Neese, F. Software Update: the ORCA Program System, Version 4.0. Wiley Interdiscip Rev Comput Mol Sci 2018, 8, 1–6. (56) Götz, A.; Clark, M.; Walker, R. An Extensible Interface for QM/MM Molecular Dynamics Simulations With AMBER. J. Comput. Chem. 2014, 35, 95–108. (57) Walker, R.; Crowley, M.; Case, D. The Implementation of a Fast and Accurate QM/MM Potential Method in Amber. J. Comput. Chem. 2008, 29, 1019–1031. (58) Weigend, F.; Ahlrichs, R. Balanced Basis Sets of Split Valence, Triple Zeta Valence and Quadruple Zeta Valence Quality for H to Rn: Design and Assessment of Accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297–3305. (59) Weigend, F. Accurate Coulomb-fitting Basis Sets for H to Rn. Phys. Chem. Chem. Phys. 2006, 8, 1057–1065. (60) Hellweg, A.; Hättig, C.; Höfener, S.; Klopper, W. Optimized Accurate Auxiliary Basis Sets for RI-MP2 and RI-CC2 Calculations for the Atoms Rb to Rn. Theoretical Chemistry Accounts 2007, 117, 587–597. (61) Meagher, K.; Redman, L.; Carlson, H. Development of Polyphosphate Parameters for Use With the AMBER Force Field. J. Comput. Chem. 2003, 24, 1016. (62) Klähn, M.; Schlitter, J.; Gerwert, K. Theoretical IR Spectroscopy Based on QM/MM Calculations Provides Changes in Charge Distribution, Bond Lengths, and Bond Angles of the GTP Ligand Induced by the Ras-Protein. Biophys. J. 2005, 88, 3829–3844. 29

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(63) Cramer, C. J. Essentials of Computational Chemistry, 2nd ed.; Wiley, 2004. (64) Jensen, F. Introduction to Computational Chemistry, 2nd ed.; Wiley, 2007.

30

ACS Paragon Plus Environment

Page 30 of 35

Page 31 of 35

Graphical TOC Entry       

``` ```

c c c c c c

Switch I

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

P-loop S SwitchSII S S S S SS

31

ACS Paragon Plus Environment

``` ```

Page 32 of 35 Gly 13 Gly 12 Gln 61 Pβ Pα



Tyr 32 c c c c c c

Switch I

1 2 3 4 5 6 7 8 9 (a)

The Journal of Physical Chemistry

P-loop

Mg2+

ACS Paragon Plus Environment Switch II

(b)

Thr 35

Gly 60

(a)WT Page 33QM/M of 35M

1 2 3 4 5 6 (c)Q61E 7 8 9 10 11 12 13 14 15 (f )Q61L 16 17 18 19 20 21 22 23

(b)WT The Journal M M of Physical Chemistry Water molecules coordi-

nated to the Mg2+ ion appear as gray sticks while transient water molecules are depicted in a colored ball-andstick representation. (d)Q61P

(e)Q61H

(g)Q61K

(h)Q61R

ACS Paragon Plus Environment

Wild-type structure from MM vs QM/MM

The Journal of Physical Chemistry Page 34 of 35 WT 2 H O (QM/MM)

H H O H H O H H O H H O O2’ O2’ H2’ C2’ H3’ O3’ O3’ C3’ C4 N3 H22 H21 N2 C2 H1 N1 O6 C6 C5 N7 H8 C8 N9 H1’ C1’ O4’ H4’ C4’ 5’’ H5’ C5’ O5’ O2a O1a Pa O3a O2b O1b Pb O3b O3g O2g Pg O1g Mg

1 2 3 4 5 6 7 8 9 10 11

with

2

WT with 1 H2O (MM)

(a)

(b)

ACS Paragon Plus Environment -0.8

-0.6

-0.4

-0.2

0

0.2

0.4

Lowdin Reduced Atomic Charges

0.6

0.8

(a) WT vs Q61E

Page 35 of 35 H H O H H O H H O H H O O2’ O2’ H2’ C2’ H3’ O3’ O3’ C3’ C4 N3 H22 H21 N2 C2 H1 N1 O6 C6 C5 N7 H8 C8 N9 H1’ C1’ O4’ H4’ C4’ 5’’ H5’ C5’ O5’ O2a O1a Pa O3a O2b O1b Pb O3b O3g O2g Pg O1g Mg

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39

H H O H H O H H O H H O O2’ O2’ H2’ C2’ H3’ O3’ O3’ C3’ C4 N3 H22 H21 N2 C2 H1 N1 O6 C6 C5 N7 H8 C8 N9 H1’ C1’ O4’ H4’ C4’ 5’’ H5’ C5’ O5’ O2a O1a Pa O3a O2b O1b Pb O3b O3g O2g Pg O1g Mg

Q61E

(i)

-0.8

-0.6

(ii)

-0.4

-0.2

0

(b) WT vs Q61P

The Journal of Physical Chemistry Wild-type

0.2

0.4

0.6

0.8

Wild-type Q61P

(i)

-0.8

-0.6

(ii)

-0.4

Lowdin Reduced Atomic Charges

-0.2

(c) WT vs Q61H

H H O H H O H H O H H O O2’ O2’ H2’ C2’ H3’ O3’ O3’ C3’ C4 N3 H22 H21 N2 C2 H1 N1 O6 C6 C5 N7 H8 C8 N9 H1’ C1’ O4’ H4’ C4’ 5’’ H5’ C5’ O5’ O2a O1a Pa O3a O2b O1b Pb O3b O3g O2g Pg O1g Mg

(i)

-0.6

(ii)

-0.4

-0.2

0

0.4

0.6

0.8

0.2

0.4

0.6

0.8

Wild-type Q61L

(i)

-0.8

-0.6

(ii)

-0.4

Lowdin Reduced Atomic Charges

-0.2

0

0.2

0.4

0.6

0.8

Lowdin Reduced Atomic Charges

(e) WT vs Q61K

H H O H H O H H O H H O O2’ O2’ H2’ C2’ H3’ O3’ O3’ C3’ C4 N3 H22 H21 N2 C2 H1 N1 O6 C6 C5 N7 H8 C8 N9 H1’ C1’ O4’ H4’ C4’ 5’’ H5’ C5’ O5’ O2a O1a Pa O3a O2b O1b Pb O3b O3g O2g Pg O1g Mg

0.2

(d) WT vs Q61L H H O H H O H H O H H O O2’ O2’ H2’ C2’ H3’ O3’ O3’ C3’ C4 N3 H22 H21 N2 C2 H1 N1 O6 C6 C5 N7 H8 C8 N9 H1’ C1’ O4’ H4’ C4’ 5’’ H5’ C5’ O5’ O2a O1a Pa O3a O2b O1b Pb O3b O3g O2g Pg O1g Mg

Wild-type Q61H

-0.8

0

Lowdin Reduced Atomic Charges

(f) WT vs Q61R H H O H H O H H O H H O O2’ O2’ H2’ C2’ H3’ O3’ O3’ C3’ C4 N3 H22 H21 N2 C2 H1 N1 O6 C6 C5 N7 H8 C8 N9 H1’ C1’ O4’ H4’ C4’ 5’’ H5’ C5’ O5’ O2a O1a Pa O3a O2b O1b Pb O3b O3g O2g Pg O1g Mg

Wild-type Q61K

(i)

(ii)

Wild-type Q61R

(i)

(ii)

ACS Paragon Plus Environment -0.8

-0.6

-0.4

-0.2

0

0.2

0.4

Lowdin Reduced Atomic Charges

0.6

0.8

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

Lowdin Reduced Atomic Charges

0.6

0.8