Efficient Empirical Valence Bond Simulations with GROMACS

We describe a protocol to perform empirical valence bond (EVB) simulations using GROMACS software. EVB is a fast and reliable method that allows one to determine the reaction free-energy profiles in complex systems, such as enzymes, by employing classical force fields to represent a chemical reaction. Therefore, running EVB simulations is basically as fast as any classical molecular dynamics simulation, and the method uses standard free-energy calculations to map the free-energy change along a given reaction path. To exemplify and validate our EVB implementation, we replicated two cases of our earlier enzyme simulations. One of these addresses the decomposition of the activation free energy into its enthalpic and entropic components, and the other is focused on calculating the overall catalytic effect of the enzyme compared to the same reaction in water. These two examples give virtually identical results to those obtained with programs that were specifically designed for EVB simulations and show that the GROMACS implementation is robust and can be used for very large systems.


I Preparing the system and building the topologies
1. Calculate the relaxed geometries and ESP charges with Gaussian09 6 for the reacting moieties which include the 3-oxovalerate substrate, nicotinamide adenine dinucleotide cofactor (NADH) and tyrosine 161 (Tyr161) in reactant state (RS) and product state (PS) .We will obtain the following files for RS: 3ov.log (substrate), nah.log (NADH) and twn.log (Tyr161); and for PS we will have: hov.log (substrate), nad.log (NADH), and twm.log (Tyr161).
2. Calculate the RESP charges with antechamber from AmberTools17. 7The command below will repeat for each output file from Gaussian, but for convenience we will show an example only for one of them.antechamber -fi gout -i 3ov.log -fo ac -o 3ov.ac -c resp -nc -1 3. Extract the pdb from the gaussian output file.antechamber -fi gout -i 3ov.log -fo pdb -o 3ov.pdb -rn 3OV we must build only the residues corresponding to RS. NOTE: For now, we must do this step manually; in future versions we will consider building .itpfiles which will allow us to skip the steps from 6 to 8. 15.Build the topologies for EVB; there will be 51 topologies for a 51 frames FEP simulation.gmx4evb.py-f 51 -r 3ov nah twn -p hov nad twm 16.Generate the position restraints files (these files will correspond to the restraints mentioned in step 13).The following example will build a file with 209.2 kJ mol -1 nm -2 (0.5 kcal mol -1 Å -2 ) on EVB atoms and 12.6 kJ mol -1 nm -2 (0.03 kcal mol -1 Å -2 ) for the rest of the protein.
genposre.py--r1 209.2 --r2 12.6 17.Generate tabulated potential files for soft-repulsion interactions.The soft-repulsion potential has the following expression: and it is passed to GROMACS as a tabulated bond of type 9.These files will contain three columns with discrete data corresponding to the following functions:   ,  −  ,  −  .

II Equilibration and production run
1. Copy the position restraints files, the tabulated potential files, topol_000.topand the modified force field directory into the equilibration folder (equil) and then run the equilibration steps sequentially.
2. Copy the last equilibrated .grofile inside the production folder (fep/rep_000) together with the tabulated potential files, all topology files, and the force field directory.First run restart.mdpfile which further equilibrates the system starting from randomized velocities -this way we can provide a different starting point to each FEP replica.Then submit the FEP frames sequentially, passing to each frame its corresponding topology (topol_000.top to fep_000.mdp,topol_001.topto fep_001.mdp and so on).

III Extracting the valence energies and building the EVB profile
1. Move all the trajectory (.trr) files to rerun folder (fep/rep_000/rerun).Copy also the initial .grofile, the force field directory, the tabulated potential files, topol_000.top,topol_050.top(fist and last topology), and evbless.top.This last topology file is a special topology where all the EVB atoms have been converted into dummies.Then rerun all the FEP frames on the coordinates from the trajectory files using the topology corresponding to RS (topol_000.top), then with the PS topology (topol_050.top) and one last time using evbless.toptopology.
2. Extract the potential energies from each rerun.To this end, you can use the bash script get_ene.shfound inside the fep folder.Inside this file, you must give the correct number to nosys and noless variables which correspond to the numbers passed to the energy tool of GROMACS to extract the potential energy from the rerun with the RS and PS topologies (nosys) and with evbless.toptopology (noless).At the end of this process, we will obtain two folders, sysA and sysB, with files containing the extracted potential energies.
3. Calculate the valence energies H11 and H22 (see main text) and write them into a format appropriate for qfep analysis.
gmx2qfep.py-f 51 4. Prepare the input file for qfep (see qfep.inp file inside fep folder).For Q users, this file is the same as for Q simulations. 9,10 ble S1.EVB parameters for the reaction of 3-oxovalerate substrate catalyzed by (R)-3-hydroxybutyrate dehydrogenase from Psychrobacter arcticus enzyme (PaHBDB).The parameters include the RESP charges for the EVB atoms in reactant (RS) and product (PS) states, soft exponential repulsion for the atoms that form or break chemical bonds and Morse potential parameters.The EVB atoms include the side chain of tyrosine 161 (TYR), the 3-carbamoyl-pyridine group of NADH (NAH) and the entire substrate (3OV).The atom names are the same as in topology and are shown in Figure S1.The residue names are the same as in topology as well.

inside aminoacids.rtp file (see 3OV residue inside aminoacids.rtp);
Reacting moieties with atom names.The atom names correspond to those in TableS1and in the topology.

Table S2 .
EVB parameters for the reaction of dihydroxyacetone phosphate (DHAP) substrate catalyzed by triosephosphate isomerase (TIM) enzyme.The parameters include the ffld_server generated charges for the EVB atoms in reactant (RS) and product (PS) states, soft exponential repulsion for the atoms that form or break chemical bonds and Morse potential parameters.The EVB atoms are shown with labels in FigureS2.The labels in this table correspond to those in FigureS2and in topology.The residue names are the same as in topology as well.Reacting moieties with atom names.The atom names correspond to those in TableS2and in the topology.