Figure
1.
Pathway of the oxidation of sulfur by H2O2 and the forming of disulfide bond.
| Citation: |
Shengwei Yuan, Guohui Li. Parameterization of Reactive Force Field for Thiol Oxidation and Disulfide Linking[J]. Chinese Journal of Chemical Physics . DOI: 10.1063/1674-0068/cjcp2408115
|
This study employed a simulated annealing algorithm to develop new reactive force field (ReaxFF) parameters specifically tailored for thiol oxidation and disulfide cross-linking reactions. The optimization strategy demonstrably enhanced the capability of ReaxFF to model these chemical processes. To mitigate potential overfitting during parameterization, we strategically expanded the training set’s conformational space relevant to the reactions through normal mode analysis. This approach effectively addressed concerns about overfitting, and the resulting parameters exhibit a degree of scalability for similar reactions. The optimized parameters achieve good agreement with experimental reaction energies. While surpassing the original parameters in accuracy, the characterization of reaction energy barriers is not ideal due to the exclusion of parameters related to O–H bond during optimization and inherent limitations within the ReaxFF framework. However, the current level of accuracy is sufficient for preliminary semi-quantitative studies. Notably, molecular dynamics simulations utilizing the new parameters successfully captured the binding mode of reactants in solvent-mediated proton transfer reactions, even for systems beyond the training set. These advancements in the ReaxFF force field offer a more efficient and reliable computational tool for exploring thiol oxidation and disulfide cross-linking reactions relevant to biological systems.
Thiol functional groups, present in essential biomolecules like cysteine residues and glutathione, are susceptible to oxidation. This process is intricately linked to enzymatic activity across various organisms. For instance, it is a crucial step in the catalytic cycle of the cytoplasmic enzyme ribonucleotide reductase [1]. Furthermore, the oxidation state of thiol groups directly impacts protein function, including photosynthetic enzymes in plants [2] and transcription factors in certain bacteria [3, 4]. Notably, thiol oxidation plays a vital role in cell signaling [5, 6]. Here, hydrogen peroxide (H2O2), a reactive oxygen species, reacts with thiols, constituting a key step in this pathway [7]. The process is postulated to involve some elementary reactions [8] (FIG. 1). Firstly, H2O2 facilitates the oxidation of sulfhydryl (RSH) to sulfoxide (RSHO). Subsequently, H2O2 directly oxidizes RSH to sulfenic acid (RSOH), then the sulfenic acid can be oxidized to sulfinic acid (RSO2H) and sulfonic acid (RSO3H) in sequence. Meanwhile, RSHO can transform to RSOH. Finally, a molecule of RSH reacts with a molecule of RSOH to form a disulfide bond (RSSR) and eliminate a water molecule. Disulfide bonds play a critical role in protein structure, significantly impacting the stability of proteins like human monoclonal antibodies [9, 10]. They also serve as a key factor influencing protein folding [11, 12] within living organisms.
Considering the importance of these reactions in biological systems, extensive investigations [6, 13−18] have been conducted by numerous researchers to elucidate the mechanism of disulfide cross-linking process between a pair of cysteine residues and H2O2. Quantum mechanics (QM)-based methods, such as density functional theory (DFT) [19, 20], offer powerful tools for describing chemical reactions at the atomic level. Extensive computational chemistry studies have been conducted using these methods to investigate the reactions mentioned above. However, such calculations are typically limited to static properties due to their high computational cost, restricting the accessible time and spatial scales. In contrast, empirical force field methods, like classical molecular dynamics (MD) simulations, are computationally more efficient and enable the study of dynamic system evolution over nanoseconds to microseconds with spatial resolution in the nanometer range. However, because the atom connectivity has been predefined, MD-based empirical force fields cannot describe chemical reactions correctly.
ReaxFF is a kind of reactive force field which can efficiently simulate the chemical reactions [21−23]. It has been successfully employed in numerous chemical reactions in aqueous solution over the past few decades, such as protonation of amino acids [24], proline-catalyzed aldol reaction [25, 26] and carbon dioxide capture reaction by phosphoenolpyruvate [27]. While ReaxFF has been successfully applied to a wide range of chemical systems, limited literature exists regarding its application to reactions involving functional groups containing sulfur atoms. The closest work reported by Müller et al. [28] focuses on disulfide mechanochemistry. However, the specific parameterization of sulfur atoms within their model may not be directly transferable to the current system and the target chemistries under investigation. Consequently, this existing force field might not accurately capture the formation and breakage of sulfur-containing chemical bonds in our study.
In this work, we introduce new ReaxFF parameters based on the original ReaxFF parameters (abbreviated as ff_origin), which are optimized for amino acids. First, we use the intrinsic reaction coordinate (IRC) calculations for relative reaction, then obtain molecular conformations of different reaction extent and the corresponding energetic data are added to the training set. Second, we employ a novel method to generate conformations based on the normal modes of molecules and subsequently expand the training set. The trained ReaxFF parameters have been successfully applied to the simulations of the disulfide crosslinking process. Due to the limitations of ReaxFF, accurately reproducing the reaction energy barriers remains challenging; however, it effectively captures the structural characteristics of the molecules during the reaction process. The optimized parameters also exhibit significant improvements in accurately describing responses outside the training set, which indicate that the parameters have a certain scalability. Consequently, the enhanced ReaxFF for disulfide bond condensation not only expands the capabilities of ReaxFF but also makes it possible to explore reactions involving sulfur-containing biomolecules such as cysteine.
In contrast to traditional classical force fields, ReaxFF offers the capability to accurately describe bond breaking and formation processes. The total energy of the system comprises two components: non-bonding energy terms and the terms based solely on bond order. In ReaxFF, atoms belonging to the same chemical element are assigned a single atom type. Moreover, rather than predefining bonding connectivity between calculated atoms, it is dynamically recalculated during the simulation using the bond order function.
The bond order function between specified atoms comprises three exponential terms, corresponding to single, double, and triple bonds respectively. Consequently, with an increasing interatomic distance, their bond order gradually diminishes from three to zero, signifying the dissociation of the bond between these atoms. Regarding the non-bond energy terms, the calculations of van der Waals and Coulomb forces are consistent with traditional classical force fields, encompassing all possible atomic pairs. The introduction of a shielding term prevents non-bonding interactions between closely positioned atoms. Lastly, in ReaxFF, atomic charges are regularly updated using the QEq charge equilibration method or an extended method during the simulation to reflect polarization effects.
The ff_origin parameter was developed by Susanna Monti and co-workers [29] based on the 2012 ReaxFF parameter for glycine [24]. Although this parameter is not the most up-to-date one optimized for biological system reactions, we still choose it as the initial parameters in our optimization process due to concerns raised in literature by Moerman et al. [30] regarding the newest parameter misrepresentation of chemical reaction pathways.
Considering the varying contributions of different parameters to our target system, we focuse our analysis on specific parameters directly associated with the reaction, such as C–S bond, S–S bond, S–O bond, and S–H bond parameters, along with their corresponding angles and dihedral angles. These parameters were re-optimized. In addition, though they related to the reaction, parameters describing amino acids or water molecules (e.g., O–H bonds) kept unchanged.
In this article, we choose methyl mercaptan as the model molecule. As the focus of our study does not encompass the conformational distribution of cysteine side chains, it is beyond the scope of our investigation. All quantum chemical calculations were performed using Gaussian 16 [31]. The M062X functional [32] and the basis set def2TZVP [33] combined D3 dispersion correction [34] were used in DFT calculations to simulate the geometrically optimized structures and get the energetic data. All molecular structures were optimized and subsequently subjected to vibrational frequency analysis, ensuring the absence of any imaginary frequencies in all reactants and products, as well as only one imaginary frequency in all transition states (TSs).
The training set includes rigid scans of related parameters for bond lengths, bond angles, and dihedral angles as they are essential for optimizing the parameters. Since our objective is to enhance the performance of the ReaxFF in accurately describing reactions such as sulfhydryl oxidation and associated disulfide bond condensation, the aforementioned geometric conformations are insufficient. Extensive literature has consistently demonstrated that without incorporating the relevant conformation of the reaction pathways into the training set, optimized parameters fail to adequately capture the correct reaction pathway. Consequently, we calculated IRC curves of four specific reactions in FIG. 1.
During the training process of ReaxFF parameters associated with proline-catalyzed aldol condensation, Hubin et al. [26] discovered that incorporating the conformations along IRC curves of a specific reaction into the training set significantly enhances the performance of said reaction. However, this may potentially decrease parameter scalability, necessitating the training of many sets of parameters to describe different reaction stages. Therefore, it is imperative to employ a strategic approach during parameter training by incorporating additional conformations closely associated with the reaction in the training set.
In this study, we used the IRC curves depicted in FIG. 1 as a central element to establish a systematic workflow for expanding the conformational space of the training set. A more comprehensive description is provided below.
(1) Select 10 to 15 points on each IRC curve in the training set. The distribution of these points is non-uniform, with a denser selection near the TS of the reaction, and at both ends encompassing the reactants and products.
(2) Calculate the frequencies for each selected point. After excluding the imaginary frequencies corresponding to the direction of the IRC curve, record the frequencies and direction vectors of all normal vibration modes.
(3) Since there is the abundance of regular vibration modes, and some higher frequency vibration modes, such as those involving methyl carbon–hydrogen bond, are not of interest to us, we exclude vibration modes with wavenumbers exceeding 2000 cm–1. However, visual inspection is still required to prevent missing conformations that are closely related to the parameters to be optimized.
(4) Then generate new conformations from vibrational modes along the respective normal mode direction vectors.
(5) Subsequently, eliminate high-energy structures resulting from potential atomic overlap, certain conformations are discarded after quantum chemical calculations.
(6) Finally, a selection of conformations from the IRC curves of reactions 1–4 in FIG. 1 was added to the training set.
Many studies have introduced the most stable cluster structures in the training set for reactions involving multiple reactant molecules. However, these parameters have limited ability for the reaction. Instead, we incorporated a ‘molecular pair’ structure with restricted interatomic distances in the training set, and our subsequent results demonstrate that this addition significantly improves the ability to accurately characterize the reaction TSs.
In order to determine the optimal parameters, we minimize the following objective functions.
| Error=n∑i=1[(ref_rffi−ref_qmi)σi]2 |
Here, ref_qm represents the energy difference between the conformation in the training set and the reference conformation calculated using quantum chemistry. Similarly, ref_rff denotes the energy difference calculated using ReaxFF calculation. Additionally, σi signifies the weight assigned to different conformations.
Weight selection is employed to ensure approximate equality in objective function values calculated for each conformation type. Additionally, for conformations generated through rigid scanning, their weights are reduced if ref_qm exceeds 100 kcal/mol. Moreover, corresponding conformations near the TS on the IRC curve are assigned artificially increased weights. Lastly, for conformations obtained via normal vibration modes, decreasing weights are assigned as they deviate further from the conformation nearby IRC curve.
Due to the large number of parameters needed in the optimization process, we have chosen a combination of simulated annealing [35] (SA) and one-parameter search technique [36] to minimize the use of artificially selected parameters during parameter optimization, both of them are widely used in optimization of ReaxFF parameters [22, 27, 36, 37]. To improve our efficiency of the optimization process, we gradually expand the training set and incrementally adjust the upper and lower limits of the parameters under consideration.
Although certain combinations may exhibit superior performance on the objective function, they fail to support stable MD simulations. Hence, after the error function converges, a short MD simulation of 2 ps with a timestep of 0.25 fs and temperature set at 300 K is performed here with the selected initial conformation. Subsequently, the trajectory is visually inspected to eliminate parameter combinations that produce unreasonable structure. The process is iterated until the overall performance of the new parameters achieved acceptable errors, ensuring that the key conformation matched the training data. All MD simulations in this article are performed using LAMMPS [38] and PLUMED [39]. The trajectories are analyzed by VMD [40] and MDTraj [41]. The specific workflow is shown in FIG. 2. The scripts we utilized are written in Python and Bash.
To evaluate the performance of optimized parameters, we conduct extensive MD simulations. Except for the reactions in FIG. 1, the processes with different number of solvents assisted with varying numbers of explicit solvent molecules are also included.
The reactant molecules are confined within a cubic water box with sides measuring 20 Å. The energy minimization of each system was conducted at T = 0 K. Subsequently, the temperature of the systems was gradually increased from 0 K to 298.15 K over 25 ps. Once reaching 298.15 K, NPT MD simulations are performed for various simulation systems with total time ranging from 25 ps to 50 ps. Temperature was maintained using a Nosé-Hoover thermostat with a damping constant of 25 fs, and pressure was maintained using a Nosé-Hoover barostat with a damping parameter of 250 fs. To ensure the reliability of the results, each MD simulation was performed multiple times in parallel. In order to observe the reaction occurring within the simulated time, it is imperative to apply additional potential on atoms associated with the reaction. Simultaneously, to minimize artificial constraints, we choose linear combination of bond lengths in reaction as collective variables (CVs).
In the following section, we denote the newly trained force field ff_new, which will be compared to the original force field denoted as ff_origin. The performance of ff_new is comprehensively evaluated by conducting a detailed analysis of the energetics associated with various scans outlined in Section II.
The first step involved a rigid scan of the dissociation processes occurring in various chemical bonds. The results depicted in the aforementioned figure (FIG. 3) exhibit a good agreement. For chemical bonds that can be reasonably described by ff_origin, such as S–S bond, although the optimized parameters cannot fully capture the energy curve calculated by QM, they exhibit varying degrees of improvement. In the case of the S–S bond, ff_origin overestimates the energy of the molecule after dissociation, while the optimized parameters effectively addressed this issue.
The ff_origin fails to accurately describe the S–O bond, whereas our newly optimized parameter effectively addresses this issue. As depicted in FIG. 3, the energy calculated using the original parameter consistently decreases with increasing S–O bond length. This observation indicates that ff_origin is highly unreasonable, as it does not account for effective bonding between sulfur and oxygen atoms or an equilibrium bond length that aligns with reality, thus rendering it incapable of simulating relevant reactions. The novel parameter significantly enhances the characterization performance of the S–O bond in both two molecules, while accurately reproducing the equilibrium bond length and bond dissociation energy. We observe that ff_new exhibits a non-smooth component during the dissociation process of S–OH functional group; however, its impact on the simulation of chemical reactions is negligible. Similarly, in the dissociation process of S=O bond, the new parameter slightly overestimates the force constant near equilibrium. It is noteworthy that underestimation poses a greater risk in practical applications as it may result in non-physical atomic distances. Conversely, overestimation can generally be mitigated through energy minimization prior to simulation, thereby avoiding significant errors in system description. Nevertheless, the descriptions provided by ff_new regarding bond dissociation in equilibrium and at asymptotic stage are quite similar to the results of QM calculations, this circumstance does not lead to significant error consequences.
The optimized parameters provide a more accurate description of the key angle parameters associated with the target reaction compared to ff_origin. While ff_origin adequately captures the energetics of C–S–H and C–S–S angles, our optimization efforts focus on improving C–S–O and S–O–H angles related to the S–O bond. Although our training set includes a rigid scan of the angle in MeSHO, considering the coupling of multiple angles involved, it is more reasonable to utilize data from the rigid scan in MeSOH to represent optimized parametric performance. FIG. 4 demonstrates that the original parameters overestimate the slope of the energy curve. Given that changes in atom angles directly related to the reaction are significant, this undoubtedly diminishes our ability to simulate the reaction. These issues are resolved by employing optimized parameters which can accurately reproduce both bond angles and equilibrium positions.
The capability of the novel parameter to describe a chemical reaction is illustrated in FIG. 5. The energy profile computed by ReaxFF is obtained through fitting the structure along the IRC curve calculated by DFT. We list the reaction energy barrier, the reaction energy and the root mean squared error (RMSE) between ReaxFF and QM curves obtained from FIG. 5 in Table I.
| R | ΔE | ΔE‡ | RMSE | |||||||
| QM | ff_origin | ff_new | QM | ff_origin | ff_new | ff_origin | ff_new | |||
| 1 | −40.79 | 203.92 | −38.18 | 48.08 | 156.47 | 81.54 | 137.92 | 17.04 | ||
| 2 | −55.12 | 68.75 | −62.85 | 57.78 | 134.5 | 65.11 | 82.96 | 17.18 | ||
| 3 | −15.50 | −91.96 | −8.91 | 46.29 | 5.47 | 52.22 | 49.88 | 10.52 | ||
| 4 | −25.25 | −101.51 | −20.96 | 48.21 | 44.11 | 47.19 | 57.81 | 11.95 | ||
DownLoad:
CSV
According to Table I, it is evident that the ff_origin exhibits a mean unsigned error (MUE) of 130.32 kcal/mol in comparison with the QM calculation for the four reactions in the training set. On the other hand, after optimization, the MUE between the ff_new and QM calculation results decreases to 5.31 kcal/mol. This outcome indicates that optimization effectively addresses the inadequacy of the unoptimized parameters in describing reactants, products, and intermediates in chemical reactions. Concurrently, optimization of the force field parameters yields a marked improvement in the description of reaction energy barriers. The MUE for energy barriers decreased from 57.51 kJ/mol to 11.93 kJ/mol, signifying a substantial enhancement in the representation of key reaction TSs. However, the improvement for these high-energy species was less pronounced compared to the improvement observed for stable structures like reactants and products. This observation aligns with inherent limitations of ReaxFF, particularly regarding its tendency to overestimate energy barriers due to its reliance on bond order. Several studies have reported this overestimation issue with ReaxFF [30].
It is noteworthy that during simulations of reactions 3 and 4, ff_new, while accurately reproducing the overall reaction energy, generated several unphysical minima in the energy curve (FIG. 5(c, d)), similar to those observed with ff_origin. Notably, reactions 1, 2, and 4 involve breaking and formation of O–H bonds. However, only reaction 1 maintains the same number of O–H bonds before and after the reaction, potentially mitigating any errors introduced by the optimization process. Conversely, reactions 2 and 4 involve more complex bond-breaking and bond-forming scenarios, leading to these unsatisfactory results. This phenomenon is not exclusive to aqueous solution reaction systems [42]; it is also evident in ReaxFF-based simulations of combustion reactions [43]. Refinement of key parameters associated with O–H interactions may be necessary to address this issue. Furthermore, the inherent complexity of the ReaxFF function can introduce instabilities in the energy curve, despite the use of taper functions to mitigate abrupt changes. To address this challenge, some ongoing research [44, 45] focuses on refining the functional form of ReaxFF.
Overall, despite the presence of occasional false peaks, ff_new exhibits commendable performance in characterizing reaction energy and significantly enhances the energy barrier compared to ff_origin. This enables characterization of reactions that were previously unattainable using the original force field parameters.
We selected reactions 5 and 6 in FIG. 1 outside of the training sets for testing extensibility. The approximate reaction energy and energy barrier are determined using the same methodology as described above. It is evident that ff_origin fails to accurately describe these reactions. Although ff_new tends to overestimate the energy barrier, it is suitable for qualitative research purposes (FIG. 6). Notably, closely related parameters, such as O–S–O angle, have not been optimized yet. These tests demonstrate that the newly optimized parameters exhibit excellent performance in thiol oxidation and disulfide bond formation while also displaying potential for future expansion into reaction systems involving disulfide bond oxidation.
The tests in the previous section to evaluate the new parameters are conducted without MD simulations and in the absence of explicit water molecules. A recurrent challenge during force field parameter optimization is the potential to obtain parameters that accurately reproduce the IRC curves but yield unrealistic structures during subsequent MD simulations [30]. These structures often violate fundamental chemical principles. To assess their suitability for simulating biochemical reactions in aqueous solutions, it is imperative to examine their performance in MD simulations. Solvent-assisted proton transfer (SAPT) is a common phenomenon in aqueous systems. We simulated eight reactions about reactions 1–4 in FIG. 1. Each reaction involved one or two solvent molecules facilitating the proton transfer process. It is worth noting that there are no SAPT-related conformations in our training set, so further testing of the scalability of the parameters can be performed while evaluating their performance in MD simulations. To efficiently capture the structure of key transition states within a limited simulation time, we employed the technique of steered molecular dynamics [46] (SMD), with specific CVs listed in Supplementary materials (SM). All simulations successfully identify reaction events within 50 ps.
As a case study, we present a detailed analysis of reaction 1 in FIG. 1 assisted by a single water molecule. FIG. 7(a) depicts the initial configuration of methyl mercaptan, hydrogen peroxide, and the intervening water molecule involved in the reaction. The key interatomic distances (d1, d2, d3) serve as the CV monitored during the simulation, and their evolution over time is presented in FIG. 7(b). Choosing the sum of d1, d2 and d3 as the collective variable (CV) effectively captured the dynamic changes in these distances during the reaction. Notably, the simulation did not exhibit any unrealistic scenarios where interatomic distances became abnormally close or excessively large simultaneously. Furthermore, the observed changes in d1, d2 and d3 displayed consistent trends throughout the simulation. In contrast to the expected smooth evolution of bond lengths observed in typical QM/MM simulations, our results reveal a significant lack of coordination among the changes in d1, d2, and d3. Notably, during the later stages of the simulation, d1 exhibits a substantial decrease even when d2 and d3 remain nearly constant. Unfortunately, further refinement of the force field parameters did not alleviate this issue. This observation aligns with findings in the literature, suggesting inherent limitations within the ReaxFF regarding its ability to capture the concerted nature of bond-breaking and bond-formation during complex reactions [30].
The simulation result of the reaction in FIG. 7(a) successfully reproduces a linear configuration for the S1–O6–O8 angle. The average angle value throughout the simulation remained at 167°, faithfully reproducing the relative spatial arrangement of the reacting molecules. Encouragingly, the optimized parameters consistently capture the binding mode during proton transfer involving one or two water molecules. Snapshots of different stages of the reactions in FIG. 1 will also be listed in FIG. S1 to FIG. S8 in SM, which can be compared with the results derived from QM/MM simulations. Notably, they accurately reproduce the planar, four-membered or five-membered ring structures formed by the participating sulfur and oxygen atoms [13].
Similar to the situation shown in FIG. 7(b), it is worth mentioning that all these transition states exhibit uncoordinated proton transfer, indicating a tendency towards sequential proton transfer reactions—a departure from QM/MM simulation results. However, this simulation approach significantly enhances computational efficiency and provides valuable insights for preliminary investigations.
Our findings suggest that a training set solely comprised of IRC-related conformations can be detrimental for capturing reactions in extended MD simulations. While this approach may yield satisfactory results for static properties and short MD simulations during parameter optimization, it can lead to inaccurate reaction simulations in longer time scales. Limiting the training set to conformations within the range of 155°–180° for the S5–O7–O8 angle restricts the data available for parameter optimization and can lead to overfitting. This manifests as the appearance of unphysical minima (pseudo-minima) in non-chemical regions of the potential energy surface. These artifacts can introduce errors in reaction simulations, such as the error cleavage of carbon-sulfur bonds.
FIG. 8 exemplifies this concept for reaction 1 in FIG. 1 without water molecules. The average S5–O7–O8 angle in the early simulation, before and after the introduction of extended conformations, was 102° and 165° respectively. This discrepancy highlights the inability of the parameters trained without extended conformations to capture the correct binding mode of the reactants. Conversely, incorporating extended conformations through normal mode analysis allows for a more accurate representation, as evidenced by the close agreement with QM/MM calculations.
These observations emphasize the importance of employing a comprehensive training set. Expanding the training set conformations through normal modes is a very effective strategy. This approach effectively mitigates overfitting and improves the reliability of force field predictions.
We present an extension of the ReaxFF for its application to thiol oxidation and disulfide cross-linking reactions in biological systems. The original parameters were found to be inadequate for describing these reaction mechanisms. Therefore, we primarily optimized parameters associated with the S–O bond and its bond angle. To mitigate potential overfitting during optimization, we augmented the IRC-derived conformational space through normal mode analysis. The optimized force field exhibits good agreement between the relative energies of reactants, products, and stable intermediates in reactions involving thiols with hydrogen peroxide and disulfide bond formation. Additionally, the calculated reaction energy barriers provide a sufficient foundation for preliminary semi-quantitative studies. Furthermore, MD simulations with the new parameters accurately capture the binding mode of reacting molecules and reproduce solvent-mediated proton transfer reactions, even in systems beyond the training set. But systems with significantly different chemical environments or bonding motifs may require further targeted optimization. In summary, these new parameters facilitate efficient simulations of thiol oxidation and disulfide-related reactions.
Supplementary materials: More details about MD simulation and the final force field are available.
The authors declare no competing financial interests.
| [1] |
A. Jordan and P. Reichard, Annu. Rev. Biochem. 67, 71 (1998). doi: 10.1146/annurev.biochem.67.1.71
|
| [2] |
S. D. Dai, C. Schwendtmayer, P. Schürmann, S. Ramaswamy, and H. Eklund, Science 287, 655 (2000). doi: 10.1126/science.287.5453.655
|
| [3] |
M. Zheng, F. Åslund, and G. Storz, Science 279, 1718 (1998). doi: 10.1126/science.279.5357.1718
|
| [4] |
U. Jakob, W. Muse, M. Eser, and J. C. A. Bardwell, Cell 96, 341 (1999). doi: 10.1016/S0092-8674(00)80547-4
|
| [5] |
A. Bindoli, J. M. Fukuto, and H. J. Forman, Antioxid. Redox Signal. 10, 1549 (2008). doi: 10.1089/ars.2008.2063
|
| [6] |
C. A. Bayse, Org. Biomol. Chem. 9, 4748 (2011). doi: 10.1039/C1OB05497J
|
| [7] |
E. A. Veal, A. M. Day, and B. A. Morgan, Mol. Cell 26, 1 (2007). doi: 10.1016/j.molcel.2007.03.016
|
| [8] |
C. Jacob, A. L. Holme, and F. H. Fry, Org. Biomol. Chem. 2, 1953 (2004). doi: 10.1039/B406180B
|
| [9] |
T. Martinez, A. Guo, M. J. Allen, M. Han, D. Pace, J. Jones, R. Gillespie, R. R. Ketchem, Y. L. Zhang, and A. Balland, Biochemistry 47, 7496 (2008). doi: 10.1021/bi800576c
|
| [10] |
A. McAuley, J. Jacob, C. G. Kolvenbach, K. Westland, H. J. Lee, S. R. Brych, D. Rehder, G. R. Kleemann, D. N. Brems, and M. Matsumura, Protein Sci. 17, 95 (2008). doi: 10.1110/ps.073134408
|
| [11] |
B. S. Mamathambika and J. C. Bardwell, Annu. Rev. Cell Dev. Biol. 24, 211 (2008). doi: 10.1146/annurev.cellbio.24.110707.175333
|
| [12] |
P. Kosuri, J. Alegre-Cebollada, J. Feng, A. Kaplan, A. Inglés-Prieto, C. L. Badilla, B. R. Stockwell, J. M. Sanchez-Ruiz, A. Holmgren, and J. M. Fernández, Cell 151, 794 (2012). doi: 10.1016/j.cell.2012.09.036
|
| [13] |
M. A. Hagras, M. A. Bellucci, G. Gobbo, R. A. Marek, and B. L. Trout, J. Phys. Chem. B 124, 9840 (2020). doi: 10.1021/acs.jpcb.0c07510
|
| [14] |
M. A. Hagras, R. A. Marek, F. Hatahet, and B. L. Trout, J. Phys. Chem. B 126, 5972 (2022). doi: 10.1021/acs.jpcb.2c03588
|
| [15] |
L. A. H. van Bergen, G. Roos, and F. De Proft, J. Phys. Chem. A 118, 6078 (2014). doi: 10.1021/jp5018339
|
| [16] |
A. Zeida, R. Babbush, M. C. G. Lebrero, M. Trujillo, R. Radi, and D. A. Estrin, Chem. Res. Toxicol. 25, 741 (2012). doi: 10.1021/tx200540z
|
| [17] |
M. Khavani, M. Izadyar, and M. Reza Housaindokht, Phosphorus, Sulfur, Silicon Relat. Elem. 190, 1680 (2015). doi: 10.1080/10426507.2015.1019069
|
| [18] |
B. Cardey and M. Enescu, J. Phys. Chem. A 111, 673 (2007). doi: 10.1021/jp0658445
|
| [19] |
R. G. Parr and W. T. Yang, J. Am. Chem. Soc. 106, 4049 (1984). doi: 10.1021/ja00326a036
|
| [20] |
R. G. Parr and W. T. Yang, Annu. Rev. Phys. Chem. 46, 701 (1995). doi: 10.1146/annurev.pc.46.100195.003413
|
| [21] |
A. C. T. van Duin, S. Dasgupta, F. Lorant, and W. A. Goddard, J. Phys. Chem. A 105, 9396 (2001). doi: 10.1021/jp004368u
|
| [22] |
K. Chenoweth, A. C. T. van Duin, and W. A. Goddard, J. Phys. Chem. A 112, 1040 (2008). doi: 10.1021/jp709896w
|
| [23] |
T. P. Senftle, S. Hong, M. Islam, S. B. Kylasa, Y. X. Zheng, Y. K. Shin, C. Junkermeier, R. Engel-Herbert, M. J. Janik, H. M. Aktulga, T. Verstraelen, A. Grama, and A. C. T. van Duin, npj Comput. Mater. 2, 15011 (2016). doi: 10.1038/npjcompumats.2015.11
|
| [24] |
O. Rahaman, A. C. T. van Duin, W. A. Goddard III, and D. J. Doren, J. Phys. Chem. B 115, 249 (2011). doi: 10.1021/jp108642r
|
| [25] |
P. O. Hubin, D. Jacquemin, L. Leherte, J. M. André, A. C. T. van Duin, and D. P. Vercauteren, Theor. Chem. Acc. 131, 1261 (2012). doi: 10.1007/s00214-012-1261-4
|
| [26] |
P. O. Hubin, D. Jacquemin, L. Leherte, and D. P. Vercauteren, J. Comput. Chem. 37, 2564 (2016). doi: 10.1002/jcc.24481
|
| [27] |
Y. H. Huang, A. S. Wexler, K. J. Bein, and R. Faller, J. Phys. Chem. C 126, 9284 (2022). doi: 10.1021/acs.jpcc.2c01841
|
| [28] |
J. Müller and B. Hartke, J. Chem. Theory Comput. 12, 3913 (2016). doi: 10.1021/acs.jctc.6b00461
|
| [29] |
S. Monti, A. Corozzi, P. Fristrup, K. L. Joshi, Y. K. Shin, P. Oelschlaeger, A. C. T. van Duin, and V. Barone, Phys. Chem. Chem. Phys. 15, 15062 (2013). doi: 10.1039/c3cp51931g
|
| [30] |
E. Moerman, D. Furman, and D. J. Wales, J. Chem. Theory Comput. 17, 497 (2021). doi: 10.1021/acs.jctc.0c01043
|
| [31] |
M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E . Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Jr. Montgomery, J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman, and D. J. Fox, Gaussian 16, Revision B. 01, Wallingford: Gaussian, Inc., (2016).
|
| [32] |
Y. Zhao and D. G. Truhlar, Theor. Chem. Acc. 120, 215 (2008). doi: 10.1007/s00214-007-0310-x
|
| [33] |
F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys. 7, 3297 (2005). doi: 10.1039/b508541a
|
| [34] |
S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J. Chem. Phys. 132, 154104 (2010). doi: 10.1063/1.3382344
|
| [35] |
E. Iype, M. Hütter, A. P. J. Jansen, S. V. Nedea, and C. C. M. Rindt, J. Comput. Chem. 34, 1143 (2013). doi: 10.1002/jcc.23246
|
| [36] |
A. C. T. van Duin, J. M. A. Baas, and B. van de Graaf, J. Chem. Soc. Faraday Trans. 90, 2881 (1994). doi: 10.1039/FT9949002881
|
| [37] |
A. Rahnamoun, M. C. Kaymak, M. Manathunga, A. W. Götz, A. C. T. van Duin, K. M. Merz Jr., and H. M. Aktulga, J. Chem. Theory Comput. 16, 7645 (2020). doi: 10.1021/acs.jctc.0c00874
|
| [38] |
A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in 't Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, Comput. Phys. Commun. 271, 108171 (2022). doi: 10.1016/j.cpc.2021.108171
|
| [39] |
G. A. Tribello, M. Bonomi, D. Branduardi, C. Camilloni, and G. Bussi, Comput. Phys. Commun. 185, 604 (2014). doi: 10.1016/j.cpc.2013.09.018
|
| [40] |
W. Humphrey, A. Dalke, and K. Schulten, J. Mol. Graph. 14, 33 (1996). doi: 10.1016/0263-7855(96)00018-5
|
| [41] |
R. T. McGibbon, K. A. Beauchamp, M. P. Harrigan, C. Klein, J. M. Swails, C. X. Hernández, C. R. Schwantes, L. P. Wang, T. J. Lane, and V. S. Pande, Biophys. J. 109, 1528 (2015). doi: 10.1016/j.bpj.2015.08.015
|
| [42] |
T. Trnka, I. Tvaroška, and J. Koča, J. Chem. Theory Comput. 14, 291 (2018). doi: 10.1021/acs.jctc.7b00870
|
| [43] |
L. W. Bertels, L. B. Newcomb, M. Alaghemandi, J. R. Green, and M. Head-Gordon, J. Phys. Chem. A 124, 5631 (2020). doi: 10.1021/acs.jpca.0c02734
|
| [44] |
D. Furman and D. J. Wales, J. Phys. Chem. Lett. 10, 7215 (2019). doi: 10.1021/acs.jpclett.9b02810
|
| [45] |
D. Furman and D. J. Wales, J. Chem. Phys. 153, 021102 (2020). doi: 10.1063/5.0013906
|
| [46] |
S. Izrailev, S. Stepaniants, B. Isralewitz, D. Kosztin, H. Lu, F. Molnar, W. Wriggers, and K. Schulten, Proceeding of the 2nd International Symposium on Algorithms for Macromolecular Modelling Computational Molecular Dynamics: Challenges, Methods, Ideas, Berlin Heidelberg, (1999).
|
Figures(8) / Tables(1)
CJCP WeChat
Copyright© 2016 CPS
Address: No.96, JinZhai Road, Hefei, Anhui, 230026, P.R.China
Tel: 0551-63601122 E-mail: cjcp@ustc.edu.cn
Supported by:
Beijing Renhe Information Technology Co., Ltd.
| R | ΔE | ΔE‡ | RMSE | |||||||
| QM | ff_origin | ff_new | QM | ff_origin | ff_new | ff_origin | ff_new | |||
| 1 | −40.79 | 203.92 | −38.18 | 48.08 | 156.47 | 81.54 | 137.92 | 17.04 | ||
| 2 | −55.12 | 68.75 | −62.85 | 57.78 | 134.5 | 65.11 | 82.96 | 17.18 | ||
| 3 | −15.50 | −91.96 | −8.91 | 46.29 | 5.47 | 52.22 | 49.88 | 10.52 | ||
| 4 | −25.25 | −101.51 | −20.96 | 48.21 | 44.11 | 47.19 | 57.81 | 11.95 | ||
DownLoad:
CSV
