The Force Field Toolkit ( fg TK) Christopher G. Mayne, Emad - PowerPoint PPT Presentation

Parametrizing Small Molecules Using: The Force Field Toolkit ( fg TK) Christopher G. Mayne, Emad Tajkhorshid Beckman Institute for Advanced Science and Technology University of Illinois, Urbana-Champaign Klaus Schulten University of Illinois, Urbana-Champaign James C. (JC) Gumbart, Anna Pavlova Georgia Institute of Technology Computational Biophysics Workshop | DICP | July 11 2018

（） （） MD Simulations of Biological Systems Molecular Mechanics Force Fields U = U bonds + U angles + U dihedrals + U vdW + U coulombic bonded non-bonded The CHARMM Force Field Σ k i Σ k i Σ Σ U = U = bond bond ( r i - r 0 ) 2 ( r i - r 0 ) 2 angles k i angles k i angle angle ( θ i - θ 0 ) 2 ( θ i - θ 0 ) 2 + + + + bonds bonds Σ Σ [ 1 + cos( n i ɸ i + δ i )] [ 1 + cos( n i ɸ i + δ i )] + + dihedrals k i dihedrals k i dihedral dihedral i Σ i Σ σ ij σ ij σ ij σ ij i Σ i Σ Σ Σ Σ Σ 12 12 q i q j q i q j 6 6 + + - - j ≠ i 4 ∈ ij （ j ≠ i 4 ∈ ij （ r ij ） r ij ） r ij r ij r ij r ij j ≠ i j ≠ i

Parameter Transferability In Biopolymers Parameter set describes molecular behavior in varied 
 chemical (connectivity) and spatial (conformation) contexts Peptides and Proteins Nucleic Acids H 2 N O NH 2 O S Me O NH N Me N N N HN OH O NH 2 O O O N N H H H N N O N O N N N N N N N H H H O O O O Me O O O O O O Me Me O P O O P O O P O O O NH 2 O O O limited set of isolated 
 R Key Features: building blocks repetitive backbone unit

Parametrization as an Impasse non-standard or 
 small molecule ligands engineered amino acids N HO Me Me HN N N O N Me O S HN N COOH H 2 N COOH O OH H Me O N O S N Tiotropium (Spiriva) Imatinib (Gleevec) cofactors metal centers NH 2 H 2 N NH 2 N O O O O Me Me N S S HS O P O P O N N O N N O O H H O O Fe OH S S O OH O O O P O Coenzyme A NH 2 H 2 N O

CGenFF Parametrization Workflow General Parametrization Workflow Find Missing Parameters build init PAR PSF/PDB System Preparation update PDB Geometry Optimization (QM) Water Interaction En. (QM) Charges update PSF Charge Optimization Hessian Calculation (QM) Bonds & Angles Bond & Angle Optimization update PAR Torsion Scan (QM) PAR File Dihedrals / Torsions Torsion Optimization update PAR Calculation Action K. Vanommeslaeghe et al. , J. Comput. Chem. 2010 , 31, 671-690. C.G. Mayne et al. , J. Comput. Chem. 2013 , 34, 2757-2270.

fg TK Interface standard file dialogs tasks organized under tabs action buttons action menus

Functionality Provided by fg TK Setup & Perform 
 Multi-dimensional Optimizations Core Functions Abstraction of Gaussian I/O (QM) Assess Performance of Parameters by Visualizing Optimization Data Support Functions • Auto-detect Water Interaction Sites • Visualize Target Data in VMD • Auto-detect Charge Groups • Create Graphic Objects in VMD • Label Atoms in VMD • Auto-detect Non-redundant Torsions • Build & Update Parameter Files • Read Input Parameters from File • Browse Existing Parameter Sets • Read/Write Data From Opt. Logs • Export Plot Data to File • Write Updated Charges to PSF • Reset Opt. Input from Output • Monitor Optimization Progress

fg TK Exemplified by Charge Optimization build init PAR Find Missing Parameters PSF/PDB update PDB Geometry Optimization (QM) Water Interaction En. (QM) update PSF Charge Optimization Hessian Calculation (QM) update PAR Bond & Angle Optimization Torsion Scan (QM) PAR File update PAR Torsion Optimization Calculation Action

Generating Charge Optimization Target Data Load QM optimized geometry | Auto-detect interaction sites | Generat pyrrolidine VMD main window fg TK GUI

Generating Charge Optimization Target Data mized geometry | Auto-detect interaction sites | Generate Gaussian Input Files | Run Q Donor Acceptor VMD main window fg TK GUI

Generating Charge Optimization Target Data raction sites | Generate Gaussian Input Files | Run QM | Inspect water optimization Compute water position Optimize 
 distance & rotation fg TK GUI

Generating Charge Optimization Target Data iles | Run QM | Inspect water optimization Visually assess 
 QM-optimized 
 water position(s) fg TK GUI

Charge Optimization Setup Optimization Objective Function Load QM Target Data Prepare Optimization Σ 𝑔 (U MM -U QM ) Optimizer: wat. int. Assign Charges + Compute U MM , d MM , µ MM Σ Compute Objective Function 𝑔 (d MM -d QM ) wat. int. Return Optimized Charges + 𝑔 ( μ MM - μ QM ) Analyze Performance Write Charges to PSF

Assessing MM Water-Interaction Profiles 10 Initial Charges Intermediate Charges 8 Final Optimized Charges Literature Charges � U MM-QM (kcal/mol) 6 4 2 0 -2 -0.4 -0.2 0 0.2 0.4 � d MM-QM (Å)

Sampling MM Water-Interaction Profiles Mode: Simulated Annealing 6 4500 H H H H 5 H N N N N 4000 4 H Δ U MM-QM (kcal/mol) 3 3500 H H 2 3000 Iteration 1 2500 0 2000 -1 -2 1500 -3 1000 -4 500 -5 -6 0 -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 Δ d MM-QM ( Å )

Tuning the Optimization Objective Function 2 U MMmin -U QMmin Σ + w i U tol w i wat. int. 2 d MMmin -d QMmin Σ w d w i + d tol wat. int. w d w μ 2 2 μ MM - μ QM 𝛴 + nw μ 𝛴 tol μ tol In practice, it is impossible to fit all of these perfectly! Often we decrease w d and w µ to improve the fit to the energies

Plotting Charge Optimization Data 8 export 7 data 6 Δ E from Target (kcal/mol) 5 4 3 2 1 0 -1 -2 0 50 100 150 200 250 300 Optimization Iteration

Restrained Electrostatic Potential (RESP) fitting An alternative to water interactions for charges, commonly used in Amber The QM electrostatic potential is calculated and then fit by optimizing the MM charges Has problems with buried atoms, which may not noticeably affect the ESP https://studynights.blogspot.com/2015/03/the- single-point-energy-of-mnh2o6-and.html RESP fitting is supported by FFTK, requires downloading the resp program as part of AmberTools (free)

Fitting of Bonds and Angles Bond and angle are fit by creating a small distortion of the bond/angle and calculating the QM energy and the MM energy, then choosing the force constants to match Σ k i Σ bond ( r i - r 0 ) 2 angles k i angle ( θ i - θ 0 ) 2 + bonds

Fitting of Dihedrals Dihedrals are scanned in QM in 10-15 deg. increments Energies of each conformation are fit in MM according to the pre-determined dihedral terms included Σ [ 1 + cos( n i ɸ i + δ i )] dihedrals k i dihedral periodicity (1-6 possible) phase (always 0 or 180 deg.) Energies above a threshold (e.g., 8-10 phase (always kcal/mol) are ignored 0 or 180 deg.) periodicity

Fitting of Dihedrals QM PES: black After 7 rounds of initial: red simulated annealing, First fit: blue the fit is much better

Two Approaches to Fitting the Dihedrals Several multiplicities One multiplicity and and free phase locked phase Pros : limits incorrect behavior, Pro: very good fit of QM PES sets of force constants Cons: possible incorrect behavior Cons : fit to QM PES not always multiple sets of force constants possible In practice a tradeoff is needed!

Example of Overfitting imidazole-pyrridine moiety of antibiotic telithromycin 3 multiplicities, free phases for each dihedrals k 1 [1 + cos( φ + δ 1 )] + k 2 [1 + cos(2 φ + δ 2 )] + k 3 [1 + cos(3 φ + δ 3 )] Dihedral Fit QM PES MM Fit 10 9 8 7 E (kcal/mol) multiple dihedrals 6 scanned; plotted 5 simultaneously to 4 reference a common 3 global minimum 2 1 0 0 100 200 50 150 Configuration Pavlova, Gumbart. Parametrization of macrolide antibiotics using the force field toolkit. J. Comp. Chem . 2015, 36 , 2052 − 2063.

Example of Overfitting imidazole-pyrridine moiety of antibiotic telithromycin Using too many dihedral multiplicities can leads to distortion of a planar molecule! Pavlova, Gumbart. Parametrization of macrolide antibiotics using the force field toolkit. J. Comp. Chem . 2015, 36 , 2052 − 2063.

Example of Overfitting imidazole-pyrridine moiety of antibiotic telithromycin planar dihedrals have multiplicity 2 and phase 180 deg. k 2 [1 + cos(2 φ + π )] Dihedral Fitting QM PES MM Fit 10 9 8 The fit looks (surprisingly) 7 E (kcal/mol) better despite using fewer 6 terms (why?) 5 4 3 Fitting a lot of parameters 2 simultaneously cannot 1 always find the best fit! 0 0 50 100 150 200 Configuration Pavlova, Gumbart. Parametrization of macrolide antibiotics using the force field toolkit. J. Comp. Chem . 2015, 36 , 2052 − 2063.

Example of Overfitting imidazole-pyrridine moiety of antibiotic telithromycin planar dihedrals have multiplicity 2 and phase 180 deg. k 2 [1 + cos(2 φ + π )] Planarity is maintained! Problems persist! Eclipsed conformation of the alkane Restraining phase of CH dihedrals to 0 prevents eclipsed conformations Pavlova, Gumbart. Parametrization of macrolide antibiotics using the force field toolkit. J. Comp. Chem . 2015, 36 , 2052 − 2063.

How to know what terms to include? https://mackerell.umaryland.edu/~kenno/cgenff/faq.php