This item is only available as the following downloads:
1 C OMPUTATIONAL MODEL L ING AND SIMULATION OF PRENYLTRANSFERASES AND CORRESPONDING PRENYLATION By YUE YANG A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENT S FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2011
2 2011 Yue Yang
3 To Amelia
4 ACKNOWLEDGMENTS I would like to express my gratitude to everyone that has been in my life in the past five years. To pursue a doctoral degree is never an easy task. All of your help and support have made this road much easier to walk through. It is five years of hard work ing and a ll of your assistance that make this dissertation here today. Please accept my sincerely grateful A special thank you is due to my family especially my dearest daughter, Amelia for all the happiness, hope, faith and courage they gave me. Another special thank you should be given to my advisor Kennie thank you for giving me this precious opportunity to work with you a nd offering me help throughout my entire graduate career it has been an enjoyable and memorable five years that I will always remember. Thank you Adrian, Erik, Nicole and Joanna for joining my supervisory committee and giving m e countless help advices. An d t hank you also to Yuning, Yilin and all my friends in Merz group, QTP, Department of Chemistry and University of Florida.
5 TABLE OF CONTENTS page ACKNOWLEDGMENTS ................................ ................................ ................................ ........... 4 LIST OF TABLES ................................ ................................ ................................ ...................... 8 LIST OF FIGURES ................................ ................................ ................................ .................... 9 LIST OF ABBREVIATIONS ................................ ................................ ................................ .... 12 ABSTRACT ................................ ................................ ................................ ............................. 15 CHAPTER 1 INTRODUCTION ................................ ................................ ................................ ............. 17 1.1 Prenylation and Prenyltransferase ................................ ................................ ................. 17 1.2 Protein Prenyltransferase ................................ ................................ .............................. 18 1.3 Aromatic Prenyltransferase ................................ ................................ ........................... 22 1.3.1 Bacterial Aromatic Prenyltransfera se NphB ................................ ........................ 22 1.3.2 Fungal Indole Prenyltransferase FtmPT1 ................................ ............................ 25 1.4 Computational Chemistry/Biochemistry ................................ ................................ ....... 27 2 THEORY AND METHODS ................................ ................................ .............................. 29 2.1 Quantum Mechanics and the Born Oppenheimer Approximation ................................ .. 29 2 .2 Molecular Mechanics ................................ ................................ ................................ ... 31 2.2.1 Bond Stretching ................................ ................................ ................................ .. 33 2.2.2 Angle Bending ................................ ................................ ................................ .... 34 2.2.3 Torsion and Improper Torsion ................................ ................................ ............ 35 2.2.4 Van der Waals Interaction ................................ ................................ .................. 36 2.2.5 Electrostatic Interaction ................................ ................................ ...................... 38 2.2.6 Cross Terms ................................ ................................ ................................ ....... 39 2.2.7 Force Field Parameterization ................................ ................................ .............. 41 2.3 Energy Mini mization and Molecular Simulation ................................ ........................... 41 2.3.1 Energy Minimization ................................ ................................ .......................... 42 126.96.36.199 Steepest descent method ................................ ................................ ........... 43 188.8.131.52 Conjugated gradient method ................................ ................................ ..... 43 184.108.40.206 Newton Raphson method ................................ ................................ .......... 43 2.3.2 Molecular/Comp utational Simulation ................................ ................................ 44 220.127.116.11 Phase space ................................ ................................ ............................... 45 18.104.22.168 Monte Carlo method ................................ ................................ ................. 45 22.214.171.124 A brief introduction to molecular dynamics ................................ .............. 46 2.4 Molecular Dynamics ................................ ................................ ................................ .... 46 2.4.1 Finite Difference Method: Ver let and Leap frog Algorithm ................................ 47 2.4.2 Time Step and SHAKE ................................ ................................ ....................... 48
6 2.4.3 Molecular Dynamics at Constant Temperature ................................ .................... 50 2.4.4 Pressure Regulation in Molecular Dynamics ................................ ....................... 52 2.4.5 Solvent Effects and Water Models ................................ ................................ ...... 53 126.96.36.199 Implicit solvent models ................................ ................................ ............. 54 188.8.131.52 Explicit water models ................................ ................................ ............... 55 184.108.40.206 Periodic boundary conditions and partic le mesh Ewald ............................. 57 2.5 Free Energy Calculations in Molecular Dynamics ................................ ........................ 59 2.5.1 Steered Molecular Dynamics ................................ ................................ .............. 60 2.5.2 Umbrella Sampling ................................ ................................ ............................. 61 2.6 Simulating Chemical Reactions and QM/MM ................................ .............................. 63 2.6.1 QM/MM ................................ ................................ ................................ ............. 63 2.6.2 Semiempirical Theory in QM/MM ................................ ................................ ..... 65 2.6.3 Semiempirical DFT SCC DFTB ................................ ................................ ...... 66 2.6.4 Transition State Theory ................................ ................................ ...................... 68 2.6.5 Kinetic Isotope Effect ................................ ................................ ......................... 71 3 APPLY MOLECULAR DYNAMICS SIMULATION TO PREDI CT THE BINDING OF MAGNESIUM CATION IN FANESYLTRANSFERASE ACTIVE SITE ................... 74 3.1 Background ................................ ................................ ................................ .................. 74 3.2 Methods and Theories ................................ ................................ ................................ .. 78 3.2.1 Building Initial Model ................................ ................................ ........................ 78 3.2.2 Modelling and Simulating Mg 2+ in FTase ................................ ........................... 79 ................................ ................................ .............. 82 3.2.4 Modeling Removal of Mg 2+ from Active Site ................................ ...................... 82 3.2.5 PMF Validation ................................ ................................ ................................ .. 83 3.2.6 Modelling and Simulation of GGTase I ................................ .............................. 83 3.3 Results and Discussion ................................ ................................ ................................ 84 3.3.1 Modeling and Simulation of FTase/CVIM/FPP 3 /Mg 2+ Complexes ..................... 84 2+ /H 2 O Swap Simulation ................................ .................. 91 3.3.3 PMF Studies ................................ ................................ ................................ ....... 92 3.3.4 Modelling and Simulation of GGTase I ................................ .............................. 96 3.4 Summary ................................ ................................ ................................ ...................... 98 4 SYSTEMATICALLY MECHANISTIC STUDY OF FTASE COMPLEXED WITH DIFFERENT PEPTIDE SUBSTRATES ................................ ................................ .......... 101 4.1 Background ................................ ................................ ................................ ................ 101 4.2 Methods and Simulation Details ................................ ................................ ................. 103 4.2.1 FTase/CVIM/FPPH 2 ................................ ........................ 103 4.2.2 FTase/CVLM/FPPH 2 ................................ ................................ ....................... 105 4.2.3 FTase/CVLS/FPPH 2 ................................ ................................ ........................ 106 4.2.4 QM/MM Validation of FTase/CVIM/FPP 3 /Mg 2+ Model ................................ .. 107 4.3 Results and Discussions ................................ ................................ .............................. 108 4.3.1 A Glimpse of Enti re Product Formation Process of FTase/CVIM/FPPH 2 ......... 108 220.127.116.11 Transition state and free energy barrier ................................ ................... 108
7 18.104.22.168 Comparison of QM/MM and MM PMF of conformational transition step ................................ ................................ ................................ ................. 111 4.3.2 Effect of Mutations ................................ ................................ ........................... 112 ................................ ................................ .... 113 22.214.171.124 Interesting context dependent effect ................................ ........................ 116 SKIE results of FTase/CVLS complex ....................... 117 4.3.3 Mg 2+ in Chemical Step ................................ ................................ ..................... 122 4.4 Summary ................................ ................................ ................................ .................... 125 5 COMPUTATIONAL APPROACH TO UNDE RSTAND AROMATIC PRENYLTRANSFERASE CATALYTIC MECHANISMS ................................ ............. 127 5.1 Background ................................ ................................ ................................ ................ 127 5.2 Methods and Simulation Details ................................ ................................ ................. 130 5.2.1 QM/MM Study of NphB Prenylation ................................ ................................ 130 5.2.2 Modelling and Simulation of Proton Transfer Step in NphB Catalysis .............. 132 5.2.3 Classical and QM/MM study of FtmPT1 ................................ .......................... 134 5.3 Results and Discussions ................................ ................................ .............................. 136 5.3.1 Orientation of Substrate Orientation in QM/MM Equilibration of NphB ........... 136 5.3.2 Free Energy Surface of NphB Catalyzed Prenylation ................................ ........ 138 126.96.36.199 S1 state ................................ ................................ ................................ ... 138 188.8.131.52 S3 state and S4 state ................................ ................................ ................ 142 5.3.3 Proton Transfer Step of NphB Catalysis ................................ ............................ 145 5.3.4 Substrate Protonation State Preference of FtmPT1 ................................ ............ 148 5.3.5 A Look at FtmPT1 Catalysis ................................ ................................ ............. 150 5.4 Summary ................................ ................................ ................................ .................... 152 LIST OF REFERENCES ................................ ................................ ................................ ........ 155 BIOGRAPHICAL SKETCH ................................ ................................ ................................ ... 168
8 LIST OF TABLES Tab le page 2 1 Parameters of several water models. ................................ ................................ .............. 57 3 1 Zinc parameters used in MD simulations ................................ ................................ ....... 81 4 1 Free energy parameters obtained from QM/MM and MM PMF studies. ....................... 114
9 LIST OF FIGURES Figure page 1 1 Chemic al structure of isoprene and prenyl group. ................................ .......................... 17 1 2 Protein farnesyltransferase (FTase) catalyzed farnesylation. ................................ .......... 19 1 3 Protein gera nylgeranyltransferase (GGTase) catalyzed geranylgeranylation. .................. 19 1 4 NphB (middle up) catalyzed geranylation. ................................ ................................ ..... 22 1 5 Indole prenyl transferase FtmPT1 (middle up) catalyzed attachment of dimethylallylic group to brevianamide F (left bottom). ................................ .................. 25 2 1 Schematic representations of five components in general MM force field: ..................... 32 2 2 Torsion (or dihedral) angles in protein structure. ................................ ............................ 36 2 3 Schematic representation of some multipole examples. ................................ .................. 39 2 4 Schematic representation of angle closing coupling with bonds stretching. .................... 40 2 5 Simple interaction site water models. ................................ ................................ ............. 56 2 6 Part of 2 D periodic boundary conditions (PBC). ................................ ........................... 58 2 7 Schematic representation of Ewald summation. ................................ ............................. 59 2 8 Schematic representation of building PMF with umbrella samplings. ............................. 62 2 9 The effect of k w on umbrella sampling (US). ................................ ................................ 62 2 10 Schemetic representation of QM and MM partition. ................................ ....................... 65 2 11 Reaction coordinate diagram for the biomolecular nucleophilic substitution (S N 2) reaction between bromomethane and the hydroxide anion. ................................ ............. 69 2 12 Reaction mechanism of S N 1 reaction between S 3 chloro 3 methylhexane and iodide ion that yields both R and S 3 iodo 3 methylhexane ................................ .................... 73 2 13 Reaction mechanism of S N 2 reaction between bromomethane and hydroxide that yields ethanol and bromide. ................................ ................................ ........................... 73 3 1 Snapshot of FPP binding pocket. Two reacting centers are labeled as C1 (C 1 of FPP) and SG (S of Cys1p). ................................ ................................ ................................ ... 75 3 2 Active site snapshot of GGTase I. C 1 of GGPP and S of Cys1p are labeled as C1 and SG. ................................ ................................ ................................ .......................... 77
10 3 3 Snapshots of FTase active site for different systems ................................ ....................... 79 3 4 Distribution of d RC in the original 6 ns MD simulation of WT2 model. .......................... 90 3 5 Snapshot of active site from the extend 10 ns MD simulation of WT2 model. ................ 90 3 6 Free energy profiles of the modeled systems. ................................ ................................ 94 3 7 Comparison of FPP 3 orientation found in WT2 MD simulation and from previous FTase/CVIM/FPPH 2 (without Mg 2+ ) MD simulation. ................................ ................... 95 3 8 Distances variation in the simulation of GGTase/CVI L/GGPP 3 (top) and GGTase/CVIL/GGPPH 2 (bottom). ................................ ................................ ................ 97 4 1 2+ and diphosphate and the first isoprene unit of FPPH. ................................ .......................... 104 4 2 Free energy profile of FPPH model obtained from a set of 63 windows of umbrella sampling. Active site snapshots are given at key points. ................................ ............... 110 4 3 Snapshot of the active site of FPPH at the transition state. O 1 and C 1 of FPPH, and S of Cys1p are label ed in red. ................................ ................................ ......................... 111 4 4 Free energy profiles of all five models included in this study. ................................ ...... 113 4 5 Snapshot of the active site of Y300F mu tant at the transition state. O 1 and C 1 of FPPH, and S of Cys1p are labeled in red. ................................ ................................ .... 115 4 6 Snapshot of the active site of CVLM mutant at the transition state. O 1 and C 1 of FPPH, and S of Cys1p are labeled in red. ................................ ................................ .... 117 4 7 Free energy profile associated with the conformational transition step of CVLS model. ................................ ................................ ................................ ......................... 118 4 8 Active site snapshot of CVLS model at the transition state. O 1 and C 1 of FPPH, and S of Cys1p are labeled in red. ................................ ................................ ..................... 119 4 9 H 1 H 2 C 1 C 2 dihedral variation in the farnesylation of FTase/FPPH/CVIM (red) and FT ase/FPPH/CVLS (green) complexes ................................ ................................ ........ 120 4 10 Ramachandran plot of a 2 residue of the Ca 1 a 2 X motif throughout the farnesylation of different systems ................................ ................................ ................................ ......... 121 4 11 The active site snapshot of WT2 at the transitio n state. C 1 and O 1 of FPP, and S of Cys1p are labeled in red. ................................ ................................ .............................. 124 5 1 QM partition of NphB: GPP and 1,6 DHN. Important atom specifications are given.. .. 131
11 5 2 Ac tive site configuration of the initial structure selected for proton transfer step simulation. ................................ ................................ ................................ ................... 1 33 5 3 Active site configuration of PDB ID 3O2K. DMAPP and Brevianamide F are treated with QM in QM/MM simulation. Important atoms are labeled in red. .......................... 135 5 4 1,6 DHN orientation in S1 (left) and S3 (right) State. ................................ ................... 138 5 5 Free energy profile associated with prenylation step for both S1 and S3 model.. .......... 139 5 6 Active site configuration at the intermediate state of S1 model. O 1P of PPi, C 1P of carbocation and C 5A of 1,6 DHN are labeled in red. ................................ ..................... 140 5 7 H 1(P) H 2(P) C 1(P) C 2(P) dihedral variation in the chemical reaction of NphB catalyzed 5 prenylation (red) and FTase/FPPH/CVIM catalyzed farnesylation (blue). ................. 141 5 8 Active site conformation in S3 model at the intermediate state. O 1P of PPi, C 1P of carbocation, and C 2A of 1,6 DHN are labeled in red. ................................ .................... 144 5 9 Free energy profile associated with S4 model prenylation. ................................ ........... 145 5 10 Active site snapshot of S4 prenylation. O 1P of PPi, C 1P of carbocation and C 4A of 1,6 DHN are red labeled. ................................ ................................ ................................ ... 146 5 11 Free energy profile associated with the S1 proton transfer step. ................................ ... 147 5 12 Active site configuration after proton abstract ion from 5 prenyl 1,6 DHN to H3O (WAT in Figure 5 2). ................................ ................................ ................................ ... 148 5 13 Carbocation formed in FtmPT1 catalysis. Important atoms are labeled in red. .............. 151
12 LIST OF ABBREVI ATIONS 2(/n) D 2(/n) dimensional AA Amino Acid APTase Aromatic prenyltransferase C terminus Carboxyl terminus or COOH terminus DFT Density functional theory DHN Dihydroxynaphthalene DMAPP Dimethylallyl diphosphate DMAPP Dimethylallyl diphosphate FES Free e nergy surface FPP Farnesyl diphosphate FTase Protein prenyltransferase G proteins Guanine nucleotide binding proteins GAFF Generic AMBER force field GB Generalized Born GGA Generalized gradient approximation GGPP Geranylgeranyl diphosphate GGTase Protein g eranylgeranyltransferase GPP Geranyl diphosphate GPU Graphics Processing Unit GTPase Guanine triphosphate hydrolase HF Hartree Fock HPPi Monoprotonated diphosphate leaving group (P 2 O 7 H 3 ) IRC Intrinsic reaction coordinate KIE Kinetic isotope effect
13 KS Kohn Sham LDA Local density approximation LES Local enhanced sampling LJ Lennard Jones MC Monte Carlo MD Molecular Dynamics MEP Minimum energy path MM Molecular mechanics PB Poisson Boltzmann PBC Periodic boundary conditions PDB Protein data bank PES Potential energy surface PKIE Primary kinetic isotope effect PME Particle mesh Ewald PMF Potential of mean force PPi Diphosphate leaving group (P 2 O 7 4 ) PPTase Protein prenyltransferase PT barrel Prenyltransferase barrel QM Quantum mechanics QM/MM Quantum mechanics/ molecular mechanics RC Reaction coordinate RESP Restrained electrostatic potential RLS/RDS Rate limiting(/determining) step RMSD Root mean square deviation RMSF Root mean square fluctuation
14 SCC DFTB Self consistent charge density functional tight binding S KIE Secondary kinetic isotope effect SMD Steered MD TS Transition state TSS Transition state structure TST Transition state theory USP Umbrella sampling VDW Van der Waals VTST Variational transition state theory WHAM Weighted histogram analysis method WT W ild type ZPE Zero point energy ZPVE Zero point vibrational energy
15 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy COMPUTATIONAL MODEL L ING AND S IMULATION OF PRENYLTRANSFERASES AND CORRESPONDING PRENYLATION By Yue Yang August 2011 Chair: Kenneth M. Merz, Jr. Major: Chemistry Now adays with the power of supercomputers a lot of labor extensive work especially in scie ntific research can be accomplished efficiently. Such a good example is the application of computational simulations into chemistry and biology research, e.g., protein folding, high throughput virtual screening and so on. In fact, due to the limitation of experimen tal instruments or techniques, many aspects of enzyme and their associate enzymatic reactions have not been understood completely or even clearly, therefore, the advantage of computational simulation makes itself an excellent candidate to solve th ese problems. In our group, computational studies have been conducted on a variety of enzymes and biological systems, such as human carbonic anhydrase, lactamase, prenyltransferase and even interactive enzyme RNA systems. Among these functional important enzymes, prenyltransferase has attracted more and more attention in recent years. As an important candidate for anti cancer drug design, protein prenyl transferases play a role in the posttranslational modification of Ras superfamily of enzymes and many others that have appeared in about 30% of human cancers 1 In this dissertation, we focus on our effort of conducting computational study to obtain knowledge from various prenyltransferases and corresponding enzymat ic reactions catalyzed by them.
16 Chapter one gives an o verall introduction of terms such like prenyltransferase prenylation, as well as the prenyltransferase enzymes that will be discussed in this dissertation. The background and (pharmaceutical) significa nce of this class of enzymes will also be explained. Chapter two is the method section, where computational theories, techniques mentioned in this dissertation, along with difficulties of our approach and the strategy employed to tackle the problem will be introduced. The next three chapters describe the work we have already accomplished in this big project including our prediction of magnesium ion binding in the FTase and related validation, the catalytic mechanism of FTase catalyzed farnesylation along w ith impact of key mutation and the prenylation mechanism and the proton transfer pathway of NphB With c omputational modeling and sampling we have solve d these chemical/biological problems t hat experiments cannot explain and revealed a lot of interesting facts of how enzymes adopt different strategies based on different environments or situations.
17 CHAPTER 1 INTRODUCTION 1.1 Prenylation and Prenyltransferase Prenylation, also called isoprenylation, refers to a type of attachment of a hydrophobic function group called isoprenoid or terpenoid Isoprenoid group is the derivative of a five carbon isoprene unit. S uch a group is usually called a prenyl group (Fig. 1 1) and commonly appears in this king of lipid modification reaction s The prenylation is of grea t importance to researchers because this reaction enabl es certain enzyme s, such like Ras superfamily of enzymes, to anchor in to the inner side of plasma cellular membrane 2 as well as modifies the biological activity of certain aromatic complex es 3 5 Figure 1 1. Chemical structure of isoprene and prenyl group. Ubiquitous in natural products, isoprenoids are very important chemicals that can be found in almost every class of living things, including plants, fungi, animals and organisms 6 Among them, plant terpenoids are extensively used beca use of their aromati c qualities. A number of plant terpenoids are well known, such as citral, menthol, camphor, cannabinoids and Salvinorin A. Their contributions to nature such as the color of yellow flowers and the flavor of cinnamon and ginger, have been well re cognized a nd studied. In addition this kind of plant terpenoids play s an important role in traditional herbal remedies. Recently investigation showed that some isoprenoids have certain pharmaceutical desired characteristics such like antibacterial, antineoplastic anticancer effects, making this king of chemicals more important
18 More importantly, several members of this isoprenoid family play a crucial role in the signal transduction pathway. Examples of these terpenoids include the 15 carbon farnesyl and the 20 c arbon geranylgeranyl 2 7 The prenylation that attaches one of these two group s to a cysteine residue at or near the C terminus of the target intracellular enzyme is usually referred to as protein prenylation. Many protein s including various members of Ras superfamily require such a posttranslational attachment for certain functions. 1.2 Protein Prenyltransferase Common ly, enzymes that catalyze the prenylation reaction are called prenyltransferase. Consequently, the enzyme that catalyze s the protein prenylation is usually termed as protein prenyltransferase 2 8 14 There are three members in this family: protein farnesyltransferase (FTase) 7 15 18 protein geranylgeranyltransferase type I (GGTase I) 19 23 and type II (GGTase II) 8 23 25 Protein prenyltransfer ases can be further categorized into two subclasses: the CaaX prenyltransferase subclass that includes FT ase (Figure 1 2) and GGTase I (Figure 1 3) and Rab GGTase subclass, which contains GGTase II Suggested by the name CaaX prenyltransferases are in charge of the catalysis of attaching a farnesyl group or a geranylgeranyl group to a structurally fixed cys teine residue at the fourth position from the C terminus of the target enzyme. In this four letter term, C stands for cysteine, a refer s to aliphatic AA, while X is usually taken by alanine, serine, methionine or phenylalanine in FTase targets and occupied by either leucine or phenylalanine for GGTase I catalysis 26 27 On the other hand, GGTase II is distinct from both FTase and GGTase I due to the fact that it mainly attaches two geranylgeranyl groups to two C terminus cysteine residues in a Cy s Cys, Cys X Cys or Cys Cys X X motif in Rab family of enzymes 13 28
19 Figure 1 2. Protein farnesyltransferase (FTase) catalyzed farnesylation. Figure 1 3. Protein ge ranylgerany ltransferase (GGTase ) catalyzed geranylgeranylation.
20 Many e xperimental studies, in particular those s equence studies, has uncovered some secrets of the CaaX motif: first, Reiss and coworkers made further classification of Ca 1 a 2 X and revealed that a charged residue at a 1 position would slightly reduce the affinity while at a 2 position or X position would be associate with more drastic reduction 29 ; later Troutman and coworkers have discovered that different C a 1 a 2 X motifs would give different stimulate to the product release, which would further give FTase different target selectivity 30 ; more recently, Hougland and co workers investigated the dependency of FTase catalysis rate constant on variant Ca 1 a 2 X motifs and claimed the recognition of target substrate by FTase is not invariant but context dependent 31 In addition although FTase and GGTase I primarily target on the Ras superfamily of enzymes possessing this CaaX motif, various short peptide s with this same motif can also be recognized and catalyzed by these two prenyltran sferases 2 9 Because of their important pharmaceutical potential, many studies have been conducted to disclose the nature of these enzymes. As a result, a great deal of useful structural information of these PPTase s has been exposed FTase is a heterodimer composed of a 48 kDa 46 erestingly, it shares an identical while possesses a of 43 kDa that is slightly smaller than FTase GGTase II is distinctly different to FTase and GGTase I, although it is a heterodimer too. This enzyme possesses a which is 60 kDa and the kDa. All of three transferases are metalloenzymes, however, with different metal dependency. FTase requires a zinc ion for substrate binding and chemical reactivity and a se cond magnesiu m cation for optimal reactivity 18 32 GGTase I also requires a zinc ion for substrate binding and cata lysis activity, but does not need a magnesium ion, which is possibly due to the substitute role
21 played by Lys 311 18 20 21 GGTase II on the other hand, appears to require a single Mg 2+ for reactivity but not showing any dependency on Zn 2+ making the metal dependency of this family of enzymes more intricate and complicated. Each of three prenyltransferases b inds a corresponding isoprenoid substrate and an enzyme ( or peptide ) substrate prior to catalysis. The substrates in FTase are farnesyl diphosphate (FPP) the visual signal transduction pathway, or short peptide s with the CaaX motif. GGTase I adopts geranylgeranyl diphosphate (GGPP) as isoprenoid substrate and an enzym e substrate of either a Ras related GTPase or one of the the same isoprenoid substrate as GGTase I but targets primarily on R ab family of enzymes that possess those specific double Cys motifs. Although the binding of diphosp hate substrate does not require the binding of zinc ion in both FTase and GGTase I, such a metal binding is the prerequisite condition for either enzyme in order for the binding of enzyme or peptide substrate 19 33 Therefore, i t explains why a zinc ion is ind ispensible to the activity of CaaX prenyltransferase, despite the fact that FTa se also re quires a magnesium ion to fully restore its activity. However, from experimental results, it is still not clear whether Zn 2+ has the same function in both enzymes It is also unclear that whether the zinc ion play a solo structural role or it also directly involves in the catalysis. An assumption is made th at the zinc ion activate s the sulfur or sulfhydryl of the substrate cysteine residue in CaaX motif making it more nucleophilic, as observed in Ada, an enzyme that is known for DNA repair via a catalysis similar to the pren ylation s discussed here. 34 In this dissertation, effort will be made focusing on FTase and GGTase I. GGTase II on the other hand wi ll not be widely discussed.
22 1.3 Aromatic Prenyltransferase An incredible increase in the amount of studies focused on prenylation has been seen in the past ten to twenty years. Though a fairly large portion attributed to PPTases that catalyzed the posttran slational S prenylation required by certain enzymes in order for the later penetration into the plasma membrane, other prenyltransferase, such as those capable of the lipidic modification of various aromatic substrate s, have also been identified and invest igated 5 35 This type of prenyltransferases is known as aromatic prenyltransferases (APTase) 36 40 Examples of such APTase include NphB (also called Orf2 before rename), a 33 kDa soluble monomeric bacterial prenyltransferase with 307 residues expressed in Escherichia co li 3 and FtmPT1, a 464 residues fungal indole prenyltransferase 4 1.3.1 Bacterial Aromatic Prenyltransferase NphB Figure 1 4. NphB (middle up) catalyzed geranylation. C 5 prenylation corresponds to the major product. C 2 prenylation represents the minor product, which is in a product ration of 1:10 to the major product. C 4 prenylation refers to the extra minor product. NphB (Figure 1 4) i s obtained from Streptomyces sp. Strain CL190 3 It is originally called Orf2 that is named after being revealed as one of three new open reading frames (ORFs) in the
23 upstream region of MVA pathway genes containing gene clusters. It shows fairly high degree homologies to three other bacterial pr enyltransferases : HypSc (GenBank accession # AF939130) CloQ (GenBank accession # AF329398) and NovQ (GenBank accession # AF170880) 5 However, with regard to 3 di mensional (3 D) structure, substrate selectivity, and metal catalyst dependency it exhibits distinct difference to the other three APTases mentioned above not to mention any other P P Tases 7 18 41 45 CloQ, NovQ and HypSc display strong selectivity of five carbon DMAPP substrate but show little activity for isoprenoid diphosphate with longer carbon chain, such like ten carbon geranyl diphosphate (GPP) or fifteen carbon FPP I n addition, all three are independent to magnesium catalysis 3 5 To the contrary, NphB shows significant Mg 2+ dependency. It exhibits no re activity to dimethylallyl diphosphate ( DMAPP ) but strong activity to GPP and relative low but observable reactivity to F PP Additionally, NphB also di s plays flexible aromatic substrate s electivity. Its known substrates include flaviolin, 1,6 dihydroxynaphthalene (1,6 DHN), 2,7 DHN, apigenin and many others. Structurally, NphB consists of a single domain The 3 D structure h ighlights a novel barrel fold, which is termed as prenyltransferase barrel ( PT barrel ) consisting of ten anti sheets surround ed helices that are expo sed to solvent. Although this barrel looks similar to the TIM barrels, secondary structure of NphB shows no sharing with TIM barrels 46 49 The ternary complexes of NphB with Mg 2+ GSPP, a GPP analog and 1,6 DHN or other aromatic s ubstrate reveals a solvent accessible core and a spacious binding pocket inside provided by this PT barrel Such a voluminous binding pocket allows the binding of a large variety of aromatic substrates, thus partially explaining the flexible substrate selectivity exhibited by NphB
24 Although this spacious binding pocket is capable of binding various pr enyl diphosphates including DMAPP, GPP and FPP, NphB displays highest reactivity to GPP 3 Inside the binding pocket, GPP is encircled by several positively charged residues including Lys119, Arg228 and Lys284. These positive charge residues, along with Tyr21 6 and Asn173, tether the diphosph ate group through hydrogen bond interactions helping to anchor the isoprene substrate. Magnesium provides additional stabilization to GPP by interacting with a non di phosphate oxygen atom th r ough an octahedral coordination The coordinating ligands also include Asp62 and four water molecules. However, a (N/D)DXDD signature motif observed in many other Mg 2+ dependent isoprenoid diphosph ate containin g enzymes is absent in NphB. Unlike CloQ, NovQ and HypSc, NphB shows significant Mg 2+ dependency Whether this metal cation directly involved in the catalysis is yet know n What has not been explained well also is the regioselectivity displaye d in the system of NphB complexed with Mg 2+ GPP and 1,6 DHN. T wo products were originally isolated with the isoprenoids group attaching to different carbon atoms each neighbouring to one of the hydroxide connecting carbon atoms Under the condition that e ach substrate reaches its saturation the product ion ration is 10:1. Later, more geranylation loc ation of 1,6 DHN has been discovered although the production ratio was even lower. Evidences show ed that certain polyketide based aromatic complexes such as the anti oxidant natural product naphterpin have quite different biological activities to their isoprenoid modified derivatives 50 53 In addition, these aromatic derivatives that possess either a 5 carbon dimethylallyl, a 10 carbon geranyl or 15 carbon farnesyl display ed desired pharmaceutical properties such as anti microbial, anti oxidant, anti inflammatory, anti viral and anti cancer
25 effect. Therefore, it is of great importance to study NphB chemistry computationally to complete our knowledge of this enzyme 1.3.2 Fungal Indole Prenyltransferase FtmPT1 FtmPT1 (Figure 1 5) is another APTase that has just been purified from Aspergillus fumigatus in recent years 4 Unlike NphB or CloQ that was obtained from bacteria, FtmPT1 and several other indole prenyltransferases were characterized from fungi, thus they were usually referring to as fungal indole prenyltransferase. In general, these fungal indole prenyltransferases share very few sequence similarities with other APTases and show no divalent metal ion dependency 54 55 Th is type of prenyltransferases catalyze s the isoprenoid attachmen t of tryptophan and its derivatives. The product of such catalysis such like tryprostatin B and tryprostatin A, display high cytotoxicity for several cancer cell lines, making them selves promising anti tumor agents 56 58 Figure 1 5. Indole prenyltransferase FtmPT1 (mid dle up) catalyzed attachment of dimethylallylic group to brevianamide F (left bottom). The product is tryprostatin B (right).
26 The recently characterized 3 D structure of FtmPT1 features a PT barrel consisting 10 anti h elices, similar to that observed in NphB. The secondary 5 representation helices coupled two wo substrate, dime thylallyl diphosphate (DMAPP) and brevianamide F, sit uat e in the center of the PT barrel. Similar to NphB, the active site in FtmPT1 is also spacious and solvent accessible. Despite similar ities in the core structures the active sites of FtmPT1 and NphB are distinct FtmPT1 T he hydroxyl groups of Tyr203, Tyr296, Tyr382 and Tyr450 make up a hydrophobic circle therein binds brevianamide F and DMAPP. Although NphB also possesses an aromatic rich environment inside its bi nding po W ith the dimethylallyl group completely embedded, it seems like specifically designed to stabilize and engineer the potential carbocation through cation Comparing AA seque nce of FtmPT1 and homologue enzyme FgaPT2, which is another fungal indole prenyltransferase possessing a PT barrel and represents an excellent superposition (RMSD of C based on 402 residues measures 1.3) with FtmPT1, differences at several key residues ar e discovered and worth noting. For instance, Gly115 and His279 (in FtmPT1) are replaced by Thr102 and Arg244 in FgaPT2, respectively. These two substitutions spare more space in the binding pocket of FtmPT1, allowing the enzyme to bind larger indole substr ate, such like brevia namide F. On the other hand, Glu 102 of FtmPT1 has been found conserved in nearly all fungal indole prenyltransferases. This negative charged residue forms hydrogen bond with indole nitro gen, stabilizing the substrate. It also possibly functions to abstract the abundant proton after the prenylation step, in order to reestablish the aromatic qualify of the prenylated
27 indole substrate In addition, it is interest ing to notice that a small modification of AA sequence can result in the attac hment of dimethylallyl group from a different DMAPP carbon to a different carbon atom in indole substrate. Fungal indole prenyltransferases display interesting diversity and play an important role in biosynthesis of a variety of chemica lly important compou nds 54 55 In this dissertation, progress in the study of FtmPT1 is reported. 1.4 Computational Chemistry /Biochemis try The past 30 years have seen significant improvement in computer science and technology. The introduction of workstations and supercomputers has opened the door to huge calculations, generating a new branch o f chemistry that utilizes the assistance fr om computer science to solve chemical problem s Such branch of chemistry is usually known as computational chemistry. Computational chemistry has been well known for its capability of calculating properties such as absolute or relative energies, electronic charge distributions and vibrational frequencies, based on the application of quantum mechanics (QM). Back into the 1940s, the idea of solving complicated and labor intensive wave functions for atomi c systems has already come out. I n the 1950s the first s emi empirical atomic orbital calculation an d Hartree Fock (HF) calculation of diatomic molecule have been carried out. L ater in the 1960s the empirical methods based on Hckel method were introduced. Coming into the 1970s, methods of molecular mechanics (MM) were developed. As of today, more advanced QM calculations and more complex hybrid QM/MM simulations have been widely implemented resulting in many exciti ng results. The implementation of computational chemistry used to be strictly limited by several factors such as the size of the system due to the rapidly increasing in computational expense associated with the implementation of higher level theory or the calculation of larger system. I n present, with supercomputers and parallel computing, it is comm on to perform highly accurate
28 coupled cluster calculations on large organic molecules or carry out molecular dynamics (MD) simulations on large bigger systems. The implementation of graphics processing unit (GPU) allows 100 ns simulation per day for classical MD simulation s of a system composed of ten s of thousands atoms making in silico investigation of microsecond level or even bigge r time scale biological activity possi ble. Moreover, with super fast machines specifically optimized for MD simulation such as Anton introduced by D. E. Shaw Research Lab, a y of classical simulation of a 20,000 + atoms system has become feasible 59 More and more powerful computer techniques have pushed computational chemistry into divers e research areas such as bioinformatics. The ability to build and analyze huge databases with high efficiency has made computational appro ach a very important part of modern drug design in pharmaceutical industries. By providing an alternative and more environmental way to study the properties of dru g molecules, the correlations betwee n structures and properties, as well as the possibility of finding or even designing a more suitable substrate, computational chemistry has saved countless labors and other lab expense s making huge impact in industry. T oday, with less expensive computational expense and higher accuracy, computational chemistry has already become an i ndispensable branch of chemistry and biochemistry It is able to investigate many areas that are not accessible by traditional experimental approaches In addition, it can also make valuable predictions that can in turn gui de experiments design ing In this dissertation, computational effort has been made on mechanistic study of several important prenyltransferases. Great results and interestin g prediction s are reported.
29 CHAPTER 2 TH EORY AND METHODS Before computers bec ame popular in scientific research, molecular modelling was used to associate with the utilization of pens and paper to make simple calculation s of small molecules mimicking the ir behaviors. In present with the rapid development of computer science and tech nologies, molecular modelling becomes rarely independent to computational chemistry. Utilizing the power of supercomputers and parallel computing, molecular modelling has been able to apply on to macromolecules and even bigger systems like enzymes. T hus it becomes more and more important in the researches of che mistry, biology, drug de s ign and many other areas Today, molecular modellers have implemented all available computatio nal chemistry methods, including QM, MM, MD, minimizations, simulations, conformational analysis, protein ligand interactions, ch em/ bi oinformatics and other computer based techniques to elucidate and predict molecular behaviors. In the rest of thi s chapte r, many important concep ts and methods in molecular modelling will be discussed. 2.1 Quantum Mechanics and the Born Oppenheimer Approximation The fundament and heart of computational chemistry based molecular modellin g is QM, which was introduced to descri be the behavior of microscopic particles. Today, although molecular modelling does not necessarily mean QM, every successful model must be able to find its basis in quantum mechanics. It is well know n that in classical world, macroscopic systems obey Newto nian mechanics T hus the change of their energies is continuous. However, rules in the microscopic world are quite different: particles that form matter have a characteristic of quantization. As a result, their energies are discrete. In fact, microscopic systems possess both waveli ke and particle like
30 properties. This characteristic makes no pre existing theory, including wave mechanics that possesses the quantization characteristic ca pable to explain their behavior. Therefore a new mechanics cap able to interpret this dichotomy is required. It was under such a circumstance that QM was introduced. The base of QM is the assumption that wave functions can be applied to microscopic systems and can describe all their properties. In detail, the properties are c haracterized by applying corresponding operators onto the wave functions. E ach observable physical property is assumed to have one such operator The most important operator in QM theory is the Hamiltonian operator, H which describes and returns the syste m energy. The mathematical equation of such a relationship is (2 1) Equation 2 1 is commonly known as Schrdinger equation. Hamiltonian operator generally consists of five parts: the 1 st and 2 nd part describes the kine tic of electrons and nuclei, the 3 rd part accounts for the attraction between electrons and nuclei, the 4 th part attributes to interelectronic repulsion while the 5 th part corresponds to the internuclear repulsion. Simply speaking, by solving the Schrdinger equ ation the energy of the system will be obtained. However, in practi ce, this is difficult to accomplish since Schrdinger equation cannot be exactly s olved for any molecular systems due to the interdependency between electronic motion and nuclear motion. T herefore, another fundamental postulate is required to simplify the case, which is the famous Born Oppenheimer approximation. Considering the fact that either proton s or neutron s are more than 1000 times heavier than electronic s it is rational to assume t hat the motion of electrons is instantaneous compared to
31 nuclear motion. Thus, to calculate the electronic energies is possible since the nuclear motion and the electronic Schrdinger equation has a reform of: (2 2) Here, the electronic and nuclear, respectively. As Born Oppenheimer approximation is unconditionally applied in any case of molecular tion of the nuclear coordinate. In such a case any change of nuclei positions varies the energy of the system. The changes of energy against the all of the possible nuclear or atomic positions construct a hyper surface denoted as potential energy surface (PES). The stationary points on the PES are of particular interest to computational chemists because at those points the first derivative of energy with re spect to the coordinate is zero, t he physical meaning of which demonstrates the force s at certain poi nt s are zero. Both local minima and saddle points are stationary points however, in different types 2.2 Molecular Mechanics QM deals with electrons in molecular system, usually giving the highest accuracy among all molecular modelling methods. However, i t is not occasional for molecular modellers to find that the problems they are going to solve are too large to be considered quantum mechanically. Molecular mechanics (also called force field methods) are introduced in such a situation to With the implementation of Born Oppenheimer approximation, it is possible to compute the energy of a system as a function of nuclear coordinates with electronic motion ignored, which is the basis of MM. Comparing with QM, obviously MM uses a simpler model to mimic interactions with generally five components included, a mongst bonding stretching, angle bending
32 and bond torsion that all represent the energy penalties associated with deviation from their reference position are considered bonded interactions, whereas the other two terms, van der Waals ( VDW ) and electrostatics are categorized into non bonded interactions. Thus, the potential energy of a MM modeled system can be written in the following form: (2 3) In equat ion 2 3, the first three terms on the right are bond stretching, angle bending and bond torsions separately, while the latter two terms are VDW and electrostatics, respectively. The schematic representations of these five terms are given in Figure 2.1. Figure 2 1 Schematic representations of five components in general MM force field: bond stretching, angle bending, bond torsion, non bonded VDW and electrostatic.
33 Today more advanced MM force field such as MM4 60 62 have been introduced, however, the five contributions above nearly exis t in every known MM force field It is very important to bear in mind that all MM force fields are empirical, which means there i s absolutely no correct form for any force fields. Also, i t is worth noting that the transferability, which means whether the same set of parameters can be extend ed to successfully mod el a set of different molecules or enzymes, is the most important featur e for every force field Otherwise, although a certain force field can well model the system that it was parameterized against, another desired function of molecular modelling, the ability to make predictions (to other systems), is lost. 2.2.1 Bond Stretch ing Taking a careful examination at equation 2 3 again, the bond stretching term is in the form as shown in equation 2 4, however, it is in fact the second order truncation of Taylor expansion about the potential energy function a bond at l ength l i based on the assumption that the potential at equilibrium bond length l i0 is zero. (2 4) This truncation makes the energy function working excellent to bond stretching near the reference points but brings physically unrealistic energy when the bond is largely stretched. A practical solution is to include the cubic term: (2 5 ) Here can then counterac t the unrealistic high potential energy associated with large s tretching. However, equation 2 5 also leads to a chaotic result that the lowest potential energy represents bonds dissociation. Again, the solution is to include the quart ic term in the Taylor expansion:
34 (2 6 ) Equation 2 6 is used in the general organic force field MM3 63 65 However, including higher order terms in Taylor ex pansion increases the accuracy as well as the computational cost distinctly thus equation 2 4 In fact, there is an alternative way to provide higher accuracy of describing potential energy contribution ass ociated with bond stretching, the Morse function, (2 7 ) Here D i is the dissociation energy of bond and i is a fitting constant. However, this approach has been proved much computationally less efficient than those truncations of Taylor expansion. 2.2.2 Angle Bending The second term in equation 2 3 is to calculate the potential energy deviation associated w ith angle bending. (2 8) 8 indeed comes from another second order truncation of Taylor expansion. It is worthy noting that the force constant of angle bending is signifi cantly smaller than that of bond stretching, due to the energy required to disturb an angle away from its equilibrium pose is much less than to compress or withdraw a bond. Like any other truncation of Taylor expansion, including higher order terms a ssocia tes with greater accuracy. One such example is MM3, which includes up to sextic term in the expansion for certain angles. Certainly, the computational expense increases when more terms in the expansion are included.
35 2.2.3 Torsion and Improper Torsion The t hird term in equation 2.3 reflects the potential energy change associated with the deviation of dihedral angles or torsion angles, as shown in equation 2 9. (2 9) in equation 2 9 is the torsion angle; n is the multipli city that means how many energy minimum can be located through a 360 scan; is called the phase factor, defining the location of energy minimum; and V n is the barrier height, which refers to the re lative barrier of rotation. The torsion term might be considered as the most important term among the bonded interaction components, because a large amount of energy is required to alter the reference conformation for both bond stretching and angle bending terms while torsion term and non bonded interaction terms accounts for most variation in potential energy and system structure. Sometimes, another type of torsion is also required to computationally reproduce experimental structures. This type of torsion is always defined as improper torsion (as shown in equation 2 10 which means a torsion angle is not bonded in the order of 1 2 3 4. One good example is cyclobutanone. Experimental results exhibited that the oxygen atom is in the same plane as the cyclobut ane ring, whereas computing without such an improper torsion results in an out of plane oxygen atom. The reason for this discrepancy is simple, in real world the C C=O angle is 133 in order to have oxygen atom in the plane, however, in calculation this an gle still adopts a value close to its reference value of 120. (2 10) In proteins, dihedrals play very important role in determining their structures. There are three torsion angles in protein structure (Figure 2.2): the dihed ral of C N C dihedral of N C C C N C
36 while the other two torsions have certain conformation. Figure 2 2. Torsion (or dihedral) angles in protein structure. 2.2.4 Van der Waals Interaction The interatomic force is usually divided into two parts: the long range attractive force an d the short range repulsive force. The former is commonly known as dispersive forces, as well as quantum mechanics 66 The origin of such a dispersive force is the instantaneous electrical dipole moments developed by atoms. In general, without a permanent charge, the dipole dipole interaction resulted from the attraction between an instantaneous dipole and another dipole in the neighbor atoms induced simultaneously by the first dipole, is the strongest interaction of such a dispersive force. Schrdinger equation of this model is in the same form as that of a simple harmonic oscillator. Thus the dispersive interaction energy can be app roximately given by equation 2 11,
37 (2 11) Obviously, the dispersion energy is dependent on 1/ r 6 The dispersion force describes the attractive interaction when two atoms (or particles) are quite away from each other. When they approaching each other into a certain range, for example, within 3 , even a very small decrease in distance will result in large increase in potential energy. The quantum mechanical basis of this phenomenon is any two identical fermions may not occupy th e same set of quantum numbers, as described in Pauli exclusion principle (though some people believe it is actually attribute to nuclear repulsion) Thus the repulsive force in short range is usually known as exchange force. It is possible to calculate the dispersive force and exchange force quantum mechanically, however, certainly with a high computational cost since electron correlation must be included in the treatment. In MM, the VDW potential is modeled using a much simpler empirical expression, the Le nnard Jones 12 6 function. (2 11 ) In equation 2 11 is the collision parameter that physically corresponds to the distance o the well depth that represents the mi nimum potential energy. Clearly the sixth power term inside the square bracket is referring to the dispersive force that has been shown to have an exponentially dependence of 1/ r 6 and the other term is corresponding to the exchange repulsive force that va ries as 1/ r 12 However, it is worth noting that this twelfth exponential does not have a significant physical meaning, but instead has a favor in calculation since it is simply the square of 1/ r 6
38 There are other VDW energy functions introduced to improve on reproducing the experimental data, such as Halgren buffered 14 7 potential 67 72 However, Lennard Jones (LJ) 6 12 potential remains the most popular choice of calculating potential energy associated with VDW interactions. 2.2.5 Electrostatic Interaction The electrostatic interaction is caused by the unequal distribution of charge in a molecular system. By s acrificing the accurate electron electron interactions, MM usually adopts a so called point charges approach to represent the charge distribution of a molecule. Restricted to atomic centers, the point charges are commonly known as partial atomic charges or net atomic charges. Then it is straightforward to appl in equation 2 12, where N A and N B are the number of points charges in two molecules. (2 12) Central multipole expansion provides another way to calculate intermol ecular electrostatic interaction energy. This method considers a molecule as an entity and separate electrostatic interaction into different electric moments such as charge, dipole, quadrupole, octopole and so on (Figure 2 3) With multipole can be describ ed as a proper distribution of charges, dipole, quadrupole and octopole can be represented by two, four and eight charges correspondingly, making calculations simpler. One good example of using such an expansion is the benzene model 73 which demonstrates significant advantage of including quadrupole. However, on the other hand, thi s kind of multipole expansion c an be troublesome when it is applied to large systems and cannot be employed to describe intramolecular interactions. The implementation of fixed point charges provides a rather simple and straightforward electrostatic functi tice that the chan ge
39 of distribution of an atom or a molecule can induce dipole to neighboring atoms or molecules. Therefore, in order for better description of electrostatic interactions, it is nece ssary to include not only interactions between fixed charge distributions but also interactions caused by polarization. To describe polarization is not a straightforward task, but somewhat tricky and complex, because one must bear in mind that when a molec charge distribution of a neighbored molecule its own charge distribution is in turn be modified by the induced dipole of the neighboring molecule. Many approaches have been proposed to model polarization effect 74 77 However, due to the significant increase of computational cost, the p olarization effect is not included in all force field but instead considered only in the case that it is indispensible. Figure 2 3. Schematic representation of some multipole examples. 2.2.6 Cross Terms It has been mentioned that in electrostatic intera ction, the charge distribution of two
40 deviations of bond stretching, angle bending or bond torsion are not independent but in fact as well The nomenclatu re of Mathematically, the cross terms can be revealed by expanding the potential energy of a molecular system in a multi dimensional Taylor expansion. For example, if we consider the bonds stretching when an angle is closing, as shown in Figure 2 4. Figure 2 4. Schematic representation of angle closing coupling with bonds stretching. Such a scheme can be formula ted in the following form. (2 13) Equation 2 13 has been applied in MM force field MM2, MM3 and MM4. Besides this example, there are many others cross terms, such as stretch stretch, stretch bend, stretch torsion, bend torsion and so on. Hwang and coworkers have proposed an interesting categorizing scheme for MM force fields 78 80 : based on whether a force field includes cross terms, force fields are divided into class I that does not have cross terms integrated or class II that should have both anharmonic terms and cross terms. Later, Allinger and coworkers further proposed that more advanced class III
41 force field should also have electronegativity, hyperconjugation and other features taken into account 81 82 2.2.7 Force Field Parameterization A force field is not simply a function but a potential energy function coupled with a set of pa rameters. Developing force field parameters, or force field parameterization, is always an intricate task requiring a great deal of effort. Once the function form of a force field is chosen, the next step is to choose a series of data. The parameters will then be developed to reproduce, or fit, these data. Ideally, experimental data are preferred, however, they are not always available. Thus, QM calculations using higher level ab initio theory is re quired to provide missing data. When both force field funct ion and selected set of data are de cid ed the next step is to carry out the parameterization. Two major a pproaches squares fitting to obtain the parameters that best fit to the selected data 83 85 An important implementation of MM force field is to model large biological systems such especially when explicit solvents are included. Therefore, to parameterize against the entire system is almost prohibitiv e. In fact, most popular force fields today, such as AMBER 86 89 CHARMM 90 92 GROMOS 93 94 OPLS 95 96 and AMOEBA 97 99 are all parameterized against small organic molecules. In this manner, the transferability of such force fields is v ery important. 2.3 Energy Minimization and Molecular Simulation In molecular modelling, one of the major goals is to calculate the potential energy and reproduce the PES. With potential energy and /or PES available many molecular properties and
42 thermodynamic properties can thus be computed. There are two totally different approaches to achieve this goa l: energy minimization and molecular/computational simulation. 2.3.1 Energy Minimization Ideally, with a PES available, it is simple and straightforward to locate energy minimum of a molecular system. However, in practice, this is only working with very si mple molecules. For a large system with N atoms, there are 3 N Cartesian coordinates and (3 N 6) internal coordinates; complicated. In such a case, locating the m inimum points along the entire PES, which is referred to as global minimum, can be troublesome. Energy minimizations are designed to search the global minimum and other energy minimums along the PES. Considering the potential energy as a function F of a se t of variables r1, r2, ri, the energy minimums have the following mathematical characteristics, (2 14) Equation 2 14 means the first derivative of function F with regard to any variables equals to zero while the second d erivatives are positive. It is worthy noting that s addle points, represent ing those stationary points on the PES but not local extremum points also have their first derivatives equal to zero. Although minimization algorithms with no derivatives of potent ial energy considered, prefer to use such derivatives since they can provide many useful information such as the shape of PES. In fact, experienced molecular modeler might find the application of energy derivati ves can boost the efficacy of locating the important points. Examples of most common minimization alg orithms includes, but not limited to, the steepest descent method, conjugated gradient method and Newton Raphson method.
43 184.108.40.206 Steepest descent m ethod Th is algorithm calculates the first derivative of energy function determining the direction of force. Then the coordinates of the system are gradually mo dified to make the system move in a direction that is parallel to the net forc e to approach the minimum point In this method, the structure of the system obtained at one step is used as the starting structure for the next step. Also, it is important to note that both the gradients and the direction of two connecting steps are orthogonal. Practically, steepe st descent method is good at rapidly finding out a structure that is very close to the minimum but rather inefficient in the final approaching, thus, this method is always carried out in the beginning stage of minimization if multiple minimization methods are carried out in order to accurately locate the minimum. 220.127.116.11 Conjugated gradient m ethod This algorithm also approaches the minimum by consequently moving the system just like what steepest descent method does However, different to steepest descent m ethod, conjugated gradient method does not require the orthogonal directions for two consecutive steps. In fact, within this approach, the gradients of two successive steps are orthogonal but the directions are Generally speaking this method is better at locating the exact minimum and in some cases can achieve this goal with fewer steps than steepest descent method when a rather good structure is given. However, it is interesting to notice that the move of first step of both steepest d escent method and conjugated gradient method should be identical. 18.104.22.168 Newton Raphson m ethod Different to steepest descent or conjugated gradient method that only use the gradient (the first derivative), Newton Raphson method uses both the first and the second derivatives to search the minimum. The inclusion of the second derivatives provides more information, e.g. the
4 4 curvature of the energy function, however, on the other hand, increases computational cost by possibly a significant margin because withi n this algorithm the Hessian matrix must be calculated and inverted. Newton Raphson method is always applied to small systems due to the high computational expense. In addition, this approach might not be a good choice given a structure not even close to i ts minimum. On the other hand, application of such feature is to search for saddle points. 2.3.2 Molecular/Computa tional Simulation Although energy minimization can discover the minimum points along a PES and thus provide information for predicting properties of a system, it is not reliable that the information provided is sufficient to make predictions especially to those thermodynamic properties. This is especially true given a large biological system with thousands of atoms. In such cases, the huge number of degrees of freedom and multiple minima separated by energy barriers make it extremely difficult to gather inf ormation covering the entire energy surface. In addition, most experimentally measurable properties especially those thermodynamic properties are indeed time average (equation 2 15) which cannot be theoretically reproduced or predicted using only energy m inimizations. Thus, it is clear that a technique different to energy minimization is required to accompany those attempts to understand the nature of the PES of a system. Computational simulation, or call it molecular simulation, is a good candidate for be tter understanding of the intrinsic characteristics of energy surface of molecular systems. Integrated evolving in time and thus to calculate those properties averaged over time using certain numerical equations However, the biggest problem encountered is that it is rarely possible to
45 find the initial configuration. To solve this problem, statistical mechanics was introduced by Boltzmann and Gibbs in order to c onvert studying a time evolving system into investigating a large number of system replicas simultaneously. Whereafter in associate with the ergodic hypothesis, ensemble average properties (equation 2 16) can substitute time average properties. (2 15) (2 16) In both equations, p and r corresponds to momentum and coordinate separately, while in equation 2 16 the angle p N r N ) is the probability density. Monte Carlo (MC) method and molecular dynamics (MD) method are the two most 22.214.171.124 Phase s pace For a classical system containing N atoms, each of its conformation has 3 N coordinates and 3 N mo menta for. Hence, the state of this system can be characterized by a set of 6 N values. Such a s tate that is defined by these 6 N values is called phase space. All computational simulations sample the phase space of the given molecular system. Sampling enoug h phase space is the key to calculate important system properties especially thermodynamic properties. 126.96.36.199 Monte Carlo m ethod In MC simulation, each configuration is generated by makin g a random modification, such like mov ing an atom or rotating a bond to its predecessor, thus this kind of simulation stores no time information. For each configuration, the potential energy is calculated. By applying a certain criteria, usually Metropolis criterion 100 in today, whether the new configuration is accepted or rejected is decided
46 After sufficient simulation of conf igurations a certain property can be calculated using equation 2 17, (2 17) where M refers to the numb er of configuration simulated. It is noteworthy that traditional MC simulations are carried out in canonical ensemble. In such ensemble, the number of particles ( N ), the volume ( V ) and the temperature ( T ) remain constant. 2.3.2. 3 A b rief i ntroduction to m o lecular d ynamics Unlike MC, MD derives new configurati 18). In detail, the force of each particle is first calculated by differentiating the potential energy function, then applying the equation 2 18 will give in the acceleration and velocity information, thus a real trajectory over time is generated. (2 18) (2 19) In MD simulations, thermodynamic properties can be obtained via averaging over time, as s hown in equation 2 19. Comparing equation 2 19 and equation 2 17, clearly the only difference is that the function of a certain property in MD simulations has momentum dependence that is absent in the function of MC simulations. 2.4 Molecular Dynamics MD g motion, propagating the system in phase spaces and generat in g a trajectory varies with time. For
47 a simple system such as a n one dimensional classical harmonic oscillator or a hard sphere model moving straightly with a constant velocity 101 the position and momentum information can be obtained analytically as functions of time. However, for more realistic molecular models, the force on each atom (or particle) is not changing independently and the motion of the system can hardly be solv ed analytically. Hence, in such situation, the equation of motion must be solved using certain integration algorithms. 2. 4.1 Finite Difference Method: Verlet and Leap frog Algorithm Finite difference method is commonly used to integrate the equation of mot ion when dealing with large and complicated molecular systems. The basic idea of this technique is that the integration can be separated into different small parts with a very small time interval of and the force on each particle is assumed to be consta nt during interval. Thus, the force, coordinate, velocity and acceleration can be solved for each step and the simulation can be propagated. One pop ular integrator based on finite difference method is the Verlet algorithm 102 The central idea of this algorithm is to use the position and accelerati on information at time t and the position information from previous step ( t ) to determine the position of new step at ( t + ), as demonstrated in equation 2 20. (2 20)
48 However, there are several clear disadvantages of this algorithm: the possible precision loss associated with the use of 2 r (t) and r (t ) along with much smaller term a (t) 2 the lack of an explicit velocity term and the trouble of getting positional information of r (0 ) Leap frog algorithm is actually a variant form of Verlet method 103 The mathematical expression is shown as: (2 21) H ere a (t) can b v (t) can be calculated from (2 22) The major advantage of leap frog algorithm over Verlet algorithm is that velocity term in included explicitly, however, the coordinates and velocities are unsynchronized, making this method unable to count the kinetic energy contribution at the same step that the coordinates are defined. 2. 4.2 Time Step and SHAKE The time interval mentioned in 188.8.131.52 is usually referred to as time step. Choosing a prop er time step is clearly important for setting up a MD simulation giving the fact that an inappropriate time step choice can make the simulation either unable to sample sufficient phase spaces (when it is too small) or instable (when it is too large). Most biological systems are fairly flexible. In such systems, the highest frequency vibration corresponds to the bond stretching, especially for those heavy atom to hydrogen bond. This kind of vibration is usually at 10 fs level. Hence, the appropriate time st ep should not exceed 1 fs. With such a small time step, it is
49 large time scale, e.g. ns or ms. In order to accelerate MD simulations without making it instable, co nstraint dynamics are introduced. In such a scheme, bonds or angles are forced to adopt specific values thus the degrees of freedom associated with those terms are tot ally removed. The most popular constraint technique is SHAKE introduced by Berendsen and coworkers 104 O ppositely equal forces are a pplying on each of the two atoms of a constraint bond. Thus, a constraint ij that freezes the bond length between atom i and j can be expressed as, (2 23) Here d ij denotes the constraint distance. The force on each involving atom resulting from this const raint can be calculated base on, (2 24) Here is the Lagrange multiplier. Applying equation 2 23 into equation 2 24, it is easy to get hence (2 25) In application, the constraints are applied iteratively, consequently, the Lagrange multipliers following incorporating equation 2 23 and 2 24 into propagating function such like equation 2 20 or 2 21. Practically, the SHAKE algorithm is usually applied on those bonds involving both heavy atom and hydrogen. Removing the degr ees of freedom associated with such bonds will not affect the remaining degrees of freedom since they are only weakly coupled.
50 By doing so, the time step can be set to up to 2 fs that is twice as fast as simulation without invoking SHAKE. 2. 4. 3 Molecular D ynamics at Constant Temperature In general, MD simulations are carried out in microcanonical ensemble where the number of particles N the volume of system V and the total energy E remain constant. Ideally, thermodynamic properties obtained in NVE ensemble are able to transfer to other ensembles such as canonical (constant NVT) ensemble. However, the prerequisite condition of such transferability is that the size of system should be infinite. Therefore, it is important to apply MD in other ensembles such as constant NVT ensemble. For this reason, different thermostats for MD simulations are developed. Kinetic energy in a MD simulation is closely related to its temperature, thus temperature can be calculated base on, (2 26) Here k B is the Boltzmann constant, N is the total number of particles and n c is the degrees of freedom being constraint. Obviously, temperature can be written as a function of momentum or velocity. Therefore, rescale velocity is the simplest way to regulate tem perature 105 Multiplying th e velocities at time t by a factor of the corresponding change in temperature can be expressed in the following form: (2 27) Now, the scaling factor can be written as. (2 28)
51 Thus, the temperature of a system can be controlled. Ther e are several alternative approaches to regulate the temperature for MD simulations. One of them is the Berendsen thermostat that is acting like weakly coupling to an external heat bath 106 The mathematical expression of such algorithm is given in equation 2 29. so (2 29) H ere variable is actually functioning as a coupling constant. Bigger results in weaker coupling while smaller brings tighter coupling. Empirically, good performance is achieved by setting coupling constant to 40 times of the time step ( t in equation 2 29). The ad vantage of Berendsen thermostat is that the temperature of a system is allowed to fluctuate, however, there is no guarantee that the temperature is canonical averaged over the entire system that sometimes results in the temperature of solute is lower than that of solvent. In fact, this is the biggest disadvantage of any velocity scaling thermostats. Stochastic collisions method is one of the thermostat algorithms developed to improve the temperature distribution thus to reach the canonical ensemble. Langevi n dynamics 107 is classified into this category 108 By mimicking t he Brownian movement, the force on each particle in Langevin dynamic s is composed of three com ponents, as shown in equation 2 30. (2 30) The first component F i is resulting from interaction between particle i and other particles thereby depends on positions of particles, the second component is due to the frictional drag on particle i from solvent with a collision frequency of i while the third part A i is a random force arisen from random fluctuation of particle i due to interactions with solvents. The random force A i obeys Gaussian probability distrib ution with correlation function,
52 (2 31) Here ij is the Kronecker delta function while ( t is the Dirac delta function. With equation 2 31, clearly the equation of motion of Langevin dynamics is temperature dependen t thus it can be applied to temperature regulation. 2. 4 .4 Pressure Regulation in Molecular Dynamics In molecular modelling, it is frequently desired to simulate in isothermal isobaric ensemble in which the number of particles N the system temperature T a nd the system pressure P are reversed. This ensemble is actually very close to the real lab condition under which most experimental measurements are made. Pressure fluctuation is usually large cause by rapid change of virial, rdU(r)/dr with r The relatio n of system pressure and virial is demonstrated in equation 2 32. (2 32) U(r ij ) is the interaction energy between atom i and j Pressure regulation in MD simulation is usually accomplished by controlling system volume, because of clear dependence of pressure on volume from equation 2 32. While on the other hand, since N is fixed, regulating density can also be applied to barostat. Volume fluctuation can be quantified by introducing isothermal compressibility defined by equation 2 33. (2 33) Under isobaric condition, can be expressed in another form.
53 (2 34) t is possible to use Berendsen barostat to regulate system pressure using equation 2 35. (2 35) H ere p is the pressure coupling constant Supposing the volume is scaled by a factor of which corresponds to scaling the coord inates by a factor of 1/3 Then the mathematical solution of can be expressed as: (2 36) Constant NTP ensemble is widely used in MD simulations of biological systems. Important thermodynamic properties such as Gibbs free e nergy can only be determined under this ensemble (whereas under constant NTV canonical ensemble Helmholtz free energy is instead measured). 2. 4.5 Solvent Effects and Water Models Simulation in the gas phase is straightforward because of its simplicity. How ever, most chemical reactions and biological activities are taken place in solution. One good example is the Menschutkin reaction of ammonia and chloromethane, which will take place in solution (water) but instead show no reactivity in gas phase. Therefore it is of great importance to i ncluding solvent effects into molecular simulations. In MD simulations, the most important solvent model is the water model. In present, there are two very different types of approaches developed: the implicit solvent models that appl y a
54 so called reaction field to the simulation while the explicit water model that impends a water box to the system. 2. 4.5 .1 Implicit solvent m odels In implicit solvent models, solution water is treated as a dielectric continuum and the solvent effect is incorporated into the simulation as a reaction field, which is actually a continuous electric field. The free energy required to build the charge distribution can be calculated based on equation 2 37. (2 37) ( r ) i s the electrostatic potential and r ) is the charge density. Depending on different models, the latter can be written in the form of either continuous functions of r or discrete charges. The Poisson Boltzmann (PB) model is based on P oisson equation (equat ion 2 38), but practically the linearized PB equation (equation 2 39) is always used 109 (2 38) and (2 39) H ere denotes Debye Hckel parameter and I is the ionic strength of the electrolyte solution. It is worthy noting that the linearized PB equation cannot be solved analytically, thus, it has to be calculated iterativel y, which is time consuming Suppose a charge q is diffused uniformly to a conducting sphere surface with a radius of d then, at any point on the surface, the charge density is (2 40)
55 Since the dielectric constant of a condu ctor is infinite and at anywhere inside the conductor the electrostatic potential is zero, the electrostatic potential at r on the surface can only be determined from outside by using equation 2 41. (2 41) Thus, equation 2 37 can be written as (2 42) Since the solvation free energy G solvation is the difference of free energy in gas phase and in solution, it can thus be calculated. (2 43) Equation 2 43 is the Generalize d Born (GB) equation, which is the basis of GB model 110 114 Comparing to PB model, GB equation can be calculated analytically and thus computationally more effic ient. By extending the GB equation to polyatomic system, equation 2 43 becomes where (2 44) Variable d i is the effective Born radius of charge q i 2. 4.5 .2 Explicit water m odels Under certain circumstances, for instance, no water or other solvent molecule directly involve in interactions with solute system, implicit solvent models are popular choices. How ever, in most cases, if computational expense is affordable, u sing explicit water models is more desired, since the essential nature of chemistry lies on molecules and atoms.
56 There are several types of explicit water models developed. The basic type of mo d els adopts a rigid geometry with three to five interaction sites (Figure 2 5), thus it is also referred to as simple interaction site model. Because of rigid geometry, their potential energy function s have no bonded term. In addition, only VDW interactions between oxygen atoms are calculated. TIP3P 115 is one of the most popular water models in present. It falls into 3 site category with a 0.9572 O H distance and a 104.52 H O H angle. The 4 position with the n egative charge transferred from oxygen. One such example is TIP4P 115 which has exactly the same bond length and angle as TIP3P. ST2 116 is the most commonly used 5 site wa ter model, with negative charge of oxygen put on two lone pair sites. It is notable that two models within same number of sites do not have to have identical paramet ers, as demonstrated in Table 2 1 Figure 2 5. Simple interaction site water models. Par ameters developed for each water model should allow bulky water to reproduce their natural properties, such as to give a maximum density at 4 C. Unfortunately, the most widely used TIP3P model failed to reproduce this property, but its successors TIP4P an d TIP5P 117 are capable of doing so. The implementation of rigid geometry in these water models is able to save computational cost, however, on the other hand unable predict some properties, such as vibrational spectrum du e to the lack of internal flexibility. More advanced water models that
57 allow flexible geometry 118 or include polarization effect 77 119 have been developed to improve the simulation of water solvent effect in biological systems. 2. 4.5 .3 Periodic boundary c ondition s and p article m esh E wald When explicit solvents are employed, it is important t are one of such applications that allow measurable properties to be obtained from simulations with much smaller numbers of particles 120 The central idea of PBC is that a box is duplicated in all direction and a Figure 2 6. In practice, a two dimensional (2 D) PBC application will have 8 duplicated boxes while a 3 D application will have 26 such boxes. The duplicated box is also called image box and coordinates of a particle in the image box can be determined based on that of the corresponding in the original box. Today, there are several box choices with different shapes, amongst the rectangular and truncated octahedral box are the most widely used. The latter is especially good for simulating enzymes/proteins because its spherical shape is close to the real shape of enzymes. Table 2 1. Parameters of several water models. SPC/E 121 TIP3P BF 122 TIP4P ST2 # of sites 3 .0000 3 .0000 4 .0000 4 .0000 5 .0000 d OH () 1. 0 0 00 0.9572 0.95 00 0.9572 1.0 000 HOH () 109.47 00 104.52 00 105.7 000 104.52 00 109.47 00 q(O) 0.8476 0.834 0 0.0 000 0.0 000 0.0 000 q(H) 0.4238 0.417 0 0.49 00 0.52 00 0.2375 q(M)/q(L) 0.0 000 0.0 000 0.98 00 1.04 00 0.2375
58 In MM simulations, electrostatic decays as r 1 thus it is important to extend the calculation to distant. However, to calculate this kind of long range forces can be very time consuming. Ewa ld summation 123 was developed to sp eedup such calculations by divide the original summation into two series, as demonstrated in equation 2 45 and Figure 2 7. (2 45) Figure 2 6. Part of 2 D periodic boundary conditions (PBC). Note in 2 D 8 image boxes should surround the original box while in 3 D this number should go up to 26. By represent function f(r) by using complimentary error function (equation 2 46). (2 46) B oth terms on the right side of equation 2 45 converge faster th an the original form 1/r Therefore, the summation is expedited. More recently, several variants of Ewald summation has
59 been developed, including particle mesh Ewald (PME) 124 and particle particle particle mesh Ew ald 125 126 Amongst PME is more popular a choice at present. Figure 2 7. Schematic representation of Ewald sum mation. 2.5 Free Energy Calculations in Molecular Dynamics As one of the most important property in thermodynamics, free energy calculation is free energy in molecular simulations due to the insufficient sampling of certain phase spaces. In macromolecular of enzymatic systems, many energy minimums exist, separated by different energy barriers. MD (or MC) is usually unable to sample those barrier areas, therefor e making it impossible to obtain an accurate free energy that requires sufficient sampling through the entire phase space In particular, people are interested in the potential of mean force (PMF), which represents the free energy surface along certain int ernal coordinates. Such an internal coordinate is also
60 known as reaction coordinate (RC). PMF provides valuable information with regard to the rea ction. F or example, the location of the free energy barrier on the surface usually corresponds to the transiti on state while the barrier height can be used to calculate the rate constant, an important kinetic property. Various sampling te chniques have been developed for such purposes. Umbrella sampling (US) and steered MD (SMD) are two of them that are heavily imp lemented in studies discussed in this dissertation. 2.5.1 Steered Molecular Dynamics The basic idea of SMD is that by putting an addition al force onto the system, the simulation will be able to overcome the free energy barrier an d sample the entire phase s pace The force applied is usually in the form of a harmonic restraint (equation 2 47). (2 47) Here, k denotes the restraint force constant, while x is the RC that could be distance, angle, dihedral or more complexed linear c ombination of coordinates. By integrating the force over RC, the work, W is obtained, however, it is from a non equilibrium process Jar zynsky relationship (equation 2 48) can convert this generalized work to a n equilibrium betwee n the initial and final state. (2 48) Many factors have been proved to play important roles in successful SMD simulations. It is particularly important to notice that if applied solo, several parallel SMD simulations must be carried out using different starting structures obtained from an equilibrium simulation 127 128 In is vital to SMD. F or instance, fast
61 moving could result in higher barrier height than reality SMD has been a perfect way to set up subsequent umbrella samp lings throughout all the studies discussed in this dissertation. 2.5.2 Umbrella Sampling A nother popular choice to construct the PMF is US P The implementation of USP also lies on an external potential term that in general in quadratic form, (2 49) H owever, unlike SMD, the external potential W here functions to bias the potential energy ble to be sampled sufficiently. The biased potential energy can be expressed as equation 2 50. (2 50) In practice, the step, an appropriate force constant k w is chosen for each window to guarantee enough sampling of phase spaces through the RC covered by t he window. The k w for adjacent windows should be able to provide sufficient overlap that is required for collecting distribution information The obtained distribution is certainly a non Boltzmann distribution. Next, in the third step, the unbiased PMF as a Boltzmann average, can be extracted from the non Boltzmann distribution by using equation 2 51 129 (Figure 2 8) (2 51) The advantage of USP over SMD is that the obtained PMF obtained from a series USP improper SMD moving speed. However, the success of USP strongly relies on the choice of
62 force constant k w for each window: a too small k w usually is not able to bias the energy surface enough for sufficient sampling over all phase spaces, while on the other hand, a too large k w converge (Figure 2 9). Figure 2 8. Schematic representation of building PMF with umbrella samplings. Figure 2 9 The effect of k w on umbrella sampling (US).
63 The weighted histogra m analysis method (WHAM) 130 is routinely used to extract the non Boltzmann distribution information from USP followed by constructing the unbiased PMF. The WHAM equation i s, (2 52) Both equations in equation 2 52 must be solved self consistently, which is achieved by iterative ly solving them in turn until both are satisfied. 2.6 Simulating Chemical Reactions and QM/MM The implementation of USP or SMD in MM MD simulations enable s computational chemists to study the free energy associated with certain transition processes, such as the trans and gauche conformational transition of butane 131 132 However, the limitation of MM force field impedes the free energy studies associated with chemical reactions. While on the other hand, the large size of molec ular system always keeps pure QM calculations away from simulating such scenarios. As a result, hybrid QM/MM method arises as the promising candidate to handle the modelling of chemical reactions. 2.6.1 QM/MM In general, QM/MM simulation divides the system into two parts: a small part that is of great important to research interest will be modeled using QM methods, while the large part of the system (most times including the solvent) will be modeled using MM force field. Thus the Hamiltonian can be rewritt en as: (2 53)
64 H ere the first two parts on the right side apply to QM region and MM region separately while the third part, H QM/MM is in charge of the interactions between QM particles and MM nuclei. Since H MM is governed by po tential energy function of chosen force field only but not sensitive to QM particles at all, the energy of t he system can be represented as: (2 54) The H QM is controlled by chosen QM theory, while the H QM/MM is the most complic ated part, as demonstrated in equation 2 55. (2 55) In equation 2 55 q is the charge while the subscript i j e refer to QM nuclei, MM nuclei and QM electrons, respectively. When partition the system into QM and MM regions, it is often necessary to cut a covalent bond (although one can also cut polarized bond, it is not recommended). In such scenarios, the H QM/MM becomes even more complicated and requires special treatments. Approaches to tackle this problem include link atoms, fr ozen orbitals and others, amongst the link atoms is most closely related to this dissertation. In implementation, a hydrogen atom will be placed on the boundary QM atom at default 1.09 al ong the QM atom MM atom vector, as depicted in Figure 2 10. For QM in Figure 2 10) will have zero charge in order to avoid double counting the electrostatics. 1 2 and 1 3 VDW interactions between QM and MM atoms are set to zero that means only non bonded VDW interactions are considered, in addition, 1 4 interaction is s ometimes rescaled to reproduce physical properties. Bond stretching, angle
65 bending and bond torsions with at least one atom from each side included are calculated with MM force field. Figure 2 10. QM and MM regions are defined by cutting the covalent bo nd between (MM) and C (QM). The hydrogen in pink is the capping link atom added to C along QM atom MM atom vector. In QM/MM simulation, when the PBC is applied, the number of interactions between QM and MM will be infinite. To handle such a case, a PME approach similar to that of MM has been developed. While for non periodic system, such a worry is not necessary, but a large QM cutoff is always recommended. 2.6.2 Semiempirical Theory in QM/MM The implementation of QM/MM in molecular simulations can provide great er accuracy at chemically interested area with moderate computational cost However, even as of today, the application of ab initio level QM/MM simulation is still very expensive, thus prohibiting this approach to be extended to large molecular systems. Fo r instance, the popular density functional theory (DFT) B3LYP coupled with MM potential has been rarely reported to apply on biological systems. Under such circum stances, semiempirical QM theories have become popular choice s for QM/MM simulation of macromo lecules and proteins.
66 Comparing to ab initio theories, apparently semiempirical method ologie s have made several approximations. The first approximation is that the two electron part is not explicitly included in the Hamiltonian Instead, Hckel theory is e electrons while extended Hckel theory is applied on all other valence electrons where the core electrons are considered ignorable The second approximation made is that the development of each semiempirical method includes some param eterization s to reproduce experimental data, thus these theories partially corrects the HF energy by including some electron c orrelation effects. As a result, the semiempirical methods usually have speed advantage over those ab initio theories and can prov ide comparable accuracy to those molecular systems that are similar to which the method is parameterized against. However, it is important to bear in mind that semiempirical methods could also lead to very wrong results on those systems that are not in the dataset used for parameterization. There have been various semi empirical methodologies developed. Amongst AM1 and PM3 are perhaps the most popular methods at present and they have been tested through a large variety of calculations. However, the transfer ability and the lack of d orbital are the major shortcomings of these two theories. 2.6.3 Semiempirical DFT SCC DFTB Recently, a new semi empirical theory called self consistent charge density func tional tight binding (SCC DFTB) based on DFT has been dev eloped and shown promising future. Starting by picking up a reference density 0 the energy can be expressed as (2 56) Here H KS is the Kohn Sham (KS) Hamiltonian,
67 (2 57) 0 is considered as t he sum of neutral atomic densities, In order to expedite the computation, the i can be expanded as By applying this basis to the first term of equation 2 56, it can be rewritten as, (2 58) Here is the KS eigenvalue of neutral, unperturbed atom, v eff is the effective KS potential consist of Coulombic and exchange correlation parts, while A and B define two neutral atoms. The three center terms are neglected, thus the off diagonal elements are evaluated by the two center terms. The Hamiltonian matrix element H KS and overlap matrix element are calculated with respect to interatomic distance between A and B and tabulated, thus making DFTB calculation faster. However, with above assumption, the charge density of a system is expressed as the sum of the charge density of its component atoms, which is in contradict to the reality that polarization should occur. To address this issue, the contribution regar ding to the second order density variation should be included (2 59) And (2 60)
68 In equation 2 60, q A is the Mulliken charge while q 0 A is the number of valence electrons of neutral atom A Equation 2 59 can also be written as. (2 61) AB can be expressed by either Hubbard parameter U A or chemical hardness A With equation 2 59 or 2 61 included in the total energy term, the approach is now self consistent in charge. With the repulsive potential, E rep is defined as, (2 62) Thus, the total energy of SCC DFTB can be rewritten as, (2 63) It is worthy noting that strictly speaking SCC DFTB is not a semi empirical method because th e parameterization includes no fitting to empirical data as other semi empirical theories do. In addition, SCC DFTB employs non orthogonal basis set which is the key for the improved transferability 133 Also, it is important to include an explicit dispersion term when apply SCC DFTB, since it has been well known that most LDA and GGA classes of DFT functionals have a poor record on dispersion force. With p arameters for more and more atoms developed, S CC DFTB has been widely applied to numerous molecular modelling and sampling studies and given good results. 2.6.4 Transition State Theory Transition state theory (TST) was developed to help people qualitativel y and quantitatively understand the process of chemical reactions. In order to understand TST it is important to make clear several concepts first.
69 Consider the reaction describe d in Figure 2 11, to illustrate the reaction path of a single set of molecule s involved, the PES is a good choice. For such a S N 2 type reaction, the PES has a similar shape to the free energy surface (FES): the two minimums correspond to the reactant and product sep arately, while the saddle point implies the transition state struct ure (TSS) The paths that connect the TSS downward to the reactant or product are called the minimum energy path (MEP) or intrinsic reaction coordinate (IRC) However, in macroscopic world, the equilibria and kinetics of reacting systems are consisted of c opies of single set of reacting system, and the rules are governed by free energy instead of potential energy. The free energy path is not necessarily the MEP, while the point associated with the highest free energy on the FES, the transition state (TS), m ay not exactly be the TSS on the PES. Figure 2 11. Reaction coordinate diagram for the biomolecular nucleophilic substitution (S N 2) reaction between bromomethane and the hydroxide anion. (Source: http://en.wikipedia.org/wiki/File:Rxn_coordinate_diagram_ 5.PNG )
70 TST invokes a so called quasi equilibrium assumption that considers a quasi equilibrium is reached between the TS and the reactant even the reactants and products are not equilibrium with each other. Hence, the quasi equilibrium constant K for the reaction in Figure 2 11 is (2 64) Therefore the rate equation of product formation can be expressed as (2 65) While the rate constant of product formation can be written as (2 66) Statistical mechanics provides us a temperature dependent relationship between the equilibrium constant and free energy, which is (2 67) Since k is proportional to the vibrational frequency associated with converting TS complex to product, rate constant of product formation, k can be rewritten as (2 68) H ere denotes the transmission coefficient referring to how much vibration contributes to the formation of product. TST has several major drawbacks, such as it does not account for quantum tunneling and it fails at high temperatures. There are several variants of TST, for example, variational transition state theory (VTST) i s a popular method that employs a variationally optimizing procedure to reach more accurate rate constant.
71 2.6.5 Kinetic Isotope Effect The rate constant of product formation under certain c ircumstances can be altered. One such example is the quantum tunneling effect mentioned above, which allows the system to tunnel through the barrier even when the system does not have an energy higher than the saddle point. Hence, such an effect can increa se the reaction rate constant. Another example is the kinetic isotope effect (KIE), which refers to the difference in the observed reaction rate constants of two differently isotopically labeled reactants. There are two types of KIE. When the isotopical su bstitution is involved in a bond that will brake or form in the rate limiting step (RLS) it is referring to as primary KIE (PKIE); while in all other situation, the effect is denoted as secondary KIE (SKIE). Based on the location of substitution, SKIE can SK IE, which means the isotope connects one S KIE, which refers to isotope connecting to the neighbor atom to the reacting center and so on. On a PES, the electronic energies of stationary points are in depen dent of atomic mass, thus the difference of zero point energies (ZPE) associated with two isotopically labeled reactions is the origin of KIE. For PKIE, the ZPE is dominant by the zero point vibrational energy (ZPVE). This is easy to understand becaus e in the case of either a bond is broken or formed, there is no such a vibrational mode in the TSS. If KIE is mathematically defined as the ratio of k light and k heavy as shown in equation 2 69 (assume it is a first order reaction), (2 69) T hus, the PKIE can be expressed as, (2 70)
72 The proportionality symbol can be replaced by an equal sign if the bond is fully broken or formed in the TSS. In the case of SKIE, however, there is no such vibrationa l mode disappearing associated with bond broken for formation, thus the ZPVE difference of the TSS also plays an important role in determining the k light / k heavy making the situation more complicated. Numerically, the PKIE measurements in general are signi ficantly larger than that of SKIE, while in SKIE, the measurements for S N 1 reactions are usually larger than S N 2 reactions. The former is easy to understand because of the contribution from ZPVE difference of TSS cancelling the ZPVE difference of reactants while the latter is somewhat tricky. Consider an S N 1 reaction of 3 chloro 3 methylhexane and iodide ion, as illustrated in Figure 2 12. The RLS is the formation of the carbocation intermediate This process is associated with the hybridization of chlorid e connecting carbon change from sp 3 to sp 2 which results in a variety of vibrational modes affected, especially the in plane and ou t of plane bending. The typical SKIE measurement of an sp 3 to sp 2 hybridization change is 1.1 to 1.2, while that of an sp 2 t o sp 3 change is typically 0.8 to 0.9. For S N 2 reactions, such like bromomethane and hydroxide yielding ethanol and bromide as shown in Figure 2 11 and Figure 2 13, features a TSS with both the nucleophile (hydroxide) and exiphile (bromomethane) being prese nt. Thus, the energies with regard to the old bond breaking and the new bond forming compensate each other, maintaining the ZPE difference throughout the reaction process. Depending on how well the nucleophile attacking and the exiphile retreating are sync hronized, the magnitude of SKIE varies slightly, in general between 0.95 and 1.06. If the hydroxide entry and bromide removal take place exactly simultaneously, the measured SKIE could be close to a unity value. However, it is important to notice that a le ss than 1.07 of SKIE measurement should not be automatically linked to an S N 2 reaction. For instance,
73 in enzymatic reactions, a unity SKIE value could be implying that the RLS is other than the chemical step. Figure 2 12. Reaction mechanism of S N 1 reac tion between S 3 chloro 3 methylhexane and iodide ion that yields both R and S 3 iodo 3 methylhexane. (Source: http://en.wikipedia.org/wiki/File:SN1Stereochemistry.png ) KIE measurement provides useful insights of determining the reaction mechanism, espec ially for enzyme catalyzed reactions. Several computational packages have been developed to enable computational chemists to reproduce experimental KIE values. Figure 2 13. Reaction mechanism of S N 2 reaction between bromomethane and hydroxide that yield s ethanol and bromide. (Source: https://secure.wikimedia.org/wikipedia/en/wiki/File:BromoethaneSN2reaction small.png )
74 CHAPTER 3 APPLY MOLECULAR DYNA MICS SIMULATION TO P REDI C T THE BINDING OF MAGNESIUM CATION IN FANESYLTRANSFERASE A CTIVE SITE 3.1 Backgrou nd FTase catalyzed the posttranslational attachment of a 15 carbon farnesyl group to a cysteine residue at a conserved CaaX motif at or near the carboxyl terminus (C terminus) of certain biological systems, enabling them to be recognized by the inner side of the plasma membr ane. For many small G proteins including the Ras superfamily of enzymes, such a modification is essential to their attachment to the cell membrane prior to play their important roles in the signal transduction 134 138 Mutations in Ras enzymes are responsible for nearly 30% of human cancers Therefore, the potential of inhibiting mutant Ras activity makes FTase catalyzed farnesyla tion a popular research area. FTase is a metelloenzyme that requires a zinc ion for activity and a second magnesium ion for optimal reactivity. The location of zinc subunit near the subunits interface: the metal ion forms a tetrahedral coordination with three while either a water molecule (before target peptide or enzyme binding) or a cysteine residue (Cys1p) from the substrate biological systems makes up the fourth coordination. Results from pKa experiments reveal that Cys1p interacts with zinc via a form of thiolate (Zn S ) instead thiol (Zn SH). In fact, a weak zinc thiolate interaction is believed to a key to the catalysis cycle as such a bond can enhance the nucleophilicity of the bound thiolate 139 In contrast to its important role in chemical step, the zinc ion plays little role in the FPP binding, which is required for the later peptide binding. Thus, the substrate binding scheme for FTa se can be described as follows: zinc ion first bi nds into FTase with coordinating to 3 AA side chains and a water molecule, then FPP binds into the high ly positively charged pocket surrounded by
75 ion if pre sent, finally the peptide or enzyme substrate localize s itself with the Cys1p replacing the water molecule to fill up Zn 2+ tetrahedral coordination. However, the resolved ternary crystal structure of FTase FPP CaaX complex shows an approximate 7 distance between the thiolate and another reacting center, C 1 atom of FPP. (The active site snapshot is given in Figure 3 1.) In order to bring two reacting atoms close enough into reacting range, a conformational rearrangement features two substrate approaching e ach other is necessary. Unfortunat ely, experiments could not give further support to reproduce this procedure. Figure 3 1. Snapshot of FPP binding pocket. Two reacting centers are labeled as C1 (C 1 of FPP) and SG (S of Cys1p). Magnesium cannot restore the full reactivity of FTase without zinc. However, its existence enhance s the rate of product formation by up to 700 fold In terms of Gibbs free energy, the bound of Mg 2+ causes approximate 4 kcal/mol difference. Add itionally recent study reveals that (1) Mg 2+ phosphate group and (2) such a coordination can also contribute to the physical step that makes required conformation transition. Howev er, unlike many other Mg 2+ metalloenzymes, a rich aspartate
76 DDXXD motif is not present in FTase, making the localization of this ion somewhat difficult. In catalysis, a millimolar of Mg 2+ is enough to make the huge impact on reaction rate enhancement, resu lting in the lack of this ion in any of currently available FTase structures. Computational effort has already been made to analyze the FTase structure and illustrate its conformational transition mechanism 140 141 With the Mg 2+ ion excluded, the results show ed the monoprotonated form of FPP (FPPH 2 ) is cap able of making such a transition while the fully deprotonated FPP (FPP 3 ) needs a very high energetic cost thus is unlikely able to undergo this procedure Due to the limitation of classical MD simulation, this PMF study was carried out to move the RC ( d C1 ) from 8 to 4.5 . Such a range is sufficient to cover the transition from the resting state and also able to avoid drastic VDW clashes associated with closely placing two atoms In the light of this success, a computational modelling of Mg 2+ in FTa se active site and a series of consequent validations have been conducted and will be presented here. Also included in this chapter is the modelling and simulation of GGTase I. GGTase I and FTase share a great deal of similarities : t heir s are highly homologues, with approximately The active sites of two zinc metalloenzymes are particularly similar. In FTase, FPP possesses a binding pocket surrounded (Figure 3 2) the binding pocket of GGPP is circled by In fact, these residues are highly conserved in CaaX prenyltransferases An interesting fact is that GGTase I shows no Mg 2+ dependency for reactivity, whereas FTase requires this second metal ion for optimal activity. Experimental study of GGTase I 2+ in FTase. Consider also this place in FTase is the proposed most important Mg 2+ binding ligand, this fact becomes more
77 interesting since in such highly related enzymes this change may relate to important functional 2+ dependency, i ndicating that Mg 2+ binding site could be restored with an aspartic acid replacing the lysine. Nevertheless, in the crystal structure of GGTase I complex, this lysine is over 5 away from the nearest diphosphate oxygen, making its role in the catalysis pr ocess unclear. It is also unclear whether GGTase I binds a monoprotonated GGPP (GGPPH 2 ), as FTase does in the absence of Mg 2+ or a fully deprotonated GGPP (GGPP 3 ) due to the existence also of great interest to perform simulation with GGTase I and make comparison with FTase. Figure 3 2. Active site snapshot of GGTase I. C 1 of GGPP and S of Cys1p are labeled as C1 and SG.
78 3.2 Methods and Theories 3.2.1 Bu ilding Initial Model Two crystallographic structures of FTase complex from protein data bank (PDB), PDB code 1QBQ and 1FT2, are superimposed to give an active FTase complex composed of FTase, FPP and acetyl CVIM peptide 140 Crystal water molecules from PDB code 1QBQ are completely al protonation states under physiological pH are adopted. Histidine residues are all modeled in the neutral form and Cys1p is mode led as thiolate, as suggested based on experimental observations 142 Zinc ion is coordinated with Cys1 p, The enzyme and peptide are modeled with AMBER force field ff99SB 89 while the zinc site is modeled using the parameters developed by Cui et al. based on a bonded approach (Table 3 0 =0.00125, developed by Merz et al. are adopted 143 FPP is modeled in its fully deprotonated state (FPP 3 ). FPPH 2 is excluded from this study since the protonation of diphosphate unit leads to a loss of magnesium binding affinity 144 Force field parameters applied to FPP are obtained from generic AMBER force field (GAFF) 145 for organic molecules, the torsion parame ters for the farnesyl group developed by Merz and coworkers and the diphosphate parameters developed by Car l son lab are adopted 146 In order to obtain the atomi c charges for FPP 3 a geometry optimization of FPP 3 is first performed using B3LYP/6 31G* level of theory, then the electrostatic potential derived atomic charges are computed using HF/6 31G* level of theory, finally the atomic charges are derived using a two stage restrained electrostatic potential (RESP) fitting procedure 147 148 Within such a procedure, bot h QM calculations are conducted with Gaussian 03 149 while the fitting process is performed using ANTECHAMBER integrated in the AMBER 10 suite of programs 150 AMBER 10 is also chosen to carry out all MD simulations.
79 3.2.2 Modelling and Simulating Mg 2+ in FTase Figure 3 3 Snapsh ots of FTase active site from A ) initial structure with 10 Mg 2+ placement loci; B) WT1 model; C) WT2 model; D) WT3 model; E) MUT model and F ) SWAP model. In each model, Zn 2+ and Mg 2+ are represented as grey sp here and green sphere separately. Mg 2+ coordination is labeled in bolded dashed line. Zn 2+ coordination is shown in light dashed line. Important hydrogen bond interactions are labeled in red dashed lines. Residue from and subunit are labeled with a and b separately. By carefully examining the FPP 3 binding pocket, three aspartate residues, namely magnesium binding. The latter two have been recognized as the most possible Mg 2+ binding ligands based on availa b l e experimental evidence 151 A total of 10 guess starting stru ctures are created (Figure 3 3 A ). In each of them, Mg 2+ is placed close to negatively charged diphosphate moiety and one o f three candidate aspartic acid residues. Neighbouring crystal water molecules make up the octahedral coordination of Mg 2+ which has been identified in many Mg 2+
80 0 =0.8947 152 is applied. Unli k e Zn 2+ which parameters are derived through a bonded approach, Mg 2+ is modeled as a f ree ion since success has been achieved with such a strategy in previous MD simulations of NphB 153 Each of new generated FTase/CVIM/FPP 3 /Mg 2+ complexes is so lvated in a truncated octahedron periodic water box. Eac h side of the box is at least 10 away from the closest solute atom. The TIP3P water model is employed and counter ions are added to neutralize the system. PME method 124 is applied t o compute long range electrost atic interactions with 8 and 1 set as cutoff of the nonbonded interactions in real space and reciprocal space grid spacing separately The translational movement of the system is removed. The motions of covalent bonds with at least one hydro gen atom in volved are constrained with SHAKE algorithm. Each system is first relaxed with 5000 steps of steepest descent minimization followed by 5000 steps conjugate gradient minimization in order to remove any close contacts. Within both stages of minimization, a w eak positional restraint is applied to the enzyme, peptide and zinc, while free movements of magnesium and FPP 3 are allowed. Then the system is slowly heating up to 300 K over 100 ps time period using Langevin dynamics with a 1 ps 1 collision frequency in canonical ensemble. A 1 fs time step is used in this stage to improve the stability of MD. All restraints are gradually removed over the same period of time. Finally, a 6 ns production MD simulation is carried out in an NPT ensemble with a 2 fs time step. For every picosecond, a snapshot is saved. The last 4 ns trajectory is used for analysis, validating the corresponding Mg 2+ binding pose. One of the Mg 2+ Mg site (WT1 model in Figure 3 3B), proposed by Fierke and coworkers 151 is artificially constructed by applied
81 strong NMR restraint in the minimization stages to forc e the Mg 2+ place of a Mg 2+ water coordination. The remaining stages of simulations are carried out using the same strategy as described previously. Table 3 1. Zinc parameters used in MD simulations 140 Atom /Bond/Angle Types Atomic Mass /Force Constant ( kcal/mol 2 / kcal/mol radian 2 ) At omic Charges (e) /Equilibrium Value (/degree) ZN 65.37 00 0.6928 S1 32.06 00 0.6797 S2 32.06 00 0.6666 ZN O2 180.94 00 1.90 00 ZN NB 94.82 00 2.24 00 ZN S1 167.6 0 00 2.21 00 ZN S2 154.53 00 2.48 00 CT S1 237 .00 00 1.81 00 CT S2 237 00 00 1.81 00 ZN O2 C 88.03 00 101.53 00 ZN S1 CT 80.88 00 111.28 00 ZN S2 CT 77.34 00 124.24 00 ZN NB CV 42.29 00 112.69 00 ZN NB CR 41.83 00 133.68 00 O2 ZN S1 12.57 00 99.95 00 O2 ZN S2 10.54 00 90.73 00 S1 ZN NB 34.48 00 131.59 00 S2 ZN NB 25.15 00 97.17 00 S1 Z N S2 11.84 00 111.29 00 An extend MD simulation in addition to the 6 ns MD simulation has been performed later with the most successful model to evaluate the FPP 3 interactions and Mg 2+ coordinati on With such a purpose, canonical ensemble that allows more flexible conformational change is adopted,
82 while Langevin dynamics with a 1 ps 1 collision frequency is again selected. This simulation lasts 10 ns. The time step is 2 fs. 3.2 Based on experimental evidences, the mutation of As in the presence of Mg 2+ would result in significant loss in reactivity model in Figure 3 3 E ) is prepared based on the final snapshot of a successful Mg 2+ binding motif. All counter ions a long with water molecules that are more than 6 away from FPP 3 or Mg 2+ residue. Such a change leave the only Mg 2+ binding candidate. The resulting system is then simulat ed following the same strategy as described in Chapter 3.3.2. 3.2 .4 Modeling Removal of Mg 2+ from Active Site A good tactic to validate the Mg 2+ binding motif is to test whether the removal of Mg 2+ from FPP 3 binding site will restore the behavior of FTase /CVIM/FPP 3 complex. To achieve this goal, the most straightforward approach is to switch the position of Mg 2+ and a solute water molecule far away from the biding pocket (SWAP model in Figure 3 3 F ) The advantage of this approach over val of Mg 2+ from the system is the avoidance of rebuilding the system (topology) thus any a ssociated artificial impacts will also be avoided As a result, the optimal comparability between SWAP model and original Mg 2+ binding model is reached In practice, the coordinates and velocities of bound Mg 2+ and a water oxygen atom are exchanged. The involved water is referred to as H 2 O swap The coordinates of two hydrogen atoms of H 2 O swap are modified accordingly and the velocities are maintained. The remaining po rtions of the system are totally preserved. This approach also saved much effort such like re minimization and re heating up. However, a longer equilibration (4 ns) is required to allow the system to fully relax itself. This is followed by a 4 ns productio n MD under NPT ensemble.
83 3.2 .5 PMF Validation In order to further validate the binding motifs, PMF for each reasonable binding motif is calculated as a function of C 1 S which is referred to as RC. To compare to previous results from Cui and Merz 141 the same coverage of RC from 8.0 to 4.5 , is chosen in this study. A MD scanning process is performed first to generate a set of initial structures for US. 15 such structures are prepared with an interval of 0.25 between adjacent two windows. A harmonic potential of 20 kcal/mol 2 is applied to each window throughout 2 ns simulation with a 2 fs time step. Data are collected every 5 steps for the latter 1 ns simulation, resulting in 100,000 data per window. More windows, longer samplings, and larger force constant (100 kcal/mol 2 ) are 154 is utilized to unbias the distribution of RC a nd reconstruct the FES. 3.2.6 Modelling and Simulation of GGTase I PDB ID 1N4S represents a heterohexamer of GGTase I complexed with GGPP and a LCVIL peptide 20 A monomer is select from the original PDB file and prepared as the initial structure. Two models are generated, adopting GGPP 3 and GGPPH 2 separately. Each of the resulting model s is solvated into a truncated octahedron water box. TIP3P water model is adopt ed. AMBER force field ff99SB is applied to model the protein while GAFF is employed to model the isoprenoid substrate. The charge of GGPP is derived following the two stage RESP fitting procedure, as described in Chapter 3.2.1. Long range electrostatic pot ential is handled by PME method and SHAKE is applied to constraint hydrogen involve bond motion. The resulted structures are minimized carefully with a weak harmonic positional restraint applied. After that, each of the systems is heating to 300 K over 10 0 ps time span with 1 fs time step. Langevin dynamic s is employed as thermostat. Positional restraint is gradually removed
84 during this process. Then each system is equilibrated in NPT ensemble for 900 ps with 1 fs time step, followed by a 4 ns production r un. 3.3 Results and Discussion 3.3.1 Modeling and Simulation of FTase/CVIM/FPP 3 /Mg 2+ Complexes A set of 10 possible Mg 2+ binding motifs are generated by manually placing magnesium ion at different loci in the FPP 3+ binding pocket. All binding schemes are test ed through MD simulations, resulting in three possible Mg 2+ binding models, namely WT1, WT2 and WT3, respectively ( Figure 3 3 B,C,D 2+ coordination at least once in three models. 2+ is found to interacting with a pair of oxygen atom s diphosphate, a single carboxylic oxygen atom from water molecules. In WT2, Mg 2+ didate aspartates subunit of FPP 3 along with three oxygen atoms from three water molecul es make up the octahedral coordination with Mg 2+ 2+ coordination. In such a 2+ An diphosphate oxygen atom four water molecules fi ll up the remaining of the octahedral coordination. Each of these three binding motifs appear s rational upon visual inspection. Octahedral coordination is observed in all three models, in agreement the natural preference of such a cation in many biologica l systems. Each motif includes at least one Mg 2+ Asp interaction with an
85 exception of WT1 model which has two such interactions. Despite the absence of an aspartate rich DDXDX motif identified in many Mg 2+ metalloenzymes, a solvent accessible FPP 3 binding pockets allows multiple water molecules to make up the coordination. In fact, not only crystal water molecules, but also those from solvent box, are identified in the Mg 2+ coordination in three possible binding schemes. Recently, local ly enhanced samplin g (LES) method has become popular in determining such binding locations 155 157 However, we are not abl e to achieve success with such an approach. A rational explanation is that Mg 2+ carries a fairly strong positive charge thus it is rather difficult to move it around. All three binding schemes are validated through MD simulations. Stable octahedral coordi nation of Mg 2+ is observed in each of these simulations. Hence, at least one other important characteristic is required for further validation. Such an example is the distance of C 1 S It is experimentally determined that FTase complexes possess an over 7 distanc e between two reacting centers and require a conformational transition prior to the chemical step. Under such circumstances, d RC binding motifs, because it represents the level of difficulty to per form such a transition. The d RC is measured throughout MD simulations and an average value is calculated for each of three binding poses in order for comparison. In WT1 and WT2 models, su ch an average distance is below 6.5 . Comparing to the 7.2 of d RC measured in FTase/CVIM/FPPH 2 system from previous studies 140 141 it is reasonable to draw a conclusion that both WT1 and WT2 models represent more favorable structures than the same FTase and peptide complex with Mg 2+ in absence, since less distance in RC is required to overcome to reach the reacti ng range. On the other hand, however, in contrast to WT1 and WT2 models, the average distance of C 1 S is found to be over 10 for WT3
86 model. Such a distance represents an approximately 50% increasing in length from the FTase/CVIM/FPPH 2 complex, indicating WT3 is unlikely an active complex In addition, Mg 2+ coordination in WT3 features a Mg 2+ subunit non bridging oxygen atom of FPP, although a similar Mg 2+ scheme has been identified in aromatic prenyltransferase NphB, subunit FPP oxygen atom in the octahedral coordination is contradict to the model propose d by Fierke and coworkers 151 Therefore further effort has been more focused on WT1 and WT2 models, though the possibility of WT3 binding scheme cannot be completely excluded. In nearly every crystallographically resolv found at least responsible for been found particularly important for the magnesium chemistry in this prenyltransferase. Mutagenesis experiments has revealed an approximate 27 fold decrease in catalytic activity associate fold increase in catalytic activity associated with Mg 2+ binding 151 Based on experimental observations, the assumption has been proposed th at Mg 2+ plays an active role not only in the chemical step but also in the physical step 151 With the chemical step remaining unclear, results from this study support the presumption that including Mg 2+ in the FPP 3 bindi ng pocket makes a subunit diphosphate oxygen atoms involve in the Mg 2+ coordination in both models, consisting with the model proposed by Fierke et al 151 Furthermore, the involvement of Mg 2+ in FPP 3 binding pocket provides extra stabilization to the big negative charge on diphosphate moiety in addition
87 to the hydrogen bond interactions between diphosphate oxygen atoms and side chains of of chemical step b y stabilizing the developing 4 units of negative charge of the leaving group. In summary, it is reasonable to regard the Mg 2+ ed in both WT1 and WT2 models as our best estimate of the Mg 2+ binding motifs in the FTase/CVIM/FPP 3 /Mg 2+ complex. Similarities have been discovered in both models. The active site conformations from both WT1 and WT2 share great similitude. The key hydrog en bonds between diphosphate moiety and surrounding functionally important AA residues are well preserved, and a key hydrogen bond is identified between S of Cys1p and a Mg 2+ bound s of is found to move into FP P 3 binding pocket, localizing itself within hydrogen bond both cases (as well as one Mg 2+ bound water in WT2 model). Such a movement brings an extra stabilization t o the Mg 2+ binding site. Another metal site, the tetrahedral coordination of Zn 2+ maintains fairly well in both models. The bound of Mg 2+ introduces strong electrostatic interactions with diphosphate moiety, disrupting the salt bridge interactions between the temporarily However, in the extend MD simulations with WT2 model, these two salt bridges have been restored ( Figure 3 5 ) An estimation of approximate 5 kcal/mol 2 free energy penalty associated with breakin g these interactions has been made by Cui and Merz 141 while experimental results revealed a binding free energy of appro ximate 3 4 kcal/mol provided by these interactions. Thus, it is necessary to identify the exact roles these two positively charged residues played in the chemical step.
88 e subunit of diphosphate. In previous FTase/CVIM/FPPH 2 and FTase/CVLM/FPPH 2 studies, such a hydrogen bond was clearly identified. However, in the 6 ns MD simulation of both models, this hydr ogen bond is not with the acetyl capping group. Another positively charged residue that acts as a hydrogen bond aintains a stable salt bridge interaction with a non subunit oxygen atom throughout entire simulation for each model. It is interesting to note that these two residues change roles from previously FTase/CVIM/FPPH 2 subunit of diphosphate moiety. With such a roles exchange, the total contribution from these two residues to the system is maintained, as each of the peptide and the subunit of di phosphate moiety are stabilized by one of two residues. Comparing WT1 and WT2, obviously the main difference is the function of zinc bound 2+ binding ligand, which resembles the Mg 2+ binding scheme proposed by Fier ke and coworkers, while in WT2 model it functioning as only a zinc bound ligand. In WT1 model, the two metal centers are 3.3 apart, and the equilibrium d RC is 5.8 , which suggests this scheme might represent the energetically more favored FTase/CVIM com plex than either FTase/CVIM/FPP 2 or FTase/CVIM/FPP 3 studied previously. With two 2+ active site less flexible. Such an assumption is supported by the evidence that d RC only fluctuates in a fairly small range in WT1 simulation. In this sense, despite the short equilibrium d RC the
89 bridged bimetallic site may result in increased difficulty in performing required conformational transition to make two reacting atoms further approaching each other. In WT2 model, there i s no such an aspartic acid bridge. The two metal ions are around 5.1 separated in space. A water molecule takes th coordination of Mg 2+ The equilibrium distance of C 1 S is fluctuating around 6.1 , which is 0.3 longer than that found in WT1 model. A magnesium bound water molecule has been found to form a stable hydrogen bond with Cys1p in the peptide, resembling a substitute bridge two metal binding sites is less than that in WT1 model thus a more flexible active site in WT2 model can be expected. As expected, a larger fluctuation in d RC is observed in t he simulation of WT2 ( Figure 3 4 ), making it easier for two reacting centers to further approach each other and reach a pre organized state prior to the chemical step. In cont rast to the model proposed by Fierke and coworkers, either WT1 or WT2 model does not adopt a double Mg 2+ donates a single carboxylate oxygen atom to Mg 2+ coordination. Attempts were made to te st the possibility of Mg 2+ interacting with both carb a stable structure cannot be obtained. In the 10 ns extend ed MD simulations carried out with WT2 model, the interactions of su bunit of diphosphate moiety are restored ( Figure 3 5 ). The FPP 3 are all back into range to form hydrogen bond interactions with diphosphate. The distance between Zn 2+ and Mg 2+ is greatly relaxed to approximate 7.6 , allowing higher degree of
90 back d RC is maintained as 6.1 , in support of Mg 2+ playing active role in conforma tional transition step. Figure 3 4 Distribution of d RC in the original 6 ns MD simulation of WT2 model. Figure 3 5 Snapshot of active site from the extend 10 ns MD simulation of WT2 model.
91 2+ /H 2 O Swap Simulation An approxim ate 30 fold difference in the rate of product formation has been observed MD simulations are also conducted to study the impact of this mutation. The initial structure is taken from and modified based on a snapshot from the simulation of WT2 model. Through the n ew series of simulations, a substitute Mg 2+ binding motif is observed and to that in WT1 model. However, different to WT1, a water molecule picks up the role of is observed, marginally longer than that observed in WT1 but significantly shorter than that in sequence conserved. Several important hydrogen bond interactions observed in WT2 are well preserved in this the capping acetyl group Water molecules play important roles in stabilizing both substrates. diphosphate moiety, while a couple of Mg 2+ bound water molecules interacting with Cys1p via hydrogen bonds diphosphate, although it is still within a range to form hydrogen bond with back forming a salt away from the closest oxygen atom of FPP 3
92 The equilibrium d RC is approximate 6.5 , slightly longer than that of both WT1 and WT2 model. However, it is still clearly shorter than that observed in FTase/CVIM/FPPH 2 complex in the absence of Mg 2+ The value of d RC somehow can be interpreted as a measure of difficulty level to undergo a conformational transition required by the system prior to the chemical step. I n this sense, the predicted order of difficulty is: FTase/ CVIM/FPPH 2 > MUT > WT1/WT2. Nevertheless, this trend is not obtained based on free energy, thus it is necessary to carry out PMF studies to determine the real trend. Another series of MD simulatio ns are carried out to study the SWAP model, in which the Mg 2+ has been switched with an outside water molecule. Such a switch mimics the effect of removal Mg 2+ from the enzyme complex, resulting in the system relaxing itself. Ideally, such a relaxation sho uld give the same conformation of FTase/CVIM/FPP 3 without Mg 2+ from previous study. Through the MD simulation, a similar FPP 3 binding pocket as that observed in the resting state of FTase/CVIM/FPP 3 complex is obtained. More importantly, key interactions are all restored as expected. In t the hydrogen bond range to the nearest diphosphate oxygen. Asp away from the binding pocket accompanied with the disappearance of Mg 2+ Also moved away is Cys1p, leaving the equilibrium d RC back to approximate 7.4 , well represent its status without Mg 2+ binding. 3.3.3 PMF Studies Anot her good approach to validate these Mg 2+ binding models is the PMF study. However, it is not a straightforward task to study free energy in MD simulations due to the difficulty to sample certain phase spaces, in particular those associated with high energi es, along the RC. In order to achieve this goal, special techniques must be implemented in order to obtain distribution
93 information over all phase spaces required by free energy calculation. In this study, USP method is adopted. For each model, a set of at least 15 simulations is carried out with a harmonic bias potential placed between 4.5 and 8.0 . For WT1 model, in order to achieve better overlap between adjacent windows, a much stronger bias potential is employed and trajectories are propagated at 29 p oints along the same range of RC with an interval of 0.125 . The most important observation is that the inclusion of Mg 2+ clearly alters the free energy profile ( Figure 3 6 ) For instance, the global minimum, which represents the resting state, shifts fr om 7.4 in FTase/CVIM/FPPH 2 complex to 5.8 in WT1 model. The 1.8 decrease in the equilibrium d RC implies WT1 model energetically more favored than FTase/CVIM/FPPH 2 system, since less distance is required to bear a conformational transition. However, the free energy profile suggests that it might not be a correct Mg 2+ binding model. In fact, the free energy curve of WT1 model represents steep increases on both side of the global minimum, in accordance with the small fluctuations of d RC observed in th e MD simulation. The system, as a result, seems locked into the conformation associated with this resting state. A reasonable explanation is that the bridged bimetallic Zn 2+ Mg 2+ site is too rigid for the two substrates to get closer. Another possi bility is that the zinc parameters based on a bonded approach are not good enough. In fact, the bond parameters to describe Zn 2+ length. However, a set of newer fas hioned parameters built with MTK++ 158 that is based on bonded approach and Seminario method gives a 51.96 kcal/mol force constant and a 2.01 equilibrium bond l ength. Obviously, the latter set of parameters represents more flexibility. However, this discussion is not relevant to this study thus will not be further expanded here.
94 Figure 3 6 Free energy profiles of the modeled systems. WT2 model, on the other h and, shows a global minimum at 6.1 along the RC and a second minimum at around 5.0 corresponding to an intermediate state during conformation transition. Although there is an about 0.3 increase in equilibrium d RC the free energy curve appears more re asonable than WT1 model. The highest barrier of the transition is approximate 0.8 kcal/mol at around 5.6 , and the intermediate on the other hand is about 0.5 kcal/mol higher in free energy than the resting state. Furthermore, the curve on the left side o f the global minimum is fairly flat, indicating it is easy for the system to undergo a conformational transition to close the gap between two reacting atoms and bring them into reacting range. During the transition, FPP 3+ approaches the peptide substrate v ia the rotation of its first and second isoprene unit ( Figure 3 7 ), which is also in agreement with the isoprene rotation mechanism. Comparing to the free energy profile of FTase/CVIM/FPPH 2 complex that has a global minimum at around 7.4 and an 0.5 kcal /mol transition intermediate at around 5.0 separated by a barrier of about 1.0 kcal/mol at around 5.8 , clearly WT2 possesses a much shorter d RC at equilibrium and a
95 lower transition barrier. In addition, the free energy at 4.5 is lower than 1.0 kcal/ mol for WT2 model, which is nearly only a half of that observed for FTase/CVIM/FPPH 2 complex. However, due to the strong VDW clashes possibly at distances shorter than 4.5 , a QM/MM level PMF study is required to further validate the last point. Figu re 3 7 Comparison of FPP 3 orientation found in WT2 MD simulation and from previous FTase/CVIM/FPPH 2 (without Mg 2+ ) MD simulation. The free energy profile computed for MUT model possesses a global minimum at around 6.6 , approximately 0.5 longer than that before mutation. No second minimum has been found for this mutant, while the free energy penalty to bring two reacting centers from 6.6 to 4.5 is approximate 6 kcal/mol. From mutagenesis study, such a mutation results in a loss of reactivity by ne arly 30 fold. In terms of free energy, the difference is approximate 2 kcal/mol. Apparently this value is much lower than the 5 kcal/mol difference between WT2 and MUT found from PMF study. However, by breaking down the profile of MUT into two parts, clear ly a large portion of the free energy penalty comes with the 4.5 5.5 part, while from 5.5 to the resting state the transition requires only less than 1.0 kcal/mol free energy. This interesting fact
96 once again lifts the discussion on the quality of cu rrently employed zinc parameters and the limitation of MM force field at the edge of potential VDW clashes. Nevertheless, it appears that further determine its im pact in the following chemical step. The PMF of SWAP model finds great similarity with the profile obtained from previously studied FTase/CVIM/FPP 3 complex in the absence of Mg 2+ Both models find the global minimum at around 7.5 , and both systems requi re an approximate 3 kcal/mol free energy penalty to close the big gap between two reacting centers. The excellent agreement between these two models further confirms the validity of our Mg 2+ binding models complexed with FTase, in particular the WT2 model. 3.3.4 Modelling and Simu lation of GGTase I The independency of GGTase I reactivity on Mg 2+ is interesting if one considers the great homologous shared by GGTase I and FTase. In the absence of Mg 2+ the most important question is whether there is an AA side chain can play a substi tute role. Experimental results reveal that GGTase I 21 However, from the simulations carried out with both GGTase/CVIL/GGPP 3 and GGTase/CVIL/GGPPH 2 it is revealed that this residue will approach diphosphate after around 5 ns of simulation and forms salt bridge interaction to stabilize that moiety to asso ciates with the decrease of d C1 (see Figure 3 8). This correlation is particularly evident in GGTase/CVIL/GGPPH 2 complex. Another question of the same importance is whether GGTase I adopts a monoprotonated GGPPH 2 like the FTase does in the absence of Mg 2+ or it possesses a fully deprotonated GGPP 3 This question can not be answered based on current results from regular MD simulations at present A lthough diphosphate may prefer a more negatively charg ed diphosphate moiety, the possibility of
97 GGPPH 2 cannot yet be totally excluded Hence, t he best solution to these two questions is the PMF study, which is still on going currently. Figure 3 8. Distances variation in the simulation of GGTase/CVIL/GGPP 3 (top) and GGTase/CVIL/GGPPH 2 (bottom). X axis represents frames. Y axis represents while red line represents d C1 the diphosphate. On the other hand MD simulations provide useful information to evaluate important ligand substrate interactions. GGPP shares very similar binding pocket with FPP in FTase, the bonds be
98 simulations for both models, making it impossible to determine a more favorable model based on eracting diphosphate via hydrogen bond interaction after the conformational transition. Thus it is reasonable to assume diphosphate accompanied by system bearing conformational transition. The 8.2 equilibrium d RC represents more challenge for the conformation transition preformed by GGTase I than that conducted by FTase. A reasonable explanation of this increasing in d RC lies on GGT ase I possessing a binding pocket that locks the terminal isoprene unit of GGPP thus making the rotation of this substrate to peptide cysteine more difficult Kinetic experiments reveal that the k cat for GGTase I is (0.18 0.02) s 1 corresponding to 18 .6 kcal/mol in Gibbs free energy. This value is lower than FTase complexes without Mg 2+ but higher than the FTase/Mg 2+ complexes. It is yet clear to say whether the difference attributes to the longer es. Furthermore GGTase I can also accept FPP and catalyze the farnesylation. A QM/MM study is need ed to elucidate the nature of GGTase I chemistry. 3.4 Summary Mg 2+ binding schemes are modeled and validated through regular MD simulations and PMF studies w ith umbrella sampling. Among the three possible binding motifs, only WT2 is 2+ as proposed based on mutagenesis results, however, another proposed Mg 2+ bound ligand, 2+ binding. The Zn 2+ Mg 2+ bridged bimetallic site appears incapable of lowering the free energy barrier associated with conformational transition, thus making itself an unreasonable choice for Mg 2+ coord ination. A
99 hydrogen bond formed between Cys1p and Mg 2+ bound water molecule is detected and proposed to contribute to reduce the d RC tetrahedral Zn 2+ site and octahedral Mg 2+ site. However, difficulty in closing the distance gap between two reacting atoms is found in this model, although the bridge seem s increasing the connection between two binding sites thus shortened the d RC USP method is employed to study the PMF of both models. WT2 possesses a flat free energy curve with less than 1.0 kcal/mol difference between the resting state and 4.5 However, an approximate 6 kcal/mol free energy penalty is observed for the MUT system, making it much more difficult for FPP 3 to approach Cys1p in this case. Introducing Mg 2+ into the FPP 3 binding pocket enables the fully deprotonated FPP to the zinc bound peptide cysteine residue via isoprene rotation mechanism. It also disturbs the active site structure. For example, diphospha at least temporarily. The loss of interactions between negatively charged diphosphate and these two positively charged residues is compensated by Mg 2+ itself via coordinating to both diphosphate oxygen atom. However, with longer relaxation time, the salt bridges between Mg 2+ and Zn 2+ After removing Mg 2+ from the active site, the system relax es itself back to a configuration represents the model of FTase/CVIM/FPP 3 without Mg 2+ The equilibrium d RC is also reverted to 7.4 . Such a complete restoration strengthens our conclusion that Mg 2+ binding is the solo reason for the conformational chang e.
100 GGTase I, a highl y homologous system to FTase, has also been modeled and simulated. 2+ affinity and at least partially plays a substituting role in GGTase I chemistry. However, during MD simulations of GGTase/CVIL complexed with both GGPP 3 and GGPPH 2 the distance between this lysine and its nearest diphosphate oxygen is more than 5 . Thus it either reaches in during the conformation transition or stabilizes the diphosphate indirectly v ia interaction with water molecule in between. Both protonation states give stable MD simulations, thus PMF study seems the only way to identify the substrate preference of GGTase I.
101 CHAPTER 4 SYSTEMATICALLY MECHANISTIC STUDY OF FTASE COMPLEXED W ITH DI FFERENT PEPTIDE SUBSTRATES 4.1 Background Conformational transition step is required by FTase complexes to close the big distance gap between two reacting atoms. However, this step does not account for the highest free energy barrier in FTase chemistry. Ki netic studies reveal that the rate limiting step for FTase catalysis is the product release step 159 160 With single turnover kinetics, it is feasible to measure the free energy barrier for the product formation that covers both the conformational transition and the chemical reaction step 139 161 It is understandable the chemical step draws more at tention, sinc e reaction is the real essence of chemistry. In addition, as an important pharmaceutical target enzyme, more ef ficient FTase inhibitors are always of great interest to researchers 23 162 176 In drug discovery, more than often the enzyme inhibitor is designed to adopt or mimic the pose of transition state or intermediate state of enzyme catalyzed reactions. In this sense, it is of great importance to study the chemical step of FTase catalyzed farnesylation. Many usef ul information and interesting facts have been uncovered by experimental approaches. Long, Casey and Beese have resolved the crystal structures of FTase complex at several key steps of its catalysis process By analyzing the active site configuration of reactant complex and product complex, they proposed a transition state configuration which resembles a S N 2 reaction TS 18 This assumption is suppo rted by various experimental evidence such as the guration in FTase catalysis 177 and the nearly unity measurement SKIE by Weller and Distefano 178 However, not only S N 2 characteristics were found for FTase chemistry. Harris, Derdowski and Poulter have proposed a zinc dissociation mechanism featured the break of Cys1p an d Zn 2+ coordination before the chemical step 179 Nevertheless, this is inconsistent to the product complex struct ure resolved by Long, Casey and Beese 18 thus
102 making itself irrelevant to the discussion Also, s everal experimental groups have proposed a carbocation invoked dissociative mechanism for farnesylation 139 180 181 To explain the existence of both nucleophilic and electrophilic characters identified in FTase chemistry 139 144 182 184 another mechanism has been proposed. This new m echanism possesses an S N 2 like TS with dissociative characters. More evidence in support of this combined mechanism came from a computational approach conducted by Klein and coworkers. QM/MM MD simulations coupled with PMF study was employed and as a resul t a clear S N 2 like TS with S N 1 like dissociative characteristics was identified 185 Unfortunately, the energetic parameters obtained from their simulation are fairly off track, partially attributed to the use of FPP 3 instead of FPPH 2 as isoprenoid substrate in the absence of Mg 2+ In addition, the QM theory they adopted in the study, BLYP, has been proved problematic with transition metals especially zinc 186 Kinetic studies under single turnover condition provide valuable information regarding to the rate of product formation. One such example is the approximate up to 700 fo ld difference in rate of product formation observed in the presence or absence of Mg 2+ FTase complexed with Dns GCVLS peptide gives a k chem of 8.1 0.3 s 1 at saturating concentration of MgCl 2 and a k chem of 0.025 0.001 s 1 in the absence of MgCl 2 151 Such k chem values correspond to 16.3 kcal/mol and 19.8 kcal/mol Gibbs free energy with Mg 2+ present or absent, respectively. This kind of energetic information, coupled with TST, could efficiently guide theoretical studie s to illustrate the nature of FTase catalysis. Some recent in teresting experiment results mad e FTase catalytic mechanism more fascinating. Fierke and coworkers carried out 3 SKIE experiments to examine the reaction mechanism for FTase complexed with GC VLS and TKCVIF peptide separately. The measurement for TKCVIF peptide reads 1.154 0.004, while the GCVLS peptide gives a k H / k T
103 ratio of 1.00 0.04 184 Based on SKIE measurement for analogues enzymatic system, the first value falls well into the so called associative mechanism wit h dissociative characteristics 187 H owever, the second value, which is nearly unity is not straightforward to interpret. Fierke and coworkers proposed that such a value ind icates the switch of the RLS from the chemical step to the conformational transition step. Later, the same group has also systematically analyzed the recognition of different Ca 1 a 2 X motifs by FTase 31 and proposed an interesting context dependent a 2 and X selectivity of the substrate recognition. The impact of different a 2 and X could occur either at or before the chemical step. Co mbining these two studies a very i nteresting scheme with regard to the nature of FTase chemistry is given Despite the large amount of studies carried out on FTase, it is unfortunate to say that the stud y of farnesylation mechanism, in particular the transition state configuration, has not gone far away from proposals and hypotheses on paper. This is partially resulted from the limitation of current expe rimental techniques. However computational study at QM/MM level can provide valuable insight into the catalytic mechanism of this importan t enzyme. 4.2 Methods and Simulation Details 4.2.1 FTase/CVIM/FPPH 2 FTase/CVIM/FPPH 2 system has been well studied in our lab 141 This complex will be referred to as FPPH model in this chapter. It is reasonable to adopt one of the configurations extracted from the production simulation performed previously as the initial structure for this study. Because the system has been well equilibrated, any more classical simulation with this model is redundant. The structure is partitioned into two parts: QM region and MM r egion. In order to limit the QM region to a moderate size, only those residues, ligands, or p arts of them that functionally important to the chemical reaction are included. In this sense, the diphosphate group, the first
104 coordinating ligands, along with Zn 2+ itself are selected as QM atoms (Figure 4 1) All the side chains are cut at the covalent bon d between C and C thus there is no polarized bond cut. Figure 4 2+ and diphosphate and the first isoprene unit of FPPH. Note only labeled carbon and their hydrogen atoms in t he farnesyl group are considered in QM. Due to the existence of diphosphate, in particular the phosphor, traditional semiempirical methods such as PM3 and AM1 cannot handle the simulation properly 188 Therefore, SCC DFTB 189 190 becomes the best, if not the only, choice The system is first minimized using steepest descent method for 1000 steps to fully relax the system with hybrid SCC DFTB/ff99SB potential 191 A 100 ps Langev in dynamic with a 0.5 fs time step is employed to heat the system temperature back to 300 K in the canonical ensemble. Followed is an equilibration of 500 ps in NPT ensemble with a 1 fs time step. A cutoff of 8 is selected for both QM and MM portion of t he system. SHAKE is applied to MM part to constrain
105 the motion of bond with hydrogen involved. Long range electrostatic interaction s in MM region are computed using PME method with the default settings. In QM region, SHAKE is turned off so that the free mo tion of hydrogen involved covalent bond is allowed, and a modified PME implementation designed for QM region is employed. A couple of SMD simulations with a 5000 kcal/mol A force constant are conducted to propagate the trajectory along the RC to both 1.8 (product state) and 8.0 (reactant state). The RC is still defined as the d C1 Alt hough at present a RC in the form of d RC = d C1 S d C1 O1 is more popular, our pre liminary results showed this selection will lead to a very unreasonable dissociative pathway, in accordance with the claim given by Klein and coworkers 185 63 struct ures are extracted with an interval of 0.1 along the RC in order for the subsequent US. A set of 63 windows of USP is carried out to obtain the free energy profile. Each window is sampled in NTP ensemble for 150 ps wherein the first 50 ps is for equilibr ation purpose while the latter 100 ps is used for date collection. A harmonic force constant of 10 0 kcal/mol is applied to bias all the windows An approximate 10 ns total simulation is performed for US. employed to unbias the d istribution and construct the PMF. resulting system denoted as and treated in the same way as FPPH model, except the QM/MM study for this model does not cover the 4.5 to 8.0 range that has already been studied classically. As a result, 28 windows, instead of 63, are simulated with USP method. 4.2.2 FTase/CVLM/FPPH 2 FTase/CVLM/FPPH 2 system (denoted as CVLM model in this chapter) is another complex that has been well studied by Cui and Merz 141 The original purpose of this study is to testify the possibility of a switch of RLS from chemical step to physical step, as observed in
106 TKCVIF and GCVLS complexes. The computational results show an increase of conformational transition barrier, which is, however, too small to cause a RLS switch. Considering the context dependent selectivity of both a 2 and X position in the Ca 1 a 2 X motif, it is not surprising to see this result, since only one position of the motif is changed in the simulation thu s it is insufficient to reproduce the impact of double change studied in the experiment. On the other hand, a complete mechanism of both the chemical and physical step will help to understand the substrate dependency at a 2 position of this four unit peptid e and provide useful information to illustrate the more complicated context dependent substrate selectivity. In this sense, a QM/MM PMF study is carried out to energetically study the chemical step. In the structure preparation step, isoleucine in the pep tide is modified to leucine within the FPPH model. Thereafter all the procedures and treatments are exactly identical to that applied on 300F model. 4.2.3 FTase/CVLS/FPPH 2 The most interesting mutational effect associated with FTase chemistry seems related to the FTase/CVLS complex. In nature, GCVLS motif represents the C terminus of human Ras enzyme, making this complex more attractive. A crystal structure (PDB code: 1TN8) 22 obtained at 2.25 resolution is adopted to prepare the initial structure. This structure repr esents FTase enzyme complexed with GCVLS peptide. In order to focus our study on the context dependent of a 2 and X position, the glycine is modified to an acetyl capping group in order to be consistent with our previous study The result acetyl Ca 1 a 2 X mot if is consistent with all the models investigated in this chapter. The modified system is then solvated in an octahedral box with a distance of 8 to separate the edge of box and closest solute atom. The resulting model is denoted as CVLS in this study.
107 C VLS model is first treated classically. A cutoff of 8 is applied, while SHAKE and PME method are employed to improve the simulation. A series of minimizations with a combination of steepest descent method and conjugated gradient method is performed in th e first stage to remove possible close contacts in the system. A weak positional restraint is applied to all of the non water components. In the next stage, the system is slowly heat ing to 300 K with Langevin dynamics in NVT ensemble. During this state the restraint potential gradually removed. Then the model is allowed to fully relax in the NPT ensemble followed with a 10 ns production run. A snapshot extracted from the last 1 ns of production simulation is chosen as the starting structure for the followi ng MM PMF study and QM/MM simulations. The MM PMF study adopts the same strategy as described for WT2 model in Chapter 3.2.5. The QM region and MM region are selected following the same procedure described before. The subsequent treatment is identical to t hat applied to FPPH model. 4.2.4 QM/MM Validation of FTase/CVIM/FPP 3 /Mg 2+ Model Mg 2+ plays an important catalytic role in the FTase chemistry. Based on the available evidence this divalent metal cation does not cause the approximate 3.5 kcal/mol differe nce in free energy during the conformational transition. Thus, it is rational to suppose the difference i n rate of product formation only comes from the chemical step. An approach is made to study the chemical step of our WT2 model from chapter 3 with QM/ MM potential. Due to the lack of Zn Mg parameter, Mg coordinating ligands with the exception of FPP 3 are not modeled with SCC DFTB. This leaves the QM selection of WT2 model identical to that of FPPH except one less hydrogen from diphosphate group The s tarting structure of WT2 model is extracted from the trajectory of the extra 10 ns relaxation instead of the original 6 ns simulation. This is a reasonable move since it is found that during the extra 10 ns, the system relaxed itself more deeply. Thus the strain associated with the
108 inclusion of Mg 2+ will be minimal. The Mg 2+ and Zn 2+ are more than 6 apart, adding less extra trouble to the troublesome SCC DFTB convergence procedure. The WT2 model is pa rti ti oned and treated with QM/MM potential in the same way as applied to FPPH and other systems. 4.3 Results and Discussions 4.3.1 A Glimpse of Entire Product Formation Process of FTase/CVIM/FPPH 2 In FTase complex, two reacting atoms, C 1 of FPP and S of Cys1p are more than 7 away in the resting state 18 In order to comple te the catalysis, both a physical step that closes the distance gap and a chemical step that attaches the farnesyl moiety to the target cysteine are required. Thus, a comprehensive understanding of the catalytic mechanism requires a PMF study to cover both steps. To accomplish this goal, we carried out a couple of SMD simulations to force the RC moving from the resting state to the product state and reactant state, respectively. The obtained free energy profile reveals a resting state at approximate 7.0 , an intermediate state at about 5.3 , a second intermediate state at approximate 3.9 that is very likely represent a pre organized state of the complex, and a TS at approximate 2.6 as well as a product state at 1.8 . Although the barrier predicted by SMD simulation is frequently less accurate than that computed with USP it usually an S N 2 like TS is found, however, with a unusual long C S bond observed that characteristics. 184.108.40.206 Transition state and free energy b arrier With the structures extracted from SMD trajectories, a set of 63 windows USP is then conducted in order to refine the results. The generated free energy profile ( Figure 4 2 ) is simila r to that gathered from SMD simulations The resting state is located at 6.90 , while the first
109 intermediate state is identified at 5.32 , separated by an approximate 1.0 kcal/mol barrier. These results are in excellent agreement with the MM PMF results obtained by Cui and Merz 141 In the region that has not been covered in their previous study, a second intermediate is fou nd at 3.90 , associated with an approximate 6.5 kcal/mol free energy. It is worthy noting that at this point, O 1 and C 1 of FPPH 2 and S of Cys1p start to align linearly. Such evidence can be well explained by the pre organized state theory. The TS is located at 2.63 and the barrier height is approximate 20.6 kcal/mol. The experimental measurements are 21.1 kcal/mol for TKCVIF and 20.0 kcal/mol for GCVLS 184 Obviously, our result perfectly falls into this range. The product state is detected at 1.82 , associated with a free energy of approximate 4.5 kcal/mol higher than the resting state. Upon reaching the transition state (see Figure 4 3) the C 1 O 1 bond starts to break and the distance between them starts to explode. The S of Cys1p, C 1 of FPPH 2 and O 1 of FPPH 2 still stay linearly, however, with C 1 slightly out of plane. The H 1 H 2 C 1 C 2 (see Figure 4 1 for atom specification) dihedral is giving the highest value at this point, reading 168 that is very close to the ideal value of 180 for a typical S N 2 like TS. Moreover, such a dihedral changes sign after this point, meaning the two hydrogens bend out of plane in a different direction. The distance of C 1 S is 2.63 , approximately 0.2 longer than d C S observed in a S N 2 reaction calculation carried out at MP2/6 31+G**//MP2/6 31+G* level 192 The difference in TS C S bond length mi ght be a good indicator of the dissociate characteristic. Also, the distance of Zn S elongates to a bout 2.45 from its equilibrium value 2.38 , implying a weaker coordination of zinc. After carefully examining the structure of each window along the RC, it is conspicuous that the distance gap between two reactants is closed by FPP rotating its first and second
110 binding pocket play important roles in the catalysis pro diphosphate is stabilized via salt bridges formed between its non protonated oxygen Figure 4 2 Free energy profile of FPPH model obtained from a set of 63 windows of umbrella sampling Active site snapshots are given at key points. From right to left, (1) resting state; (2) first transition intermediate; (3) second transit ion intermediate; (4) transition state and (5) product state. Zinc coordinates the peptide cysteine throughout the entire process. An increase of approximate ly 0.3 of d Zn S from the FTase reactant complex (2.4 ) to the product complex (2.7 ) has been observed though comparing the crystallographic structures 18 Such a d ifference results in a weaker Zn S coordination, which has been proposed to increase the nucleophilicity
111 at S thus improving the reactivity. From our simulations, a change of about 0.2 has been detected. Furthermore, the Zn S starts elongating when the C 1 O 1 bond is about to break. Su ch observations seems in support to the hypothesis of weak zinc sulfur interaction improve the nucleophilicity at S 139 Figure 4 3. Sna pshot of the active site of FPPH at the transition state. O 1 and C 1 of FPPH, and S of Cys1p are labeled in red. 220.127.116.11 Com parison of QM/MM and MM PMF of conformational transition s tep The difference between our QM/MM PMF results and MM PMF results of Cui and Merz in regard to the conformational transition step is fairly marginal, exp ect at close to 4.5 of RC where the potential of VDW clashes becomes significant (due to the sum of carbon and sulfur VDW radius is about 4 ). The little discrepancy at the resting state location, 6.9 of QM/MM
112 prediction and 7.2 of MM results, as we ll as at the intermediate state, 5.3 from QM/MM and 5.0 from MM, are acceptable, since QM and MM are totally different approaches. In fact, the energy barrier separating the two states predicted by QM/MM and MM matches fairly well. The biggest differe nce is the free energy value at 4.5 , which QM/MM gives around 2.7 kcal/mol while MM shows only 1.8 kcal/mol. As mentioned before, at this point, results from MM is suspicious due to the potential of VDW clashes. Whereas QM method does not have this limit ation, WT2, the QM/MM PMF studies do not cover the 4.5 to 8.0 region, instead, MM results are used to predict the total reaction barrier bas ed on the difference between corres ponding value at 4.5 and that of CVIM at 4.5 . This is a reasonable approach, since the MM error at this point for each model will very likely cancel. Hence, a great deal of computational cost is saved. 4.3.2 Effect of Mutations Mutagenesis study plays a very important role in the mechanistic study of enzymatic catalysis. For example, i t helps to determine the function of key AA residues Hence, it is of great importance to reproduce the results from such experiments in computational simulations. diphosphate in the chemical step of FTase catalysis. The k cat measured for the product formation changes from 0.016 0.008 s 1 of wild type (WT) FTase /GCVLS complex to 0.0002 0.00004 s 1 h in the absence of Mg 2+ ). A nearly 80 fold loss of activity accompanies by the mutation, which corresponds to 2.6 kcal/mol difference in Gibbs free energy 182 The mutation in the peptide substrate plays more important role. It not only alter s the free energy barrier height, but also appears to disturb the reaction mechanism. One good example is SKIE experiment of FTase/TKCVIF complex and FTase/GCVLS complex 184 The k H / k T measurement of the former system is 1.154 in the absence of Mg 2+ while for the latter is nearly
113 unity, 1.00 0.04, under same condition. Such a phenomenon has been explained as the RLS switches from the chemical step (for TKCVIF ) to the physical step (for GCVLS). However, consider the free energy barrier of FTase/GCVLS complex in the absence of magnesium ion is approximate 20.0 kcal/mol, it is rather unusual to see such high a energy barrier accompanies with the conformational tr ansition. utation resulting system is decently minimized and then simulated in the same manner as applied to CVIM model. However, PMF study has not been applied to the 4.5 8.0 range of RC due to the reason explained before. Figure 4 4. Free energy profile s of all five models included in this study.
114 The result free energy profile (see Figure 4 4) shows a TS at 2.61 associated with a free energ y of 22.5 kcal/mol (see Table 4 1). Comparing to the 20.6 kcal/mol free energy barrier of CVIM model covering the same range of RC, the 0.7 kcal/mol discrepancy between calculated difference and experiment observation is almost within the chemical accuracy (1.0 kcal/mol) to the experimental observed difference. Table 4 1. Free energy parameters obtained from QM/MM and MM PMF studies. TS CONF (kcal/mol) 1 4.5 MM (kcal/mol) 2 4.5 QMMM (kcal/mol) 3 TS QMMM (kcal/mol) 4 TS PROD (kcal/mol) 5 FTase/FP PH/CVIM 1.0 (1.0 6 ) 1.8 2.7 6 17.9 20.6 6 FTase/FPPH/CVLM 2.7 3.0 3.9 7 19.3 22.7 FTase/FPPH/CVLS 2.5 2.8 3.7 7 18.3 22.0 Y300F/FPPH/CVIM 1.5 2.3 3.2 7 19.3 22.5 FTase/FPP/Mg/CVIM 0.8 1.0 1.9 7 16.5 18.4 1 TS CONF refers to the free energy barrier associat ed with the conformational transition taken play in 4.5 2 4.5 MM refers to the free energy at 4.5 obtained from MM PMF calculations. 3 4.5 QMMM refers to the free energy at 4.5 obtained from QM/M M PMF calculations. 4 TS QMMM refers to the free energy barrier of the chemical step, obtained from QM/MM calculations in 4.5 8.0 region. 5 TS PROD refers to the free energy barrier of the entire product formation, covering reaction coordinate of 1.8 8.0. 6 These values are directly produced from QM/MM PMF in 1.8 8.0 region. 7 4.5 QMMM 4.5 MM 4.5 FTase/FPPH/CVIM where 4.5 FTase/FPPH/CVIM 4.5 QMMM 4.5 MM = 2.7 1.8 = 0.9 kcal/mol. A nalyzing the TS structure (see Figure 4 5), it is not surprise to see similarities to that of ibility of FPPH 2 and later developed diphosphate leaving group (PPi). The shorter length of the forming C 1 S bond at TS seems resulting from this extra flexibility, so does the longer distance between C 1 and O 1 lay a role in the catalysis by stabilizing the
115 farnesyl group via cation 1 C 2 or C 3 atom of farnesyl (where based on the dissociative characteristic the delocalization of positive charge should take plac transition from TS to the product state. However, more analysis such as higher level QM calculation is needed to validate this assumption. Figure 4 5 Snapshot of the active site of Y300F mutant at the transition state. O 1 and C 1 of FPPH, and S of Cys1p are labeled in red. Note the position and orientation of mutant moiety via cation residue cause the increase of free energy barrier. The 1.9 kcal/mol difference between the WT and mutant is divided into two parts, approximate 30% (0.5 kcal/mol) in th e physical step and
116 70% (1.4 kcal/mol) in the chemical step. With experimental measurements of other enzyme residue mutation avaialbe, it will be interesting to find out what is the impact associated with those residue mutations. 18.104.22.168 Interesting c ont ext dependent e ffect The most interesting mutational effect of FTase complex might be the so called context dependent selectivity of the a 2 and X in the Ca 1 a 2 X motif 31 Such an effect does not only impact the free energy barrier height but also alters the reaction pathway 184 The most solid evidence SKIE ex periments of FTase/TKCVIF and FTase/GCVLS complexes. The most straightforward effort based on our previous study would be to simply substitute CVIM with CVLM. CVLM model is prepared and studied with the same strategy as applied to The result free energy profile features a TS at 2.68 and a barrier height of 22.0 kcal/mol. The experimental value with regard to this motif in the absence of Mg 2+ is yet available, while a peptide with this motif in the presence of Mg 2+ gives a Gibbs free energy of 16.6 kcal/mol, 0.3 kcal/mol larger than that measured for the CVIM motif (16.3 kcal/mol) under identical experimental condition. Thus the difference of calculated free energy difference between CVLM and CVIM model, 1.4 kcal/mol (22.0 20.6 ), is again i n excellent agreement with the experimental evidence and almost reaches the chemical accuracy requirement. The TS active site configuration (see Figure 4 6) once again represents great similarity to the CVIM model. The biggest difference this time is the longer C 1 S bond. The 2.68 bond length corresponds to a 0.05 increase, possibly resulting from the increased hydrophobic interaction between the side chain of a 2 residue and the farnesyl group of FPPH 2 In addition, the mutation of Ile to Leu results in bulkier side chain terminal, thus possibly also increasing the steric effect of the residue to the surrounding enzyme AA side chains.
117 In spite of increased barrier height and d C1 S at the TS, the RLS remains unchanged. Hence, it seems a single change at a 2 posit ion is not sufficient to result in mechanism change. In order to further illustrate this fascinating context dependent effect, another mutant, CVLS model, is then prepared and studied. Figure 4 6. Snapshot of the active site of CVLM mutant at the transi tion state. O 1 and C 1 of FPPH, and S of Cys1p are labeled in red. 4. SKIE results of FTase/CVLS c omplex The double mutation at both a 2 and X position of the peptide substrate accompanies bigger change in peptide conformation, therefore simply modifying the existing CVIM model might not be as good as utilizing a crystal structure of FTase/GCVLS complex. The preparation and simulation detail has already been described in 4.2.3. The free energy curve of MM PMF study
118 clearly locates the resting state at approximate 7.65 (Figure 4 7), which is significantly longer than the WT of either CVIM or CVLM, and a transition intermediate state at around 5.0 . The free energy barrier associated with the transition is around 2.75 kcal/mol. Although this value is inarguably higher than either WT CVIM or WT CVLM, it is definitely unable to cause the RLS SKIE result observed on this system. An explanation of such observations might be the chemi cal reaction mechanism alters from an associative mechanism with dissociative characteristics to a N 2 mechanism. Figure 4 7. Free energy profile associated with the conformational transition step of CVLS model. The QM/MM PMF study re sults provide more support to this conclusion. The free energy profile represents a 2.50 of TS C 1 S distance and a 22. 7 kcal/mol free energy barrier. Comparing to CVIM model, markedly the d C1 S is 0.13 shorter and more close to the results
119 of S N 2 reaction computationally studied at MP2/6 31+G**//MP2/6 31+G* level 192 Moreover, the S of Cys1p, C 1 of FPPH 2 and O 1 of FPPH 2 are almost linearly aligned, despite C 1 is slightly out of the plane (see Figure 4 8). Unlike in CVIM model, where the H 1 H 2 C 1 C 2 dihedral remains fluctuating around 160 for a period of time corresponding to the RC change from 2.65 to 2.50 , the same dihedral in CVLS model reach es a peak at TS, which reads 178 .6 (see Figure 4 9) and changes rapidly on both sides of the peak. All of this structural information strongly pushes toward the conclusion of a more in dependent S N 2 reaction mechanism is adopted by FTase/CVLS complex. Figure 4 8 Active site snapshot of CVLS model at the transition state. O 1 and C 1 of FPPH, and S of Cys1p are labeled in red.
120 On the other hand, the experimentally observed free energy difference for FTase complexed with CVLS (16.8 kcal/mol) and CVIM (16.3 kcal/mol) in the presence of Mg 2+ is 0.5 kcal/mol, while the barrier for FTase/GCVLS without M g 2+ is about 20.0 kcal/mol, both of which are in excellent agreement with our results. Before reaching the TS, there is another intermediate state with regard to the conformational transition is identified at approximate 3.7 . Like in CVIM model, the thr ee atoms involved in bond formation and broken, namely O 1 of FPPH 2 C 1 of FPPH 2 and S of Cys1p, appear to align linearly upon this point, suggesting this intermediate state can be referred to as the pre organized state. Figure 4 9. H 1 H 2 C 1 C 2 dihedral variation in the farnesylation of FTase/FPPH/CVIM (red) and FTase/FPPH/CVLS (green) complexes. X axis of the figure represents number of frames, Y axis of the figure represents the dihedral in degree.
121 Comparison between CVLS and CVIM reveals structural information of interactions between peptide substrate and enzyme active site residues. In CVIM, t he sulfur of methionine backbone of serine forms water mediate hydrogen bond inte observations are supported by crystal structures of FTase/KKKSKTKCVIM complex (PDB 1D8D) 193 and FTase/CVLS complex (PDB 1TN8). Thus, despite less bul ky in terminal residue size, the peptide of CVLS model shows the same level of stabilization as CVIM model. Figure 4 10 Ramachandran plot of a 2 residue of the Ca 1 a 2 X motif throughout the farnesylation of FTase/FPPH/CVIM (red), FTase/FPPH/CVLS (blue) and FTase/FPPH/CVLM (green) complexes.
122 The Ramachandran plots 194 of the peptide a 2 residue of CVIM, CVLM and CVLS model during the chemical step (see Figure 4 10) provide more interesting facts about the context dependent a 2 X selectivity. In CVIM or CVLM model, either the isoleucine or the leucine stays fluctuating in the transition area conn implies tha t the peptide CVLS might sacrifice stability to help two reactants approaching. Furthermore, the closest contact between leucine side chain and farnesyl moiety in CVLM model is around 3.5 while nearly remains over 3.9 in CVLS. Such a difference might r esult in an extra flexibility for FPPH 2 to rotate its farnesyl group due to the less steric stain from peptide leucine side chain. SKIE measurement of FTase/GCVLS complex, calculations to address the ZPE difference at TS and resting state of both CVLS and CVIM model are undergoing, since ZPE difference dominates the free energy differen ce resulted from isotope substitution. Moreover, what role that the residue before Ca 1 a 2 X plays in the catalysis and what impact the length of the peptide does to the reaction pathway are of great interests to us and needs to be clarified in future studies 4.3.3 Mg 2+ in Chemical Step In Chapter 3, a reasonable Mg 2+ binding model of FTase/CVIM complex has been proposed. Although results from MM level simulations match experimental observations fairly well, it is necessary to carry out QM/MM level studies, i n particular the PMF study, to further testify this binding scheme and illustrate how this metal ion alters the free energy barrier. The preparation and simulation strategy is identical to that applied to CVIM model. The original simulation fails because of the difficulty of convergence possibly resulting from the complex electrostatic interactions associated with the close distance between two metal ions. However, freezing water molecules at least 10 away from the active site ligands seems
123 improving the electrostatic interaction calculation and thus solving the convergence p roblems. This approximation is not unusual to see in many other successful QM/MM simulations. Thus, the restraint has been applied throughout all simulations with regard to this model The QM/MM PMF study of WT2 model gives a free energy profi le featuring an approximate 18.3 kcal/mol at approximate 2.69 along the RC. It is interesting to notice the TS d C1 S increases about 0.06 from an identical complex without Mg 2+ present. Consider the residue, it seems the inclusion of Mg 2+ coordination introduces extra degree of inflexibility as well as extra degree of stabilization to the diphosphate moiety. The TS active site configuration looks somewhat different to that of the CVIM TS ( Figure 4 11 ). The major difference is of course the Mg 2+ binding throughout the reaction process. Mg 2+ coordinates to two diphosphate oxygen atoms from each subunit, a carboxylate oxygen atom of apping acetyl group in MM PMF study, however, in this study it appears to form a hydrogen bond with the peptide cysteine stabilizing the 4 charged PPi leavin g group via water mediated hydrogen bond interaction. Mg 2+ and Zn 2+ are more than 7 apart throughout the process, however, the two metal sites are correlated via a The computed free energy barrier is 18.4 kcal/mol. It is approximately 2.2 kcal/mol lower than the barrier height of CVIM model but about 2 .1 kcal/mol higher than the experimental
124 measurement. This discrepancy is significantly bigger than any other model studied in this series. A rational explanation to this increased discrepancy might be the implementation of freezing water 10 away from th e substrates. However, unless an improved set of SCC DFTB Mg 2+ parameters including the most important Zn Mg interaction parameters is introduced, to simulate this system with SCC DFTB/AMBER level QM/MM theory would remain challenging thus the freeing outs ide water approach would still be necessary. Other approaches, such as ab initio level of QM/MM simulation, should also be able to contribute to the mechanistic study of FTase chemistry in the presence of Mg 2+ nevertheless, the computational cost almost p rohibits the implementation of such methods coupled with USP approach. Figure 4 11 The active site snapshot of WT2 at the transition state. C 1 and O 1 of FPP, and S of Cys1p are labeled in red.
125 4.4 Summary In this series of studies the chemical step of FTase catalyzed farnesylation has been systematically studied. The resulting TS active configuration of CVIM model with Mg 2+ absent is in perfect agreement with expe rimental evidence both structurally and energetically. In particular, the 20.6 kcal/mol free energy barrier is impressively falling into the range of 20.0 kcal/mol for FTase/GCVLS complex and 21.1 kcal/mol for FTase/TKCVIF complex with zero concentration of MgCl 2 On the other hand, the 2.63 TS d C1 is significantly longer than the computed TS d C S of a S N 2 reaction in gas phase at MP2/6 31+G**//MP2/6 31+G* level, indicating this reaction adopts an associative mechanism with dissociative characteristics. In addition, the QM/MM PMF study of the c onformational transition step finds its results and the results from MM PMF study match excellently. and CVLS mo del in the absence of Mg 2+ all give excellent free energy pred ictions. Additionally, the difference in TS distance of C 1 S varies among these models. Although the exact reason remains puzzle, the trend seems related to the level of flexibility possessed by FPPH 2 and/or the peptide substrate. The recognition of pept ide residue as a 2 and X position in the Ca 1 a 2 X motif exhibits fascinating context dependent selectivity. Amongst, the FTase/CVLS model is particularly important since it seems making impact not only on the free energy barrier height, but also on disturbing SKIE experiments reveal a nearly unity measurement for FTase/GCVLS complex while a value of 1.154 for FTase/TKCVIF complex. The latter obviously falls into the range of an associative mechanism with dissociative characteristi cs while the former has been explained as the RLS switches from the chemical step to the physical step. However, our results lead to an alternative interpretation of the nearly unity 3 SKIE result of CVLS model that the chemical step remains the RLS bu t the mechanism
126 N 2) mechanism. This assumption is SKIE calculations at MP2/6 31+G** level. The Ramachandran plot of a 2 residue shows great difference between CVL S and CVIM/CVLM stably remains in the region. Finally, the WT2 model obtained fro m studies described in Chapter 3 is validated through QM/MM PMF simulations. The result free energy barrier feat ures a free energy barrier of 18.4 kcal/mol corresponding a TS d C1 S e discovered to play different roles in the presence of Mg 2+ than in the absence of such a metal ion. With an approximation of freezing water molecules 10 outside any substrates applied, the resulted free energy barrier shows the biggest discrepancy amon g all the models discussed in this chapter. Therefore, an improved simulation of this model is still required.
127 CHAPTER 5 COMPUTATIONAL APPROA CH TO UNDERSTAND ARO MATIC PRENYLTRANSFERASE CA TALYTIC MECHANISMS 5.1 Background Prenyltransferases include not o nly three members of protein prenyltransferases that are famous for their important roles in G proteins regulated signal transduction pathway but also aromatic prenyltransferase such as bacterial prenyltransferase NphB 3 ( or Orf2 before rename) purified from S treptomyces sp. Strain CL190 and fungal indole prenyltransferase FtmPT1 4 extracted from Aspergillus fumigatus This kind of aromatic prenyltransferase makes great co ntribution to the biosynthesis of terpenoid in addition to a variety of terpenoid synthase 195 The products of such a prenylation, the terpenoids, are very widely distributed in a large variety of natural products, such as in fungi, animals and plants. Amongst these terpenoids products, plant isoprenoids have been well known for their application into traditional herbal remedies because of their aromatic qua lities 6 Some isoprenoids and their derivatives have been revealed to possess certain pharmaceutically desired characteristics such as antiviral, antioxidant as well as anticancer effect 196 205 In addition, the low cellular toxicity and excellent ability to penetrate cell membrane associated with this kind of compounds have made them desired drug templates The biosynthesis of these terpenoids has been of great interest and importance to biochemistry and medical chemistry researches 195 NphB, identified recently from Streptomyces possesses a rather interesting 3 D structure barrel. Inside the barrel a spacious bi nding pocket is provided to two reacting substrates, the isoprenoid substrate GPP and the aromatic substrate 1,6 DHN GPP is stabilized via interaction between negatively charged diphosphate and several AA side chains, including Lys119, Lys169, Arg228, Tyr216 and Lys284, in addition to a Mg 2+
128 diphosphate ox ygen atom and four water molecules, in spite of the missing of a DDXDX motif usually found in Mg 2+ metalloenzymes. Asp110 situates within 3.5 to a Mg 2+ bound water molecule, provide extra stabilization to the metal site. NphB exhibits interesting substr ate selectivity and product regioselectivity 3 153 206 In the catalysis of geranylation of 1,6 DHN, at least 3 products have been identified. The major product, a trans 5 geranyl 1,6 DHN, and the minor product, a 2 geranyl 1,6 DHN, are originally characterized, wh ile recently a extra minor product has been purified and identified as a 4 geranyl 1,6 DHN. Previous computational study has revealed that the different orientation of 1,6 DHN cause s the product diversity. The major product is defined as S1 state, the mino r product is defined as S3 state, while an intermediate state referred to as S2 state connects them in the FES. Unfortunately, no binding orientation has been found responsible for the extra minor state. We will call this extra minor state the S4 state in this study. Indole prenyltransferase FtmPT1 has been only recently identified from Aspergillus fumigatus This family of indole prenyltransferases catalyzes the attachment of isoprenoid moiety to tryptophan and its derivatives. The product of FtmPT1 compl exed with Brevianamide F, the Tryprostatin B, shows high cytotoxicity to a number of cancer cells, making this enzyme a promising anticancer target 56 58 207 210 The 3 helices. However, this barrel is FtmPT1 prefers DMAPP as isoprenoids substrate. The active site consists of a highly positive charged binding pocket for DMAPP and a tyrosine rich environment. The diphosphate is anchored by the salt bridge
129 interactions between itself and two positively ch arged residues, namely Arg113 and Lys294. The side chains of Tyr203 and Tyr450 help stabilize the diphosphate via hydrogen bond interactions. These two tyrosine residues, along with Tyr296 and Tyr382, circle a ring DMAPP binding site provide stabilization to the dimethylallylic cation via cation onally, the reaction between DMAPP and brevianamide F very likely includes a rotation of the dimethylallylic cation due to the large separate between two reacting atoms. Therefore, the on. FtmPT1 also exhibits interesting regioselectivity. However, it is activated by a single mutation at Gly115. Mutating it to Thr115 results in the dimethylallylic cation connecting to the C 3 (C 3A ) instead of C 2 (C 2A ) position of the bre vianamide F via its own C 3 (C 3P ) atom instead of C 1 (C 1P ) atom. However, G115A mutant does not show the same regioselectivity. The reaction mechanism of NphB has been proposed to be an associative pathway featured an S N 2 like TS or a dissociative mechanism featured a car bocation invoked electrophilic capture. Although the dissociative mechanism is more preferred, the possibility of an associative mechanism cannot be totally excluded without further study. On the other hand, the catalysis of FtmPT1 has been assumed to adop t a dissociative pathway. Three steps mechanism has been proposed: in the first step the hydrolysis takes place to form a PPi leaving group and a carbocation, next in the second step the carbocation rotates then attaches either C 1P or C 3P (in G115T mutant) to the indole nucleus forming a new C C bond consequently the proton transfer occurs in the third step to restore aromatic quality of the indole substrate. KIE experiment has determined that both the hydrolysis and bond forming steps are partially rate l imiting Consider
130 the rate constant of product formation ( k cat ) is 5.57 s 1 under single turnover condition, the KIE result is reasonable since the time scale for such a rotation is usually not in this order of magnitude. 5.2 Methods and Simulation Detail s 5.2.1 QM/MM Study of NphB Prenylation The initial structure of NphB is obtained from the QM/MM X ray refinement conducted in the previous study 153 The S1 state structure is selected and then solvated in a truncate d octahedron box filled with TIP3P water models. The edge of the box is at least 8 away from the nearest solute atom. The enzyme is modeled using AMBER force field ff99SB, while the substrates are modeled using GAFF. The atomic charges of two substrates from previous study are adopted here. PME is applied to describe the long range electrostatic interactions. An 8 cutoff of the real space nonbond interaction and an approximately 1 for the reciprocal space grid spacing are used. All hydrogen involved b onds are constrained with SHAKE algorithm. The resulting system is first fully minimized with a weak, harmonic restraint applied to remove any possible close contacts. The temperature of the system is slowly heated to 300 K using Langevin dynamics during a 100 ps time period. During this heating process, the positional restraints are gradually removed. Consequently, the system is simulated in the NPT ensemble for 900 ps with a 1 fs time step to further equilibrate the density. After that a 4 ns NPT simula tion with a 2 fs time step is performed to obtain a descent starting point for QM/MM simulations. A snapshot of the final 1 ns simulation is taken to serve as the initial structure for QM/MM simulation. The system is then partitioned into QM region and MM region. Only two substrates are included in the QM region since preliminary results showed this selection produces the most stable simulation while does not sacrifice accuracy (Figure 5 1) The system is quickly minimized with the SC C DFTB/ff99SB hybrid p otential followed by re heating to 300 K in
131 NVT ensemble during a 50 ps period with a 0.5 fs time step. The pressure regulation is then turned on to have the system equilibrated in the NPT ensemble for 450 ps with a 1 fs time step. Figure 5 1. QM parti tion of NphB: GPP and 1,6 DHN. Important atom specification s are given. O1p in figure denotes O 1P in text (which means O 1 from diphosphate substrate ), and C5a denotes C 5A in text (which means C 5 in aromatic substrate). The resulting system is then used for a SMD scan to access a quick view of the FES shape. The RC is defined as the distance between C 1 (C 1P ) of GPP and C 5 (C 5A ) of 1,6 DHN. The linear combination of RC, which defined as d RC = d C1P O1P d C1P C5A appears brin g ing instability into the PMF stud y, possibly due to the existence of multiple reacting capable atoms in 1,6 DHN. While a 2 d RC 1 = d C1P O1P and d RC 2 = d C1P C5A increases the computational cost significantly but does not provide additional accuracy. A set of 26 starting structures are extracted from the trajectory of SMD covering the RC from 1.5 (product state) to 4.0 (resting state), with a 0.1 interval. Another 14 windows are prepared to cover the middle points of two originally generated windo ws for the reaction important area, namely 1.6
132 to 3.0 . A set of resulting 40 windows of USP is then carried out with ea ch windows simulated for 250 ps, thus a total of 10 ns QM/MM simulation is accomplished for this study. The last 150 ps simulation of each window is used for data collecting. The results are then processed with Dr. The S3 state is reached by applying the following NMR restraints (C 1P C 5A : 5.9 , 10 kcal/mol ; C 1P C 2A : 4.0 , 10 kcal/mol ) to an extended QM/MM simulation after the already mentioned 4 ns simulation. A SMD scan is then carried out with the new generated S3 state structure with the new RC defined as d C1P C2A The same strategy as described for the S1 state is adopted into the QM/MM PMF study of S3 state coupled with USP method. WHAM cod e is utilized again to generate the free energy profile. The S4 state is also studied. The initial structure is obtained from a snapshot of QM/MM equilibration o f S3 state. The reason of such a selection will be explained later. A SMD scan is carried out with the resulting structure to propagate the trajectory over the RC, where RC is defined as d RC = d C1P C3A A following QM/MM PMF study coupled with USP techniqu e, however, has not yet been performed. 5.2.2 Modelling and Simulation of Proton Transfer Step in NphB Catalysis After the prenylation, a proton transfer step is required to remove the redundant hydrogen from the prenylated 1,6 DHN, thus restoring its arom atic quality. The proton transfer is usually much more complicated than other reactions due to the highly possible quantum tunneling effect resulting from the small atomic mass of hydrogen atom. In this study, we still adopt the SMD US strategy to study th is reaction with QM/MM simulation. However, ab initio level calculation will be conducted later to study the tunneling effect. The initial structure (see Figure 5 2) is a snapshot taken from the S1 QM/MM USP at d RC = 1.5 that corresponds to the product s tate. After carefully examine the resulted structure, two
133 candidates to perform such a proton abstraction have isolated. One of them is the Tyr216 and another one is a crystal water molecule. Figure 5 2. Active site configuration of the initial structu re selected for proton transfer step simulation. WAT refers to the water molecule that is included in QM region at this stage. Tyr216 stabilize the PPi diphosphate, and it situates within 4.5 away from C 5A of 1,6 DHN. A reasonable tyrosine mediate proton transfer process would be the diphosphate first abstracts the hydroxyl oxygen and then consequently the tyrosine compensates itself by capturing the extra proton from geranylated 1,6 DHN. However, our preliminary resul ts show the first step associates with an unreasonably large free energy barrier (45+ kcal/mol from SMD scan). Therefore, this hypothesis has not been further studied.
134 On the other hand, a crystal water molecule has been observed approaching 1,6 DHN when the new C C bond starts to form. This water molecule i nteracts with the hydroxyl group at C 6A thus staying close to the redundant proton at C 5A An equilibration of 250 ps is first performed to validate that this water molecule remains its position during longer simulation. Within this equilibration, this water molecule is added into QM region. Then a snapshot is extracted from the trajectory resulted from the last 50 ps of such equilibration. The resulting structure is scanned with a SMD simulation with R C defined as d RC = d OW H5A A following QM/MM USP study has not been performed yet. 5.2.3 Classical and QM/MM study of FtmPT1 The FtmPT1 model is prepared based on the structure of PDB 3O2K (Figure 5 3) The system is first gone through MolProbity (at http ://molprobity.biochem.duke.edu/ ) to improve the structure. Several flip are taken place at a couple of glutamines and an asparagine residue. The resulting system is then prepared with LeaP in the AMBER suite of programs. AMBER force field ff99SB and GAFF a re used to model the enzyme and substrates, respectively. The atomic charge of DMAPP and brevianamide F and derived following a two stage RESP fitting procedure. After that, the system is solvated into a truncated octahedron water box with TIP3P water mode l applied. The edge of the box situates at least 8 away from the closest solute atom. Long range electrostatic interactions are computed using PME method. SHAKE algorithm is implemented to apply constraints to all bonds with hydrogen atoms. This model is termed as DPP in the following text. Another system is prepared in the same way for FtmPT1 complexed with monoprotonated DMAPPH 2 in order to testify the substrate preference of this enzyme. This new model will be called DPH in this study.
135 Figure 5 3. Active site configuration of PDB ID 3O2K. DMAPP and Brevianamide F are treated with QM in QM/MM simulation. Important atoms are labeled in red. Two resulting sy stems are first fully minimized. A weak, harmon ic force is applied to restrain the coordinates of protein and ligand ato ms. Then the systems are heated to 300 K during a 100 ps period with a 1 fs time step in canonical ensemble. The positional restraints are gradually removed during the heating processes. The resulting systems are equilibrated in th e NPT ensemble for 900 ps with a 1 fs time step, followed by a 4 ns simulation for each system. A snapshot of the final 1 ns simulation of DPP is extracted for each system and served as the initial structure for the following QM/MM study. T he QM/MM study of DPH has not yet been conducted. Similar to NphB, only two reactants are selected in the QM partition of the system, while all remaining system is modeled with MM potential. The resulting two systems are then minimized with hybrid QM/MM potential and re heating to 300 K during a 50 ps time
136 period in NVT ensemble with a 0.5 fs time step. Then both systems are equilibrated with QM/MM potential for 450 ps with a 1 fs time step. A snapshot from the last 50 ps simulation for each system is chosen as the initia l structure for the PMF study. Due to the large separation between C 1P and C 2A a 2 D PMF study would be very computationally consuming. Th us a two stage 1 D PMF study is appli ed to each system. In the first stage, the RC is defined as d C1P O1P while in t he second state d C1P C2A denotes the RC. A first stage SMD scan for each model is conducted to simulate the hydrolysis of DMAPP. The resulting structures are then undergoing another SMD scan at the second stage to study the bond forming process. A two stag e study of US P method coupled PMF is also required Such a study is currently undergoing in our lab. After obtaining knowledge of the RC variation from two stage 1 D PMF studies, it is possible to carry out 2 s that only cover th ose important areas while leaving physically unfavorable area s out of study. In such a study, the two d RC 1 = d C1P C2A and d RC 2 = d C1P O1P This study would provide additional support to the conclusion drawn from 1 D PMF study. 5.3 Results and Discussions 5.3.1 Orientation of Substrate Orientation in QM/MM Equilibration of NphB A QM/MM X ray refined structure of NphB S1 state has been adopted as the initial structure and solvated into a water box in truncated octahedro n shape. The resulting system is first fully relaxed through MM equilibration. During this procedure, d C1P C2A and d C1P C5A fluctuates as the system is transforming between the major product of S1 and the minor product state of S3. A snapshot from the last 1 ns MM simulation that best represents the S1 state is extracted to serve as the initial structure of the consequent QM/MM simulations.
137 During QM/MM simulation, GPP and its binding pocket are fairly stable. The salt bridge interactions between negatively charged side chains of Lys119, Arg228 and Lys284 help stabilizing the isoprenoids substrate. In addition, Tyr216 functions as a hydrogen bond donor to diphosphate. Tyr121 situates on the opposite side of geranyl group to the aromatic substrate thus it looks like the lipid group being sandwiched by the two aromatic units. Such a pose would be essential to the catalysis if the reaction adopts a dissociative pathway, since two aromatic groups could stabilize the developin g carbocation through cation 211 218 Carefully examining the relative positions of Tyr121 and 1,6 DHN, it is found that the aromatic ring of the tyros ine residue points to the C 3P position of the geranyl moiety while the aromatic group of 1,6 DHN (in the S1 state orientation) is in a perfect position to interact with C 1P Consider ing the fact that the positive charge will be delocalized between C 1P and C 3P this staggered arrangement of two aromatic units is particular interesting and important. The orientation of 1,6 DHN fluctuates between S1 and S3 state through the QM/MM equilibration (see Figure 5 4) The side chains of Ser214 and Tyr288 play import ant role to help stabilize the S3 state orientation via hydrogen bond interactions. A configuration of 1,6 DHN that leads to the C 4A closer to C 1P of GPP 3 than either C 5A or C 2A has been observed during the orientation transition, however, is fairly unsta ble as it nearly instantly transforms to other orientations. That finding leads to an assumption that the 4 prenyl 1,6 DHN is an associated product of either S1 state or S3 state orientation. This assumption makes sense if one considers the fact that at th e equilibrium conformation of S3 state C 4A is quite accessible to the C 1P since the distance between C 1P and C 4A is fluctuating between approximately 1.0 shorter to 1.2 longer than that d C1P C2 A
138 Figure 5 4. 1,6 DHN orientation in S1 (left) and S3 ( right) State. 5.3.2 Free Energy Surface of NphB Catalyzed Prenylation 22.214.171.124 S1 s tate QM/MM PMF studies are carri ed out to study the FES associated with different product formation processes of NphB catalysis. A snapshot from the final 50 ps QM/MM equili bration of S1 state that best represents the resting state of such an orientation ( d C1P C5A d C1P C2A is prepared as the initial structure for the future PMF study. A SMD scan is performed with the resulting structure. 40 structures ar e consequently extracted from the SMD trajectory to serve as starting structures for the following US. After unbias ing the distribution of each windows using WHAM code, a free energy profile is generated. An intermediate is identified on such a FES at appr oximate 2.10 along the RC. Two free energy barriers, at 2.55 and 2.00 respectively, separate the intermediate state to either the reactant state or the product state. The barrier at 2.55 (TS1) is the highest barrier, corresponding to 12.35 kcal/mol free energy (see Figure 5 5 ), while the second barrier at 2.00 (TS2) is fairly small (about 5.55 kcal/mol) comparing to the intermediate state (approximately 5.10 kcal/mol).
139 Figure 5 5 Free energy profile associated with prenylation step for both S1 and S3 model. The profile of S3 does not start from 0 kcal/mol, instead, it starts at 2.3 kcal/mol, which is the free energy difference between S1 and S3 orientation. The intermediate features a carbocation (Figure 5 5). The H 1P H 2P C 1P C 2P (see Figure 5 1 for atom specifications) dihedral reads nearly 180 upon the forming of such a carbocation, which strongly indicating the generation an sp 2 carbon at C 1P Moreover, this dihedral stays above 160 during the RC moving from approximately 2.6 to 2.3 and fluctuates above 150 until the intermediate state is passed. Comparing to the fluctuation of such a dihedral in the farnesylation of FTase/FPPH/CVIM complex (see Figure 5 7), the characteristic of a carbocation in this geranylation is more evident. Ther efore, we propose this NphB catalyzed geranylation adopts a typical S N 1 reaction mechanism. Tyr121 sandwiches the developing carbocation along with 1,6 DHN via cation interactions. The staggered arrangement of two aromatic groups guarantees that the posi tive charge is stabilized by at least one cation charge. In fact, after the C 1P O 1P
140 On the other hand, Tyr216, which forms diphosphate, is within 4.0 to 4.8 away from the C 1P atom after the carbocation forms. In spite of the distance significantly longer than that between the c arbocation and either Tyr121 or 1,6 DHN, it is still within the interaction range of such a cation shaped cation interaction pocket is formed and possibly steers the carbocation to electrophilic att ack the C 5A of 1,6 DHN. Figure 5 6 Active site configuration at the intermediate state of S1 model. O 1P of PPi, C 1P of carbocation and C 5A of 1,6 DHN are labeled in red. Note Tyr121 and 1,6 DHN are staggere The diphospha te leaving group, with a big 4 units negative charge, is stabilized by Lys119, Asn173, Arg228, Lys284 and Tyr216 via hydrogen bond interactions. Lys169 might also
141 diphosphate moiety Upon passing the TS2, the C 1P C 5A bond starts forming. At the same time, a crystal water molecule reaches in to within reacting range to the 1,6 DHN. Such a water molecule might function to abstract the extra proton at C 5A Figure 5 7 H 1(P) H 2(P) C 1 (P) C 2(P) dihedral variation in the chemical reaction of NphB catalyzed 5 prenylation (red) and FTase/FPPH/CVIM catalyzed farnesylation (blue). The computed free energy barrier is 12.35 kcal/mol, while the experimentally measured k cat for this reaction is 4200 200 s 1 corresponding to 12.6 kcal/mol in barrier height associated with the product formation. In addition to the excellent agreement between our calculation and experimental observation, this value suggests that the hydrolysis step is the RLS for this NphB catalyzed geranylation. Although the possibility of the proton abstraction step being the RLS cannot be totally excluded, it is very unlikely that the free energy barrier associated with such a
142 proton transfer will be larger than 12 kcal/mol, le ave alone the highly possible accompanied quantum tunneling effect due to the small mass of hydrogen atom. 126.96.36.199 S3 state and S4 s tate The initial structure of S3 state for the QM/MM simulation is obtained by artificially building such a configuration w ith NMR restraints applied. One of the snapshots during last 10 ps of such a restraint simulation is then extracted and prepared for the QM/MM PMF study. Like S1 study, a SMD scan is first conducted to propagate the trajectory along the RC from 4.0 (rea ctant state) to 1.5 (product state) A set of 40 structures are then extracted from the trajectory and prepared for the consequent US After 40 windows of simulations finish, data is collect ed from each window and processed with WHAM code to generate the unbiased free energy profile. The resulting profile once again implied a dissociative mechanism ( Figure 5 5 ) Similar to S1, TS1 is identified at approximate 2.5 w ith a barrier height of about 13. 50 kcal/mol, and TS2 is found at approximate 1.95 with a barrier height of around 11.3 5 kcal/mol. Two TS are connected by an intermediate state at approximately 2.10 with a free energy approximately 10.4 0 kcal/mol higher than the S1 resting state. Comparing to the profile previously obtained for S1. An inte resting note is the difference at the RLS barrier height. S1 gives a 12.35 kcal/mol barrier while S3 produces a barrier of 13.50 kcal/mol, which is 1.15 kcal/mol higher Experimental results show an approximate 4 fold difference in their corresponding k cat that results in a 0.8 kcal/mol difference (S3 higher than S1). Such a value is perfectly within chemical accuracy comparing to the experimental measurement. In previous MM PMF study, S1 state is found approximately 2.3 kcal/mol energetically more favorab le than S3 state. Put this energetic difference into consideration a 11.15 kcal/mol barrier height for S3 chemical step is obtained.
143 S3 clearly has a much higher free energy than S1 with regard to the second barrier and the intermediate state. Stru ctural analysis provides usually information to help understanding this difference. In their corresponding resting state, 1,6 for the hydrophobic interactions with geranyl group. The resulted dispersion might be the source of an easier hydrolysis, or C 1P O 1P bond broken. However, the orientation of 1,6 DHN in the product state is more similar to the orientation in S1 resting stat e. One possible explanation of this orientation change is in its original orientation the repulsion between aromatic groups and the attaching geranyl group is too strong thus 1,6 DHN has to adopt a pose similar to S1. Such a repulsion force or the reorient ation of 1,6 DHN possibly causes the energy increase associated with the intermediate state and TS2. The active site configuration of the intermediate state (see Figure 5 8 ) represents a very similar look to S1 state with the exception of a different orien tation of 1,6 DHN. The carbon atoms in the first isoprene unit almost situate in the same plane, while the H 1P H 2P C 1P C 2P dihedral is nearly 180. The carbocation is sandwiched between the aromatic moiety of Tyr121 DHN. C 1P C 2P and C 3P where the charge delocalization could take place, all have close aromatic moiety available to provide cation side chains of Lys119, Lys169, Asn173, Tyr216 Arg228 and Lys284 help stabilized the PPi via hydrogen bond interactio ns, in addition to the stable Mg 2+ coordination. Tyr288 forms a hydrogen bond with the hydroxyl group at C 6P of 1,6 DHN, by serving as either a hydrogen bond donor or acceptor. When the C 1P and C 2A reach into the reacting range, a couple of crystal water m olecules, one of which is exactly the possible proton acceptor in S1, approach to the hydroxyl group of 1,6 DHN at C 1A and possibly serve as the proton abstractor.
144 Figure 5 8 Active site conformation in S3 model at the intermediate state. O 1P of PPi, C 1 P of carbocation, and C 2A of 1,6 DHN are labeled in red. Notice 1,6 DHN starts to rotate counter clockwise in order to reduce steric repulsion for the following attachment. The free energy profile obtained from SMD simulation (see Figure 5 9 ) also support s a dissociative mechanism for this model. The active site configuration after carbocation forming looks very similar to that of S3 state ( Figure 5 10 ). 1,6 able to stabilize the carbocation via cation ar 1P C 2P and C 3P positions. Tyr121, although still cooperates with 1,6 DHN sandwiching the carbocation, seems moving away from C 3P Thus it is unclear whether it still helps stabilizing the carbocation through cation ractions. The free energy barrier associated with this product is more than 15.0 kcal/mol. Qualitatively, this value is in agreement with the experiment evidence that S4 associated with a free energy higher than both S1 and S3. Quantitatively, however, thi s value is
145 not conclusive enough since several other important factors still need to be considered, such like SMD usually produces a free energy barrier than USP as well as the free energy difference associated with the orientation of 1,6 DHN between S1 or S3 and S4 remains unknown. Figure 5 9 Free energy profile associated with S4 model prenylation. This profile starts at 2.3 kcal/mol, which is the free energy difference between S1 and any observed state orientation state of 1,6 DHN. 5.3.3 Proton Trans fer Step of NphB Catalysis The NphB chemistry includes both a prenylation step and a proton transfer step that restores aromatic quality of prenylated 1,6 DHN. Proton transfer is usually a tricky reaction in computational simulation. The quantum tunneling effect is usually strong in this kind of reactions because hydrogen possesses such a small atomic mass. Therefore, QM calculation is required to fully understand such a scenario. In this study, we mainly focus on whether or not the deprotonation of gerany lated 1,6 DHN is the RLS of NphB catalyzed prenylation reaction. Regarding this question, our previous analysis of the prenylation step has strongly implied that the deprotonation step is not the RLS. In addition, KIE experiment results of indole prenyltra nsferase system also disaffirm the deprotonation step
146 in a similar catalysis as the RLS. However, a more persuasive conclusion still requires QM/MM simulation coupled with PMF study. Figure 5 10 Active site snapshot of S4 prenylation. O 1P of PPi C 1P o f carbocation and C 4A of 1,6 DHN are red labeled. The initial structure is first equilibrated with QM/MM potential for 250 ps in NPT ensemble. During the last 50 ps simulation, the water oxygen (O w ) remains in the middle of diphosphate group and geranyl 1 ,6 DHN (Figure 5 2) The distance between O w and H 5A fluctuates around 3.0 , while the closest distance between water hydrogen (H w ) and O 1P stays around 3.5 . In addition, this water molecules forms hydrogen bond with both the hydroxyl group of geranyl 1 ,6 DHN and a non O 1P diphosphate oxygen atom, thus making itself able to maintain its position and serve as a proton capturer as well.
147 Figure 5 11 Free energy profile associated with the S1 proton transfer step. This profile is obtained from a SMD scan. Note the startin g point (2.6 ) corresponds to the end point (1.5 ) in Figure 5 5 A SMD scan is then taken to propagate the trajectory along the RC, d OW H5A from 3.0 to 0.95 , which is the default O H bond length of TIP3P water model though this water is no longer a TIP3P water model. The resulting free energy profile gives an approximate 13.6 kcal/mol free energy barrier (see Figure 5 11 ). However, consider (1) the start point that is the end point of prenylation step (of S1 state) corresponds to a 6.0 kcal/mol fr ee energy difference to the resting state and (2) the SMD free energy barrier height is usually higher than the USP result and (3) the existence of quantum tunneling effect will significantly drops the real barrier height, it seems the prenylation step is indeed the RLS of the NphB catalyzed geranylation. The resulting product structure, 5 prenyl 1,6 DHN, displays a relaxed configuration (Figure 5 12). Tyr288 help stabilize this product by donating a hydrogen bond to one of the carboxyl groups of 5 prenyl 1 ,6 DHN. The H 3 O resulted from the proton abstraction forms hydrogen bond with the hydroxyl group of Tyr216 as a hydrogen bond donor, while the tyrosine
1 48 serves as a hydrogen bond donor to interact with the PPi leaving group. Such a configuration implies the re might be a context proton transfer scheme between these three residues, however, it is out of our discussion. Figure 5 12 Active site configuration after proton abstraction from 5 prenyl 1,6 DHN to H3O (WAT in Figure 5 2). O 1P of PPi, C 1P of geranyl moiety and C 5A of aromatic product are labeled in red. 5.3.4 Substrate Protonation State Preference of FtmPT1 The MM simulation of FtmPT1 complexed with brevianamide F and DMAPP has been simulated with two models with two different protonation state of DM APP adopted. In the first model, fully deprotonated DMAPP (DMAPP 3 ) is selected while in the second model the monoprotonated form of such a n isoprenoid diphosphate (DMAPPH 2 ) is chosen.
149 Both models have been equilibrated with MM potential in NPT ensemble for several nano seconds. The resulting trajectories are examined and compared carefully. RMSD (root mean square deviation) and RMSF (root mean square fluctuation) information of two models in the last 2 ns of simulation are extracted. Both simulations app ear stable over this period of time, though the fluctuation of DPP appears slightly bigger than DPH. The same trend has been observed in the RMSF plot, whereas the important areas appear stable in both models. The distance s between two reacting atoms, C 1 P of DMAPP and C 2A of brevianamide F, have also been pulled out. The result shows that DPH adopts a much shorter RC than DPP In previous computational study of FTase 14 1 the monoprotonated form the FPP seems the only form able to steer the conformational transition in the absence the magnesium. In FtmPT1 chemistry however, Mg 2+ dependency is not observed, thus whether the enzyme prefers DMAPPH 2 or its deprotonated f orm, DMAPP 3 is of great interest to us. The active site configurations of two models are examined separately. In DPH model, side chains of positively charged residues such as Arg113, Lys201 and Lys294 anchor the negatively charged diphosphate moiety via salt bridges. In addition th e side chains of Tyr203, Gl n380 Thr448 and Tyr450 help stabilize the diphosphate group via hydrogen bond interactions. Four tyrosine residues, including Tyr203, Tyr296, Tyr382 and Tyr450 (see Figure 5 3) surround the binding pocket of DMAPPH 2 F itself, form an aromatic rich environment, which is proposed to function as a steer to engineer the rotation of developing carbocation required by the new bond forming. Two reacting atoms are approximately 4.5 apart, while C 3P of DMAPPH 2 is also within 4.5 to the C 3A atom of brevianamide F, where an alternative reaction will take place upon the G115T mutation.
150 On the other hand, in DPP model, the diphosphate moiety is al so anchored by Arg113, Lys201 and Lys294 via salt bridge interactions, and the side chains of Tyr203, Gln308 and Tyr450 still form hydrogen bond interactions with diphosphate, providing additional support. The distance between C 1P and C 2A fluctuates around 5.0 while the C 3P C 3A distance is approximate 4.0 . The tyrosine shield and aromatic environment are stable throughout the simulation as well. Different to FTase, which requires a conformational transition step prior to the chemical step, no evidence s hows that FtmPT1 require s such a physical step, although the two reacting atoms are at least 4.5 apart (in DPH model more than 5 away in DPP model ). Therefore, a PMF study at MM level might be unnecessary and useless. As a result, QM/MM PMF study is t he only way to testify the substrate preference of this enzyme. 5.3.5 A Look at FtmPT1 Catalysis The chemical reaction of FtmPT1 catalysis is studied by dividing the PMF into two parts. The first part is the hydrolysis step that adopts d C1P O1P as the RC, while the second part is the bond formation step that a RC of d C1P C2A is employed. Such a strategy is a compromise of a 1 D PMF study with a RC of either d C1P O1P or d RC = d C1P O1P d C1P C2A and a 2 D PMF study with d RC 1 = d C1P O1P and d RC 2 = d C1P C2A The former case fails to give a reasonable FES because the simulations tend to move the DMAPP close to brevianamide F before breaks the C 1P O 1P bond thus the free energy barrier is mistakenly raised due to the huge steric effect associated with the movemen t of DMAPP. On the other hand, the 2 D PMF will be a very expensive and time consuming approach without gaining any knowledge about the reaction pathway. The adopted strategy is rational because within each step the free energy is described as a function o f the most important property.
151 Figure 5 13 Carbocation formed in FtmPT1 catalysis. Important atoms are labeled in red. H1p corresponds to H 1P C 2A A QM/MM SMD scan is conducted to study the hydrolysis step first. The RC is moved from 1.5 (reactant state) to 3.5 where C 1P and O 1P are sufficiently apart thus a carbocation should have already formed. The resulted free energy profile features a steadily increase before d RC = 2.7 and a flat curve thereafter. Such a FES implies the formation of a carbocation after the two bonding atoms are forced sufficiently apart. The active site configuration extracted from the trajectory provides addition s upport to the existence of a carbocation as five carbon atoms of dimethylallylic group nearly situate in the same plane, while the dihedral of H 1P H 2P C 1P C 2P remains nearly 180 after 2.7 point (Figure 5 13). The resulted dimethylallylic cation is sandw iched by two aromatic units: brevianamide F and a tyrosine residue involved in the
152 tyrosine shield. Strong cation carbocation. In addition, the shape of the tyrosine shield enables the steer of th e carbocation rotation in order to approach another reactant, since continuous cation at least one of four tyrosine residues regardless the orientation of dimethylallylic cation inside the pocket. An initial structure of the seco nd step is extracted from the SMD trajectory with d RC = 2.8 upon where a carbocation is already generated. The resulted structure is then undergoing another QM/MM SMD scan with a RC of d C1P C2A In the configuration of such an initial structure, C 1P and C 2A are 5.8 apart, thus the trajectory is propagated to move the RC from 5.8 to 1.5 . The resulted free energy curve (data not shown) indicates the existence of an intermediate state after the carbocation rotates its orientation slightly. Another free energy barrier is observed at d RC = 2.1 , due to the strong repulsion generated in the approaching of two carbon atoms. The hydrogen at C 2A position distinctly bends out of the aromatic plane, indicating a proton removal step is required to restore the a romatic quality of the aromatic substrate Despite the qualitatively understanding of FtmPT1 catalyzed prenylation, the quantitative properties, such as free energy barrier height or the exact location of important state, obtained from SMD simulation are n ot sufficiently conclusive. In fact, it is not unusual to find distinguishable difference between SMD generated free energy profile and USP resulted free energy curve. Therefore, a set of QM/MM USP is still required to complete our understanding of this ca talysis. In addition, effort should be also made to study the proton abstraction step and the impact of key residue mutation such as G115T. These works are currently undergoing in our lab. 5.4 Summary In this study, effort is made to illustrate the reactio n mechanism catalyzed by two aromatic prenyltransferases: NphB and FtmPT1.
153 NphB is determined to catalyze the geranylation via a dissociative pathway. A stable carbocation is observed for each of three models leading to there different products. The free energy param eters obtained from QM/MM USP (12.35 kcal/mol for S1 and 13.50 kcal/mol for S3) are in excellent agreement with results from kinetic experiments (12.6 kcal/mol for S1 and 13.3 kcal/mol for S3). Although the free energy barrier obtained for S4 t hought a QM/MM SMD simulation qualitatively matches the experimentally observed trend, quantitatively it is not convincing enough. The proton transfer step is studied with S1 model the given free energy profile suggests a barrier height lower than the hyd rolysis of GPP, therefore, the RLS is determined to be the hydrolysis step in the prenylation. Indole prenyltransferase FtmPT1 is also studied. In the MM simulation of both FtmPT1 complexed with DMAPPH 2 and FtmPT1 complexed DMAPP 3 stable simulation s ar e reached. The distance between two reacting atoms is slightly smaller in DPH model than in DPP model. Except that, everything looks similar in two active site configurations. QM/MM PMF is required to identify which protonation state of DMAPP is more favor able to FtmPT1. However, due to the large separation between two reactants, traditional 1 D PMF study fails because the system prefers to move the reactants into a close range before forcing bond cleavage. A 2 D PMF without preliminary knowledge of reactio n pathway requires a huge number of sampling windows, thus making itself not an ideal implementation. As a compromise, a two stage 1 D PMF strategy is adopted and gives exciting results. The 2 stage 1 D SMD scan reveals a stable carbocation and an intermed iate state associated with the rotation of such a carbocation. Tyrosine rich environment in DMAPP biding pocket plays an important role to stabilize the lipid group and furthermore engineer its rotation via cation aromatic
154 The aromatic rich environment in the active sites of aromatic prenyltransferases makes huge impact on the mechanism selection of this type of enzymes. Their ability to stabilize the developing carbocation via cation is essential to the S N 1 like mechanism. Such findings make an interesting argument about the target dependent mechanism sele ctivity of prenyltransferases, in addition to the mechanism switch associated with the fascinating context dependent substrate sele ctivity across different FTase/peptide complexes.
155 LIST OF REFERENCES (1) Cox, A. D. Drugs 2001 61 723 32. (2) Zhang, F. L.; Casey, P. J. Annu. Rev. Biochem. 1996 65 241 6 9. (3) Kuzuyama, T.; Noel, J. P.; Richard, S. B. Nature 2005 435 983 7. (4) Jost, M.; Zocher, G.; Tarcz, S.; Matuschek, M.; Xie, X. L.; Li, S. M.; Stehle, T. J. Am. Chem. Soc. 2010 132 17849 58. (5) Pojer, F.; Wemakor, E.; Kammerer, B.; Chen, H. W.; Walsh, C. T.; Li, S. M.; Heide, L. Proc. Natl. Acad. Sci. U. S. A. 2003 100 2316 21. (6) Piironen, V.; Lindsay, D. G.; Miettinen, T. A.; Toivo, J.; Lampi, A. M. J. Sci. Food Agric. 2000 80 939 66. (7) Park, H. W.; Boduluri, S. R.; Moomaw, J. F.; Cas ey, P. J.; Beese, L. S. Science 1997 275 1800 4. (8) Casey, P. J.; Seabra, M. C. J. Biol. Chem. 1996 271 5289 92. (9) Clarke, S. Annu. Rev. Biochem. 1992 61 355 86. (10) Caldwell, G. A.; Naider, F.; Becker, J. M. Microbiological Reviews 1995 59 406 22. (11) Schafer, W. R.; Rine, J. Annu. Rev. Genet. 1992 26 209 37. (12) Omer, C. A.; Gibbs, J. B. Mol. Microbiol. 1994 11 219 25. (13) Glomset, J. A.; Farnsworth, C. C. Annual Review of Cell Biology 1994 10 181 205. (14) Newman, C. M. H.; Ma gee, A. I. Biochim. Biophys. Acta 1993 1155 79 96. (15) Park, H. W.; Beese, L. S. Curr. Opin. Struct. Biol. 1997 7 873 80. (16) Fu, H. W.; Beese, L. S.; Casey, P. J. Biochemistry 1998 37 4465 72. (17) Long, S. B.; Casey, P. J.; Beese, L. S. Bioche mistry 1998 37 9612 8. (18) Long, S. B.; Casey, P. J.; Beese, L. S. Nature 2002 419 645 50. (19) Zhang, F. L.; Casey, P. J. Biochem. J. 1996 320 925 32. (20) Taylor, J. S.; Reid, T. S.; Terry, K. L.; Casey, P. J.; Beese, L. S. EMBO J. 2003 22 59 63 74.
156 (21) Hartman, H. L.; Bowers, K. E.; Fierke, C. A. J. Biol. Chem. 2004 279 30546 53. (22) Reid, T. S.; Terry, K. L.; Casey, P. J.; Beese, L. S. J. Mol. Biol. 2004 343 417 33. (23) Maurer Stroh, S.; Washietl, S.; Eisenhaber, F. Genome Biology 2 003 4 (24) Pfeffer, S. R.; Diracsvejstrup, A. B.; Soldati, T. J. Biol. Chem. 1995 270 17057 9. (25) Nuoffer, C.; Balch, W. E. Annu. Rev. Biochem. 1994 63 949 90. (26) Reiss, Y.; Goldstein, J. L.; Seabra, M. C.; Casey, P. J.; Brown, M. S. Cell 1 990 62 81 8. (27) Seabra, M. C.; Reiss, Y.; Casey, P. J.; Brown, M. S.; Goldstein, J. L. Cell 1991 65 429 34. (28) Seabra, M. C.; Goldstein, J. L.; Sudhof, T. C.; Brown, M. S. J. Biol. Chem. 1992 267 14497 503. (29) Reiss, Y.; Stradley, S. J.; Gie rasch, L. M.; Brown, M. S.; Goldstein, J. L. Proc. Natl. Acad. Sci. U. S. A. 1991 88 732 6. (30) Troutman, J. M.; Andres, D. A.; Spielmann, H. P. Biochemistry 2007 46 11299 309. (31) Hougland, J. L.; Lamphear, C. L.; Scott, S. A.; Gibbs, R. A.; Fierk e, C. A. Biochemistry 2009 48 1691 701. (32) Zhang, F. L.; Fu, H. W.; Casey, P. J.; Bishop, W. R. Biochemistry 1996 35 8166 71. (33) Reiss, Y.; Brown, M. S.; Goldstein, J. L. J. Biol. Chem. 1992 267 6403 8. (34) Myers, L. C.; Jackow, F.; Verdine, G. L. J. Biol. Chem. 1995 270 6664 70. (35) Edwards, D. J.; Gerwick, W. H. J. Am. Chem. Soc. 2004 126 11432 3. (36) Sacchettini, J. C.; Poulter, C. D. Science 1997 277 1788 9. (37) Bohlmann, J.; Meyer Gauen, G.; Croteau, R. Proc. Natl. Acad. Sci. U. S. A. 1998 95 4126 33. (38) Ashby, M. N.; Kutsunai, S. Y.; Ackerman, S.; Tzagoloff, A.; Edwards, P. A. J. Biol. Chem. 1992 267 4128 36. (39) Suvarna, K.; Stevenson, D.; Meganathan, R.; Hudspeth, M. E. S. J. Bacteriol. 1998 180 2782 7.
157 (40) Yaza ki, K.; Sasaki, K.; Tsurumaru, Y. Phytochemistry 2009 70 1739 45. (41) Koyama, T.; Tajima, M.; Sano, H.; Doi, T.; KoikeTakeshita, A.; Obata, S.; Nishino, T.; Ogura, K. Biochemistry 1996 35 9533 8. (42) Kharel, Y.; Koyama, T. Nat. Prod. Rep. 2003 20 111 8. (43) Liang, P. H.; Ko, T. P.; Wang, A. H. J. Eur. J. Biochem. 2002 269 3339 54. (44) Tarshis, L. C.; Yan, M. J.; Poulter, C. D.; Sacchettini, J. C. Biochemistry 1994 33 10871 7. (45) Cane, D. E. Isoprenoids, Including Carotenoids and Steroid s, in Comprehensive Natural Products Chemistry ; Elsevier: London, 1998. (46) Gerlt, J. A.; Raushel, F. M. Curr. Opin. Chem. Biol. 2003 7 252 64. (47) Sacchettini, J. C.; Gordon, J. I.; Banaszak, L. J. J. Mol. Biol. 1989 208 327 39. (48) Xu, Z. H.; B ernlohr, D. A.; Banaszak, L. J. J. Biol. Chem. 1993 268 7874 84. (49) Branden, C. I.; Tooze, J. Introduction to Protein Structure ; 2nd ed. Garland, New York, 1999. (50) Shinya, K.; Imai, S.; Furihata, K.; Hayakawa, Y.; Kato, Y.; Vanduyne, G. D.; Clardy J.; Seto, H. J. Antibiot. 1990 43 444 7. (51) Kazuo, S. Y.; Furihata, K.; Hayakawa, Y.; Seto, H. Tetrahedron Lett. 1990 31 6025 6. (52) Seto, H.; Watanabe, H.; Furihata, K. Tetrahedron Lett. 1996 37 7979 82. (53) Botta, B.; Vitali, A.; Menendez, P.; Misiti, D.; Delle Monache, G. Curr. Med. Chem. 2005 12 713 39. (54) Steffan, N.; Grundmann, A.; Yin, W. B.; Kremer, A.; Li, S. M. Curr. Med. Chem. 2009 16 218 31. (55) Li, S. M. Nat. Prod. Rep. 2010 27 57 78. (56) Zhao, S.; Smith, K. S.; Deve au, A. M.; Dieckhaus, C. M.; Johnson, M. A.; Macdonald, T. L.; Cook, J. M. J. Med. Chem. 2002 45 1559 62. (57) Woehlecke, H.; Osada, H.; Herrmann, A.; Lage, H. Int. J. Cancer 2003 107 721 8. (58) Jain, H. D.; Zhang, C. C.; Zhou, S.; Zhou, H.; Ma, J. X.; Liu, X. B.; Liao, X.; Deveau, A. M.; Dieckhaus, C. M.; Johnson, M. A.; Smith, K. S.; Macdonald, T. L.; Kakeya, H.; Osada, H.; Cook, J. M. Bioorg. Med. Chem. 2008 16 4626 51.
158 (59) Shaw, D. E.; Dror, R. O.; Salmon, J. K.; Grossman, J. P.; Mackenzie, K M.; Bank, J. A.; Young, C.; Deneroff, M. M.; Batson, B.; Bowers, K. J.; Chow, E.; Eastwood, M. P.; Ierardi, D. J.; Klepeis, J. L.; Kuskin, J. S.; Larson, R. H.; Lindorff Larsen, K.; Maragakis, P.; Moraes, M. A.; Piana, S.; Shan, Y.; Towles, B. In Proceed ings of the Conference on High Performance Computing Networking, Storage and Analysis ; ACM: Portland, Oregon, 2009, p 1 11. (60) Nevins, N.; Chen, K. S.; Allinger, N. L. Journal of Computational Chemistry 1996 17 669 94. (61) Nevins, N.; Lii, J. H.; Al linger, N. L. Journal of Computational Chemistry 1996 17 695 729. (62) Nevins, N.; Allinger, N. L. Journal of Computational Chemistry 1996 17 730 46. (63) Allinger, N. L.; Yuh, Y. H.; Lii, J. H. J. Am. Chem. Soc. 1989 111 8551 66. (64) Lii, J. H.; Allinger, N. L. J. Am. Chem. Soc. 1989 111 8566 75. (65) Lii, J. H.; Allinger, N. L. J. Am. Chem. Soc. 1989 111 8576 82. (66) London, F. Zeitschrift Fur Physik 1930 63 245 79. (67) Halgren, T. A. J. Am. Chem. Soc. 1992 114 7827 43. (68) Halgre n, T. A. Journal of Computational Chemistry 1996 17 490 519. (69) Halgren, T. A. Journal of Computational Chemistry 1996 17 520 52. (70) Halgren, T. A. Journal of Computational Chemistry 1996 17 553 86. (71) Halgren, T. A.; Nachbar, R. B. Journal of Computational Chemistry 1996 17 587 615. (72) Halgren, T. A. Journal of Computational Chemistry 1996 17 616 41. (73) Claessens, M.; Ferrario, M.; Ryckaert, J. P. Mol. Phys. 1983 50 217 27. (74) Rick, S. W.; Berne, B. J. J. Am. Chem. Soc. 1996 118 672 9. (75) Sprik, M.; Klein, M. L. J. Chem. Phys. 1988 89 7556 60. (76) Dang, L. X.; Rice, J. E.; Caldwell, J.; Kollman, P. A. J. Am. Chem. Soc. 1991 113 2481 6. (77) Rick, S. W.; Stuart, S. J.; Berne, B. J. J. Chem. Phys. 1994 101 6141 56.
159 (78) Maple, J. R.; Hwang, M. J.; Stockfisch, T. P.; Dinur, U.; Waldman, M.; Ewig, C. S.; Hagler, A. T. J. Comput. Chem. 1994 15 162 82. (79) Hwang, M. J.; Stockfisch, T. P.; Hagler, A. T. J. Am. Chem. Soc. 1994 116 2515 25. (80) Maple, J. R.; Hwang, M. J.; Stockfisch, T. P.; Hagler, A. T. Isr. J. Chem. 1994 34 195 231. (81) Allinger, N. L.; Chen, K. S.; Lii, J. H. Journal of Computational Chemistry 1996 17 642 68. (82) Allinger, N. L.; Chen, K. S.; Katzenellenbogen, J. A.; Wilson, S. R.; Anstea d, G. M. Journal of Computational Chemistry 1996 17 747 55. (83) Lifson, S.; Warshel, A. J. Chem. Phys. 1968 49 5116 &. (84) Hagler, A. T.; Huler, E.; Lifson, S. J. Am. Chem. Soc. 1974 96 5319 27. (85) Hagler, A. T.; Lifson, S. J. Am. Chem. Soc. 1 974 96 5327 35. (86) Weiner, S. J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.; Alagona, G.; Profeta, S.; Weiner, P. J. Am. Chem. Soc. 1984 106 765 84. (87) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M .; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995 117 5179 97. (88) Wang, J. M.; Cieplak, P.; Kollman, P. A. Journal of Computational Chemistry 2000 21 1049 74. (89) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Proteins Structure Function and Bioinformatics 2006 65 712 25. (90) MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph McCarthy, D.; Kuc hnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz Kuczera, J.; Yin, D.; Karplus, M. J. Phys. Chem. B 1998 102 3586 616. (91) Mackerell, A. D.; Feig, M.; Brooks, C. L. Journal of Computational Chemistry 2004 25 1400 15. (92) MacKerell, A. D.; Banavali, N.; Foloppe, N. Biopolymers 2000 56 257 65. (93) Scott, W. R. P.; Hunenberger, P. H.; Tironi, I G.; Mark, A. E.; Billeter, S. R.; Fennen, J.; Torda, A. E.; Huber, T.; Kruger, P.; van Gunsteren, W. F. J. Phys. Chem. A 1999 103 3596 607.
160 (94) Christen, M.; Hunenberger, P. H.; Bakowies, D.; Baron, R.; Burgi, R.; Geerke, D. P.; Heinz, T. N.; Kastenh olz, M. A.; Krautler, V.; Oostenbrink, C.; Peter, C.; Trzesniak, D.; Van Gunsteren, W. F. Journal of Computational Chemistry 2005 26 1719 51. (95) Jorgensen, W. L.; Maxwell, D. S.; TiradoRives, J. J. Am. Chem. Soc. 1996 118 11225 36. (96) Jorgensen, W. L.; Tiradorives, J. J. Am. Chem. Soc. 1988 110 1657 66. (97) Ren, P. Y.; Ponder, J. W. Journal of Computational Chemistry 2002 23 1497 506. (98) Ren, P. Y.; Ponder, J. W. J. Phys. Chem. B 2004 108 13427 37. (99) Ren, P. Y.; Ponder, J. W. J. Phy s. Chem. B 2003 107 5933 47. (100) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. J. Chem. Phys. 1953 21 1087 92. (101) Alder, B. J.; Wainwright, T. E. J. Chem. Phys. 1957 27 1208 9. (102) Verlet, L. Physical Revie w 1967 159 98 &. (103) Hockney, R. W. Methods Computation Phys. 1970 9 135. (104) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977 23 327 41. (105) Woodcock, L. V. Chem. Phys. Lett. 1971 10 257 &. (106) Berendsen, H. J. C.; Postma, J. P. M.; Vangunsteren, W. F.; Dinola, A.; Haak, J. R. J. Chem. Phys. 1984 81 3684 90. (107) Langevin, P. Comptes Rendus Hebdomadaires Des Seances De L Academie Des Sciences 1908 146 530 3. (108) Leach, A. R. Molecular Modelling: Principl es and Applications ; 2nd ed.; Prentice Hall: New York, 2001. (109) Nicholls, A.; Honig, B. Journal of Computational Chemistry 1991 12 435 45. (110) Hawkins, G. D.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. 1996 100 19824 39. (111) Hawkins, G. D.; Cramer, C. J.; Truhlar, D. G. Chem. Phys. Lett. 1995 246 122 9. (112) Bashford, D.; Case, D. A. Annu. Rev. Phys. Chem. 2000 51 129 52. (113) Onufriev, A.; Bashford, D.; Case, D. A. J. Phys. Chem. B 2000 104 3712 20.
161 (114) Onufriev, A.; Bashford, D.; Case, D. A. Proteins Structure Function and Bioinformatics 2004 55 383 94. (115) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983 79 926 35. (116) Stilling.Fh; Rahman, A. J. Chem. Phys. 1974 60 1545 57. (117) Mahoney, M. W.; Jorgensen, W. L. J. Chem. Phys. 2000 112 8910 22. (118) Ferguson, D. M. Journal of Computational Chemistry 1995 16 501 11. (119) Barnes, P.; Finney, J. L.; Nicholas, J. D.; Quinn, J. E. Nature 1979 282 459 64. (120) Allen, M. P.; Tildesley, D. J. Computer simulation of liquids ; Clarendon Press ; Oxford University Press: Oxford [England]: New York, 1987. (121) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987 91 6269 71. (122) Bernal, J. D. ; Fowler, R. H. J. Chem. Phys. 1933 1 515 48. (123) Ewald, P. P. Annalen Der Physik 1921 64 253 87. (124) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993 98 10089 92. (125) Luty, B. A.; Davis, M. E.; Tironi, I. G.; Vangunsteren, W. F. Mol. Simul. 1994 14 11 20. (126) Luty, B. A.; Tironi, I. G.; Vangunsteren, W. F. J. Chem. Phys. 1995 103 3014 21. (127) Crespo, A.; Marti, M. A.; Estrin, D. A.; Roitberg, A. E. J. Am. Chem. Soc. 2005 127 6940 1. (128) Xiong, H.; Crespo, A.; Marti, M.; Estrin, D.; Roitberg, A. E. Theor. Chem. Acc. 2006 116 338 46. (129) Torrie, G. M.; Valleau, J. P. J. Comput. Phys. 1977 23 187 99. (130) Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, J. M. J. Comput. Chem. 1992 13 1011 21. ( 131) Ryckaert, J. P.; Bellemans, A. Faraday Discuss. 1978 95 106. (132) Jorgensen, W. L.; Gao, J.; Ravimohan, C. J. Phys. Chem. 1985 89 3470 3. (133) Goringe, C. M.; Bowler, D. R.; Hernandez, E. Reports on Progress in Physics 1997 60 1447 512.
162 (134 ) Roskoski, R. Biochem. Biophys. Res. Commun. 2003 303 1 7. (135) Resh, M. D. Cell. Signal. 1996 8 403 12. (136) Sinensky, M. Biochimica Et Biophysica Acta Molecular and Cell Biology of Lipids 2000 1484 93 106. (137) Seabra, M. C. Cell. Signal. 19 98 10 167 72. (138) Adjei, A. A. Journal of the National Cancer Institute 2001 93 1062 74. (139) Huang, C. C.; Hightower, K. E.; Fierke, C. A. Biochemistry 2000 39 2593 602. (140) Cui, G. L.; Wang, B.; Merz, K. M. Biochemistry 2005 44 16513 23. (141) Cui, G.; Merz, K. M., Biochemistry 2007 46 12375 81. (142) Hightower, K. E.; Huang, C. C.; Casey, P. J.; Fierke, C. A. Biochemistry 1998 37 15555 62. (143) Merz, K. M.; Murcko, M. A.; Kollman, P. A. J. Am. Chem. Soc. 1991 113 4484 90. (14 4) Saderholm, M. J.; Hightower, K. E.; Fierke, C. A. Biochemistry 2000 39 12398 405. (145) Wang, J. M.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. J. Comput. Chem. 2004 25 1157 74. (146) Meagher, K. L.; Redman, L. T.; Carlson, H. A. J. Comput. Chem. 2003 24 1016 25. (147) Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. J. Phys. Chem. 1993 97 10269 80. (148) Cieplak, P.; Cornell, W. D.; Bayly, C.; Kollman, P. A. J. Comput. Chem. 1995 16 1357 77.
163 (149) Frisch, M. J. T., G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, Jr., J. A.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Str atmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al Laham, M. A.; Peng, C. Y.; Nanayakka ra, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; and Pople, J. A.; Gaussian, Inc.: Wallingford CT, 2004. (150) Case, D. A.; Darden, T. A.; T.E. Cheatham, I.; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Crow ley, M.; Walker, R. C.; Zhang, W.; Merz, K. M.; B.Wang; Hayik, S.; Roitberg, A.; Seabra, G.; Kolossvry, I.; K.F.Wong; Paesani, F.; Vanicek, J.; X.Wu; Brozell, S. R.; Steinbrecher, T.; Gohlke, H.; Yang, L.; Tan, C.; Mongan, J.; Hornak, V.; Cui, G.; Mathews D. H.; Seetin, M. G.; Sagui, C.; Babin, V.; Kollman, P. A.; University of California, San Francisco: 2008. (151) Pickett, J. S.; Bowers, K. E.; Fierke, C. A. J. Biol. Chem. 2003 278 51243 50. (152) Aqvist, J. J. Phys. Chem. 1990 94 8021 4. (153) C ui, G. L.; Li, X.; Merz, K. M. Biochemistry 2007 46 1303 11. (154) Grossfield, A.; 1.6d ed. (155) Ulitsky, A.; Elber, R. J. Chem. Phys. 1993 98 3380 8. (156) Simmerling, C.; Fox, T.; Kollman, P. A. J. Am. Chem. Soc. 1998 120 5771 82. (157) Wang, X. S.; Roitberg, A. E.; Richards, N. G. J. Biochemistry 2009 48 12272 82. (158) Peters, M. B.; Yang, Y.; Wang, B.; Fusti Molnar, L.; Weaver, M. N.; Merz, K. M. J. Chem. Theory Comput. 2010 6 2935 47. (159) Pais, J. E.; Bowers, K. E.; Stoddard, A. K.; Fierke, C. A. Anal. Biochem. 2005 345 302 11. (160) Furfine, E. S.; Leban, J. J.; Landavazo, A.; Moomaw, J. F.; Casey, P. J. Biochemistry 1995 34 6857 62. (161) Mathis, J. R.; Poulter, C. D. Biochemistry 1997 36 6367 76.
164 (162) Gibbs, J. B.; Oliff A.; Kohl, N. E. Cell 1994 77 175 8. (163) Graham, S. L.; Desolms, S. J.; Giuliani, E. A.; Kohl, N. E.; Mosser, S. D.; Oliff, A. I.; Pompliano, D. L.; Rands, E.; Breslin, M. J.; Deana, A. A.; Garsky, V. M.; Scholz, T. H.; Gibbs, J. B.; Smith, R. L. J. Med. Chem. 1994 37 725 32. (164) Garcia, A. M.; Rowell, C.; Ackermann, K.; Kowalczyk, J. J.; Lewis, M. D. J. Biol. Chem. 1993 268 18415 8. (165) Nigam, M.; Seong, C. M.; Qian, Y. M.; Hamilton, A. D.; Sebti, S. M. J. Biol. Chem. 1993 268 20695 8. ( 166) James, G. L.; Goldstein, J. L.; Brown, M. S.; Rawson, T. E.; Somers, T. C.; Mcdowell, R. S.; Crowley, C. W.; Lucas, B. K.; Levinson, A. D.; Marsters, J. C. Science 1993 260 1937 42. (167) Patel, D. V.; Schmidt, R. J.; Biller, S. A.; Gordon, E. M.; Robinson, S. S.; Manne, V. J. Med. Chem. 1995 38 2906 21. (168) Manne, V.; Yan, N.; Carboni, J. M.; Tuomari, A. V.; Ricca, C. S.; Brown, J. G.; Andahazy, M. L.; Schmidt, R. J.; Patel, D.; Zahler, R.; Weinmann, R.; Der, C. J.; Cox, A. D.; Hunt, J. T.; Go rdon, E. M.; Barbacid, M.; Seizinger, B. R. Oncogene 1995 10 1763 79. (169) Manne, V.; Ricca, C. S.; Brown, J. G.; Tuomari, A. V.; Yan, N.; Patel, D.; Schmidt, R.; Lynch, M. J.; Ciosek, C. P.; Carboni, J. M.; Robinson, S.; Gordon, E. M.; Barbacid, M.; S eizinger, B. R.; Biller, S. A. Drug Dev. Res. 1995 34 121 37. (170) Patel, D. V.; Gordon, E. M.; Schmidt, R. J.; Weller, H. N.; Young, M. G.; Zahler, R.; Barbacid, M.; Carboni, J. M.; Gullobrown, J. L.; Hunihan, L.; Ricca, C.; Robinson, S.; Seizinger, B R.; Tuomari, A. V.; Manne, V. J. Med. Chem. 1995 38 435 42. (171) Tamanoi, F. Trends Biochem. Sci. 1993 18 349 53. (172) Tamanoi, F.; Esson, K.; Gomez, R.; Wood, D.; Uh, M.; Hara, M.; Akasaka, K.; Okabe, M.; Nakano, H. FASEB J. 1993 7 A1049 A. ( 173) Bishop, W. R.; Bond, R.; Petrin, J.; Wang, L.; Patton, R.; Doll, R.; Njoroge, G.; Catino, J.; Schwartz, J.; Windsor, W.; Syto, R.; Schwartz, J.; Carr, D.; James, L.; Kirschmeier, P. J. Biol. Chem. 1995 270 30611 8. (174) Chen, Z.; Sun, J. Z.; Pradi nes, A.; Favre, G.; Adnane, J.; Sebti, S. M. J. Biol. Chem. 2000 275 17974 8. (175) Maurer Stroh, S.; Washietl, S.; Eisenhaber, F. Biol. Chem. 2003 384 977 89.
165 (176) Gelb, M. H.; Van Voorhis, W. C.; Buckner, F. S.; Yokoyama, K.; Eastman, R.; Carpente r, E. P.; Panethymitaki, C.; Brown, K. A.; Smith, D. F. Mol. Biochem. Parasitol. 2003 126 155 63. (177) Mu, Y. Q.; Omer, C. A.; Gibbs, R. A. J. Am. Chem. Soc. 1996 118 1817 23. (178) Weller, V. A.; Distefano, M. D. J. Am. Chem. Soc. 1998 120 7975 6 (179) Harris, C. M.; Derdowski, A. M.; Poulter, C. D. Biochemistry 2002 41 10554 62. (180) Dolence, J. M.; Poulter, C. D. Proc. Natl. Acad. Sci. U. S. A. 1995 92 5008 11. (181) Cassidy, P. B.; Poulter, C. D. J. Am. Chem. Soc. 1996 118 8761 2. ( 182) Pickett, J. S.; Bowers, K. E.; Hartman, H. L.; Fu, H. W.; Embry, A. C.; Casey, P. J.; Fierke, C. A. Biochemistry 2003 42 9741 8. (183) Wu, Z.; Demma, M.; Strickland, C. L.; Radisky, E. S.; Poulter, C. D.; Le, H. V.; Windsor, W. T. Biochemistry 1999 38 11239 49. (184) Pais, J. E.; Bowers, K. E.; Fierke, C. A. J. Am. Chem. Soc. 2006 128 15086 7. (185) Ho, M. H.; De Vivo, M.; Dal Peraro, M.; Klein, M. L. J. Chem. Theory Comput. 2009 5 1657 66. (186) Yang, Y.; Weaver, M. N.; Merz, K. M. J. Phys Chem. A 2009 113 9843 51. (187) Scheuring, J.; Berti, P. J.; Schramm, V. L. Biochemistry 1998 37 2748 58. (188) Yang, Y.; Yu, H. B.; York, D.; Elstner, M.; Cui, Q. J. Chem. Theory Comput. 2008 4 2067 84. (189) Elstner, M.; Porezag, D.; Jungnicke l, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Physical Review B 1998 58 7260 8. (190) Elstner, M.; Cui, Q.; Munih, P.; Kaxiras, E.; Frauenheim, T.; Karplus, M. J. Comput. Chem. 2003 24 565 81. (191) Seabra, G. D.; Walker, R. C. ; Elstner, M.; Case, D. A.; Roitberg, A. E. J. Phys. Chem. A 2007 111 5655 64. (192) Gronert, S.; Lee, J. M. J. Org. Chem. 1995 60 6731 6. (193) Long, S. B.; Casey, P. J.; Beese, L. S. Structure 2000 8 209 22. (194) Ramachandran, G. N.; Ramakrishn an, C.; Sasisekharan, V. J. Mol. Biol. 1963 7 95 &. (195) Christianson, D. W. Chem. Rev. (Washington, DC, U. S.) 2006 106 3412 42.
166 (196) Santos, F. A.; Rao, V. S. N. Phytomedicine 1998 5 115 9. (197) Blanco Colio, L. M.; Tunon, J.; Martin Ventura, J. L.; Egido, J. Kidney Int. 2003 63 12 23. (198) Erdmann, K.; Berndt, G.; Grosser, N.; Schroder, H. Naunyn Schmiedebergs Archives of Pharmacology 2004 369 R14 R. (199) Grosser, N.; Oberle, S.; Berndt, G.; Erdmann, K.; Hemmerle, A.; Schroder, H. Bio chem. Biophys. Res. Commun. 2004 314 351 5. (200) Chowdhury, S. A.; Kishino, K.; Satoh, R.; Hashimoto, K.; Kikuchi, H.; Nishikawa, H.; Shirataka, Y.; Sakagami, H. Anticancer Res. 2005 25 2055 63. (201) Jahangir, T.; Khan, T. H.; Prasad, L.; Sultana, S. Redox Rep. 2005 10 303 10. (202) Soria Mercado, I. E.; Prieto Davo, A.; Jensen, P. R.; Fenical, W. J. Nat. Prod. 2005 68 904 10. (203) Zhou, Y. D.; Kim, Y. P.; Mohammed, K. A.; Jones, D. K.; Muhammad, I.; Dunbar, D. C.; Nagle, D. G. J. Nat. Prod. 2005 68 947 50. (204) Boucher, K.; Chad, S. S.; Sharma, P.; Hauschka, P. V.; Solomon, K. R. Microvasc. Res. 2006 71 91 102. (205) Jahangir, T.; Khan, T. H.; Prasad, L.; Sultana, S. Hum. Exp. Toxicol. 2006 25 235 42. (206) Kumano, T.; Richard, S. B .; Noel, J. P.; Nishiyama, M.; Kuzuyama, T. Bioorg. Med. Chem. 2008 16 8117 26. (207) Rabindran, S. K.; He, H. Y.; Singh, M.; Brown, E.; Collins, K. I.; Annable, T.; Greenberger, L. M. Cancer Res. 1998 58 5850 8. (208) Rabindran, S. K.; Ross, D. D.; Doyle, L. A.; Yang, W. D.; Greenberger, L. M. Cancer Res. 2000 60 47 50. (209) van Loevezijn, A.; Allen, J. D.; Schinkel, A. H.; Koomen, G. J. Bioorg. Med. Chem. Lett. 2001 11 29 32. (210) Allen, J. D.; van Loevezijn, A.; Lakhai, J. M.; van der Valk, M.; van Tellingen, O.; Reid, G.; Schellens, J. H. M.; Koomen, G. J.; Schinkel, A. H. Mol. Cancer Ther. 2002 1 417 25. (211) Ma, J. C.; Dougherty, D. A. Chem. Rev. (Washington, DC, U. S.) 1997 97 1303 24. (212) Dougherty, D. A. J. Nutr. 2007 137 15 04S 8S.
167 (213) Zacharias, N.; Dougherty, D. A. Trends Pharmacol. Sci. 2002 23 281 7. (214) Gallivan, J. P.; Dougherty, D. A. J. Am. Chem. Soc. 2000 122 870 4. (215) Gallivan, J. P.; Dougherty, D. A. Proc. Natl. Acad. Sci. U. S. A. 1999 96 9459 64. (216) Dougherty, D. A. Science 1996 271 163 8. (217) Kumpf, R. A.; Dougherty, D. A. Science 1993 261 1708 10. (218) Mccurdy, A.; Jimenez, L.; Stauffer, D. A.; Dougherty, D. A. J. Am. Chem. Soc. 1992 114 10314 21.
168 BIOGRAPHICAL SKETCH Yang, Yue was born in Beijing, C hina in the end of 1979. He graduated from University of Science and Technology of China in 2003 with a Bache lor of Science degree major in chemical p hysics. After college, he worked in a lab in the Technique Ins titute of Physics and Chemistry, the Chinese Academy of Science of a couple of years before coming to University of Florida and Joining the Kenneth Merz group. During his five years Ph.D. career, he has been focusing on the modeling and simulation of important enzymes and biolo gical systems.