rational 
 enzyme design has a wide scope ranging from general industrial 
 applications to medicine . in fact , designing 
 of an enzyme with a novel function can be considered as the best manifestation 
 of the understanding of enzyme catalysis and enzyme evolution . however , 
 the present generation of designers enzymes are much less efficient 
 than naturally evolved enzymes . the problems with current 
 efforts of rational design is most likely due to an incomplete modeling 
 of the transition state ( ts ) in the enzyme active site , and in particular 
 to the limited awareness to the key role of the reorganization energy . thus , an effective enzyme design methodology 
 should be judged by its ability to determine the activation free energy , 
 along with firm understanding of the factors governing the change 
 in the ts energy in directed evolution experiments . the challenges 
 of modeling enzymatic transition states is far from 
 trivial as it requires both , extensive sampling and reliable potential 
 surfaces . here perhaps the most effective option is the use of the 
 empirical valence bond ( evb ) . the evb is a semiempirical quantum mechanics / molecular 
 mechanics ( qm / mm ) approach , where the qm part is 
 represented by empirical approximations of the relevant valence bond 
 integrals . the evb has been successfully 
 used in reproducing and predicting mutational effects , as well as in quantitative screening of design proposals 
 and in reproducing observed effect of directed evolution refinement 
 of kemp eliminases . this type of approach is in principal 
 effective , but at present it involves major difficulties in obtaining 
 reliable free energies by sampling the surfaces obtained with high 
 level ab initio methods . some effective options like paradynamics method can help in this respect . in considering the evb as an effective tool for computer - aided 
 enzyme design , it is useful to note that this approach has reproduced 
 reliably the observed activation barriers for different mutants of 
 trypsin , dihydrofolate reductase and kemp eliminase . nevertheless , it is important to further validate the evb approach 
 with newer sets of designed enzyme and different types of active sites . 
 in this work we will focus on a designed mononuclear zinc metalloenzyme , 
 which catalyzes the hydrolysis of a model organophosphate . the design of this metalloenzyme started from 
 adenosine deaminase with was manipulated by a denovo methodology with the aim of generating an enzyme that can 
 catalyze the hydrolysis of an organophosphate . as in other previous cases , the most effective steps in the refinement 
 were achieved by directed evolution experiments that mimic natural 
 evolution by selecting mutations which are beneficial to the overall 
 catalytic activity of an enzyme . thus , 
 studies of this designed enzyme give us both an opportunity to validate 
 our approach on metalloenzymes , and provide ( at least in principle ) 
 the opportunity to study an evolutionary trajectory where enzyme evolves 
 to perform a completely new function . as stated above , the enzyme 
 chosen for this study is a designed mononuclear zinc metalloenzyme , 
 which catalyzes hydrolysis of diethyl 7-hydroxycoumarinyl phosphate 
 ( decp ) ( figure 1a ) ( mimicking organophosphate 
 nerve agents ) . this enzyme was designed 
 from adenosine deaminase which is a mononuclear zinc metalloenzyme , 
 where metal ion is thought to be primarily acting as an activating 
 agent for a hydroxyl ion nucleophile . the first variant that was found to show detectable 
 activity ( kcat / km ) contains 
 eight mutations ( designated as pt3 ) . three other variants , pt3.1 , 
 pt3.2 , and pt3.3 , in the evolutionary trajectory were found to have 
 activities of ( kcat / km , m s ) of 4 , 154 , 959 , and 9750 , respectively , and kcat ( 10 s ) of 
 5  10 , 0.2 , 4 , 47 , and 351 , respectively . ( a ) . chemical 
 structure of diethyl 7-hydroxycoumarinyl phosphate 
 ( decp ) . evolutionary trajectory of the decp hydrolysis activity . in order to verify our ability 
 to reproduce the results of the 
 directed evolution experiments , we have simulated the activation barriers 
 for the hydrolysis of decp by adenosine deaminase and its four variants 
 ( pt3 , pt3.1 , pt3.2 and pt3.3 ) ( figure 1b ) . the calculations used as starting points of the crystal structures 
 of adenosine deaminase ( pdb i d : 1a4l ) and pt3.1 ( pdb the structures of pt3 , pt3.2 and pt3.3 were created by generating 
 the corresponding mutations in the crystal structure of pt3.1 ( pdb 
 i d : 3t1 g ) . dealing with a designed enzyme , which performs a 
 new function , 
 it is imperative to examine the binding mode and reaction mechanism 
 of the given substrate . the structure of the substrate  enzyme 
 complex was created with the help of the model of the enzyme  ts 
 complex used in a previous study and with 
 the help of docking using autodock 4 . of course , the configurations generated by autodock were subsequently 
 subjected to extensive relaxation in the evb calculations . figure 2 depicts both the enzyme ( with the bound substrate ) 
 and a representative schematic of the binding modes of the ligands 
 around the zn metal . as seen from the figure the activated hydroxyl 
 and phosphoryl oxygen both coordinate with the zn metal and occupy 
 two out of six coordinating positions in octahedral structure around 
 the metal center . all five systems were subjected to relaxation through 
 100 ps of molecular dynamics simulation prior to the evb calculations . the first step in establishing the mechanism of the phosphotriester 
 hydrolysis requires information about the free energy surface of the 
 reference solution reaction . such information is best obtained by ab initio calculations , but 
 performing such calculations in a fully consistent way can be very 
 challenging ( see for example ) . thus , 
 considering the fact that our main interest is in mutational effect 
 rather than the whole catalytic effect , we only performed qualitative ab initio mapping of the relevant surface . this was done 
 using the gaussian03 software package , and modeling for simplicity a phosphotriester , which has only methoxy 
 groups attached to the phosphorus center . all structural optimizations 
 and energy evaluations of the ab initio potential 
 energy surface were performed using the 6 - 31++g * * basis set with the 
 b3lyp hybrid density functional . the solvent was treated with the 
 cosmo implicit solvent model and the 
 resulting free energy surface is given in figure 3 . obviously , more systematic calculations are needed , including 
 careful qm / mm studies such as those preformed in our recent studies and a systematic examination of the effect of 
 the zn ( that would require to include the ligands in the ab initio 
 calculations ) . however , our main point here is that the present surface 
 and other studies , as well as studies 
 of related systems , indicate that we 
 have a high - energy reaction intermediate or a plateau that leads to 
 a mechanism of the type described in figure 4 . this mechanism justifies the use of three states evb surface , which 
 will be discussed below . at any rate , since our effort is dedicated 
 to exploring mutational effects , we left further exploration of the 
 solution surface and its modification by the zn ion to subsequent 
 studies . ab initio energy surface for the hydrolysis of 
 dimethyl coumarinphosphate in solution ( onuc and olg designate , respectively , the nucleophilic oxygen and the 
 leaving group oxygen ) . schematic description of the energetics of stepwise hydrolysis 
 of decp . as in previous 
 studies , we performed our simulations 
 using the evb method , which is described in great details elsewhere . in setting out the evb potential , we note that the stepwise hydrolysis 
 of decp can be studied by using a three state evb description of the 
 system ( see figure 5 ) . this behavior can be described by an analytical expression1where g is the actual 
 activation barrier , max ( x;y ) is 
 the maximum of the two variables x and y and g2 = 
 g23 + g12 . although we evaluate the activation barriers , gij , by the full evb calculations we note that they can 
 also be estimated by our linear free energy relationship ( lfer ) expression2where ij , gij , 
 and hij are reorganization 
 energy , 
 reaction energy , and the off diagonal mixing term , respectively . the 
 effect of the specific environment is incorporated by taking into 
 account the changes in the corresponding by the reorganization energies 
 and/or by changing the value of gij . schematic description 
 of the three state evb model used to describe 
 the hydrolysis of decp . in calibrating the evb potential , we did not try to use the 
 observed 
 energetics of the reaction of hydroxide attack in water since it does not include the effect of the 
 zn ion . furthermore , here we are interested in the mutational effects , 
 rather than the catalytic effect relative to water . thus , we calibrated 
 the surface taking the reaction in 1a4l as the reference reaction , where we fixed 
 the rate determining barrier around 27.5 kcal / mol , while assuming 
 that in 1a4l this barrier is g2 ( this assumption 
 is based on the group contribution reported below ) . we also selected 
 a value of 20 kcal / mol for g12 in 
 1a4l ( see below ) and took for g23 in 1a4l a 
 value of 6 kcal / mol . we note in this respect that the results 
 do not depend strongly on the values of g12 and g23 . that is , as 
 can be seen by using eq 2 the trend in the rate - determining 
 barrier for different mutants does not depend the corresponding g23 ( with a proper adjustment of h23 ) and this is also partially true with regards to g12 , since we can adjust h12 to obtain the same trend . the real uncertainty is in determining 
 whether the second barrier is rate determining and at what point the 
 first barrier starts to be rate limiting ( the change in the lfer ) . 
 resolving this issue requires lfer experiments or very careful pd 
 calculations . thus , the decision on the point of change in the lfer 
 is somewhat arbitrary at the present case . at any rate , the evb calculations were performed with the molaris program in conjunction with enzymix force field . the evb activation barriers were estimated at 
 configurations selected by the same free energy perturbation umbrella 
 sampling ( fep / us ) approach described extensively elsewhere . the simulation systems were solvated by the surface constrained 
 all atom solvent ( scaas ) model , with 
 a water sphere of 18  radius around the substrate and surrounded 
 by 2  grid of langevin dipoles followed by a bulk solvent . the 
 long - range electrostatic effects were treated by the local reaction 
 field ( lrf ) method . the fep mapping 
 was evaluated by 21 frames of 20 ps each for moving along the reaction 
 coordinate using scaas model . all the simulations were performed at 
 300 k with a time step of 1 fs for integration . in order to obtain the 
 contributions from each residue to the activation barrier ( the group 
 contributions ) were estimated by calculating the effect of change 
 of substrate charges ( from rs to ts ) on the electrostatic contribution 
 of each protein residue . as discussed in our previous studies ( e.g. , 
 ref ( 6 ) ) , the electrostatic 
 contributions of all the protein residues to the activation barrier 
 can be estimated by the following expression:3here the 332 factor is the conversion to kcal / mol , qj are the 
 residual charges of the protein atoms in atomic units ( j runs over the protein residues and k runs 
 over the atoms of the jth residues and i over the substrate atoms ) , ri , k(j ) is the distance in a between the kth atom of the jth group and the ith atom of the substrate , 
 ij is the effective dielectric 
 constant for the specific interaction , and qi are the changes in the substrate charges 
 upon going from the rs to ts . decomposing this expression to the individual 
 group contributions allows one to explore the approximated 
 effect of mutating ionized or polar residues . accurate estimation of the catalytic 
 effects of the different enzyme 
 construct / mutants can be considered as the most basic requirement 
 for the effective enzyme design or understanding to evolutionary mechanism . 
 therefore , we started with systematic evaluations of the activation 
 barriers for our systems . our typical procedure of obtaining activation 
 barrier involved average over 5 free energy profiles , for each enzyme 
 variant ( mutant ) . the details of the calculations are summarized in 
 table s1 ( supporting information ) and the 
 estimated barriers are summarized in table 1 and figure 6 ) . the correlation between the calculated and observed 
 activation 
 barriers ( table 1 and figure 6 ) suggests that change in activity is driven by the change 
 in transition state binding and not by some other elusive factors 
 ( such as substrate binding or dynamics ) . the successful demonstration 
 of our ability to estimate accurate activation energies also indicates 
 that the binding mode of substrate and the reaction mechanism used 
 are reasonable . it should be noted that this is a designed enzyme , 
 and therefore , no concrete prior information about the binding mode 
 or reaction mechanism is available . we believe that rational 
 enzyme designing procedure can be improved 
 if we can quantify the contribution of each residue to the transition 
 state binding . considering the fact that the electrostatic interaction 
 is by far the most important factor in transition state stabilization 
 and therefore enzyme catalysis , we have calculated the electrostatic 
 group contributions of the protein residues . this was done , as discussed 
 in section ii.4 , by using eq 3 and collecting the contribution of each residue to the overall 
 sum ( namely the electrostatic contribution for the energy of moving 
 from the reactant to transition state ) . specifically , we have ( artificially ) 
 changed the charge of protein residues of 1a4l ( the  wild type  ) 
 from 0 to 1 , and then calculated the change in corresponding 
 group contribution upon change of the residual charges of the reacting 
 substrate . as can be seen from figure 7b , the 
 contributions of residues19 and 296 to the rate limiting c  olg bond dissociation step,g2 , are positive 
 ( note as is clear from the supporting information that figure 7a is for a barrier that does 
 not correspond to the rate limiting step ) . thus , changing the charges 
 of the corresponding residues from 1 to 0 should lead to a 
 reduction in g2. this is consistent with the 
 finding that removing the charges of d19 
 and d296 ( the d19s and d296a mutations ) in 1a4l is necessary for effective hydrolysis 
 of decp . we focus here on these two mutations since they are well - defined 
 experimentally observed electrostatic mutations . in principle we can 
 use the group contributions for further predictions but this is not 
 the purpose of the present work , since these contributions are much 
 less reliable than those obtained from evb calculations when they 
 involve residues near the substrate . the group 
 contributions should be , however , very useful for the small contributions 
 of distanced ionized residues , and exploring this point is left to 
 subsequent studies . group contributions ( in kcal / mol ) for ( a ) the nucleophilic 
 attack 
 and ( b ) the bond dissociation steps in 1a4l . the group contributions 
 reflect the interactions between the changes in the charge of protein 
 residues from 0 to 1 , with the charge change of substrate 
 upon moving from rs to ts1 and ts2 . the relatively large positive 
 contributions provide a rough guide for the optimal sites for effective 
 mutations that would enhance the catalytic effect . since the second 
 step is rate limiting in 1a4l , the corresponding group contributions are those that 
 should be compared to the observed results . the ability to accurately 
 estimate the activation energy of different 
 variant enzyme of an enzyme can significantly improve the effectiveness 
 of enzyme design efforts . at present , most enzyme design methods rely 
 on directed evolution experiments to refine and increase the activity 
 of the designed enzyme . in principle , in silico procedures 
 can help in increasing the activity of designers enzymes by accurately 
 estimating the effect of proposed mutations on the rate determining 
 activation energies . gas phase calculations or calculations which 
 explicitly focus on the electrostatic interaction between the protein 
 residues and the ts are very unlikely to have success in estimating 
 the activation barriers as they do not consider the surrounding environment 
 and its reorganization during the reaction . in principle however , the difficulties of obtaining converging free energy calculations 
 make it hard to use such methods in accurately estimating mutational 
 effects . on the other hand , the evb has been shown to be capable 
 of estimating 
 the effect of mutational change on activation as early as 1986 , where computer - aided mutations were proposed 
 for rat trypsin . as far as enzyme design is concerned , we like to 
 point out that evb has been shown to be capable of reproducing the 
 effect of mutations observed in directed evolution of kemp eliminases . however , more studies are clearly needed and 
 thus we have extended here the validation of the evb to a study of 
 the effects of multiple mutations on the activity of a designed zn 
 metalloenzyme . in doing so we note that the relatively high reactivity 
 of metalloenzyme , coupled with the wide range of reactions carried 
 out by them , makes them very attractive starting points for introducing 
 new activities . at any rate , in the present study , we have successfully 
 estimated the activities of different variants of the designed metalloenzyme 
 and have reproduced the evolutionary trajectory leading to a new catalytic 
 function ( hydrolysis of decp ) . while determining the effect 
 of different mutations on activation 
 energies is the key to effective rational design , it would be useful 
 to have a qualitative guide to propose mutations which can decrease 
 the activation energy and therefore can increase the catalytic activity . here we provide indications that the electrostatic group contributions 
 can provide an important lead for mutations , which can improve the 
 activity of an enzyme . in particular the group contributions in 1a4l reproduced the experimental 
 trend that mutations that remove the negative charges at position 
 asp19 and asp296 increase the activity . directed evolution has 
 emerged as a powerful approach that can 
 offers an effective way of optimizing enzyme activity . however , at 
 present such approach has not achieved the same impressive catalytic 
 power on enzymes that evolved by natural evolution . overcoming this 
 limitation will require exploration of mutational trajectories beyond 
 what has been suggested by directed evolution . the evb can be very 
 useful in advancing such studies . despite the encouraging results 
 of the present study it is important 
 to mention that we did not performed a sufficiently careful study 
 of the reference solution reaction or the effect of the zn ion and 
 its ligands and used relatively tentative estimates in estimating 
 the reference surface in 1a4l . to further advance in this direction it would be 
 crucial to preform ab initio qm / mm ( qm(ai)/mm ) free 
 energy calculations for the solution reaction with and without the 
 zn ion as well as pd calculations of the reaction in the enzyme . it 
 will also be very useful to have lfer experiments for the reaction 
 in the enzyme . this will help in reducing the uncertainties about 
 the rate determining barriers and increases the ability to make quantitative 
 predictions of mutational effects .
one 
 of the greatest challenges in biotechnology and in biochemistry 
 is the ability to design efficient enzymes . 
 in fact , such an ability 
 would be one of the most convincing manifestations of a full understanding 
 of the origin of enzyme catalysis . despite some progress on this front , 
 most of the advances have been made by placing the reacting fragments 
 in the proper places rather than by optimizing the preorganization 
 of the environment , which is the key factor in enzyme catalysis . 
 a 
 rational improvement of the preorganization and a consistent assessment 
 of the effectiveness of different design options require approaches 
 capable of evaluating reliably the actual catalytic effect . in this 
 work 
 we examine the ability of the empirical valence bond ( evb ) to 
 reproduce the results of directed evolution improvements of the catalysis 
 of diethyl 7-hydroxycoumarinyl by a designed mononuclear zinc metalloenzyme . 
 encouragingly , our study reproduced the catalytic effect obtained 
 by directed evolution and offers a good start for further studies 
 of this system .