Structural Articulation of Biochemical Reactions ... - ACS Publications

Article Cite This: J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

pubs.acs.org/jcim

Structural Articulation of Biochemical Reactions Using Restrained Geometries and Topology Switching Swati R. Manjari† and Nilesh K. Banavali*,†,‡ †

Laboratory of Computational and Structural Biology, Division of Genetics, NYS Department of Health, CMS 2008, Biggs Laboratory, Wadsworth Center, Empire State Plaza, Albany, New York 12201-0509, United States ‡ Department of Biomedical Sciences, School of Public Health, State University of New York at Albany, Albany, New York 12222, United States S Supporting Information *

ABSTRACT: A strategy named “restrained geometries and topology switching” (RGATS) is presented to obtain detailed trajectories for complex biochemical reactions using molecular mechanics (MM) methods. It enables prediction of realistic dynamical pathways for chemical reactions, especially for accurately characterizing the structural adjustments of highly complex environments to any proximal biochemical reaction. It can be used to generate reactive conformations, model stepwise or concerted reactions in complex environments, and probe the influence of changes in the environment. Its ability to take reactively nonoptimal conformations and generate favorable starting conformations for a biochemical reaction is illustrated for a proton transfer between two model compounds. Its ability to study concerted reactions in explicit solvent is illustrated using proton transfers between an ammonium ion and two conserved histidines in an ammonia transporter channel embedded in a lipid membrane. Its ability to characterize the changes induced by subtle differences in the active site environment is illustrated using nucleotide addition by a DNA polymerase in the presence of two versus three Mg2+ ions. RGATS can be employed within any MM program and requires no additional software implementation. This allows the full assortment of computational methods implemented in all available MM programs to be used to tackle virtually any question about biochemical reactions that is answerable without using a quantum mechanical (QM) model. It can also be applied to generate reasonable starting structures for more detailed and expensive QM or QM/MM methods. In particular, this strategy enables rapid prediction of reactant, intermediary, or product state structures in any macromolecular context, with the only requirement being that the structure in any one of these states is either known or can be accurately modeled.

■

INTRODUCTION Dynamic trajectories of biochemical reactions catalyzed by large macromolecular complexes are difficult to characterize experimentally in atomic detail.1−7 Predicting them using quantum mechanics/molecular mechanics (QM/MM) approaches provides accurate insights but requires significant computational resources.8−12 If the electronic structure details of the change in covalent bonding are not significant to the structural or dynamic aspects of the system being considered, MM energy functions can provide a sufficiently accurate model. This is also an underlying principle of QM/MM methods, which treat the nonreactive surrounding environment using an MM model, while using a QM description for the reactive regions.13,14 The use of MM methods to study biochemical reactions also has a long history, with the empirical valence bond (EVB) method, which represents the reaction as a superposition of resonance forms, being the earliest15,16 and most widely used approach.17,18 Multiple other methods, such as SEAM,19 RFF,20 Q2MM,21 multiconfigurational MM (MCMM),22 REAXXF,23 and ACE24 have been implemented © XXXX American Chemical Society

to inexpensively model chemical reactions on large scales or in complex systems. Such approaches typically utilize anharmonic bond potentials compatible with bond cleavage, which are not standard in the macromolecular MM energy functions commonly used to study macromolecular dynamics.25−27 They also usually require specific parametrization of such bond potentials that take into account the energetics of the bond formation or cleavage reaction.28 All MM programs, even those that do not have EVB or another reactive MM method implemented, allow application of distance, angular, or torsion restraints that can mimic the presence of a covalent bond. Such simple geometric restraints can be imposed and manipulated at the user-level of the MM program. It is possible to apply them to guide the reactive region through structural changes that mimic the chemical path. While this is a highly controlled structural change for the reactive region with no representation of its internal energetics, Received: December 6, 2017 Published: January 22, 2018 A

DOI: 10.1021/acs.jcim.7b00699 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Information and Modeling

Figure 1. RGATS strategy and its application in predicting structural states accompanying biochemical reactions using MM simulations. (A) Scheme showing the RGATS strategy for a typical reaction involving cleavage of one covalent bond and formation of another. (B) Implementation of RGATS in any MM program only requires the five simple steps shown.

reactant-like state is one where r10 of bond 1 is equal to a covalent equilibrium bond length, and r20 of bond 2 is equal to the distance between the two separate atoms in the reactant state. Similarly, in a canonical product-like state, r10 is equal to the distance between the two separated atoms in the product state, and r20 is equal to a covalent equilibrium bond length. Appropriate manipulation of the r10 and r20 parameters, along with all relevant angular, chiral, or torsion parameters, can thus yield quasistatic transitions between the reactant and product states in either direction. MM energy functions have distinct topologies for the reactant and product states, so to interconvert between these states and the intermediary state, a user-level patch that switches topological states is required. All calculations were performed using the c35b3 version of the CHARMM program29 or version 2.9 of the NAMD program,30 the CHARMM22 force field for proteins,31 and CHARMM27 force field for nucleic acids32,33 or lipids.34 The TIP3P model35 was used for water molecules, and previous models were used for the cations and anions.36 Since RGATS is implemented through topology and parameter file modifications, and userlevel restraints, any other macromolecular simulation software or force field can be used. Intermediary and Product State Models. No new parameters are needed for RGATS simulations, but it may be preferable to model a partial charge rearrangement in the intermediary state to make a smooth transition between reactant and product states. For this, either previous parameters for analogous entities could simply be transferred, or a partial charge distribution between the reactant and product states could be used. If greater accuracy is desired, model compounds representative of the real transition states could be generated, and a standard parametrization protocol37,38 used to optimize the intermediary state partial charges. Another choice regards the introduction or removal of the orientational restraints (angles, dihedrals, chirality) after the topology switches. These are switched fully at the user-level in the present study, but could be phased in more gradually using approaches like bias annealing,39 if deemed necessary. Adjusting their introduction

the surrounding complex molecular environment is nevertheless expected to respond to the reactive change realistically. In this study, we present a strategy based on this expectation called restrained geometries and topology switching (RGATS) to obtain atomic detail trajectories for biochemical reactions. Since it relies on simple geometric restraints, it can be used within any MM program to model biochemical reactions, even by users with just a basic familiarity with the program. We illustrate the ability of the RGATS strategy to start from nonoptimal structures and transition to reactively optimal structures by applying it to the simple model system of the transfer of a proton from one methyl imidazole group to another. To show that RGATS can be applied in a complex biomolecular explicit solvent environment, a concerted proton transfer reaction involving two histidines in the AmtB ammonia transporter embedded in a lipid environment is also modeled. The effect of small but critical changes in the active site environments is probed by the presence of either two or three Mg2+ ions during nucleotide incorporation by the DNA polymerase μ. The applicability and drawbacks of this strategy in modeling the dynamics of biochemical reactions are also discussed.

■

MATERIALS AND METHODS Restrained Geometries And Topology Switching. The main idea behind the RGATS strategy is using standard geometric restraints as drivers of the structural pathway of the biochemical reaction to be modeled. The general principle is illustrated in Figure 1A. In a typical reaction where one bond breaks and another bond forms, RGATS simulations primarily sample different versions of an “intermediary state” topology, in which both bonds are absent. Two harmonic distance restraints of the functional form: K1(r1 − r10)2 and K2(r2 − r20)2 are imposed, with force constants K1 and K2, interatomic distances r1 and r2, and minima r10 and r20. The r10 and r20 minima are gradually changed to convert a reactant-like state to a productlike state, or vice versa. If the bond being broken corresponds to r1 and the one being formed corresponds to r2, a canonical B


Article


Figure 2. Proton transfer between two methylimidazole model compounds. (A) Overlay of 50 random starting conformations of the donor methylimidazole with the exchanging proton shown as a pink sphere. (B) Overlay of 50 reactive conformations generated using RGATS. (C) Location of predicted reactive conformations (black squares) in the 2D energy landscape along the two torsions governing relative orientation of the two methylimidazoles, as shown in part D.

the X, Y, and Z axes with the magnitude and direction of these rotations based on a new random number generated for each rotation (Figure 2). A harmonic restraint was applied on the HD1−ND1 distance mediating the hydrogen bond between the two imidazoles with a force constant of 453 kcal/mol·Å2 and a minimum length of 1.0 Å (mimicking the actual ND1− HD1 covalent bond). Each of the 50 starting orientations was then optimized using 5000 SD steps and 5000 ABNR steps of minimization with an energy tolerance of 0.0001 kcal/mol, followed by 5000 1-fs timesteps of Langevin dynamics with a friction coefficient of 5.0 ps−1 applied to all non-hydrogen atoms at a temperature of 300 K, followed by another 5000 SD steps and 5000 ABNR steps of minimization with an energy tolerance of 0.0001 kcal/mol. Concerted Proton Transfer in Amtb. A previous published solvated and equilibrated Amtb system embedded in a DMPC membrane43 was provided by the Lamoureux group. It consisted of the Amtb protein, 185 DMPC molecules, an ammonium ion, 33 potassium ions, 35 chloride ions, and 13 057 water molecules in a box with dimensions 108.3 Å × 97.2 Å × 93.2 Å, with the neutral His168 protonated at the δ nitrogen (ND1) and the neutral His318 protonated at the ϵ nitrogen (NE2). The Amtb system RGATS simulations were performed using the NAMD version 2.9 with the CHARMM22 protein and CHARMM27 lipid force fields. Periodic boundary conditions were used with image centering for all molecules, and the long-range electrostatic interactions were calculated using the particle mesh Ewald method.44 Local interactions were calculated at each 1 fs time step, and the full electrostatic calculation was performed every 2 fs using a multiple time-step scheme.45 Covalent bonds to hydrogens were constrained using the RATTLE46 algorithm, except for water, where they were constrained using the SETTLE47 algorithm. The nonbonded interaction cutoff for van der Waals interactions was initiated at 14 Å and completed at 16 Å, and nonbond lists were generated every 20 fs with a 17 Å cutoff. The simulations were run in the

in this fashion can further limit the already small structural perturbations associated with the topology switches. The partial charges and geometric restraints used for the reactions modeled in the present study are based on chemical intuition. The only exceptions are the geometric restraints used for the proton transfer between histidines, which are informed by the ab initio data shown in the Figure S1 of the Supporting Information. The topologies and parameters used for the RGATS simulations in this study are included in Tables S3− S15 of the Supporting Information. These include the topologies of a protonated methylimidazole (Table S3); a neutral methyimidazole (Table S4); a patch for generating an intermediary state for the proton transfer between two methylimidazoles (Table S5); an ammonium ion (Table S6); a patch for the intermediary state (Table S7) and the product state (Table S8) for a concerted proton transfer between an ammonium ion and a pair of histidine residues; a patch to generate a deoxyribo-triphosphate (Table S9); a patch to generate diphosphate (Table S10); a patch to deprotonate the 3′-hydroxyl (Table S11); and a patch to generate the intermediary state (Table S12) and the product state (Table S13) for nucleotide addition. Methylimidazole Proton Transfer. To generate the 2D energy landscape of relative orientations between two methyl imidazoles, the initial structure was obtained from the two histidines (His168 and His318) in the AmtB ammonia transporter (PDB ID: 1U7G40). An intermediary state conformation was generated and two symmetric torsions (ND1-HD1-ND1-CG) were sampled at 10° increments by changing them in internal coordinate representation, fixing the atoms in these torsions, and then minimizing the rest of the system using 5000 steps of steepest descent (SD)41 and 5000 steps of adopted basis Newton−Raphson (ABNR)42 minimization with an energy tolerance of 0.0001 kcal/mol using CHARMM. Fifty random orientations of the protonated methyl imidazole were generated by serial rotations around C


Article


Figure 3. Concerted proton transfer step in the Amtb ammonia transporter modeled by RGATS. (A) Top and side views of the explicit solvent system used for the modeling.43 (B) Depiction of the reactant, intermediary, and product states for the ammonium ion and the conserved histidine pair modeled using an RGATS simulation. (C) Partial atomic charges for the ammonia molecule in the reactant, intermediary, and product states. (D) Cumulative distribution of protein interaction energies (Coulombic + van der Waals) with the NH3 component of the ammonium ion for the reactant, intermediary, and product states obtained by pooling five RGATS simulations. Averages are shown as vertical lines.

NPT (constant number of particles, pressure, and temperature) ensemble, with the pressure controlled using a Nose−Hoover Langevin piston48 with a target pressure of 1 atm, a period of 100 fs, and a damping time scale of 50 fs. The Langevin thermostat used for temperature control had a damping constant of 5.0 ps−1. For each of the 5 RGATS simulations, the system was minimized for 1000 steps, followed by 10 000 timesteps of dynamics at a temperature of 300 K in the reactant state. A

topology patch was applied to delete the bonds for the exchanging hydrogens on the ammonium ion and the His168 residue, and harmonic restraints with a force constant of 100.0 kcal/mol·Å2 were applied instead. These restraints were then modulated to gradually convert the resulting reactant-like intermediary state to the product-like intermediary state. The reactant-like intermediary state had the following distance restraint minima: 1.0 Å between the ammonium ion NZ and HZ2 atoms (r1), 1.8 Å between the ammonium ion HZ2 atom D


Article

Journal of Chemical Information and Modeling and the His168 NE2 atom (r2), 1.0 Å between the His168 ND1 and HD1 atoms (r3), and 1.0 Å between the His168 HD1 atom and the His318 ND1 atom (r4). The reactant-like intermediary state was first optimized using 1000 steps of minimization followed by 2000 timesteps of dynamics, and then its transition to the product-like state was enforced by changing (r1, r2, r3, r4) in 0.1 Å increments, with 2000 timesteps of dynamics at each increment, first from (1.8, 1.0, 1.0, 1.8 Å) to (1.3, 1.0, 1.0, 1.3 Å), then to (1.3, 1.3, 1.3, 1.3 Å), then to (1.0, 1.3, 1.3, 1.0 Å), and finally to (1.0, 1.8, 1.8, 1.0 Å). A patch was then applied to form the bonds between the ammonium ion HZ2 atom and the His168 NE2 atom, and the His168 HD1 atom and the His318 ND1 atom, to form the product state. This state was then minimized for 1000 steps, followed by 10 000 timesteps of dynamics. The entire RGATS simulation required only 3,000 steps of minimization and 51 000 timesteps of dynamics, which took less than 5 h on two 2.26 GHz Quad-Core Intel Xeon processors. The reactant and product states were also separately equilibrated further using an additional 1 000 000 dynamics timesteps for each RGATS simulation. Interaction energies were computed by setting an infinite cutoff and calculating Coulombic and van der Waals interactions between the ammonia molecule composed of its four constant atoms (N, H1, H3, H4) and all protein atoms (Figure 3). DNA Polymerase μ Nucleotide Addition. The DNA polymerase μ complex was generated using the PDBreader module of CHARMM-GUI49 from the crystal structure with the PDB identifier: 4M04.50 All crystallographic water molecules were maintained, and the coordinates for the third Mg ion were obtained from the extra Mg ion in the crystal structure with the PDB identifier: 4M0A.50 The structure of the ternary complex with either 2 Mg ions or 3 Mg ions was solvated using a preequilibrated cubic TIP3P water box whose each dimension was 86 Å, and all newly introduced water molecules within 2.2 Å of any non-hydrogen atom were deleted. To represent a 0.15 M salt concentration, and to neutralize the system, 58 sodium and 45 chloride ions were added in random positions and their locations were optimized using 100 Monte Carlo steps. The final system contained 63 884 atoms, with 19 332 water molecules, 389 of which were from the crystal structure. The RGATS simulation details were the same as those described above for the Amtb system, except as noted below. The nucleotide addition reaction involves formation of a bond between the O3′ atom of the last adenine nucleotide (Ade4) on the primer strand and the α phosphorus atom (PA) on the incoming dTTP molecule, and cleavage of the bond between the PA and O3A atoms on dTTP, which links the diphosphate moiety to the rest of dTTP (Figure 4). For each of the 5 RGATS simulations, the reactant state was minimized for 1000 steps, followed by 10 000 timesteps of dynamics at a temperature of 300 K. A topology patch was applied to delete the PA−O3A bond, and a harmonic restraint with a force constant of 100.0 kcal/mol·Å2 and a minimum of 1.6 Å (r1) was applied instead. The same harmonic restraint was also applied to the O3′−PA distance with a minimum of 3.0 Å (r2). This reactant-like intermediary state was first optimized using 1000 steps of minimization followed by 10 000 timesteps of dynamics. Its transition to the product-like state was enforced by changing (r1, r2) in 0.1 Å increments, with 10 000 timesteps of dynamics at each increment, first from (1.6, 3.0 Å) to (1.6, 1.8 Å), then to (1.8, 1.8 Å), then to (1.8, 1.6 Å), and finally to

Figure 4. Depiction of the reactant, intermediary, and product states of a nucleotide addition step catalyzed by DNA polymerase μ, modeled using RGATS in the presence of either two or three Mg2+ ions.

(3.0, 1.6 Å). For all intermediary state simulations, constant angular restraints were also imposed using harmonic restraint functions with a force constant of 1000 kcal/mol·Å2 to structurally guide the expected transition pathway. The angular minima for the O3′−PA−O1A, O3′−PA−O2A, O3′−PA−O5′, O3A−PA−O1A, O3A−PA−O2A, and O3A−PA−O5′ angles were set to 90°, and for the O5′−PA−O1A, O1A−PA−O2A, and O2A−PA−O5′, angles were set to 120°. A patch was then applied to form the bond between the O3′ and PA atoms to form the product state. This state was then minimized for 1000 steps, followed by 10 000 timesteps of dynamics. The entire RGATS simulation required only 3000 steps of minimization and 310 000 timesteps of dynamics, which took less than 27 h on two 2.26 GHz Quad-Core Intel Xeon processors. Separate equilibrations for the reactant and product states were carried out using 1 000 000 dynamics timesteps for each RGATS simulation. Similar to the Amtb system, interaction energies were computed by setting an infinite cutoff and calculating Coulombic and van der Waals interactions between all protein atoms and (top) the F1 fragment (base + sugar + α phosphate), (middle) the F2 fragment (diphosphate), or (bottom) the entire dTTP molecule (Figure 5). E


Article


Figure 5. Distributions for the interaction energy (Coulombic + van der Waals) between the DNA polymerase μ protein and the incoming deoxyribo-thymine trinucleotide (dTTP) residue. Each panel shows the cumulated distribution of protein interaction energies (for fragments or the whole dTTP molecule) for the reactant, intermediary, and product states obtained from five separate RGATS simulations. The left column of panels shows the results when only two Mg ions are present, the right column of panels shows the results when three Mg ions are present. As depicted on the right, the top row of panels is for the F1 fragment with the base, sugar, and α phosphate that gets transferred to the primer strand; the middle row of panels is for the F2 diphosphate fragment that dissociates; and the bottom row of panels is for the whole dTTP molecule.

■

RESULTS AND DISCUSSION RGATS Strategy. Even highly complex or concerted biochemical reactions can be deconstructed into individual components that involve the cleavage of one covalent bond and the formation of another covalent bond. In the RGATS strategy, illustrated in Figure 1A, each such component can be modeled using two switches in atomic bonding. The specific implementation of this strategy in any MM program is very straightforward and involves the five steps shown in Figure 1B. The reactant, intermediary, and product states are the three distinct MM topological states, which can be connected by the two topology switches. The first topology switch is the breakage of the first covalent bond, which enables the user-level manipulation of the structural changes in the reactive region through simple geometric restraints. The user can then guide the reactive region from a reactant-like state that resembles the reactant except for the presence of the first covalent bond, to an intermediary state that is between the reactant and product states, to a product-like state, which resembles the product, except for the presence of the second covalent bond. The second topology switch is the formation of the second covalent bond, which completes the transformation to the product state. The minimal geometric restraints that would be required would be ones changing the distances corresponding to the bonds being formed and broken. The optimal geometric restraints would also include multiple angular and torsion restraints that help mimic the expected transition pathway for the reactive region and facilitate a gradual transition from the reactant to the product state. These restraints ensure that the topology switches do not structurally perturb the system, since they are

gradually altered to assume a value compatible with the next topological state prior to each switch. Optimizing Reactant State Conformations. A frequent problem, encountered in modeling biochemical reactions that involve conformational changes associated with catalysis, is that the crystallized conformation may not be in an optimal state for the reaction to occur.51−55 To see if the RGATS strategy can inexpensively deal with this problem of bad local starting conformations, a simple test system modeling the transfer of a proton between two methylimidazole molecules was considered. This system mimics the proton transfer between two histidine sidechains. Fifty random orientations of one doubly protonated methylimidazole were generated with respect to a fixed conformation of the neutral methylimidazole, as shown in Figure 2A. The exchanging proton (pink sphere) in some of these conformations is in an orientation not conducive to transfer to the neutral methylimidazole. When these 50 random structures were minimized using the RGATS strategy, we obtain the 50 reactive structures as shown in Figure 2B, which show the proton poised for transfer. These structures were then analyzed to locate them on the potential energy surface shown in Figure 2C as a function of the two torsion angles defined in Figure 2D. The black squares in the plot in Figure 2C are the 50 RGATS minimized structures, and they appear to fall in mainly four lower energy (blue) regions. This localization of the predicted structures in reasonable lowest energy regions suggests that nonoptimal starting conformations for the reactant state can be optimized to reactive conformations by the RGATS strategy. The two methylimidazoliums are not expected to be coplanar in optimal orientations for proton transfer (Supporting Information Figure S1), and this feature is F


Article


its diffusion down the pore, which is a slower process. This process is likely to occur in an effective mean field of a tautomeric mixture of the two histidines exchanging a proton with each other.43 The ability to capture this feature demonstrates that RGATS simulations can be used to obtain preliminary insights into reaction mechanisms, which can then be refined using further QM/MM simulations or experimental approaches. Probing the Effects of Specific Differences in Reactive Environments. DNA polymerases catalyze DNA replication by helping to extend a primer DNA strand paired to a longer template DNA strand. They catalyze addition of a deoxyribonucleotide to the primer strand that is complementary to the next unpaired base on the template strand.61 The reaction involves nucleophilic attack by the 3′-oxygen of the last nucleotide of the primer strand on the α-phosphorus atom of the incoming trinucleotide. The bond between this αphosphorus and the oxygen that connects it to the other two phosphate groups gets cleaved.62 This reaction requires Mg2+ ions, with the established understanding being that two Mg2+ ions are required.63 More recent studies point to the possibility of three Mg2+ ions being required,64−67 which is similar to the revision of the alkaline phosphatase reaction mechanism.68 From a purely structural perspective, if the 3′-oxygen has already been deprotonated, this reaction also involves only one bond O−P formation and one P−O bond cleavage. This reaction was modeled using five separate RGATS simulations for both the cases: in the presence of two as well as three Mg2+ ions near the active site of the human DNA polymerase μ.50,69 Figure 4 shows the reactant, intermediary, and product state conformations obtained for the two different metal ion environments (examples of full RGATS simulations shown in Supplementary Movies S2 and S3). There is a rearrangement of partial charges within the fragmenting incoming nucleotide during the nucleotide transfer reaction, in concert with the bond cleavage and formation. These changes can be represented at the MM level and are expected to alter the interaction of the polymerase with the incoming deoxyribo-thymine triphosphate (dTTP). The distributions of protein interaction energy with dTTP, and its two fragments separated during the reaction, are shown in Figure 5. Since the RGATS simulations used to obtain these distributions are run in the presence of either two or three Mg2+ ions near the active site, the effect of this environmental difference can be assessed through these distributions. For the overall dTTP molecule, the interaction with DNA polymerase μ is weaker when there are 3 Mg2+ ions in all three topological states, with the largest difference seen for the product state (see also Supporting Information Table S2). For the fragment F1 that gets added to the primer strand (composed of the thymine base, sugar, and 1 phosphate), the product state shows lowered interaction with the protein than the reactant or intermediary states. When two Mg2+ ions are present, this interaction is repulsive to a greater degree than when three Mg2+ ions are present. For the fragment F2 (composed of the diphosphate), the product state interaction with the protein is more favorable than the reactant or intermediary states when two Mg2+ ions are present, and less favorable than the reactant or intermediary states when three Mg2+ ions are present. The interaction between the F1 fragment and the protein shifts from weakly favorable to repulsive upon completion of the reaction, which presages the motion of the primer strand incorporated nucleotide away from the active site. The interaction of the

also captured by the RGATS generated structures. The time scale accessible to MM simulations keeps increasing,56,57 and enhanced sampling methods are improved regularly.58 Large conformational changes in nonreactive regions that occur in response to a chemical reaction are therefore also likely to be amenable to the RGATS strategy. Preliminary Insights for Concerted Reactions in Complex Environments. The transfer of protons occurs in both aqueous and membrane environments to facilitate macromolecular functions.59 One example is the transport of an ammonium ion in the ammonium transporter proteins of the Amtb and Rhesus families.40 The ammonium transport has been predicted to involve the following steps: a temporary conversion of the ammonium ion to an ammonia molecule by the transfer of a proton to ϵ nitrogen of a histidine side chain (His168); proton transfer between the δ nitrogens of a pair of histidine sidechains (His168 and His318) in the channel; facile transport of the neutral ammonia molecule across the hydrophobic cell membrane;60 ammonia regaining an ϵ proton from the second histidine side chain (His318) as it exits the pore to reform the ammonium ion.43 The first two proton transfers, ammonium ion to His168, and His168 to His318, can occur either in a serial or concerted fashion. The latter possibility was modeled by five RGATS explicit solvent simulations for Amtb embedded in a lipid membrane (Figure 3A). Figure 3B depicts the reactant, intermediary, and product states for this concerted reaction (extracted from one of the RGATS simulations shown in Supplementary Movie S1). In these simulations, the intermediary state is where the two relevant N−H distances for the two simultaneous proton transfers are 1.3 Å each (see the Supporting Information text and Figure S1 for the rationale behind this choice for the intermediary distance). The RGATS approach is primarily suitable for gradually facilitating the structural transition accompanying a biochemical reaction. The topology switch between the reactant and intermediary states minimally involves deletion of the bond being broken, and elimination of the associated bonded terms. To obtain a more accurate structural response to the reaction, the partial atomic charges in the intermediary state can be set to those expected for the actual transition state. They can also be smoothly modified from the reactant state charges to the product state charges during the structural transition within the intermediary state. In modeling the concerted proton transfer in Amtb, the partial charges for the ammonium ion were altered between the reactant, intermediary, and product states, as shown in Figure 3C. For example, the partial charge on the nitrogen atom in the ammonium ion changes from −0.32 in the reactant state to −0.65 in the intermediary state to −0.99 in the product state (i.e., the neutral ammonia molecule). The other associated partial charge changes occur in the two histidine sidechains, His168 and His318, that also participate in the concerted proton transfer. If the direct electrostatic and van der Waals interaction between the Amtb protein and the ammonia molecule is assessed for just the N, H1, H2, H3 ammonia atoms that remain bonded together through the reaction, the interaction energy distributions are very different for the product state than either the reactant or intermediary states, which are similar to each other (Figure 3D). This shows that the rearrangement of partial atomic charges that occurs in transitioning from the intermediary state to the product state greatly reduces the affinity of the unfragmented part of the ammonium ion for its first binding site. This is expected to ease G


Article


multiple different reaction mechanisms in a rapid fashion. A good example would be a scenario where multiple protein residues can play the role of the proton acceptor, and each of the postulated proton transfer reactions could be simulated and compared using RGATS. These simulations could either provide answers directly, or could provide reasonable starting structures for further QM/MM simulations. In conclusion, conformational features of chemical reaction pathways can be modeled using standard macromolecular MM energy functions and programs through the RGATS strategy. This strategy yields reasonable biochemical reaction trajectories within complex macromolecular environments, at a computational cost that is orders of magnitude lower than QM/MM simulations. It also does not require parametrization of transition states beyond imposition of reasonable geometric restraints for approximating hypothesized transition pathways. It can rapidly predict structures for reactant, intermediate, or product states starting from an experimental or modeled structure in any one of these states. Covalent bonds provide a strong restraint for the conformational variability of substrates and products78 and could be a basis for reasonable structure prediction.79,80 For example, starting from just a product state conformation of the enzyme, a product molecule could be docked into the active site, and the reverse biochemical reaction could be modeled using RGATS to obtain a reasonable structural model for the substrate−enzyme complex reactant state. As the RGATS strategy is implemented through simple topology file modifications and standard geometric restraints in existing MM simulation programs, the human learning curve costs for modeling chemical reaction pathways are significantly lowered. It also allows for all methodological advances implemented in these MM program packages, such as specific implicit solvent models,81 constant pH simulations,82−84 polarizability,85 enhanced sampling methods,86 or optimized parallel processing algorithms and architectures,87 to be directly utilized to tackle a range of problems related to biochemical reactions.

F2 diphosphate fragment remains favorable, but to a lesser extent when three Mg2+ ions are present, pointing to the possibility that the presence of the third Mg2+ ion could enhance the dissociation of this byproduct to increase the processivity of the polymerase. These observations show how the effect of an environmental change as small as the presence of a single additional divalent cation is also amenable to RGATS analysis.

■

CONCLUSIONS In an active site environment poised for a chemical reaction to occur, the transitions between the reactant and product states are likely to be faster than the dynamics of the surrounding environment. The RGATS strategy can enforce the chemical transitions to occur much more gradually, thereby allowing the environment to adjust to each minuscule internal structural change accompanying the chemical transition. This can be seen clearly in the movies included in the Supporting Information, where the environment can fluctuate while the chemical transition occurs very gradually, thereby allowing greater sampling for each intermediary conformation. This is distinct from application of additional umbrella potentials in reactive simulations70 in that the internal energetics of the bond undergoing transformation is not represented. While this prevents assessment of free energy profiles for the intrinsic reaction itself, as long as the geometric details of the chemical change are represented accurately, the response of the environment to the chemical reaction is likely to be captured comprehensively. Similarly, the effect of any compositional change in the environment, such as ligand binding or mutations, can also be preliminarily assessed using RGATS, which makes it a useful supplementary tool for the design of allosteric inhibitors or proteins. Since RGATS uses MM models, and each reaction step involves relatively small atomic rearrangements, reasonable trajectories for complex biochemical reactions can be generated at a very small computational cost. The reliance on simple geometrical restraints and topology file patches makes the RGATS strategy straightforward to apply using any standard biomacromolecular MM software package such as CHARMM, 29,71 AMBER, 72 GROMACS, 73 Desmond, 74 NAMD,30 or OpenMM.75 MM energy functions are not expected to have enough terms to represent the intermediary states accurately, therefore an initial hypothesis for the internal geometrical changes accompanying the transition from reactant to product state is desirable. Application of angular, torsion, and chirality restraints based on this hypothesis can maintain the internal geometrical features of the intermediary state, while obtaining the response of the environment in an unrestrained fashion. This also allows multiple hypotheses about the reaction pathway, e.g. associative versus dissociative or concerted versus stepwise, to be modeled explicitly. The user-level control also makes it very straightforward to develop customized or automated protocols to explore any multistep biochemical reaction pathway. In the present work, the topology switches are made directly, but in the future, these could be implemented through methods analogous to thermodynamic integration,76 free energy perturbation,77 or bias annealing.39 The low cost of the approach also allows iterative sampling of forward and backward RGATS trajectories, which allows testing for convergence of the structural response of the environment. It also becomes feasible to visualize and assess the likelihood of

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jcim.7b00699. Further details, figure, and tables as described in the text (PDF) Video showing concerted proton transfer between an ammonium ion and two conserved histidines in the AmtB ammonia transporter (MPG) Video showing the addition of a nucleotide to the primer strand in DNA polymerase μ in the presence of two Mg2+ ions near the active site (MPG) Video showing and the addition of a nucleotide to the primer strand in DNA polymerase μ in the presence of three Mg2+ ions near the active site (MPG)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Mailing address: CMS 2008, Biggs Laboratory, Wadsworth Center, New York State Department of Health Empire State Plaza, PO Box 509, Albany, New York, 12201-0509. Tel.: 518-474-0569. Fax: 518-4743181. Website: http://www.wadsworth.org/senior-staff/nileshbanavali. H


Article

Journal of Chemical Information and Modeling ORCID

(18) Warshel, A.; Bora, R. P. Simulating Enzyme Reactivity: Computational Methods in Enzyme Catalysis. The Royal Society of Chemistry 2016, 1−30. (19) Jensen, F. Locating minima on seams of intersecting potential energy surfaces. An application to transition structure modeling. J. Am. Chem. Soc. 1992, 114, 1596−1603. (20) Rappe, A.; Pietsch, M.; Wiser, D.; Hart, J.; Bormann-Rochotte, L.; Skiff, W. RFF, conceptual development of a full periodic table force field for studying reaction potential surfaces. Mol. Eng. 1997, 7, 385− 400. (21) Norrby, P.-O.; Liljefors, T. Automated molecular mechanics parameterization with simultaneous utilization of experimental and quantum mechanical data. J. Comput. Chem. 1998, 19, 1146−1166. (22) Kim, Y.; Corchado, J.; Villa, J.; Xing, J.; Truhlar, D. Multiconfiguration molecular mechanics algorithm for potential energy surfaces of chemical reactions. J. Chem. Phys. 2000, 112, 2718. (23) Van Duin, A.; Dasgupta, S.; Lorant, F.; Goddard, W., III ReaxFF: a reactive force field for hydrocarbons. J. Phys. Chem. A 2001, 105, 9396−9409. (24) Weill, N.; Corbeil, C. R.; De Schutter, J. W.; Moitessier, N. Toward a computational tool predicting the stereochemical outcome of asymmetric reactions: Development of the molecular mechanicsbased program ACE and application to asymmetric epoxidation reactions. J. Comput. Chem. 2011, 32, 2878−2889. (25) Ponder, J.; Case, D. Force fields for protein simulations. Adv. Protein Chem. 2003, 66, 27−85. (26) Mackerell, A. D. Empirical force fields for biological macromolecules: overview and issues. J. Comput. Chem. 2004, 25, 1584− 1604. (27) Lemkul, J. A.; Huang, J.; Roux, B.; MacKerell, A. D., Jr An empirical polarizable force field based on the classical drude oscillator model: development history and recent applications. Chem. Rev. 2016, 116, 4983−5013. (28) Marelius, J.; Kolmodin, K.; Feierberg, I.; Åqvist, J. Q: a molecular dynamics program for free energy calculations and empirical valence bond simulations in biomolecular systems. J. Mol. Graphics Modell. 1998, 16, 213−225. (29) Brooks, B. R.; Brooks, C. L., III; MacKerell, A. D., Jr.; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; et al. CHARMM: the biomolecular simulation program. J. Comput. Chem. 2009, 30, 1545−1614. (30) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. Scalable molecular dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781− 1802. (31) MacKerell, A., Jr.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; et al. All-atom empirical potential for molecular modeling and dynamics studies of proteins. J. Phys. Chem. B 1998, 102, 3586−3616. (32) Foloppe, N.; MacKerell, A., Jr All-atom empirical force field for nucleic acids: I. Parameter optimization based on small molecule and condensed phase macromolecular target data. J. Comput. Chem. 2000, 21, 86−104. (33) MacKerell, A.; Banavali, N. All-atom empirical force field for nucleic acids: II. Application to molecular dynamics simulations of DNA and RNA in solution. J. Comput. Chem. 2000, 21, 105−120. (34) Klauda, J. B.; Brooks, B. R.; MacKerell, A. D.; Venable, R. M.; Pastor, R. W. An ab initio study on the torsional surface of alkanes and its effect on molecular simulations of alkanes and a DPPC bilayer. J. Phys. Chem. B 2005, 109, 5300−5311. (35) Jorgensen, W.; Chandrasekhar, J.; Madura, J.; Impey, R.; Klein, M. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926. (36) Beglov, D.; Roux, B. Finite representation of an infinite bulk system: solvent boundary potential for computer simulations. J. Chem. Phys. 1994, 100, 9050−9063. (37) MacKerell, A., Jr; Banavali, N.; Foloppe, N. Development and current status of the CHARMM force field for nucleic acids. Biopolymers 2000, 56, 257−265.

Nilesh K. Banavali: 0000-0003-2206-7049 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by a SUNY Research Foundation FRAP-A award (Grant Number 1123811-68578) and used computing resources provided by an Extreme Science and Engineering Discovery Environment (XSEDE) allocation, which is supported by National Science Foundation grant number ACI-1548562. We are thankful to Shihao Wang and Guillaume Lamoureux for providing us with the equilibrated AmtB structural model.

■

REFERENCES

(1) Šrajer, V.; Teng, T.; Ursby, T.; Pradervand, C.; Ren, Z.; Adachi, S.; Schildkamp, W.; Bourgeois, D.; Wulff, M.; Moffat, K. Photolysis of the carbon monoxide complex of myoglobin: nanosecond timeresolved crystallography. Science 1996, 274, 1726−1729. (2) Schotte, F.; Lim, M.; Jackson, T.; Smirnov, A.; Soman, J.; Olson, J.; Phillips, G., Jr; Wulff, M.; Anfinrud, P. Watching a protein as it functions with 150-ps time-resolved X-ray crystallography. Science 2003, 300, 1944−1947. (3) Schmidt, M.; Pahl, R.; Srajer, V.; Anderson, S.; Ren, Z.; Ihee, H.; Rajagopal, S.; Moffat, K. Protein kinetics: Structures of intermediates and reaction mechanism from time-resolved x-ray data. Proc. Natl. Acad. Sci. U. S. A. 2004, 101, 4799. (4) Haumann, M.; Liebisch, P.; Müller, C.; Barra, M.; Grabolle, M.; Dau, H. Photosynthetic O2 formation tracked by time-resolved X-ray experiments. Science 2005, 310, 1019−1021. (5) Yang, X.; Ren, Z.; Kuk, J.; Moffat, K. Temperature-scan cryocrystallography reveals reaction intermediates in bacteriophytochrome. Nature 2011, 479, 428−432. (6) Fu, Z.; Kaledhonkar, S.; Borg, A.; Sun, M.; Chen, B.; Grassucci, R. A.; Ehrenberg, M.; Frank, J. Key intermediates in ribosome recycling visualized by time-resolved cryoelectron microscopy. Structure 2016, 24, 2092−2101. (7) Srajer, V.; Schmidt, M. Watching proteins function with timeresolved x-ray crystallography. J. Phys. D: Appl. Phys. 2017, 50, 373001. (8) Friesner, R.; Guallar, V. Ab initio quantum chemical and mixed quantum mechanics/molecular mechanics (QM/MM) methods for studying enzymatic catalysis. Annu. Rev. Phys. Chem. 2005, 56, 389− 427. (9) Lin, H.; Truhlar, D. QM/MM: what have we learned, where are we, and where do we go from here? Theor. Chem. Acc. 2007, 117, 185− 199. (10) Bell, A.; Head-Gordon, M. Quantum Mechanical Modeling of Catalytic Processes. Annu. Rev. Chem. Biomol. Eng. 2011, 2, 453−477. (11) van der Kamp, M. W.; Mulholland, A. J. Combined quantum mechanics/molecular mechanics (QM/MM) methods in computational enzymology. Biochemistry 2013, 52, 2708−2728. (12) Olsson, M. A.; Ryde, U. Comparison of QM/MM methods to obtain ligand-binding free energies. J. Chem. Theory Comput. 2017, 13, 2245−2253. (13) Field, M. Simulating enzyme reactions: challenges and perspectives. J. Comput. Chem. 2002, 23, 48−58. (14) Senn, H.; Thiel, W. Atomistic Approaches in Modern Biology; Springer, 2007; pp 173−290. (15) Warshel, A.; Weiss, R. An empirical valence bond approach for comparing reactions in solutions and in enzymes. J. Am. Chem. Soc. 1980, 102, 6218−6226. (16) Warshel, A.; Weiss, R. Empirical valence bond calculations of enzyme catalysis. Ann. N. Y. Acad. Sci. 1981, 367, 370−382. (17) Kamerlin, S.; Warshel, A. The EVB as a quantitative tool for formulating simulations and analyzing biological and chemical reactions. Faraday Discuss. 2010, 145, 71−106. I


Article

Journal of Chemical Information and Modeling (38) MacKerell, A., Jr; Feig, M.; Brooks, C., III Improved treatment of the protein backbone in empirical force fields. J. Am. Chem. Soc. 2004, 126, 698−699. (39) Hu, J.; Ma, A.; Dinner, A. Bias annealing: A method for obtaining transition paths de novo. J. Chem. Phys. 2006, 125, 114101. (40) Khademi, S.; O’Connell, J., III; Remis, J.; Robles-Colmenares, Y.; Miercke, L.; Stroud, R. Mechanism of ammonia transport by Amt/ MEP/Rh: structure of AmtB at 1.35 A. Science 2004, 305, 1587−1594. (41) Curry, H. The method of steepest descent for nonlinear minimization problems. Q. Appl. Math. 1944, 2, 258−261. (42) Chu, J.; Trout, B.; Brooks, B. A super-linear minimization scheme for the nudged elastic band method. J. Chem. Phys. 2003, 119, 12708. (43) Wang, S.; Orabi, E. A.; Baday, S.; Berneche, S.; Lamoureux, G. Ammonium transporters achieve charge transfer by fragmenting their substrate. J. Am. Chem. Soc. 2012, 134, 10419−10427. (44) Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: An N log (N) method for Ewald sums in large systems. J. Chem. Phys. 1993, 98, 10089. (45) Batcho, P. F.; Case, D. A.; Schlick, T. Optimized particle-mesh Ewald/multiple-time step integration for molecular dynamics simulations. J. Chem. Phys. 2001, 115, 4003−4018. (46) Andersen, H. C. Rattle: A ”velocity” version of the shake algorithm for molecular dynamics calculations. J. Comput. Phys. 1983, 52, 24−34. (47) Miyamoto, S.; Kollman, P. A. SETTLE: an analytical version of the SHAKE and RATTLE algorithm for rigid water models. J. Comput. Chem. 1992, 13, 952−962. (48) Martyna, G. J.; Tobias, D. J.; Klein, M. L. Constant pressure molecular dynamics algorithms. J. Chem. Phys. 1994, 101, 4177−4189. (49) Jo, S.; Kim, T.; Iyer, V. G.; Im, W. CHARMM-GUI: a webbased graphical user interface for CHARMM. J. Comput. Chem. 2008, 29, 1859−1865. (50) Moon, A. F.; Pryor, J. M.; Ramsden, D. A.; Kunkel, T. A.; Bebenek, K.; Pedersen, L. C. Sustained active site rigidity during synthesis by human DNA polymerase μ. Nat. Struct. Mol. Biol. 2014, 21, 253−260. (51) Waksman, G.; Krishna, T. S.; Williams, C. H.; Kuriyan, J. Crystal Structure Of Escherichia Coli Thioredoxin Reductase Refined at 2 Å Resolution: Implication for a Large Conformational Change during Catalysis. J. Mol. Biol. 1994, 236, 800−816. (52) Huang, H.; Chopra, R.; Verdine, G. L.; Harrison, S. C. Structure of a covalently trapped catalytic complex of HIV-1 reverse transcriptase: implications for drug resistance. Science 1998, 282, 1669− 1675. (53) Ramakrishnan, B.; Qasba, P. K. Crystal structure of lactose synthase reveals a large conformational change in its catalytic component, the β1, 4-galactosyltransferase-I. J. Mol. Biol. 2001, 310, 205−218. (54) Gutteridge, A.; Thornton, J. Conformational changes observed in enzyme crystal structures upon substrate binding. J. Mol. Biol. 2005, 346, 21−28. (55) Berndsen, C. E.; Wolberger, C. New insights into ubiquitin E3 ligase mechanism. Nat. Struct. Mol. Biol. 2014, 21, 301−307. (56) Shaw, D. E.; Maragakis, P.; Lindorff-Larsen, K.; Piana, S.; Dror, R. O.; Eastwood, M. P.; Bank, J. A.; Jumper, J. M.; Salmon, J. K.; Shan, Y.; Wriggers, W. Atomic-level characterization of the structural dynamics of proteins. Science 2010, 330, 341−346. (57) Lane, T. J.; Shukla, D.; Beauchamp, K. A.; Pande, V. S. To milliseconds and beyond: challenges in the simulation of protein folding. Curr. Opin. Struct. Biol. 2013, 23, 58−65. (58) Bernardi, R. C.; Melo, M. C.; Schulten, K. Enhanced sampling techniques in molecular dynamics simulations of biological systems. Biochim. Biophys. Acta, Gen. Subj. 2015, 1850, 872−877. (59) Ishikita, H.; Saito, K. Proton transfer reactions and hydrogenbond networks in protein environments. J. R. Soc., Interface 2014, 11, 20130518.

(60) Baday, S.; Orabi, E. A.; Wang, S.; Lamoureux, G.; Bernèche, S. Mechanism of NH4+ recruitment and NH3 transport in Rh proteins. Structure 2015, 23, 1550−1557. (61) Sawaya, M. R.; Pelletier, H.; Kumar, A.; Wilson, S. H.; Kraut, J. Crystal structure of rat DNA polymerase beta: evidence for a common polymerase mechanism. Science 1994, 264, 1930−1936. (62) Steitz, T.; Smerdon, S.; Jager, J.; Joyce, C. A unified polymerase mechanism for nonhomologous DNA and RNA polymerases. Science 1994, 266, 2022−2025. (63) Steitz, T. DNA polymerases: structural diversity and common mechanisms. J. Biol. Chem. 1999, 274, 17395−17398. (64) Nakamura, T.; Zhao, Y.; Yamagata, Y.; Hua, Y.-j.; Yang, W. Watching DNA polymerase η make a phosphodiester bond. Nature 2012, 487, 196−201. (65) Perera, L.; Freudenthal, B. D.; Beard, W. A.; Shock, D. D.; Pedersen, L. G.; Wilson, S. H. Requirement for transient metal ions revealed through computational analysis for DNA polymerase going in reverse. Proc. Natl. Acad. Sci. U. S. A. 2015, 112, E5228−E5236. (66) Gao, Y.; Yang, W. Capture of a third Mg2+ is essential for catalyzing DNA synthesis. Science 2016, 352, 1334−1337. (67) Yoon, H.; Warshel, A. Simulating the fidelity and the three Mg mechanism of pol eta and clarifying the validity of transition state theory in enzyme catalysis. Proteins: Struct., Funct., Genet. 2017, 85, 1446−1453. (68) Stec, B.; Holtz, K. M.; Kantrowitz, E. R. A revised mechanism for the alkaline phosphatase reaction involving three metal ions. J. Mol. Biol. 2000, 299, 1303−1311. (69) Jamsen, J. A.; Beard, W. A.; Pedersen, L. C.; Shock, D. D.; Moon, A. F.; Krahn, J. M.; Bebenek, K.; Kunkel, T. A.; Wilson, S. H. Time-lapse crystallography snapshots of a double-strand break repair polymerase in action. Nat. Commun. 2017, 8, 253. (70) Kästner, J.; Senn, H.; Thiel, S.; Otte, N.; Thiel, W. QM/MM free-energy perturbation compared to thermodynamic integration and umbrella sampling: Application to an enzymatic reaction. J. Chem. Theory Comput. 2006, 2, 452−461. (71) Brooks, B.; Bruccoleri, R.; Olafson, B.; Swaminathan, S.; Karplus, M.; States, D. CHARMM: A program for macromolecular energy, minimization, and dynamics calculations. J. Comput. Chem. 1983, 4, 187−217. (72) Case, D. A.; Cheatham, T. E.; Darden, T.; Gohlke, H.; Luo, R.; Merz, K. M.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J. The Amber biomolecular simulation programs. J. Comput. Chem. 2005, 26, 1668−1688. (73) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (74) Bowers, K. J.; Chow, E.; Xu, H.; Dror, R. O.; Eastwood, M. P.; Gregersen, B. A.; Klepeis, J. L.; Kolossvary, I.; Moraes, M. A.; Sacerdoti, F. D.; Salmon, J. K.; Shan, Y.; Shaw, D. E. Scalable algorithms for molecular dynamics simulations on commodity clusters. Proc. ACM/IEEE 2006, 43−43. (75) Eastman, P.; Friedrichs, M. S.; Chodera, J. D.; Radmer, R. J.; Bruns, C. M.; Ku, J. P.; Beauchamp, K. A.; Lane, T. J.; Wang, L.-P.; Shukla, D.; et al. OpenMM 4: a reusable, extensible, hardware independent library for high performance molecular simulation. J. Chem. Theory Comput. 2013, 9, 461−469. (76) Bruckner, S.; Boresch, S. Efficiency of alchemical free energy simulations. II. Improvements for thermodynamic integration. J. Comput. Chem. 2011, 32, 1320−1333. (77) Ullmann, R.; Ullmann, G. A Generalized Free Energy Perturbation Theory Accounting for End States with Differing Configuration Space Volume. J. Phys. Chem. B 2011, 115, 507−521. (78) Xu, Y.; Villa, A.; Nilsson, L. The free energy of locking a ring: Changing a deoxyribonucleoside to a locked nucleic acid. J. Comput. Chem. 2017, 38, 1147−1157. (79) Yassin, A. S.; Haque, M. E.; Datta, P. P.; Elmore, K.; Banavali, N. K.; Spremulli, L. L.; Agrawal, R. K. Insertion domain within mammalian mitochondrial translation initiation factor 2 serves the J


Article

Journal of Chemical Information and Modeling role of eubacterial initiation factor 1. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 3918−3923. (80) Yassin, A. S.; Agrawal, R. K.; Banavali, N. K. Computational exploration of structural hypotheses for an additional sequence in a mammalian mitochondrial protein. PLoS One 2011, 6, e21871. (81) Kleinjung, J.; Fraternali, F. Design and application of implicit solvent models in biomolecular simulations. Curr. Opin. Struct. Biol. 2014, 25, 126−134. (82) Goh, G. B.; Knight, J. L.; Brooks, C. L. Constant pH molecular dynamics simulations of nucleic acids in explicit solvent. J. Chem. Theory Comput. 2012, 8, 36−46. (83) Wallace, J. A.; Shen, J. K. Continuous constant pH molecular dynamics in explicit solvent with pH-based replica exchange. J. Chem. Theory Comput. 2011, 7, 2617−2629. (84) Donnini, S.; Tegeler, F.; Groenhof, G.; Grubmuller, H. Constant pH molecular dynamics in explicit solvent with λ-dynamics. J. Chem. Theory Comput. 2011, 7, 1962−1978. (85) Anisimov, V. M.; Lamoureux, G.; Vorobyov, I. V.; Huang, N.; Roux, B.; MacKerell, A. D. Determination of electrostatic parameters for a polarizable force field based on the classical Drude oscillator. J. Chem. Theory Comput. 2005, 1, 153−168. (86) Liwo, A.; Czaplewski, C.; Ołdziej, S.; Scheraga, H. A. Computational techniques for efficient conformational sampling of proteins. Curr. Opin. Struct. Biol. 2008, 18, 134−139. (87) Shaw, D. E.; Deneroff, M. M.; Dror, R. O.; Kuskin, J. S.; Larson, R. H.; Salmon, J. K.; Young, C.; Batson, B.; Bowers, K. J.; Chao, J. K.; et al. Anton, a special-purpose machine for molecular dynamics simulation. Commun. ACM 2008, 51, 91−97.

K


Structural Articulation of Biochemical Reactions ... - ACS Publications

Recommend Documents