A coupled hybrid potential QMMMMD simulations in Amber

A coupled (hybrid) potential QM/MM/MD simulations in Amber Dr. Vladislav Vasilyev Supercomputer Facility The Australian National University

Hybrids are popular from the ancient time

Why Do We Need a Hybrid QM/MM Approach? Quantum chemical methods are generally applicable and allow the calculation of ground and excited state properties (molecular energies and structures, energies and structures of transition states, atomic charges, reaction pathways etc. ) Molecular Mechanical methods are restricted to the classes of molecule it have been designed for and their success strongly depends on the careful calibration of a large number of parameters.

Why Do We Need a Hybrid QM/MM Approach? The main bottleneck of quantum chemical methods is that they are CPU and memory hungry. For example, for small peptide of 126 atoms one energy evaluation requires: CPU Time Method Memory Seconds Time units KB Memory units Quantum chemical* 273. 00 1820 4889 85 Molecular Mechanical 0. 15 1 58 1 *Semi-empirical PM 3 method In general, CPU and memory requirements: Molecular Mechanical methods ~ N 2 Semiempirical Quantum Chemical methods ~ N 2 Ab initio Quantum Chemical methods ~ N 4

A Hybrid QM/MM Approach The development of hybrid QM/MM approaches is guided by the general idea that large chemical systems may be partitioned into an electronically important region which requires a quantum chemical treatment and a remainder which only acts in a perturbative fashion and thus admits a classical description.

The Simplest Hybrid QM/MM Model Hamiltonian for the molecular system in the Born-Oppenheimer approximation: “Standard” QM hamiltonian The main drawbacks of this simple QM/MM model are: it is impossible to optimize the position of the QM part relative to the external charges because QM nuclei will collapse on the negatively charged external charges. some MM atoms possess no charge and so would be invisible to the QM atoms the van der Waals terms on the MM atoms often provide the only difference in the interactions of one atom type versus another, i. e. chloride and bromide ions both have unit negative charge and only differ in their van der Waals terms.

A Hybrid QM/MM Model So, it is quite reasonable to attribute the van der Waals parameters (as it is in the MM method) to every QM atom and the Hamiltonian describing the interaction between the QM and MM atoms can have a form: The van der Waals term models also electronic repulsion and dispersion interactions, which do not exist between QM and MM atoms because MM atoms possess no explicit electrons. A. Warshel, M. Levitt // Theoretical Studies of Enzymic Reactions: Dielectric, Electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme. // J. Mol. Biol. 103(1976), 227 -49

The Hybrid QM/MM Model Now we can construct a “real” hybrid QM/MM Hamiltonian: A “standard” MM force field can be used to determine the MM energy. For example, AMBER-like force field has a form:

Choice of QM method. . . is a compromise between computational efficiency and practicality and the desired chemical accuracy. The main advantage of semi-empirical QM methods is that their computational efficiency is orders of magnitude greater than either the density functional or ab initio methods

QM Methods in Amber 9 • Available semi- empirical Hamiltonians are PM 3, AM 1, MNDO, PDDG/PM 3, PDDG/MNDO. • They can be used for gas phase, generalized Born and PME periodic simulations.

QM Methods in Amber 9 • Support is also available, on a functionally limited basis: 1. The Density Functional Theory-based-tightbinding (DFTB) Hamiltonian 2. The Self-Consistent-Charge version, SCCDFTB The DFTB/SCC-DFTB implementation does not currently support generalized Born, PME or Ewald calculations,

The elements supported by QM methods in Amber 9 • • MNDO: H, Li, Be, B, C, N, O, F, Al, Si, P, S, Cl, Zn, Ge, Br, Sn, I, Hg, Pb AM 1: H, C, N, O, F, Al, Si, P, S, Cl, Zn, Ge, Br, I, Hg PM 3: H, Be, C, N, O, F, Mg, Al, Si, P, S, Cl, Zn, Ga, Ge, As, Se, Br, Cd, In, Sb, Te, I, Hg, Tl, Pb, Bi PDDG/PM 3: H, C, N, O, F, Si, P, S, Cl, Br, I PDDG/MNDO: H, C, N, O, F, Cl, Br, I PM 3 CARB 1: H, C, O DFTB/SCC-DFTB: H, C, N, O, S, Zn

Calibration of the QM/MM potential • Crucial aspect is how the interaction between QM and MM parts is determined. • In choosing the appropriate form, it is required that the balance between attractive and repulsive forces must be preserved and the QM/MM interactions must be of the correct magnitude with respect to the separate QM and MM contributions

Calibration of the QM/MM potential: Parameterizations (1) (2) • 1) Modification of the one-electron terms arising from interaction of the electron cloud of the QM fragment with the point charge of an MM atom. • 2) By varying the radii in the van der Waals terms. • 3) By varying 1)+2)

Calibration of the QM/MM potential • 1) By hand, to find the optimum values of the parameters by calculating interaction curves for charge/ion systems and comparing them with the MP 2/6 -311++G** ab initio results. M. J. Field, P. A. Bash, M. Karplus, J. Comp. Chem. , 11(1990), 700 -733. • 2) Fitting calculated H-bond energies to experimental data on ion-molecular complexes in the gas phase. V. V. Vasilyev, A. A. Bliznyuk, A. A. Voityuk, Int. J. Quant. Chem. 44(1992), 897930.

Calibration of the QM/MM potential: • 3) Optimizing van der Waals parameters on QM atoms to reproduce the 6 -31 G(d) interaction energies for H-bonded complexes in the gas phase. P. A. Bash, L. Lawrence, A. D. Mac. Kerell, Jr. , D. Levine, P. Hallstrom, PNAS USA, 93(1996), 3698 -703. • 4) Optimizing van der Waals parameters on QM atoms to reproduce the MP 2/6 -31 G(dp) interaction energies for H-bonded complexes in the gas phase. J. Gao // Toward a molecular Orbital Derived Empirical Potential for Liquid Simulations // J. Phys. Chem. B 101(1997), 657 -63 • 5) By varying the radii in the van der Waals terms to reproduce experimental free energies of solvation using MD simulations. P. L. Cummins, J. E. Gready, J. Comp. Chem. , 18(1997), 1496 -512.

Dividing Covalent Bonds across the QM and MM Regions In many simulations it is necessary to have the QM/MM boundary cut covalent bonds, and a number of additional approximations have to be made.

Dividing Covalent Bonds across the QM and MM Regions Using a hybrid orbital on the frontier MM atom A. Warshel, M. Levitt // Theoretical Studies of Enzymic Reactions: Dielectric, Electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme. // J. Mol. Biol. 103 (1976), 227 -249 V. Thery, D. Rinaldi, J. -L. Rivail, B. Maigret, G. G. Ferenczy, J. Comp. Chem. 15 (1995), 269

Dividing Covalent Bonds across the QM and MM Regions Using “link” atoms “Link” atoms are used to gracefully cap the electron density. This approach is used in Amber

Implementation of “link” atom Approach in Amber 9 The link atom is placed along the bond vector joining the QM and MM atom The default link atom type is hydrogen It interacts with MM region only electrostatically (no VDW term). Wd. V interaction between QM and MM atoms which form 1 -2 and 1 -3 “bonded” pairs is not calculated. Bond stretching, angle bending, and torsion interactions between QM and MM regions are calculated as those in MM if 1 -2, 1 -2 -3, or 1 -2 -3 -4 terms contain at least one MM atom

Example of Application of the QM/MM Method to the Enzyme Catalysis

Tetrahedral Intermediate Formation in the Acylation Step of Acetylcholinesterases • Acetylcholinesterases (ACh. E) are the serine protease enzymes which hydrolyze the neurotransmitter acetylcholine (Ach) CH 3 COOCH 2 N+(CH 3)3 + Ach. E CH 3 CO-Ach. E + HOCH 2 N+(CH 3)3 CH 3 COO- + H+ + ACh. E and which function at a rate approaching that of a diffusioncontrolled reaction. • Remark: Human Cathepsin G is also a serine protease enzyme V. V. Vasilyev, J. Mol. Struct. (Theochem), 304(1994), 129.

Steps along the reaction pathway of serine protease catalyzed bond cleavage “Catalytic triad” Rate limiting step of the reaction

Partition into the QM and MM Parts Þ QM Region (PM 3) – CH 3 OH (as Ser-200), methylimidazole (as His-440), CH 3 CH 2 COOH (as Glu -327), and CH 3 COOCH 3 (as a substrate) MM Region – 5161 enzyme atoms

Activation of the Serine Residue Critical Step of the Reaction is activation of the Serine residue In Gas Phase: PA(CH 3 O-) = ~ -350 kcal/mole PA(Imidazole) = ~ -220 kcal/mole where PA is a Proton Affinity

Activation of the Serine Residue A proton transfer from the Serine to Histidine residue is very unfavorable in the gas phase (~34. 6 kcal/mol). A proton transfer is less unfavorable in the enzyme (energy barrier is about 21. 5 kcal/mol).

Activation of the Serine Residue Enzyme environment creates electrostatic potential which favours the proton transfer from the Serine to Histidine residue.

Tetrahedral Intermediate Formation in the Acylation Step of Acetylcholinesterases Þ Incorporation of enzyme environment in simulation changes drastically the flow of the reaction Þ A decrease in the activation energy of TI formation in the enzyme environments versus the uncatalysed gas-phase reaction is about 27 kcal/mol (experimental estimation of the reduction in activation energy for the whole enzyme reaction versus uncatalysed neutral, basic, and acidic hydrolysis of Ach. E is 14, 11, and 18 kcal/mol, respectively. )

Hints for running QM/MM calculations Choosing the QM region • There are no good universal rules here • One might want to have as large a QM region as possible • However, having more than 80 -100 atoms in the QM region will lead to simulations that are very expensive.

Hints for running QM/MM calculations Choosing the QM region • for many features of conformational analysis, a good MM force field may be better than a semi-empirical or DFTB quantum description.

Hints for running QM/MM calculations Choosing the QM region

Hints for running QM/MM calculations Parallel Simulations • At present all parts of the QM simulation are parallel except the density matrix build and the matrix diagonalisation. • For small QM systems these two operations do not take a large percentage of time and so acceptable scaling can be seen to around 8 cpus. • However, for large QM systems the matrix diagonalization time will dominate and so the scaling will not be as good.

Amber 9 QM/MM Example

Resume • Amber 9 features new and significantly improved QM/MM support • The QM/MM facility supports gas phase, implicit solvent (GB) and periodic boundary (PME) simulations • Compared to earlier versions, the QM/MM implementation offers improved accuracy, energy conservation, and performance.