Single Mutations Reshape the Structural Correlation Network of the

Feb 8, 2017 - Rich Dynamics Underlying Solution Reactions Revealed by Sampling and Data Mining of Reactive Trajectories. ACS Central Science. Zhang, Z...
0 downloads 6 Views 6MB Size
Article pubs.acs.org/JPCB

Single Mutations Reshape the Structural Correlation Network of the DMXAA−Human STING Complex Xing Che,†,§ Xiao-Xia Du,†,§ Xiaoxia Cai,† Jun Zhang,† Wen Jun Xie,† Zhuoran Long,† Zhao-Yang Ye,‡ Heng Zhang,‡ Lijiang Yang,† Xiao-Dong Su,*,‡ and Yi Qin Gao*,† †

Institute of Theoretical and Computational Chemistry, College of Chemistry and Molecular Engineering, and Biodynamic Optical Imaging Center and ‡State Key Laboratory of Protein and Plant Gene Research, and Biodynamic Optical Imaging Center, School of Life Sciences, Peking University, Beijing 100871, China S Supporting Information *

ABSTRACT: Subtle changes in protein sequences are able to alter ligand−protein interactions. Unraveling the mechanism of such phenomena is important for understanding ligand− protein interactions, including the DMXAA−STING interaction. DMXAA specifically binds to mouse STING instead of human STING. However, the S162A mutation and a newly discovered E260I mutation endow human STINGAQ with DMXAA sensitivity. Through molecular dynamics simulations, we revealed how these single mutations alter the DMXAA− STING interaction. Compared to mutated systems, structural correlations in the interaction of STINGAQ with DMXAA are stronger, and the correlations are cross-protomers in the dimeric protein. Analyses on correlation coefficients lead to the identification of two key interactions that mediate the strong cross-protomer correlation in the DMXAA−STINGAQ interaction network: DMXAA−267T−162S* and 238R−260E*. These two interactions are partially and totally interrupted by the S162A and E260I mutations, respectively. Moreover, a smaller number of water molecules are displaced upon DMXAA binding to STINGAQ than that on binding to its mutants, leading to a larger entropic penalty for the former. Considering the sensitivity of STINGAQ and two of its mutants to DMXAA, a strong structural correlation appears to discourage DMXAA−STING binding. Such an observation suggests that DMXAA derivatives, which are deprived of hydrogen-bond interaction with both 162S* and 267T, are potential agonists of human STING.



non-small-cell lung cancer (NSCLC).12 Further investigation revealed that mouse STING (mSTING), but not human STING (hSTING), shows a response to DMXAA,13 which gives an explanation to the phase-III clinical trial failure. Recent study discovered that hSTINGS162A, a single point mutation in hSITNGH232 and hSTINGR232, can be activated by DMXAA and induces IFN-β production.8 Later, G230I and Q266I have been identified to render hSTING with DMXAA sensitivity through sequence comparison between mSTING and hSTING.14 According to the crystal structure of DMXAA−hSTINGG230I (pdb code: 4QXP), 162S locates deeply in the binding pocket and is speculated not to interfere directly with DMXAA binding, given its size and relative position to DMXAA. Therefore, the mechanism of how such a single mutation S162A alters the sensitivity of human STING to DMXAA is intriguing. In fact, it is very common for proteins similar in sequence and structure to have divergent functions. Either single or multiple point mutations have the capacity to alter protein flexibility, dynamics,

INTRODUCTION STING (stimulator of interferon genes, also known as TMEM173, MPYS, MITA, and ERIS) is a vital protein involved in immune defenses to cytosolic double-stranded DNA.1 Upon binding with endogenous cyclic dinucleotides,2 STING translocates and activates TANK binding kinase (TBK1) and IkB kinase (IKK). TBK1 and IKK then phosphorylate and activate interferon regulatory factor 3 (IRF3) and nuclear factor kappa B (NF-κB), respectively, which then trigger a type-I interferon (IFN) response.3−6 Through dimerization, the carboxyl-terminal domain of STING, consisting of about 250 residues, constitutes a binding pocket. The homodimer of CTD STING takes the shape of a pair of wings. Upon binding, the dimer transforms from an open to a closed state. The distance between the apexes of the two wings is shortened by about 20 Å.7−9 In addition to binding with endogenous second messengers, mouse STING also responds to drug-like compounds10carboxymethyl-9-acridanone (CMA)10 and 5,6-dimethylxanthenone-4-acetic acid (DMXAA, Vadimezan, AS404). The xanthenone derivative compound, DMXAA, shows an antivascular effect to tumors and induces cytokine production in mice.11 However, combined with first-line platinum-based chemotherapy, DMXAA failed in the clinical trial of advanced © 2017 American Chemical Society

Received: December 11, 2016 Revised: February 4, 2017 Published: February 8, 2017 2073

DOI: 10.1021/acs.jpcb.6b12472 J. Phys. Chem. B 2017, 121, 2073−2082

Article

The Journal of Physical Chemistry B

parameters were obtained by ANTECHAMBER, using the general atom force-field and RESP procedure to fit the point charges.27 STING AQ and its two mutants, STING AQ‑S162A and STINGAQ‑E260I, were simulated in this work. The simulation models of DMXAA−STINGAQ and DMXAA−STINGAQ‑E260I complexes were built by replacing the corresponding residues of the DMXAA−STINGAQ‑S162A crystal structure with the targeted ones. Amber99SB force-field parameters were used for proteins. Counter ions were used to neutralize the systems. All systems were solvated in a cubic box of TIP3P water molecules with a buffer of 8 Å. The size of the final simulation systems is about 70.3 × 88.3 × 83.5 Å3, containing about 56 400 atoms. A cutoff of 10.0 Å was applied for the calculations of pairwise interactions (van der Waals and direct Coulomb). The particlemesh Ewald method was employed for long-range electrostatic interaction calculations. The SHAKE algorithm was used to constrain all bonds involving hydrogen. In total, 1000 steps of steepest descent minimization and 1500 steps of conjugate gradient minimization were performed to relieve any structural clash in the solvated systems, followed by a 500 ps heatingequilibration process from 0 to 330 K in the NVT ensemble and a 500 ps cooling-equilibration process from 330 to 300 K in the NPT ensemble. The temperature was regulated by Langevin dynamics to be at 300 K. The time interval between simulation steps was 2 fs, and all data were collected every 1.0 ps. Unless specified otherwise, the STING protein in this article refers to human STING. The simulation time of the production run is about 260 ns for each system and three independent trajectories of each system were obtained. Free-Energy Decomposition. Free-energy calculations together with MD simulations yield ligand−protein binding energies quantitatively. Decomposition of binding energies further enables us to identify key residues in ligand−protein interactions. In our free-energy decomposition analysis, we used the scheme “per-residue decomposition” provided by AMBERTOOLS, MMPBSA.py28 to analyze key residues in DMXAA− STING interactions. Free energies were decomposed into specific residue contributions, using implicit solvent models.29,30 The first 100 ns of MD trajectories were discarded in the calculation of binding free energy decomposition and snapshots used in the calculation were taken every 20 ps. Correlation Analyses. Correlations in the DMXAA− STING interaction network are described by a series of distance correlations. The correlation coefficients were calculated as

and, ultimately, function. It has been recognized that dynamics is important in populating available protein conformational states; rates and probabilities of states redistribution upon ligand binding15 and proteins shift conformational populations through cooperative motions, while enabling their functions.16 For example, in cyclophilin A enzyme-catalyzed reaction, rates of substrate turnover strongly correlate with the rates of conformational dynamics of the enzyme.17,18 It is also known that single mutation could exert an influence on the ligand−protein interaction network. For proteins such as GPCRs, two structural elements, rigid knobs and flexible regions, compose the allosteric communication network. Such a structure allows the protein to function like a machine, with information transmitted among flexible domains through the rigid knobs.19 A disease-associated mutation that softens the rigid knobs leads to flabby dynamics and disrupts dynamic allosteric regulation.20 The fact that mutations located far away from the binding pocket exert an influence on ligand binding indicates that conformational transitions are communicated between the mutation and binding sites, resulting in changes in the volume and shape of the pocket.21 In addition, the role of water in ligand binding has long been recognized. Studies have shown that structural water directly mediates ligand−protein interactions, and water-bridged ligand−protein hydrogen-bond interactions have been found important in a number of systems.22 Moreover, waters around the hydrophobic residues in the binding pocket may lead to a large desolvation penalty.23 Water-mediated interactions are also involved in sculpting protein functional landscapes.24 Therefore, deciphering the influence of mutated residues on the protein− ligand interaction network, including the surrounding water molecules, will shed light on understanding the response of protein to ligand binding. The complicated energy landscape of a ligand−protein complex in conjunction with solvent interactions makes it difficult to observe protein dynamics at atomic resolution in real time. Crystal structures provide useful but limited information on protein dynamics. Alternatively, molecular dynamics (MD) simulation becomes a widely used approach to investigate protein dynamics at the atomic level. MD simulations provide details of motions as a function of time, which are difficult to access experimentally but critical for understanding protein functions.25 In this study, to investigate the influence of a single mutation on STING dynamics, we performed MD simulations on DMXAA−hSTINGAQ and DMXAA−hSTINGAQ‑S162A complexes. The analyses of the simulation results allowed us to understand the detailed interactions between DMXAA and STING as well as the intramolecular interactions. On the basis of such understandings, we designed a new mutation, E260I, which was subsequently shown to endow hSTINGAQ with DMXAA sensitivity. We then investigated in detail the differences in interactions of STINGAQ and its two mutants through binding free energy calculations and structural correlation analyses. In addition, using the two-phase thermodynamic (2PT) model, we further compared entropic penalty differences of pocket water in the three systems. Such molecular-level understandings are expected to facilitate the design of effective DMXAA derivatives.

Cij =

⟨(ri − ri ̅ )(rj − rj̅ )⟩ ⟨(ri − ri ̅ )2 ⟩⟨(rj − rj̅ )2 ⟩

(1)

where ri represents the distance between atoms in the interaction network. In this work, the first 100 ns of MD trajectories were not used in the construction of the covariance matrix, and the snapshots were taken every 5 ps. Hydrogen-Bond Criteria. Hydrogen bond is considered formed when the heavy atom distance of donor and acceptor is less than 3.2 Å and the angle cutoff is 135°. Entropy Calculations. We used the 2PT model to calculate the absolute entropy of water in the binding pocket of STING. The 2PT method, developed by Lin et al., has been successfully applied to calculate the entropy of water and other fluids.31−34 The model is based on the hypothesis that the liquid partition function can be decomposed into solid-like and gas-like components and the entropy is a weighted summation of the



METHODS MD Simulations. The structure of DMXAA was extracted from the crystal structure of the DMXAA−STINGAQ‑S162A complex. Geometry optimization was performed using Gaussian26 on the 6-311G* basis set using B3LYP. HF/6-31G* level of theory was used in partial charge calculation. The force-field 2074

DOI: 10.1021/acs.jpcb.6b12472 J. Phys. Chem. B 2017, 121, 2073−2082

Article

The Journal of Physical Chemistry B

Figure 1. Key residues in the DMXAA−hSTINGAQ‑S162A interaction, where a and b represent different chains of the STING dimer.

two components. The gas-like component models the anharmonic low-frequency modes of the fluid, and the solidlike component describes the harmonic properties of the fluid. The entropy of water molecules is calculated as

∫0

Sil = kB[



dν g is(ν)W is(ν) +

∫0



dν g ig(ν)Wig(ν)]

2Δ−9/2 f 15/2 − 6Δ−3f 5 − Δ−3/2 f 7/2 + 6Δ−3/2 f 5/2 + 2f − 2

where the dimensionless diffusivity constant Δ is defined by

(2)

Δ=

where i represents the translational, rotational, and vibrational components, kB is the Boltzmann constant, W(ν) is the weighting factor, ν is the frequency, and superscripts l, s, and g refer to liquid, solid-like, and gas-like components, respectively. In this study, the liquid DOS function (gl(ν)) is calculated by the Fourier transform of the velocity autocorrelation function (VACF, C(t)) g l (ν ) =

2 lim kBT τ →∞

g s ( ν ) = g l (ν ) − g g (ν )

(3)

N

∑ m⟨υi (t )·υi (0)⟩ i=1

N

3

∑ ∑ ⟨Ijωij(t )ωij(0)⟩ i=1 j=1

(5)

where i indicates the ith water molecule, j indicates the jth principal moment of inertia of a water molecule, N is the total number of water molecules in the system, and m is the mass of a water molecule. The DOS of the gas-like component is represented by that of hard spheres g g (ν ) =

g0 ⎡ πνg ⎤2 1 + ⎢⎣ 6Nf0 ⎥⎦

(9)

W s(ν) =

βhυ − ln[1 − exp( −βhυ)] exp(βhυ) − 1

(10)

W g (ν ) =

S HS 3k

(11)

where h is the Planck constant; β equals 1/kBT; SHS is the hardsphere entropy. As the bond vibrations of the water molecules were constrained in the simulations, vibrational entropy was neglected in the 2PT calculation. The integration time step is 1 fs. Velocities and coordinates were saved every 4 fs. Trajectories of 400 ps simulations were used in the 2PT calculation. The water entropy calculated in other biological systems shows that 400 ps is enough for the convergence in the 2PT calculation.35−37 Errors were estimated as standard deviations of a series of entropy values calculated from the first 300 ps trajectories to the total 400 ps trajectories, accumulated by 20 ps in each calculation. To calculate bulk water entropy, we used a water box containing 2557 TIP3P water molecules and produced a 100 ps trajectory. Other simulation details are the same as described earlier. Luciferase Reporter Assays. Reporter assay was performed as previously described.38 Briefly, HEK293T cells were seeded in 24-well plates and transiently transfected with an IFN-β luciferase reporter plasmid and indicated expression plasmids. After 12 h, cells were stimulated with medium containing DMXAA. After an additional 12 h, cells were lysed, and the

(4)

whereas for rotational VACF, C(t) is the summation of inertiaweighted angular velocities (ω(τ)) Crot(t ) =

(8)

Furthermore, the weighting factors for the solid-like and gas-like contributions of entropy are determined by eqs 10 and 11, respectively

where T is the temperature of the simulation system. For translational and vibrational VACF, C(t) is determined from the summation over mass-weighted center-of-mass velocities and intramolecular vibration velocities (υ(τ)), respectively C trans/vib(t ) =

2g0l ⎛ πkBT ⎞1/2 ⎛ N ⎞1/3⎛ 6 ⎞2/3 ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ 9N ⎝ m ⎠ ⎝ V ⎠ ⎝ π ⎠

In eq 8, gl0 represents the DOS of liquid at zero frequency and V is the volume of the system. As the translational, rotational, and vibrational DOS of water molecules are partitioned into gas-like and solid-like components, gs(ν) is calculated as

t

∫−t C(t ) e−i2πυt dt

(7)

=0

(6)

where g0 is the DOS intensity of the real system at zero frequency. The fluidity factor f is obtained from 2075

DOI: 10.1021/acs.jpcb.6b12472 J. Phys. Chem. B 2017, 121, 2073−2082

Article

The Journal of Physical Chemistry B reporter activity was analyzed using the Dual-Luciferase Reporter Assay System (Promega).



RESULTS Designed Single Mutation E260I Endows hSTING with DMXAA Response. To unravel why the S162A mutant is able to respond to DMXAA, we first analyzed how STINGAQ‑S162A interacts with DMXAA. Here, we used the free-energy decomposition calculation to identify the key residues in the DMXAA−STINGAQ‑S162A interaction. The per-residue freeenergy decomposition scheme provides insights into the contribution of each residue in the ligand−protein interaction. Residues that contribute more than 1.5 kcal/mol to the binding free energy are considered significant. Our analysis revealed 19 key residues, and they are listed in Figure 1. These residues locate in the binding pocket (like 263T and 267T) and the lid region (such as 238R and 232R). Among them, two 238R residues play the most prominent role. This result indicates that the electrostatic interaction is critical for DMXAA binding. The failed response of STINGAQ‑S162A to DMXAA with the R238A mutation also supports the important role of the electrostatic interaction between 238R and DMXAA (Figure S1). Several other key residues, for example, 263T and 267T, are expected to act as anchors in DMXAA binding. Moreover, an earlier study showed that G230I in the lid region and Q266I in the pocket region enable STING to respond to DMXAA,14 showing the critical roles of hydrophobic interactions. These observations shed light on the discovery of new mutations that would endow STING with DMXAA sensitivity. To further enhance the attractive electrostatic interactions or hydrophobic interactions between DMXAA and STING, a number of new mutations were designed in this study (more details in Table S1). The responses of these STINGAQ mutants to DMXAA were examined by IFN-β activation through luciferase reporter assays. One mutation type, E260I, is shown to respond to DMXAA and activate IFN-β production, as illustrated in Figure 2.

Figure 3. IFN-β expression of STINGAQ with mutations to 260E.

DMXAA than when it is mutated to polar, uncharged threonine and glutamine. Binding Free Energy Comparison among STINGAQ and Its Mutants. To understand how S162A or E260I exerts influences on the DMXAA−STING interaction, we compared the binding free energy differences among the three systems. Again, the contribution of each residue to the binding free energy was calculated. Residues for which the contribution to the binding free energy is greater than 1.5 kcal/mol are listed in Tables 1 and S2 (residues with positive values are excluded). Table 1. Differences in the Per-Residue Free-Energy Contribution system number of main contributors contribution of 238R (kcal/mol) same key residues different key residuesa

STINGAQ

STINGAQ‑S162A

STINGAQ‑E260I

11

15

20

−7.7/−3.5

−17.5/−9.8

−16.2/−15.0

238R, 263T, 267T, 166G, 264P 162S

240Y, 169R, 167Y, 239Y, 232R, 178R

178R, 169R, 232R, 167Y, 163Y, 162S, 236K, 165I

a

The residues are sequenced by their contribution magnitude to binding free energy.

More residues contribute to the interaction with DMXAA significantly in mutated systems than those in STINGAQ, as the numbers of key contributors in STINGAQ, STINGAQ‑S162A, and STINGAQ‑E260I are 11, 15, and 20, respectively (Table 1). More specifically, compared to STINGAQ, 232R, 178R, 169R, and 167Y in STINGAQ‑S162A and STINGAQ‑E260I contribute significantly to the DMXAA−STING interaction. Another difference is that although the two 238R residues are dominant in DMXAA binding, their contributions are quite different among the three systems, increasing in the order STINGAQ, STINGAQ‑S162A, and STINGAQ‑E260I on average (Table 1). The value of per-residue free energy is a summation of four components: the electrostatic interaction, van der Waals interaction, polar and nonpolar solvation free energy.29 Therefore, we further compared the differences of the four components of 238R. Table 2 lists the contributions from electrostatic interaction and polar solvation free energy (the contributions from the other two parts are

Figure 2. IFN-β expression of the designed mutations.

Furthermore, we mutated 260E to other 16 kinds of amino acids to evaluate the effect of mutation in this location. Luciferase reporter assays showed that E260Q, E260T, and E260V also endow STINGAQ with DMXAA sensitivity but with lower IFN-β readouts than E260I (Figure 3). From another perspective, the mutation of the negatively charged glutamic acid to apolar isoleucine and valine leads to a better STING response to 2076

DOI: 10.1021/acs.jpcb.6b12472 J. Phys. Chem. B 2017, 121, 2073−2082

Article

The Journal of Physical Chemistry B

Moreover, distances between the two pairs of helices also have a strong correlation with those in regions iii and iv, which represent intra- and intermolecular interactions. Unlike STINGAQ, there is no obvious cross-protomer correlation in STINGAQ‑S162A, and regions iii and iv show a reduced structural correlation both in the magnitude and range; whereas in STINGAQ‑E260I, all correlations are significantly weakened. We further examined how the strong cross-protomer structural correlation arises in STINGAQ. We built a correlation network according to the strength of the correlation coefficients mentioned earlier and the relative positions of the involved residues. Distances for which the magnitude of correlation coefficients is greater than 0.5 and the corresponding residues with direct interactions, such as hydrogen-bond interactions, and forming salt bridges are considered directly coupled. The correlation network and the relative positions of the corresponding residues in the complex structure are depicted in Figure 5. As illustrated in this figure, 267T and 162S* (* is used to represent the residue on the other protomer) at the bottom of the binding pocket interact with DMXAA through hydrogen bonds. When the carbonyl group of DMXAA forms a hydrogen bond with 162S*, distances between 267T/263T and DMXAA increase. At the same time, DMXAA also moves toward 238R and 240Y and tends to form hydrogen bonds with them. The competitive hydrogen-bond interactions between DMXAA−267T and DMXAA−162S* serve to switch the interaction of DMXAA with the top of the binding pocket and with the lid region. In such a manner, the perturbations at the bottom of the binding pocket propagate to the top and lid region. Moreover, when DMXAA forms a hydrogen bond with 238R or 240Y, the hydrogen bond between 238R and 260E* or 240Y and 166G is interrupted. Because 166G and 162S locate on α1-helix, 263T and 267T locate on α3-helix, respectively, structural correlation is established between the pairs of α1-helices and α3-helices. Therefore, there are two pathways of cross-protomer interactions in DMXAA−STINGAQ: DMXAA−267T−162S* and 238R− 260E*. The S162A and E260I mutations directly eliminate the DMXAA−162S*−267T and 238R−260E* pathways, respectively. Moreover, in DMXAA−STINGAQ‑E260I, E260I also interrupts the DMXAA−162S*−267T pathway. Therefore, compared to the wild type, the structural correlation in DMXAA−STINGAQ‑S162A is weak and that in DMXAA−

Table 2. Decomposition of Free-Energy Contribution from 238R (kcal/mol) 238Ra

238Rb

system

Epolar

Eele

Etot

Epolar

Eele

Etot

STINGAQ STINGAQ‑S162A STINGAQ‑E260I

63.3 55.9 47.1

−70.5 −73.1 −62.9

−7.2 −17.2 −15.8

51.4 43.5 38.5

−54.3 −52.0 −52.9

−2.9 −8.5 −14.4

relatively small). It can be seen from Table 2 that the polar solvation free energies of the two 238R residues decrease in the order STINGAQ, STINGAQ‑S162A, and STINGAQ‑E260I. Strong Structural Correlation of STINGAQ Interacting with DMXAA. The above analyses reveal that the interaction of 238R with DMXAA is significantly affected by the S162A and E260I mutations. Interestingly, 162A and 260E are about 11.9 Å (162A Cα-238R Cζ) and 8.2 Å (260E Cδ-238R Cζ) away from 238R in the crystal structure of DMXAA−STINGAQ‑S162A. Moreover, 260E is also distant from the binding pocket. An interesting question then arises: how does the distant single mutation on 162S or 260E affect the interaction of 238R with DMXAA? In many biological systems, it has been proved that such a large distance effect often relates to long-range structural correlations.39−41 Therefore, we examined in detail the structural correlations in the three systems. To identify structural changes that relate to the mutant site, Pearson correlation coefficients of 34 pairs of distances of residues were calculated. The selected residues include those above-mentioned key residues and additional ones that locate on the pairs of α1- and α3-helices. The residues on the helices are included because they form a bundle that constitutes the U-shape binding pocket. According to the structural motif of dimeric STING, the 34 distances are categorized into four groups: (i) distances between the pair of α1-helices (nos. 1−6), (ii) distances between the pair of α3helices (nos. 7−13), (iii) distances between residues forming hydrogen bonds or salt bridges (nos. 14−19), and (iv) distances between specified residues and DMXAA (nos. 20−34) (more details in the Table S3). The correlation map (Figure 4) illustrates that the structural correlations in DMXAA−STINGAQ are much stronger than those involved in STINGAQ‑S162A or STINGAQ‑E260I. In STINGAQ, strong correlations exist between the two monomers, as the variation of distances between two pairs of α1- and α3-helices (regions i and ii) is correlated.

Figure 4. Structural correlation maps: (A) DMXAA−STINGAQ, (B) DMXAA−STINGAQ‑S162A, and (C) STINGAQ‑E260I. 2077

DOI: 10.1021/acs.jpcb.6b12472 J. Phys. Chem. B 2017, 121, 2073−2082

Article

The Journal of Physical Chemistry B

Figure 5. (A) Interaction network. “a” and “b” represent protomer a and protomer b of the STING dimer, respectively. Circles with two peripheries represent cross-protomer interactions; (B) location of the corresponding residues in the protein structure.

STINGAQ‑E260I becomes even weaker. As shown in Table 3, in STINGAQ, DMXAA, 162S*, and 267T all can form hydrogen

pocket of the wild-type STING are more integrated into structural change upon ligand binding than the other two. We also investigated changes in the thermodynamic properties of pocket water molecules. The enthalpy was estimated by the interaction energy through MD simulations. The 2PT model was used to calculate water entropy.32,33,46 In these calculations, water molecules within 5 Å from both 167TYR and 262ALA on the two chains are considered as pocket waters. As listed in Table 4, the average interaction energies per pocket water molecule of STINGAQ, STINGAQ‑S162A, and STINGAQ‑E260I are −8.5, −7.1, and −6.5 kcal/mol, respectively. The interaction energy of a water molecule in bulk water is 3.4 kcal/mol, which is in good agreement with experimental measurements47 and ab initio calculations.48,49 The pocket waters in all three systems have lower entropies (per water molecule) than bulk water (Table 4). This is not surprising because in the bulk phase, water molecules translate with a less steric hindrance compared to those surrounded by protein residues and reorient freely through concerted large-amplitude angular jumps.50,51 As listed in Table 4, the entropic penalty of individual water molecules decreases in the order STINGAQ, STINGAQ‑S162A, and STINGAQ‑E260I, which is 3.3, 3.1, and 2.9 kcal/mol, respectively. Such a result reveals that the water molecules are constrained to varied degrees in the different binding pockets. In addition, both STING mutants in the bound state accommodate fewer pocket waters than STINGAQ, whereas in the apo state, the wild type contains the fewest pocket waters. The contributions of the enthalpic and entropic changes to the binding free energy are estimated through eqs 12 and 13, respectively

Table 3. Percentage of Hydrogen-Bond Formation among DMXAA, 162S*, and 267T system

interactive pairs

STINGAQ

STINGAQ‑E260I

interaction group A

DMXa−162Sb DMXa−267Ta 267Ta−162Sb DMXb−162Sa DMXb−267Tb 267Tb−162Sa

0.30 0.50 0.33 0.08 0.89 0.01

0.94 0.84 0.02 0.89

interaction group B

bonds with each other, whereas these types of hydrogen-bond interactions do not form in STINGAQ‑E260I. In the latter, strong interactions of DMXAA with residues 238R and 232R in the lid region and at top of the binding pocket 169R, 167Y, and 163Y (Table 1) deprive the hydrogen-bond interaction between DMXAA, 162S*, and 267T at the bottom of the binding pocket. Changes Induced by Mutations to Pocket Water Molecules. Water molecules in the binding pocket have significant influences on protein−ligand binding.42−44 Even subtle structural and dynamical changes in protein are able to alter the water network in the binding pocket and thus influence the response of a protein to its ligands.45 In this section, we investigate how pocket water molecules are influenced by mutations. First, we examined whether mutations would alter the interactions between water molecules and STING. We calculated correlation functions between the number of nearby waters of residues (more details in Table S3) and the distance between pairs of residues mentioned earlier. In these analyses, water molecules within 5 Å of a residue are considered. As illustrated in Figure 6, in STINGAQ, the numbers of water molecules around 260Ea and 166Gb correlate relatively strongly with a number of residue−residue and residue−ligand distances. In STINGAQ‑S162A, only water molecules around 260Eb showed weak correlations, whereas in STINGAQ‑E260I, no distinct correlation can be observed. In this sense, the water molecules in the binding

ΔH = (Napo − Nbound)*(E bulk − Epro)

(12)

−T ΔS = −(Napo − Nbound)*T (S bulk − Spro)

(13)

where N is the number of pocket water molecules, E is the interaction energy per water molecule, and S is the entropy per water molecule. The subscripts pro and bulk are used to indicate water molecules in the binding pocket and in bulk water, respectively. The enthalpy changes (ΔH) of water molecules upon ligand binding for STINGAQ, STINGAQ‑S162A, and 2078

DOI: 10.1021/acs.jpcb.6b12472 J. Phys. Chem. B 2017, 121, 2073−2082

Article

The Journal of Physical Chemistry B

STINGAQ‑E260I are 20.4, 62.9, and 62.0 kcal/mol, respectively. The corresponding entropy changes (−TΔS) are −4.8, −23.8, and −32.0 kcal/mol, respectively.



DISCUSSION Guided by free-energy decomposition calculations on DMXAA− STINGAQ‑S162A, we identified a new mutation, E260I, endowing STINGAQ with DMXAA sensitivity. It was also found that two crucial residues, arginine 238, dominate DMXAA binding, which means the electrostatic interactions are essential for DMXAA binding. In addition, the contributions of 238R to the interaction energy are distinctively different in the interaction of DMXAA with STINGAQ, STINGAQ‑S162A, and STINGAQ‑E260I, although there is only one residue different in the sequence of the three complexes. On the basis of structural correlation analysis, we identified pathways that affect this long-distance effect. DMXAA−STINGAQ showed a stronger structural correlation than the mutant systems. In STINGAQ, not only distances between residues that directly interact with DMXAA but also distances between helices constituting the binding pocket are correlated. Two pathways mediate this cross-protomer correlation: (i) the hydrogen-bond interaction among DMXAA, 162S*, and 267T; (ii) salt bridge and hydrogen bond between 238R* and 260E. The S162A mutation eliminates the DMXAA− 162S*−267T pathway directly and DMXAA−STINGAQ‑S162A shows a weak structural correlation. The E260I mutation eliminates the 238R*−260E pathway and disrupts the second pathway, DMXAA−162S*−267T, as DMXAA only forms hydrogen bond with 162S* but not with 267T (Table 3). Therefore, there is no distinct structural correlation in DMXAA− STINGAQ‑E260I. That is how a single mutation is able to weaken the structural correlations of the ligand−substrate interaction on a large scale. Water molecules in the binding pocket of STING and its two mutants also present different thermodynamic properties. The translational and rotational properties of water molecules, characterized by entropies, can be very different in various environments. For example, the entropy penalty is ∼0.1 kcal/mol for a water molecule on the lignin surface35 and ∼0.7 kcal/mol for a water at the interface of the lipid bilayer and di-mannose,36 and for major and minor groove of DNA waters, the values are ∼0.5 and ∼0.9 kcal/mol, respectively.37 In our systems, the corresponding entropy penalty is 3.3, 3.1, and 2.9 kcal/mol in the binding pocket of STINGAQ, STINGAQ‑S162A, and STINGAQ‑E260I, respectively. From Table 4, one sees that the change in number of pocket water molecules upon ligand binding is deterministic of the total entropic change. Although entropy penalties of individual pocket water molecule are similar among the three types of STING, our simulation revealed that the two mutants expel larger numbers of pocket water than STINGAQ upon binding, which are 4, 17, and 20 for STINGAQ, STINGAQ‑S162A, and STINGAQ‑E260I, respectively, on average. Consequently, the total entropic penalties of mutants are much larger than that of STINGAQ, greater than ∼20 kcal/mol. Meanwhile, the total enthalpic changes compensate for these entropic penalties. Numerous physical mechanisms have been proposed to explain enthalpy−entropy compensation in ligand−protein interactions.52 WaterMap53 and other statistical mechanical models54,55 have been used to demonstrate that solvent reorganization is one of the origins of this compensation. Here, we showed that the change in number of pocket water molecules upon ligand binding is a factor that contributes to enthalpy−entropy compensation.

Figure 6. Correlation map of water numbers (#35−67) and distances: (A) DMXAA−STINGAQ, (B) DMXAA−STINGAQ‑S162A, and (C) STINGAQ‑E260I.

2079

DOI: 10.1021/acs.jpcb.6b12472 J. Phys. Chem. B 2017, 121, 2073−2082

Article

The Journal of Physical Chemistry B Table 4. Thermodynamic Properties of Pocket Water and Bulk Water at 300 K Napo Nbound E (kcal/mol) ΔH (kcal/mol) TS (kcal/mol) −TΔS (kcal/mol)

STINGAQ

STINGAQ‑S162A

STINGAQ‑E260I

24 ± 2 20 ± 0 −8.5 ± 0.0 20.4 3.3 ± 0.1 −4.8

33 ± 3 16 ± 1 −7.1 ± 0.0 62.9 3.1 ± 0.1 −23.8

29 ± 4 9±0 −6.5 ± 0.1 62.0 2.9 ± 0.0 −32.0

4.5 ± 0.0

Yi Qin Gao: 0000-0002-4309-9376 Author Contributions §

X.C. and X.-X.D. contributed equally to this work.

Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS The authors thank the National Natural Science Foundation of China for support (21573006, U1430237, and 21373016). We also thank Peking University for providing the computational resources.



REFERENCES

(1) Burdette, D. L.; Vance, R. E. STING and the Innate Immune Response to Nucleic Acids in the Cytosol. Nat. Immunol. 2013, 14, 19− 26. (2) Burdette, D. L.; Monroe, K. M.; Sotelo-Troha, K.; Iwig, J. S.; Eckert, B.; Hyodo, M.; Hayakawa, Y.; Vance, R. E. STING is a Direct Innate Immune Sensor of Cyclic Di-GMP. Nature 2011, 478, 515−518. (3) Ishikawa, H.; Ma, Z.; Barber, G. N. STING Regulates Intracellular DNA-Mediated, Type I Interferon-Dependent Innate Immunity. Nature 2009, 461, 788−792. (4) Zhong, B.; Yang, Y.; Li, S.; Wang, Y. Y.; Li, Y.; Diao, F. C.; Lei, C. Q.; He, X.; Zhang, L.; Tien, P.; et al. The Adaptor Protein MITA Links Virus-Sensing Receptors to IRF3 Transcription Factor Activation. Immunity 2008, 29, 538−550. (5) Keating, S. E.; Baran, M.; Bowie, A. G. Cytosolic DNA Sensors Regulating Type I Interferon Induction. Trends Immunol. 2011, 32, 574−581. (6) Paludan, S. R.; Bowie, A. G. Immune Sensing of DNA. Immunity 2013, 38, 870−880. (7) Huang, Y.-H.; Liu, X.-Y.; Du, X.-X.; Jiang, Z.-F.; Su, X.-D. The Structural Basis for the Sensing and Binding of Cyclic Di-GMP by STING. Nat. Struct. Mol. Biol. 2012, 19, 728−730. (8) Gao, P.; Ascano, M.; Zillinger, T.; Wang, W. Y.; Dai, P. H.; Serganov, A. A.; Gaffney, B. L.; Shuman, S.; Jones, R. A.; Deng, L.; et al. Structure-Function Analysis of STING Activation by c[G(2′,5′) pA(3′,5′)p] and Targeting by Antiviral DMXAA. Cell 2013, 154, 748−762. (9) Zhang, X.; Shi, H. P.; Wu, J. X.; Zhang, X. W.; Sun, L. J.; Chen, C.; Chen, Z. J. J. Cyclic GMP-AMP Containing Mixed Phosphodiester Linkages is an Endogenous High-Affinity Ligand for STING. Mol. Cell 2013, 51, 226−235. (10) Cavlar, T.; Deimling, T.; Ablasser, A.; Hopfner, K. P.; Hornung, V. Species-Specific Detection of the Antiviral Small-Molecule Compound CMA by STING. EMBO J. 2013, 32, 1440−1450. (11) Baguley, B. C.; Ching, L. M. DMXAA: An Antivascular Agent with Multiple Host Responses. Int. J. Radiat. Oncol., Biol., Phys. 2002, 54, 1503−1511. (12) Lara, P. N., Jr.; Douillard, J. Y.; Nakagawa, K.; von Pawel, J.; McKeage, M. J.; Albert, I.; Losonczy, G.; Reck, M.; Heo, D. S.; Fan, X. L.; et al. Randomized Phase III Placebo-Controlled Trial of Carboplatin and Paclitaxel With or Without the Vascular Disrupting Agent Vadimezan (ASA404) in Advanced Non-Small-Cell Lung Cancer. J. Clin. Oncol. 2011, 29, 2965−2971.



CONCLUSIONS The rational design of mutants based on computational analysis led to the discovery of a new mutation, E260I, that endows STING with DMXAA sensitivity. Further detailed analysis facilitated an understanding of how subtle changes in the protein structure influence the DMXAA−STING interaction. On the one hand, single mutations reshape the structural correlation network involved in the DMXAA−STING interaction, thereby exerting an influence on distant interactions. On the other hand, single mutations disturb the coordinated motions of water molecules and change how much water is expelled upon ligand binding. Moreover, molecular-level analyses shed light on drug optimization. We expect that DMXAA derivatives with no interaction with 162S* and 267T are candidate compounds for binding human STING. ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.6b12472. Details of all designed mutations, IFN-β expression of STINGAQ double mutations, the values of contributions of key residues to binding free energy, and collective variables used in structural correlation analyses (PDF)



−3.4 ± 0.0

ORCID

Because both mutants possessing DMXAA responses show weaker structural correlations, it might be reasonable to speculate that the strong structural correlations in STINGAQ weaken the interaction between DMXAA and 238R and consequently discourage ligand binding. 260E and DMXAA competitively interact with 238R. Mutating 260E to isoleucine enhances the interaction between DMXAA and 238R. In DMXAA− STINGAQ, the correlation coefficient of distance between DMXAA and 162S* and that between DMXAA* and 238R* is −0.7. The interaction of DMXAA and 162S* weakens the interaction between DMXAA* and 238R*. Therefore, we hypothesize that the strong structural correlation of the complex hinders protein structure adjustment, especially the adjustment of 238R, upon ligand binding. On the basis of such an argument, cross-protomer interactions, such as hydrogen-bond interactions among DMXAA, 162S*, and 267T that intermediate the structural correlation, are unfavorable for DMXAA binding. Such a hypothesis is consistent with experimental observations: the double mutation STINGAQ‑S162A/T267A induces even higher IFN-β expressions than STINGAQ‑S162A, as shown by the luciferase assay (Figure S1).



bulk water

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected] (X.-D.S.). *E-mail: [email protected] (Y.Q.G.). 2080

DOI: 10.1021/acs.jpcb.6b12472 J. Phys. Chem. B 2017, 121, 2073−2082

Article

The Journal of Physical Chemistry B

with the Two-Phase Thermodynamic Model. Phys. Chem. Chem. Phys. 2012, 14, 15206−15213. (34) Pascal, T. A.; Goddard, W. A. Interfacial Thermodynamics of Water and Six Other Liquid Solvents. J. Phys. Chem. B 2014, 118, 5943− 5956. (35) Petridis, L.; Schulz, R.; Smith, J. C. Simulation Analysis of the Temperature Dependence of Lignin Structure and Dynamics. J. Am. Chem. Soc. 2011, 133, 20277−20287. (36) Tian, J.; Sethi, A.; Swanson, B. I.; Goldstein, B.; Gnanakaran, S. Taste of Sugar at the Membrane: Thermodynamics and Kinetics of the Interaction of a Disaccharide with Lipid Bilayers. Biophys. J. 2013, 104, 622−632. (37) Jana, B.; Pal, S.; Maiti, P. K.; Lin, S. T.; Hynes, J. T.; Bagchi, B. Entropy of Water in the Hydration Layer of Major and Minor Grooves of DNA. J. Phys. Chem. B 2006, 110, 19611−19618. (38) Sun, W.; Li, Y.; Chen, L.; Chen, H. H.; You, F. P.; Zhou, X.; Zhou, Y.; Zhai, Z. H.; Chen, D. Y.; Jiang, Z. F. ERIS, an Endoplasmic Reticulum IFN Stimulator, Activates Innate Immune Signaling through Dimerization. Proc. Natl. Acad. Sci. U.S.A. 2009, 106, 8653−8658. (39) Foda, Z. H.; Shan, Y. B.; Kim, E. T.; Shaw, D. E.; Seeliger, M. A. A Dynamically Coupled Allosteric Network Underlies Binding Cooperativity in Src Kinase. Nat. Commun. 2015, 6, No. 5939. (40) Gasper, P. M.; Fuglestad, B.; Komives, E. A.; Markwick, P. R. L.; McCammon, J. A. Allosteric Networks in Thrombin Distinguish Procoagulant vs. Anticoagulant Activities. Proc. Natl. Acad. Sci. U.S.A. 2012, 109, 21216−21222. (41) Allain, A.; de Beauchene, I. C.; Langenfeld, F.; Guarracino, Y.; Laine, E.; Tchertanov, L. Allosteric Pathway Identification through Network Analysis: from Molecular Dynamics Simulations to Interactive 2D and 3D Graphs. Faraday Discuss. 2014, 169, 303−321. (42) Lemmon, G.; Meiler, J. Towards Ligand Docking Including Explicit Interface Water Molecules. PLoS One 2013, 8, No. e67536. (43) Mikulskis, P.; Genheden, S.; Ryde, U. Effect of Explicit Water Molecules on Ligand-Binding Affinities Calculated with the MM/GBSA Approach. J. Mol. Model. 2014, 20, 2273. (44) Mancera, R. L. De Novo Ligand Design with Explicit Water Molecules: an Application to Bacterial Neuraminidase. J. Comput.-Aided Mol. Des. 2002, 16, 479−499. (45) Davidson, R. M.; Lauritzen, A.; Seneff, S. Biological Water Dynamics and Entropy: A Biophysical Origin of Cancer and other Diseases. Entropy 2013, 15, 3822−3876. (46) Debnath, A.; Mukherjee, B.; Ayappa, K. G.; Maiti, P. K.; Lin, S. -T Entropy and Dynamics of Water in Hydration Layers of a Bilayer. J. Chem. Phys. 2010, 133, No. 174704. (47) Rocher-Casterline, B. E.; Ch’ng, L. C.; Mollner, A. K.; Reisler, H. Communication: Determination of the Bond Dissociation Energy (D0) of the Water Dimer, (H2O)2, by Velocity Map Imaging. J. Chem. Phys. 2011, 134, No. 211101. (48) Shank, A.; Wang, Y.; Kaledin, A.; Braams, B. J.; Bowman, J. M. Accurate Ab Initio and “Hybrid” Potential Energy Surfaces, Intramolecular Vibrational Energies, and Classical IR Spectrum of the Water Dimer. J. Chem. Phys. 2009, 130, No. 144314. (49) Leforestier, C.; Szalewicz, K.; van der Avoird, A. Spectra of Water Dimer from a New Ab Initio Potential with Flexible Monomers. J. Chem. Phys. 2012, 137, No. 014305. (50) Laage, D.; Hynes, J. T. A Molecular Jump Mechanism of Water Reorientation. Science 2006, 311, 832−835. (51) Young, T.; Abel, R.; Kim, B.; Berne, B. J.; Friesner, R. A. Motifs for Molecular Recognition Exploiting Hydrophobic Enclosure in ProteinLigand Binding. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 808−813. (52) Chodera, J. D.; Mobley, D. L. Entropy-Enthalpy Compensation: Role and Ramifications in Biomolecular Ligand Recognition and Design. Annu. Rev. Biophys. 2013, 42, 121−142. (53) Breiten, B.; Lockett, M. R.; Sherman, W.; Fujita, S.; Al-Sayah, M.; Lange, H.; Bowers, C. M.; Heroux, A.; Krilov, G.; Whitesides, G. M. Water Networks Contribute to Enthalpy/Entropy Compensation in Protein-Ligand Binding. J. Am. Chem. Soc. 2013, 135, 15579−15584.

(13) Conlon, J.; Burdette, D. L.; Sharma, S.; Bhat, N.; Thompson, M.; Jiang, Z. Z.; Rathinam, V. A. K.; Monks, B.; Jin, T. C.; Xiao, T. S.; et al. Mouse, but not Human STING, Binds and Signals in Response to the Vascular Disrupting Agent 5,6-Dimethylxanthenone-4-Acetic Acid. J. Immunol. 2013, 190, 5216−5225. (14) Gao, P.; Zillinger, T.; Wang, W. Y.; Ascano, M.; Dai, P. H.; Hartmann, G.; Tuschl, T.; Deng, L.; Barchet, W.; Patel, D. J. BindingPocket and Lid-Region Substitutions Render Human STING Sensitive to the Species-Specific Drug DMXAA. Cell Rep. 2014, 8, 1668−1676. (15) Henzler-Wildman, K.; Kern, D. Dynamic Personalities of Proteins. Nature 2007, 450, 964−972. (16) Bahar, I.; Lezon, T. R.; Yang, L. W.; Eyal, E. Global Dynamics of Proteins: Bridging Between Structure and Function. Annu. Rev. Biophys. 2010, 39, 23−42. (17) Eisenmesser, E. Z.; Bosco, D. A.; Akke, M.; Kern, D. Enzyme Dynamics during Catalysis. Science 2002, 295, 1520−1523. (18) McGowan, L. C.; Hamelberg, D. Conformational Plasticity of an Enzyme during Catalysis: Intricate Coupling between Cyclophilin A Dynamics and Substrate Turnover. Biophys. J. 2013, 104, 216−226. (19) Bhattacharya, S.; Vaidehi, N. Differences in Allosteric Communication Pipelines in the Inactive and Active States of a GPCR. Biophys. J. 2014, 107, 422−434. (20) Kumar, A.; Glembo, T. J.; Ozkan, S. B. The Role of Conformational Dynamics and Allostery in the Disease Development of Human Ferritin. Biophys. J. 2015, 109, 1273−1281. (21) Rose, R. B.; Craik, C. S.; Stroud, R. M. Domain Flexibility in Retroviral Proteases: Structural Implications for Drug Resistant Mutations. Biochemistry 1998, 37, 2607−2621. (22) Ladbury, J. E. Just Add Water! The Effect of Water on the Specificity of Protein-Ligand Binding Sites and Its Potential Application to Drug Design. Chem. Biol. 1996, 3, 973−980. (23) Papoian, G. A.; Ulander, J.; Wolynes, P. G. Role of Water Mediated Interactions in Protein-Protein Recognition Landscapes. J. Am. Chem. Soc. 2003, 125, 9170−9178. (24) Materese, C. K.; Goldmon, C. C.; Papoian, G. A. Hierarchical Organization of Eglin C Native State Dynamics is Shaped by Competing Direct and Water-Mediated Interactions. Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 10659−10664. (25) Karplus, M.; McCammon, J. A. Molecular Dynamics Simulations of Biomolecules. Nat. Struct. Biol. 2002, 9, 646−652. (26) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 09; Gaussian, Inc.: Wallingford, CT, 2009. (27) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157−1174. (28) Miller, B. R., 3rd; McGee, T. D., Jr.; Swails, J. M.; Homeyer, N.; Gohlke, H.; Roitberg, A. E. MMPBSA.py: An Efficient Program for EndState Free Energy Calculations. J. Chem. Theory Comput. 2012, 8, 3314− 3321. (29) Gohlke, H.; Kiel, C.; Case, D. A. Insights into Protein-Protein Binding by Binding Free Energy Calculation and Free Energy Decomposition for the Ras-Raf and Ras-RaIGDS Complexes. J. Mol. Biol. 2003, 330, 891−913. (30) Metz, A.; Pfleger, C.; Kopitz, H.; Pfeiffer-Marek, S.; Baringhaus, K. H.; Gohlke, H. Hot Spots and Transient Pockets: Predicting the Determinants of Small-Molecule Binding to a Protein-Protein Interface. J. Chem. Inf. Model. 2012, 52, 120−133. (31) Lin, S. T.; Blanco, M.; Goddard, W. A. The Two-Phase Model for Calculating Thermodynamic Properties of Liquids from Molecular Dynamics: Validation for the Phase Diagram of Lennard-Jones Fluids. J. Chem. Phys. 2003, 119, 11792−11805. (32) Lin, S. T.; Maiti, P. K.; Goddard, W. A. Two-Phase Thermodynamic Model for Efficient and Accurate Absolute Entropy of Water from Molecular Dynamics Simulations. J. Phys. Chem. B 2010, 114, 8191−8198. (33) Lai, P. K.; Hsieh, C. M.; Lin, S. T. Rapid Determination of Entropy and Free Energy of Mixtures from Molecular Dynamics Simulations 2081

DOI: 10.1021/acs.jpcb.6b12472 J. Phys. Chem. B 2017, 121, 2073−2082

Article

The Journal of Physical Chemistry B (54) Grunwald, E.; Steel, C. Solvent Reorganization and Thermodynamic Enthalpy-Entropy Compensation. J. Am. Chem. Soc. 1995, 117, 5687−5692. (55) Lee, B.; Graziano, G. A Two-State Model of Hydrophobic Hydration that Produces Compensating Enthalpy and Entropy Changes. J. Am. Chem. Soc. 1996, 118, 5163−5168.

2082

DOI: 10.1021/acs.jpcb.6b12472 J. Phys. Chem. B 2017, 121, 2073−2082