Computational Investigation of the Interplay of Substrate Positioning and Reactivity in Catechol O-Methyltransferase

Catechol O-methyltransferase (COMT) is a SAM- and Mg2+-dependent methyltransferase that regulates neurotransmitters through methylation. Simulations and experiments have identified divergent catecholamine substrate orientations in the COMT active site: molecular dynamics simulations have favored a monodentate coordination of catecholate substrates to the active site Mg2+, and crystal structures instead preserve bidentate coordination along with short (2.65 Å) methyl donor-acceptor distances. We carry out longer dynamics (up to 350 ns) to quantify interconversion between bidentate and monodentate binding poses. We provide a systematic determination of the relative free energy of the monodentate and bidentate structures in order to identify whether structural differences alter the nature of the methyl transfer mechanism and source of enzymatic rate enhancement. We demonstrate that the bidentate and monodentate binding modes are close in energy but separated by a 7 kcal/mol free energy barrier. Analysis of interactions in the two binding modes reveals that the driving force for monodentate catecholate orientations in classical molecular dynamics simulations is derived from stronger electrostatic stabilization afforded by alternate Mg2+ coordination with strongly charged active site carboxylates. Mixed semi-empirical-classical (SQM/MM) substrate C-O distances (2.7 Å) for the bidentate case are in excellent agreement with COMT X-ray crystal structures, as long as charge transfer between the substrates, Mg2+, and surrounding ligands is permitted. SQM/MM free energy barriers for methyl transfer from bidentate and monodentate catecholate configurations are comparable at around 21–22 kcal/mol, in good agreement with experiment (18–19 kcal/mol). Overall, the work suggests that both binding poses are viable for methyl transfer, and accurate descriptions of charge transfer and electrostatics are needed to provide balanced relative barriers when multiple binding poses are accessible, for example in other transferases.

simulated trajectories of system 1 and system 9. Both bond distances are shorter when using set 1 Mg 2+ ion parameters compared to set 9 parameters. Results from parameter sets 2-8 were intermediate between these two cases and not shown. Using set 1 parameters ensures no exchange of water molecules over the course of simulation, commensurate with experimental exchange rates for water molecules coordinating Mg 2+ .

II.B. Equilibration methods:
Before any MD simulations were carried out, the prepared proteins were minimized.
During the first stage of minimization, a positional restraint (k=200 kcal/(mol . Å 2 )) was applied to everything except the charge neutralizing ions and the water molecules. This stage consisted of 5,000 steepest descent minimization steps followed by 20,000 conjugate gradient steps. A second minimization stage consisted of releasing the positional restraints on the protein but keeping restraints with the same force constant on the substrates and Mg 2+ for an additional 2000 steepest descent and 20,000 conjugate gradient minimization steps.
In equilibration scheme 1, the system is heated to T = 300 K in the NVT ensemble while positional restraints (k=200 kcal/(mol . Å 2 )) remain on everything but the neutralizing ions and water for 20 ps. Following the quick heating stage, sequential 0.1 ns NPT equilibration runs were carried out in which the positional restraints were reduced stepwise (200, 100, 50, 25, and 20 kcal/(mol . Å 2 )) in each stage until all restraints were removed. Once all positional restraints were removed, a 100 ns NPT (T = 300 K, p = 1 bar) production run is carried out.
In equilibration scheme 2, the same initial 20 ps NVT heating stage with positional restraints is carried out as in scheme 1. Following the quick heating stage, sequential 0.1 ns NPT equilibration runs were carried out in which the positional restraints were reduced stepwise (200, 100, 50, 25, and 20 kcal/(mol . Å 2 )) in each stage until the restraints were removed. At this stage, a 20 ns equilibration run was carried out in which harmonic restraints were only applied to the Mg 2+ -O-and Mg 2+ -OH distances for the substrate-Mg 2+ coordination (d target =2.2 Å, k = 50 kcal/(mol . Å 2 )) as well as the C-C-O-H dihedral (see Figure 1) of the substrate (Ð target =-5.5°, k = 6.5 kcal/(mol . rad 2 )). Dihedral restraints were chosen by trial and error to produce a full-width half maximum of the sampled C-C-O-H dihedral of 24° centered around -5.5°. Once the Mg 2+ -O -/OH distance restraints were removed, the system was equilibrated with a 100 ns NPT (T = 300 K, p = 1 bar) run with only the dihedral restraint. Finally, a 100 ns NPT production run was carried out with no restraints. Comparison of outcomes between these two equilibration schemes is provided in Sec. 3 of the main text.
Below we describe in more detail different simulation protocols that we employed. The simulation parameters (as we describe in the manuscript) are same for all the protocols.

II.B.1. Protocol 1
Step 1: Minimization 1. The system was minimized for 5,000 steps using the steepest descent method followed by 20,000 steps of the conjugate gradient method. During minimization, positional restraints were applied on the protein, SAM, inhibitor, oxygen atom of the crystallographic water, and magnesium ion with a force constant of 200 kcal mol -1 Å -2 .
Step 2: Minimization 2. The system was minimized for 2,000 steps using the steepest descent method followed by 20,000 steps of the conjugate gradient method. During minimization, positional restraints were applied on the SAM, inhibitor, oxygen atom of the crystallographic water, and magnesium ion with a force constant of 200 kcal mol -1 Å -2 .
Step 3: Heating. The system was heated to T=300 K for 20 ps using the NVT ensemble.
Positional restraints were applied on the protein, SAM, inhibitor, one oxygen atom of the crystallographic water, and magnesium ion with a force constant of 200 kcal mol -1 Å -2 .
Step 4: Equilibration. All positional restraints were removed and the system was equilibrated for 100 ns at temperature T=300 K in the NPT ensemble.
Results: In this protocol, the hydroxyl of catecholate reorients to form an intramolecular hydrogen bond with the oxygen anion, increasing the Mg-O(H) distance nearly instantaneously.

II.B.2. Protocol 2
Step 1: Minimization 1. The system was minimized for 5,000 steps using the steepest descent method followed by 20,000 steps of the conjugate gradient method. During minimization, positional restraints were applied on the protein, SAM, inhibitor, oxygen atom of the crystallographic water, and magnesium ion with a force constant of 200 kcal mol -1 Å -2 .
Step 2: Minimization 2. The system was minimized for 2,000 steps using the steepest descent method followed by 20,000 steps of the conjugate gradient method. During minimization, positional restraints were applied on the SAM, inhibitor, oxygen atom of the crystallographic water, and magnesium ion with a force constant of 200 kcal mol -1 Å -2 .
Step 3: Heating. The system was heated to T=300 K for 20 ps in the NVT ensemble. Positional restraints were applied on the protein, SAM, inhibitor, one oxygen atom of the crystallographic water, and magnesium ion with a force constant of 200 kcal mol -1 Å -2 .
Next, we simulated the system with the C-C-O-H angle restraint during the production run in order to confirm that the substrate remained in a bidentate configuration with the Mg 2+ ion. In order to reproduce the dihedral angle distribution where no restraints were used, we ran several simulations with different force constants with the C-C-O-H target angle=-5.5°. As we can see in Table S2, the FWHM is 24° for both the systems with the force constant 5 and 6.5 kcal . mol -1. rad -2 . We also considered representative bond distances in the active site. Short nonbonded C-O distances are best preserved in the simulation with a 6.5 kcal . mol -1. rad -2 force constant.

II.D. Umbrella sampling restraint details
The Mg 2+ -OH distance was sampled from 2.0 to 4.5 A in 0.1 Å increments, and the C-C-O-H dihedral (see Figure S3 inset) was sampled from 0° to 180° in 5° increments. Initial configurations were generated for each window from a multistep process. First, an initial configuration was selected from the long direct MD production run. Next, a series of seven sequential restrained minimizations and short dynamics runs were carried out. The snapshot from the production run MD was used as a starting point to generate the first configuration (Ð CCOH =0°, d Mg-OH =2.2 Å) through a multi-step process. In this equilibration, a 2000 step steepest descent minimization with positional restraints on the protein and substrates (k=200 kcal/(mol . Å 2 )) was carried out followed by NPT (T = 300 K, p = 1 bar) equilibration for 300 ps using a 200 kcal/(mol . rad 2 ) force constant on the dihedral and 500 kcal/(mol . Å 2 ) on the bond distance. The resulting structure was then visually inspected to verify that the target distance and dihedral were satisfied without distorting the surrounding protein and substrate environment. The outcome of this simulation was also used as the starting point for the next coarse minimization and equilibration window (Ð CCOH =30°, d Mg-OH =2.2 Å). The minimization, equilibration, and inspection procedure was repeated sequentially for the other five windows (Ð CCOH =60°, d Mg-OH =2.2 Å; Ð CCOH =90°, 120°, 150°, 170° at d Mg-OH =3.5 Å).
Unbiased free energies were obtained from the sampled distributions using the weighted histogram analysis method (WHAM) 3-4 using the Grossfield lab WHAM software package 5 .
The Amber code does not multiply force constants for distances by ½ while the WHAM package does, so force constants were rescaled accordingly for analysis. Distributions were visually inspected to confirm that distributions were overlapping between windows, and the force constants for dihedrals and distances were obtained by tuning the balance between overlapping distributions and sufficiently sampling the target distance and dihedral.