Addressing Polarization Phenomena in Molecular Machines

Subscriber access provided by AUT Library

Molecular Mechanics

Addressing Polarization Phenomena in Molecular Machines Containing Transition Metal Ions with an Additive Force Field Shuangli Du, Haohao Fu, Xueguang Shao, Christophe Chipot, and Wensheng Cai J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00972 • Publication Date (Web): 24 Jan 2019 Downloaded from http://pubs.acs.org on January 29, 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 26 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

Journal of Chemical Theory and Computation

Addressing Polarization Phenomena in Molecular Machines Containing Transition Metal Ions with an Additive Force Field Shuangli Du, †,○ Haohao Fu, †,○ Xueguang Shao, †,‡,§, ‖ Christophe Chipot, ┴, #,∇ and Wensheng Cai*,†,‡, ‖ †

Research Center for Analytical Sciences, College of Chemistry, Nankai University, Tianjin

300071, China ‡

Tianjin Key Laboratory of Biosensing and Molecular Recognition, Tianjin 300071, China

§

State Key Laboratory of Medicinal Chemical Biology, Tianjin 300071, China

‖

Collaborative Innovation Center of Chemical Science and Engineering, Tianjin, 300071, China

┴ #

œuvre-les-Nancy F-54500, France y

L

Vandœuvre-les-Nancy F-54506, France ∇D

P y

y

W

G

Urbana, Illinois 61801, United States ○

Contributed equally to this work

Corresponding Author * E-mail: [email protected]

ACS Paragon Plus Environment

1

Journal of Chemical Theory and Computation 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 2 of 26

ABSTRACT: Modeling transition metals in supramolecular assemblies, in general, is extremely challenging due to polarization and charge transfer. In this paper, we demonstrate that the inherent shortcomings of additive force fields in modeling Cu+–ether-O and Cu+–olefin-C interactions are rooted in the Lorentz-Berthelot rules. A general method for investigating transition-metal-containing molecular assays using classical force fields is, therefore, proposed. In this strategy, QM/MM calculations have been performed to determine the potential of mean force (PMF) describing the interaction of a cation and a specific functional group. van der Waals parameters for the corresponding pairs of particles have then been optimized using the NBFIX feature of the CHARMM force field to fit the QM/MM PMF. This method has been applied to decipher the mechanism underlying the ―dialing‖ of a molecular machine controlled by Li+ and Cu+ cations, indicating that the process is controlled by the competition between cation–ether-O and cation–olefin-C interactions.


2

Page 3 of 26 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


1. INTRODUCTION Molecular dynamics (MD) simulations have been playing a crucial role in the exploration of a host of chemical and biochemical processes at the atomic level.1, 2 Metal elements, occupying a sizable component of the periodic table, take part in a huge number of chemical reactions.3-5 Not too surprisingly, about one-third of the molecular assemblies in the RCSB protein data bank (https://www.rcsb.org/) contain metal ions.6 Modeling with utmost accuracy metal ions, especially transition metals, using classical molecular force fields, however, represents a longstanding challenge due to the polarization, charge-transfer and covalent-bonding effects between cations and nucleophiles. So far, only a handful of methods have been proposed7-9 to model metal ions at different levels of accuracy in classical-mechanics-based computations. The most common model, described by Lennard-Jones (LJ) plus Columbic potentials, remains widely employed in classical force fields, e.g., CHARMM,10-12 Amber,13-15 and OPLS-AA,16-18 and relies on combination rules to calculate the LJ parameters associated to pairs of atoms. This model is, however, believed to be inaccurate in describing transition metal-based complexes because it grossly ignores polarization, charge-transfer and covalent-bonding effects. Recently, Li and Merz found that an LJ non-bonded model could not reproduce both the correct hydration free energy and ion-to-water distance of divalent metal ions in TIPnP water.19 They further put forth an improved LJ model, namely a 12-6-4 potential,19,

20

accounting for the

ion-induced dipole interaction by introducing an r-4 term, which was originally put forth by Minoux and Chipot in the context of cation-π

.21 This ad-hoc model successfully

describes metal ion-water interactions, though, in principle, it is unable to handle charge-transfer and covalent-bonding effects.


3


Page 4 of 26

Polarizable force fields can address some of the shortcomings of standard 12-6 models.22 For example, both the CHARMM Drude model23 and AMOEBA force field24 take into account induction effects and mimic a fraction of charge transfer.25 However, in stark contrast with additive models like CGENFF,26 GAFF27 and OPLS-AA,16 which are well-suited for the parametrization of most organic compounds, polarizable force fields are difficult to parameterize, though noteworthy efforts like those of Huang et al.28 and Wu et al.,29 have provided general strategies,

namely

GAAMP

(http://gaamp.lcrc.anl.gov/)

and

Poltype

(https://github.com/pren/poltype), for the parameterization of small, possibly drug-like molecules for the Drude and the AMOEBA force field, respectively. Aside from additive and non-additive force fields, one may turn to bonded models, wherein the metal center is explicitly bonded to its adjacent atoms. Automatic parameterization can be achieved, among others, by the Metal Center Parameter Builder in Ambertools.30 Bonded models have proven suitable to describe the fluctuation of the cation position at distinct binding sites.6 Using such a model, however, precludes the coordination number of the ion to change, which constitutes a major drawback. Dummy-atom model,31 involving fictitious particles connected to the metal center that is not directly trapped by the protein, may be regarded as a compromise between bonded and non-bonded potentials. In light of the above discussion about different types of force fields,28, 29 but also the recent developments in the AMOEBA24, 25, 29 and SIBFA32, 33 force fields, we have developed in this paper a generalized method that can implicitly describe polarization and possibly mimic part of charge transfer and covalent-bonding for pairs of particles involving a metal ion, using a standard 12-6 potential. The corresponding potential terms will be used as corrections to standard force fields, also referred to as NBFIX (Non-Bonded FIX) corrections.34,

35

Progress in the


4

Page 5 of 26 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


NBFIX corrections to the AMBER and CHARMM force fields has been made,36-40 improving description of nucleic acid41, 42 and protein folding.43 The NBFIX approach only optimizes the LJ parameters for specific atom pairs without modifying the other interactions. Generally, these LJ parameters, also termed NBFIX parameters, are optimized to reproduce the osmotic pressure data for binary solutions.35 In the present work, we have adopted an alternative strategy – the potential of mean force (PMF) between the ion and the specific group determined by QM/MM calculations is used as a reference in the optimization procedure. Because a PMF provides more atomic-level information about the interacting particles than a single experimental measurement, this strategy is extremely powerful to determine the atom-pair-based force-field parameters. As an example, the above method was used in the parameterization of the Cu+ ion. A cation-responsive rotaxane reported by Baggi and Loeb44 was then studied using the optimized NBFIX parameters. In experiment, the crown-ether macrocycle containing an olefinic group was observed to rotate about the axle, forming different spatial arrangements upon addition of Li+ and Cu+. Because the Cu+–olefin and Cu+–ether interactions cannot be well described by standard force fields, this molecular machine was selected as a test case. Our quantum-mechanical calculations demonstrate that pronounced polarization, charge-transfer and covalent-bonding effects (see the Supporting Information) coexist, in particular, between Cu+ and olefin. The LJ parameters of Cu+–olefin-C and Cu+–ether-O were optimized the strategy put forth herein. The rotation process in the rotaxane was then investigated applying a standard 12-6 model,45 a 12-6-4 potential19, 45 and a standard model augmented by NBFIX parameters, to examine the reliability of our strategy, and to further decipher the mechanism underlying the dialing process. 2. METHODS


5


Page 6 of 26

In a classical force field, the metal ion interacts with adjacent atoms through standard 12-6 and Columbic potentials,19 ∑ *

((

)

(

) )+,

where i is the metal ion, j, the adjacent atoms, kC, the electrostatic constant, q, the partial charge, rij, the distance between metal ion i and atom j, and ε and rmin, the LJ parameters. The εij and rmin,ij are usually given, employing combination rules. In the CHARMM and Amber force fields, the Lorentz-Berthelot rules are used, √

, .

In most cases, εi and rmin,i of a metal ion in classical force fields are parameterized from simulations in water.46-49 With the LJ parameters optimized based on hydration free energy and ion-to-oxygen distance, alkali metal ions, usually, can be described with the desired accuracy, but for other cations interacting with molecules other than water, a number of properties, in particular, their binding affinity with respect to Lewis bases, cannot be correctly reproduced.35 To this end, the general idea of our strategy is to x

y

εij and rmin,ij for specific atom

pairs that cannot be suitably described by the traditional combining rules, while keeping the parameters for other atomic pairs unchanged. This result can be achieved using the NBFIX35, 40, 42

feature of the CHARMM force field. The PMF characterizing the interaction of a metal ion

with a specific functional group is accurately determined by QM/MM simulations. The parameters (εij and rmin,ij) for the specific ion-atom pairs are then optimized from a fit of the PMF using the standard force field plus NBFIX corrections to the reference QM/MM result. In principle, the optimized parameters can, therefore, approximately describe the average


6

Page 7 of 26 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


polarization effects due to the specific atom pairs. It is worth noting that although the partial charges borne by the metal ions and the atoms interacting with them are not reoptimized, and, thus, remain constant, part of the charge-transfer effect that may occur is implicitly handled by the optimized NBFIX parameters. The workflow of our strategy is summarized in Figure 1A. As a test case, the LJ parameters describing the Cu+–olefin-carbon and Cu+–ether-oxygen interactions were optimized following this strategy, and were gathered in Table 1. 2.1 Molecular Models. Two simplified molecular models (Cu+–olefin and Cu+–ether) were constructed, as described in Figure 1. The two molecular assemblies were then immersed independently in a box of TIP3P50 water. The initial dimensions of the simulation cells were 30.8  30.6  30.6 Å3, featuring 1031 water molecules, and 32.5  30.9  31.9 Å3, featuring 1109 water molecules, respectively. 2.2 Free-energy Calculations. In our QM/MM calculations, electrical embedding was adopted to treat Coulombic interactions at the QM/MM interface.51,

52

All the solutes were

treated using quantum-chemical theory without applying any constraint, while all the solvent molecules were modeled by means of a standard, pairwise additive force field (see Figure S3). The PM7 semi-empirical method53 was employed in all QM/MM simulations. A comparison of the PM7 accuracy against the DFT method was made, as shown in Figure S4 of the Supporting Information. MOPAC54 and NAMD55 were used as the QM and MD engines, respectively. The recently developed meta-eABF method56 available in Colvars (the collective-variable module associated to NAMD)57 was employed in all the PMF calculations. The transition coordinate, , was defined as the distance between Cu+ and the center of mass of the ethane carbon atoms or the ether oxygen atom. To increase the efficiency of the calculations, the transition pathway, extending from 1.8 to 12 Å, was broken down into five windows. Instantaneous values of the


7


Page 8 of 26

force were accrued in bins 0.1 Å wide. A 50-ns QM/MM simulation for each molecular assembly was carried out with a 0.5-fs time step. To probe the convergence of the free-energy calculation, the time evolutions of the root-mean-square deviation over the gradients and the free-energy profiles at different simulation times are shown in Figure S5 and S6, respectively, suggesting suitable convergence of the QM/MM PMF calculations. Additional methodological description for the classical MD simulations and free-energy calculations of the rotaxanes can be found in the Supporting Information.

Figure 1. (A) Workflow of the parameterization strategy. (B) Comparison between the free-energy profiles of the Cu+–olefin and Cu+–ether interactions obtained by simulations using


8

Page 9 of 26 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


QM/MM potential, standard 12-6 potential from Li/Merz,45 12-6 potential plus NBFIX corrections derived from this work and 12-6-4 potential from Li/Merz.45

3. RESULTS AND DISCUSSION 3.1. LJ Parameters for the Interactions of Cu+–Olefin-Carbon and Cu+–Ether-Oxygen. As shown in Figure 1B, neither the Li/Merz 12-6 nor 12-6-4 model45 can suitably describe the binding affinities of Cu+ with olefin and ether. In sharp contrast, the PMFs obtained by the standard 12-6 model45 plus the NBFIX corrections can reproduce both the positions and the depths of the global minima of the reference QM/MM PMFs, which chiefly determine the length and strength of the coordination bonds. The fit between the two PMFs for the Cu+–ether interaction is suboptimal at Cu+–O distances ranging between 3 and 6 Å, due to polarization in the dehydration of the first solvation shell of Cu+. Since electrostatic embedding was adopted in this work, the quantum atoms (Cu+, olefin and ether moiety, herein) cou

“

”

borne by the classical ones (the surrounding water molecules), and, thus, were polarized by these charges. However, in the NBFIX approach, we surgically corrected the pair-specific Lennard-Jones (LJ) parameters for the solutes without modifying the water model. The polarization contribution to the first-shell dehydration is, therefore, not considered in our PMF calculations with the classical force field, even with the optimized Cu+–ether-O parameters, thus, leading to the discrepancy between the free-energy profiles obtained from QM/MM and classical MD simulations. However, the influence of such polarization on the binding affinity of the ion with respect to ether/alkene, as well as on the most stable arrangements of the atomic pairs is almost negligible. In other words, the strong coordination bond coupled with polarization and


9


Page 10 of 26

charge transfer in molecular objects featuring Cu+–olefin and Cu+–ether interactions can be suitably described using the present approximate strategy. Table 1. Lennard-Jones parameters εmin (kcal/mol) and rmin (Å) for relevant pairs modified by NBFIX. ion-atom pair

εmin

rmin

Cu+–olefin-C

-24.10000

2.000

Cu+–ether-O

-21.50000

2.400

3.2. Investigation of a Cation-Controlled Rotaxane using the NBFIX Correction Strategy. To further test the applicability of the optimized NBFIX parameters and the effectiveness of our parameterization strategy, a cation-controlled rotaxane44 formed by an H-shaped axle containing a bidentate N,N'-chelate and a 24-membered crown-ether macrocycle possessing six ether O-atoms and one olefinic group was chosen as an example. It is reported in the experiment44 that the rotaxane can dial-up different donor sets for complexation to metal ions by rotating the macrocycle about the axle after addition of the different ions (Li+ or Cu+), as w

F u

2 ‒ . Cu+–olefin and Cu+–ether interactions are considered to be extremely

important in the dialing process. We, therefore, try to use the standard force field combined with NBFIX corrections, using the parameters of Table 1 to simulate this unique movement, and further explain the mechanism that underlies dialing. In our preliminary simulations of the rotaxane, two main types of movements were observed, namely, the rotation of the macrocycle about the axle and the relative rotation of the two T-shaped bulky stoppers. Two coarse variables, θ1 and θ2, describing these rotary movements (as depicted in Figure 2D) were, therefore, chosen to form the transition coordinate.


10

Page 11 of 26 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


Figure 2. Conformations, namely, [2]rotaxane without any ions (A), Li+–containing (B) and Cu+–containing (C) complexes, obtained by coordinating to a metal ion via rotation of the y

w

. (D) D ﬁ

. Tw

y

(θ1, θ2) were

utilized to explore the putative transition pathways. θ1 denotes the rotational angle of the macrocyclic ring around the axle, and θ2, the rotation of one T-shaped bulky stopper relative to the other one around the axle. Figure 3 shows the two-dimensional free-energy landscapes delineating the rotation of the macrocycle coupled with the relative rotation of the two T-shaped bulky stoppers in three different environments, namely, in the absence of a metal ion, and in the presence of Li+ and of Cu+. As shown in Figure 3A, five local minima are identified and can be categorized into two classes (1 and 2) according to their free energies in the absence of metal ion. The structures of class 1 (site 1, site 2 and site 5) are more stable than those of class 2 (site 3 and site 4). The difference in free energy in each class is almost negligible due to the symmetric structure of the [2]rotaxane (structural detail provided in Figure S7). Only the structures at site 1 and site 4 are, therefore, discussed below. In the structure at site 1, the axle is almost centrosymmetric, wherein two hydrogen bonds form between the two imino groups and the first and fifth ether oxygen


11


Page 12 of 26

atoms, in agreement with the reported X-ray crystal structure.44 Besides, in the structure corresponding to the local minimum at site 4, the axle is axisymmetric, with a hydrogen bond formed between the axle and the ring. The free-energy profile along the least free-energy pathway determined using the LFEP algorithm58 is shown in Figure S8. The free-energy differences between the class 1 (stable states) and the class 2 (metastable states) are on the order of kBT. Moreover, the free-energy barrier against rotation from site 1 to site 5 amounts to 4.2 kcal/mol, which is sufficiently small to allow the ring to rotate about the axle. The states corresponding to the five local minima may, therefore, coexist. In summary, the standard parameter sets of the 12-6 model can properly describe the stable structure of the [2]rotaxane in the absence of a metal ion.

Figure 3. Free-energy landscapes characterizing the configurational change of (A) [2]rotaxane without any ions, (B) Li+–containing and (C) Cu+–containing complexes in CHCl3, obtained from 1.21-, 3.15-, and 4.61-μ meta-eABF simulations, respectively. The first two cases using the standard parameter sets of 12-6 model. The Cu+-encapsulated [2]-rotaxane using parameters with NBFIX terms obtained from our method.


12

Page 13 of 26 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


Upon addition of Li+, two pronounced local minima, i.e., site 6 and site 7, of nearly equal free energy, appear at (-87.5º, -42.5º) and (-72.5º, 17.5º), respectively. The corresponding structures are provided in Figure S9A and B of the Supporting Information. As can be seen, the T-shaped bulky stopper rotates by about 160º, compared with the structure at site 1, forming an axisymmetric axle. For the structure corresponding to site 6, the Li+ is chelated by two donor oxygen atoms of the macrocycle and two nitrogen atoms from the axle, forming a distorted tetrahedral geometry. In the structure at site 7, three donor oxygen atoms of the macrocycle and two nitrogen atoms of the axle are coordinated to Li+, thereby forming a distorted trigonal bipyramidal geometry structure. Both stable states are observed in the experiment by Baggi et al.44 The corresponding structures and Li-O bond lengths are consistent with those of the X-ray crystal (see Table S3 in the Supporting Information), which indicates that the complex containing Li+ ions can be described correctly using the standard parameter sets. In addition, it can be seen that the macrocycle rotates by nearly 15ºabout the axle from site 1 to site 7, “

”

( s shown in Figure 4B).

As mentioned above, due to polarization and charge-transfer effects in the supramolecular arrangements featuring Cu+–olefin and Cu+–ether interactions, dialing may not be well described using the standard parameter sets. For comparison purposes, two different parameter sets describing these two specific atomic pairs, namely, the standard parameter sets and the optimized NBFIX parameters in Table 1, were used to determine the corresponding free-energy landscapes. The result obtained using the standard parameter sets is provided in the Supporting Information (Figure S10). Comparing Figure S10 and Figure 3B, it is apparent that the difference in the free-energy landscapes between the complexes containing Cu+ and Li+ ions is subtle, which is


13


Page 14 of 26

also reflected by the representative structures shown in Figure S10B and C. This result demonstrates that Cu+ and Li+ have the same binding site, in other words,

“

”

u ,

which is inconsistent with experiment. It can, therefore, be concluded that the standard parameter sets cannot be used alone to model complexes containing Cu+ ions.

Figure 4. The dialing process for the [2]rotaxane: (A) [2]rotaxane without any ions, (B) after addition of Li+, (C) upon an addition of Cu+. The corresponding positions on the PMF shown in Figure 3 are site 1, site 7 and site 8. Only H-atoms involved in benzimidazole are shown. However, a significant difference exists between the free-energy landscapes determined using the standard force field (Figure S10) and that augmented by NBFIX parameters (Figure 3C). In Figure 3C, one obvious lowest-energy point appears at θ1= 82.5º, θ2= 17.5º, indicating a sharp rotation of the macrocycle. The representative structure is shown in Figure S9C, wherein the Cu+ cation is coordinated to the olefin moiety of the macrocycle, and chelated by the nitrogen atoms of the axle, adopting a distorted tetrahedral geometry. Comparing the average distances of the Cu+–olefin-C and Cu+–ether-O pairs in the lowest-energy structure with those in the corresponding X-ray crystal structure (summarized in Table S4), one can find that they are


14

Page 15 of 26 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


almost identical,44, 59 thus, demonstrating cogently the effectiveness and accuracy of our strategy. Besides, by comparing these representative structures, we can see that the macrocycle continues to rotate by nearly 150ºabout the axle from site 7 to site 8, which constitutes the second step of “

” of the rotaxane (shown in Figure 4C). In order to investigate the effect of the optimized Cu+–olefin-C and Cu+–ether-O parameters

on the global structure of the complex containing Cu+ ions, a representative spatial arrangement of the rotaxane corresponding to the global minimum (site 8, see Figure S9C) is compared with the crystallographic structure. Our results bring to light some discrepancies, which we ascribe to the effect of the solution. To mimic the environment of the crystal structure, we extend our investigation to the gas phase. Figure S11 depicts the free-energy landscape for the rotation of the ring determined in vacuum and the corresponding most stable structure superimposed with the crystal one. As can be seen, both the relative position of the axle and that the macrocycle are nearly identical to those of the X-ray structure, demonstrating that the NBFIX parameters optimized with our strategy can describe the structure of the molecular machine precisely. Our free-energy calculations, moreover, rationalize the mechanism underlying the relative motion in the nano-architecture observed experimentally. Detailed information is provided in the Supporting Information. 4. CONCLUSION In summary, a new strategy based on specific atom-pair adjustments to the van der Waals parameters was put forth. This strategy has proven effective to describe implicitly induction phenomena and mimic part of charge transfer and covalent bonding for metal-involved particle pairs. A cation-responsive molecular machine was further investigated using our new parameters.


15


Page 16 of 26

Due to a different affinity to olefin and ether, this mechanically interlocked supramolecular object adopts distinct conformations by rotating the macrocycle about the axle upon addition of Li+ and Cu+, inducing a “dialing” movement. Our theoretical investigation demonstrates that in stark contrast with a standard 12-6 model, the corrected van der Waals parameters for the specific atom pairs using the NBFIX feature by fitting the QM/MM PMF can describe the dialing process accurately. In principle, the approach proposed herein, using NBFIX corrections with PMF-calibrated parameters, can be applied to virtually any molecular assay involving metal-ion-induced polarization to describe with suitable accuracy the binding affinity of a metal ion toward specific functional groups. It is worth noting that the NFBIX corrections are only applied to specific atom pairs to address polarization effects, the rest of the standard parameter sets remaining untouched. Thus, this correction strategy will not alter the other interactions involved in the molecular assay, such as the metal-ion-water interaction. In addition, MD simulations using NBFIX option do not introduce any additional computational effort. From the test cases investigated in the present study, it can be seen that small alterations in the force field result in significant improvement. Put together, the theoretical framework presented here provides a cost-effective, general approach to model polarization effects in the context of an additive force field.

ASSOCIATED CONTENT Supporting Information. The Supporting Information is available free of charge on the ACS Publications website at DOI: xxxxx/xxxxxx


16

Page 17 of 26 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


Methods, simulation details, LJ 12-6 and 12-6-4 parameters used, charge transfer between Cu+– ether, Cu+–ethene, Li+–ether and Li+–ethene, covalent bonding between Cu+ and alkene, illustration of the QM and MM part in the QM/MM simulations, accuracy of PM7 with respect to the DFT method, time evolution of the root-mean-square deviation over the gradients of the free-energy profiles for the five independent windows of the simulation characterizing the association of Cu+–ethene, free-energy profiles of the Cu+–ethene interaction at different times using QM/MM potential, free-energy landscapes and the representative structures in the absence of metal ion and in the presence of Cu+ using the standard parameter sets, comparing of the representative structures using the standard 12-6 parameters with the X-ray structures in the absence of metal ion and in the presence of Li+, representative structures of the complex in site 6-site 8, comparing of the average distance of Cu+–olefin and Cu+–ether in different environments, free-energy landscape of the Cu+–containing complex in vacuum using the new strategy and the representative structure (PDF). All the files needed to perform the QM/MM simulations (ZIP).

Notes The authors declare no competing financial interests. ACKNOWLEDGMENT This study was supported by the National Natural Science Foundation of China (21773125), the Natural Science Foundation of Tianjin, China (18JCYBJC20500), and the Special Program for Applied Research on Super Computation of the NSFC-Guangdong Joint Fund (Second Phase) under Grant U1501501. C. C. is indebted to the Centre National de la Recherche Scientifique for u

j

(P

)w

P

’

u

.


17


Page 18 of 26

REFERENCES (1) Ma, P.; Cardenas, A. E.; Chaudhari, M. I.; Elber, R.; Rempe, S. B. The Impact of Protonation on Early Translocation of Anthrax Lethal Factor: Kinetics from Molecular Dynamics Simulations and Milestoning Theory. J. Am. Chem. Soc. 2017, 139, 14837-14840. (2) Stelzl, L. S.; Erlenbach, N.; Heinz, M.; Prisner, T. F.; Hummer, G. Resolving the Conformational Dynamics of DNA with Ångstrom Resolution by Pulsed Electron–Electron Double Resonance and Molecular Dynamics. J. Am. Chem. Soc. 2017, 139, 11674-11677. (3) Ene, C. D.; Maxim, C.; Rouzières, M.; Clérac, R.; Avarvari, N.; Andruh, M. Enantiopure versus Racemic Mixture in Reversible, Two-Step, Single-Crystal-to-Single-Crystal Transformations of Copper(II) Complexes. Chem. Eur. J. 2007, 13, 8164-8173. (4) Cooper, H. J.; Case, M. A.; McLendon, G. L.; Marshall, A. G. Electrospray Ionization Fourier Transform Ion Cyclotron Resonance Mass Spectrometric Analysis of Metal-Ion Selected Dynamic Protein Libraries. J. Am. Chem. Soc. 2003, 125, 5331-5339. (5) Pappalardo, S.; Villari, V.; Slovak, S.; Cohen, Y.; Gattuso,G.; Notti, A.; Pappalardo, A.; Pisagatti, I.; Parisi, M. F. Counterion-Dependent Proton-Driven Self-Assembly of Linear Supramolecular Oligomers Based on Amino-alix[5]arene Building Blocks. Chem. Eur. J. 2007, 13, 8164-8173. (6) Peters, M. B.; Yang, Y.; Wang, B.; Füsti-Molnár, L.; Weaver, M. N.; Merz, K. M. Structural Survey of Zinc-Containing Proteins and Development of the Zinc AMBER Force Field (ZAFF). J. Chem. Theory Comput. 2010, 6, 2935-2947.


18

Page 19 of 26 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


(7) Wang, Q.; Rackers, J. A.; He, C.; Qi, R.; Narth, C.; Lagardere, L.; Gresh, N.; Ponder, J. W.; Piquemal, J.-P.; Ren, P. General Model for Treating Short-Range Electrostatic Penetration in a Molecular Mechanics Force Field. J. Chem. Theory Comput. 2015, 11, 2609-2618. (8) Sakharov, D. V.; Lim, C. Zn Protein Simulations Including Charge Transfer and Local Polarization Effects. J. Am. Chem. Soc. 2005, 127, 4921-4929. (9) Dudev, T.; Lim, C. Competition among Metal Ions for Protein Binding Sites: Determinants of Metal Ion Selectivity in Proteins. Chem. Rev. 2014, 114, 538-556. (10) MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiórkiewicz-Kuczera, J.; Yin, D.; Karplus, M. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 1998, 102, 3586-3616. (11) MacKerell, A. D.; Feig, M.; Brooks, C. L. Improved Treatment of the Protein Backbone in Empirical Force Fields. J. Am. Chem. Soc. 2004, 126, 698-699. (12) Vanommeslaeghe, K.; Raman, E. P.; MacKerell, A. D. Automation of the CHARMM General Force Field (CGenFF) II: Assignment of Bonded Parameters and Partial Atomic Charges. J. Chem. Inf. Model. 2012, 52, 3155-3168. (13) Weiner, S. J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.; Alagona, G.; Profeta, S.; Weiner, P. A New Force Field for Molecular Mechanical Simulation of Nucleic Acids and Proteins. J. Am. Chem. Soc. 1984, 106, 765-784. (14) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A Second Generation Force


19


Page 20 of 26

Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. J. Am. Chem. Soc. 1995, 117, 5179-5197. (15) Cerutti, D. S.; Swope, W. C.; Rice, J. E.; Case, D. A. ff14ipq: A Self-Consistent Force Field for Condensed-Phase Simulations of Proteins. J. Chem. Theory Comput. 2014, 10, 4515-4534. (16) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225-11236. (17) Rizzo, R. C.; Jorgensen, W. L. OPLS All-Atom Model for Amines: Resolution of the Amine Hydration Problem. J. Am. Chem. Soc. 1999, 121, 4827-4836. (18) Robertson, M. J.; Tirado-Rives, J.; Jorgensen, W. L. Improved Peptide and Protein Torsional Energetics with the OPLS-AA Force Field. J. Chem. Theory Comput. 2015, 11, 3499-3509. (19) Li, P.; Merz, K. M. Taking into Account the Ion-Induced Dipole Interaction in the Nonbonded Model of Ions. J. Chem. Theory Comput. 2014, 10, 289-297. (20) Li, P.; Song, L. F.; Merz, K. M. Parameterization of Highly Charged Metal Ions Using the 12-6-4 LJ-Type Nonbonded Model in Explicit Water. J. Phys. Chem. B 2015, 119, 883-895. (21) Minoux, H.; Chipot, C. Cation−π Interactions in Proteins: Can Simple Models Provide an Accurate Description? J. Am. Chem. Soc. 1999, 121, 10366-10372. (22) Li, P.; Merz, K. M. Metal Ion Modeling Using Classical Mechanics. Chem. Rev. 2017, 117, 1564-1686. (23) Yu, H.; Whitfield, T. W.; Harder, E.; Lamoureux, G.; Vorobyov, I.; Anisimov, V. M.; MacKerell, A. D.; Roux, B. Simulating Monovalent and Divalent Ions in Aqueous Solution Using a Drude Polarizable Force Field. J. Chem. Theory Comput. 2010, 6, 774-786.


20

Page 21 of 26 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


(24) Shi, Y.; Xia, Z.; Zhang, J.; Best, R.; Wu, C.; Ponder, J. W.; Ren, P. Polarizable Atomic Multipole-Based AMOEBA Force Field for Proteins. J. Chem. Theory Comput. 2013, 9, 4046-4063. (25) Wu, J. C.; Piquemal, J.-P.; Chaudret, R.; Reinhardt, P.; Ren, P. Polarizable Molecular Dynamics Simulation of Zn(II) in Water Using the AMOEBA Force Field. J. Chem. Theory Comput. 2010, 6, 2059-2070. (26) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; Mackerell JR, A. D. CHARMM General Force Field: A Force Field for Drug‐Like Molecules Compatible with the CHARMM All‐atom Additive Biological Force Fields. J. Comput. Chem. 2010, 31, 671-690. (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) Huang, L.; Roux, B. Automated Force Field Parameterization for Nonpolarizable and Polarizable Atomic Models Based on Ab Initio Target Data. J. Chem. Theory Comput. 2013, 9, 3543-3556. (29) Wu, J. C.; Chattree, G.; Ren, P. Automation of AMOEBA polarizable force field parameterization for small molecules. Theor. Chem. Acc. 2012, 131, 1-11. (30) Li, P.; Merz, K. M. MCPB.py: A Python Based Metal Center Parameter Builder. J. Chem. Inf. Model. 2016, 56, 599-604. (31) Jiang, Y.; Zhang, H.; Tan, T. Rational Design of Methodology-Independent Metal Parameters Using a Nonbonded Dummy Model. J. Chem. Theory Comput. 2016, 12, 3250-3260.


21


Page 22 of 26

(32) Gresh, N.; Cisneros, G. A.; Darden, T. A.; Piquemal, J.-P. Anisotropic, Polarizable Molecular

Mechanics

Studies

of

Inter-

and

Intramolecular

Interactions

and

Ligand−Macromolecule Complexes. A Bottom-Up Strategy. J. Chem. Theory Comput. 2007, 3, 1960-1986. (33) Devereux, M.; Gresh, N.; Piquemal, J.-P.; Meuwly, M. A Supervised Fitting Approach to Force Field Parametrization with Application to the SIBFA Polarizable Force Field. J. Comput. Chem. 2014, 35, 1577-1591. (34) Venable, R. M.; Luo, Y.; Gawrisch, K.; Roux, B.; Pastor, R. W. Simulations of Anionic Lipid Membranes: Development of Interaction-Specific Ion Parameters and Validation Using NMR Data. J. Phys. Chem. B 2013, 117, 10183-10192. (35) Yoo, J.; Aksimentiev, A. New Tricks for Old Dogs: Improving the Accuracy of Biomolecular Force Fields by Pair-specific Corrections to Non-bonded Interactions. Phys. Chem. Chem. Phys. 2018, 20, 8432-8449. (36) Huang, J.; Rauscher, S.; Nawrocki, G.; Ran, T.; Feig, M.; de Groot, B. L.; Grubmüller, H.; MacKerell Jr, A. D. CHARMM36m: an Improved Force Field for Folded and Intrinsically Disordered Proteins. Nat. Methods 2017, 14, 71-73. (37) Miller, M. S.; Lay, W. K.; Li, S.; Hacker, W. C.; An, J.; Ren, J.; Elcock, A. H. Reparametrization of Protein Force Field Nonbonded Interactions Guided by Osmotic Coefficient Measurements from Molecular Dynamics Simulations. J. Chem. Theory Comput. 2017, 13, 1812-1826. (38) Lay, W. K.; Miller, M. S.; Elcock, A. H. Reparameterization of Solute—Solute Interactions for Amino Acid–Sugar Systems Using Isopiestic Osmotic Pressure Molecular Dynamics Simulations. J. Chem. Theory Comput. 2017, 13, 1874-1882.


22

Page 23 of 26 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


(39) Bernèche, S.; Roux, B. Energetics of ion conduction through the K+ channel. Nature 2001, 414, 73-77. (40) Luo, Y.; Roux, B. Simulation of Osmotic Pressure in Concentrated Aqueous Salt Solutions. J. Phys. Chem. Lett. 2010, 1, 183-189. (41) Yoo, J.; Aksimentiev, A. Improved Parametrization of Li+, Na+, K+, and Mg2+ Ions for All-Atom Molecular Dynamics Simulations of Nucleic Acid Systems. J. Phys. Chem. Lett. 2012, 3, 45-50. (42) Yoo, J.; Aksimentiev, A. Improved Parameterization of Amine–Carboxylate and Amine– Phosphate Interactions for Molecular Dynamics Simulations Using the CHARMM and AMBER Force Fields. J. Chem. Theory Comput. 2016, 12, 430-443. (43) Yoo, J.; Aksimentiev, A. Refined Parameterization of Nonbonded Interactions Improves Conformational Sampling and Kinetics of Protein Folding Simulations. J. Phys. Chem. Lett. 2016, 7, 3812-3818. (44) Baggi, G.; Loeb, S. J. Rotationally Active Ligands: Dialing-Up the Co-conformations of a [2]Rotaxane for Metal Ion Binding. Angew. Chem. Int. Ed. 2016, 55, 12533–12537. (45) Li, P.; Song, L. F.; Merz, K. M. Systematic Parameterization of Monovalent Ions Employing the Nonbonded Model. J. Chem. Theory Comput. 2015, 11, 1645-1657. (46) Ȧq

J.

W

P

D

F

E

y P

u

Simulations. J. Phys. Chem. 1990, 94, 8021-8024. (47) Beglov, D.; Roux, B. Finite Representation of an Infinite Bulk System: Solvent Boundary Potential for Computer Simulations. J. Chem. Phys. 1994, 100, 9050-9063.


23


Page 24 of 26

(48) Joung, I. S.; Cheatham, T. E. Determination of Alkali and Halide Monovalent Ion Parameters for Use in Explicitly Solvated Biomolecular Simulations. J. Phys. Chem. B 2008, 112, 9020-9041. (49) Won, Y. Force Field for Monovalent, Divalent, and Trivalent Cations Developed under the Solvent Boundary Potential. J. Phys. Chem. A 2012, 116, 11763-11767. (50) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926-935. (51) Melo, M. C. R.; Bernardi, R. C.; Rudack, T.; Scheurer, M.; Riplinger, C.; Phillips, J. C.; Maia, J. D. C.; Rocha, G. B.; Ribeiro, J. V.; Stone, J. E.; Neese, F.; Schulten, K.; Luthey-Schulten, Z. NAMD goes quantum: an integrative suite for hybrid simulations. Nat. Methods 2018, 15, 351. (52) Lin, H.; Truhlar, D. G. QM/MM: What Have We Learned, Where are We, and Where do We Go From Here? Theoretical Chemistry Accounts 2006, 117, 185. (53) Stewart, J. J. P. Optimization of Parameters for Semiempirical Methods VI: more Modifications to the NDDO Approximations and Re-optimization of Parameters. J. Mol. Model. 2013, 19, 1-32. (54) Stewart, J. J. P. MOPAC: A Semiempirical Molecular Orbital Program. J. Comput. Aided Mol. Des. 1990, 4, 1-103. (55) Phillips, J. C.; Braun, B.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kalé, L.; Schulten, K. Scalable Molecular Dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781-1802.


24

Page 25 of 26 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


(56) Fu, H.; Zhang, H.; Chen, H.; Shao, X.; Chipot, C.; Cai, W. Zooming across the Free-Energy Landscape: Shaving Barriers, and Flooding Valleys. J. Phys. Chem. Lett. 2018, 9, 4738-4745. (57) Fiorin, G.; Klein, M. L.; Hénin, J. Using Collective Variables to Drive Molecular Dynamics Simulations. Mol. Phys. 2013, 111, 3345-3362. (58) Ensing, B.; Laio, A.; Parrinello M.; Klein, M. L. A Recipe for the Computation of the Free Energy Barrier and the Lowest Free Energy Path of Concerted Reactions. J. Phys. Chem. B 2005, 109, 6676-6687. (59) Thompson, J. S.; Harlow, R. L.; Whitney, J. F. Copper(I)-olefin Complexes. Support for the Proposed Role of Copper in the Ethylene Effect in Plants. J. Am. Chem. Soc. 1983, 105, 3522-3527.


25


Page 26 of 26

for Table of Contents use only

Addressing Polarization Phenomena in Molecular Machines Containing Transition Metal Ions with an Additive Force Field

Shuangli Du, Haohao Fu, Xueguang Shao, Christophe Chipot, and Wensheng Cai


26

Addressing Polarization Phenomena in Molecular Machines

Recommend Documents