qm / mm 
 methods ( for reviews , see , e.g. , refs ( 17 ) ) have been used extensively in studies of enzymatic reactions . in 
 order to obtain reliable results , it is crucial to perform extensive 
 sampling and also to use a reliable qm 
 model for either calibrating a semiempirical model ( e.g. , the evb 
 model ) or actual determination of ab initio 
 ( ai ) qm ( ai)/mm free energy surfaces . in considering the reliability of such methods , one of the obvious 
 issues is the size of the qm region .
the size problem is frequently 
 attributed to the treatment of the boundaries between the qm and mm 
 regions , but this issue is not so serious if one uses the same models 
 in the reference water reaction and in the enzyme active site .
nevertheless , 
 one would like to know the limitations imposed by the size of the 
 qm region .
here we noted that the great advances in computer power 
 and the use of specialized computers ( e.g. , gpu ) allows one to use 
 very large qm regions and the results of the corresponding calculations 
 can be used to imply that it is essential to use large qm regions 
 to obtain reliable results .
thus , it is important 
 to conduct careful analysis of the size dependence of qm / mm calculations .
the issue of the size dependence in qm clusters has been a subject 
 of significant interest .
similarly , the size dependence of the qm region in qm / mm models 
 has been the subject of recent studies .
some studies explored the 
 size dependence of qm / mm calculations of spectroscopic properties .
more relevant works ( from the perspective of this paper ) examined 
 the energy difference in enzymatic reactions . however , as much as enzymatic reactions are concerned , 
 the most important issue is the reliability of the calculated activation 
 free energy , that should reflect an extensive averaging , that has 
 not been performed in most of the previous studies .
one exception 
 is the study of ref ( 26 ) , which evaluated the free energy of proton transfer ( pt ) in lysozyme 
 for qm regions of different sizes .
this work found that for a qm with 
 a radius of less than 6  we can have errors of about 3 kcal / mol . 

here we note that pt reactions involve a rather small change in charge 
 separation and one can expect a larger dependence on the size of the 
 qm region in reactions that involve a larger change in charge distribution . 

thus , we focus in this work on the dependence of the calculated activation 
 free energy on the size of the qm region in the reaction of the enzyme 
 comt ( catechol - o - methyl transferase ) that catalyzes 
 the transfer of a methyl group from sam ( s - adenosyl - l - methionine ) to a catecholate ion via an sn2 type 
 reaction , as shown in scheme 1 . the sam unit is truncated to the methyl group for clarification .
the issue of the size of the qm region in comt 
 was highlighted recently in a study which implied 
 that the inclusion of large protein residues ( up to 500 atoms ) in 
 the qm region is needed in order to obtain convergent results . while 
 ref ( 11 ) only looked 
 at the convergence of ground state properties including the donor 
acceptor 
 distance , ref ( 12 ) examined 
 the trend in activation energy , but this was done by evaluating the 
 single point energy change upon moving to the transition state ( rather 
 than the activation free energy ) .
nevertheless , such results may discourage 
 further use of qm / mm calculations with medium qm size and thus require 
 further systematic studies that will be performed in this work . in choosing comt as our benchmark
it belongs to the general class of methyltransferases 
 that plays an important role in human physiology and several diseases 
 such as cancer , and genetic diseases .
comt degrades catecholamines 
 such as dopamine and has thereby attracted a major interest in its 
 probable role in dopamine related pathologies , in particular parkinson s 
 disease .
v158 m comt polymorphism is speculated to be responsible for 
 the violent behavior of schizophrenic patients .
it is therefore imperative to study such enzymes and understand 
 their function at the molecular level .
it should also be noted that 
 the comt catalyzed decomposition of dopamine follows a parallel pathway 
 to the reaction catalyzed by mao b , which was recently studied using 
 the evb method .
hopefully , the present 
 study can also be useful by shedding additional light on the results 
 and in some cases limitations of previous qm / mm studies of comt .
our study focuses on the dependence of the activation free 
 energy of the reaction of comt on the size of the qm region .
for very 
 small sizes where , for example , the mg is treated classically , 
 we find an obvious deviation from the converging results ( although 
 this problem can be fixed by semiempirical approaches with proper 
 parametrization of the qm / mm interactions ) . for larger systems ,
our study used 
 the semiempirical pm6 method ( evaluated by the mopac2009 package ) to treat the qm region , whereas the rest of the protein was treated 
 classically , using the enzymix force field .
the qm / mm interface between the qm and mm surfaces was implemented 
 through the corresponding module of the molaris package .
the initial coordinates for the qm / mm dynamics 
 were taken from the crystal structure of the human , soluble form of 
 comt ( s - comt ) .
the first step for each 
 system involved a relaxation by slowly heating the protein from 30 
 to 300 k for at least 45 ps .
the subsequent md simulations involved 
 time steps of 1 fs at 300 k. the simulation was divided into 31 frames 
 with each frame simulated for 10 ps .
( o  c ) , was constrained 
 using a force constant of 100 kcal / mol  and varied 
 from 1.0 to 2.2  .
the pmf was obtained by combining 
 the individual simulation windows , using the weighted histogram analysis 
 method ( wham ) approach .
further calculations 
 for a few systems ( i , ii , and iii ) were done for a longer time , by dividing the simulation into 51 
 frames and running each window for 30 ps .
this yielded similar results 
 to the shorter runs , and therefore , for all other systems , we only 
 carried out shorter runs .
the boundary conditions for the simulations 
 were introduced by immersing the protein in an 18  sphere of 
 water molecules using the surface - constraint all - atom solvent ( scaas ) 
 boundary conditions .
the local reaction 
 field method ( lrf ) which is similar to having no cutoff for electrostatic 
 interactions is used to treat the long - range effects .
the geometric center of the reacting system was taken as 
 the center of the simulation sphere .
the positions of all atoms beyond 
 20  from the center of the systems were fixed as found in the 
 initial crystal structure , and their nonbonded interactions with the 
 atoms within the simulation sphere were turned off . to determine the 
 protonation state , we used our coarse grained ( cg ) model .
subsequently , all residues that lie within 12 
  from the center of the sam and catecholate anion were explicitly 
 ionized during the simulation . a large charge 
charge dielectric 
 constant was used to treat the effect of the more distanced ionizable 
 residues .
three different starting configurations were generated for 
 the simulations by taking structures from a long trajectory at a time 
 difference of 45 ps .
the resulting structures were used in the qm / mm 
 simulation , and the final result was taken as the average of the three 
 runs .
obviously , more correct barriers should be 
 obtained with high level qm methods , but the issue here is the size 
 dependence of the barrier and not obtaining an accurate barrier .
of 
 course , we hope that the size dependence of the pm6 is related to 
 the size dependence in more rigorous models .
our study focused on the activation barrier for the methyl transfer 
 step , which is depicted in scheme 1 .
the different residues considered as a part of the 
 qm region are shown in figure 1 .
the active site consists of mg ion coordinated 
 to the catecholate ion , two aspartate residues , asparagine residue , 
 and a water molecule .
the transition state with the mg ion and the ligands in the primary shell is shown in figure 1c for a representative case .
( b ) a schematic representation of the active 
 site of comt enzyme with some important residues .
( c ) the transition 
 state for the methyl transfer reaction between sam and catecholate 
 ion for a representative case ( distances are in  ) . in order to explore the size dependence of the 
 qm region
the smallest region consisted of the substrates alone , 
 while the largest region included the metal and seven residues , which 
 were selected on the basis of the distance from the substrates .
the 
 first qm region investigated ( system i ) is described 
 in figure 2 .
this system 
 includes only the two substrates , sam and the catecholate ion , and 
 consists of 26 atoms in the qm region including the linked atoms . 

the mg metal ion was treated classically as a point charge 
 with a repulsive van der waals center .
system i with only the substrates 
 ( catecholate and sam ) as the qm region and system iii with the substrates , mg ion , and the surrounding ligands 
 as the qm region .
the qm part is shown as sticks and part of the mm 
 region as wires ( the rest of the protein is not shown for clarity ) .
the next system considered ( system ii ) also included only the substrates in the qm region .
however , 
 the mg was treated using the dummy atom approach developed 
 by qvist and warshel .
this approach , 
 which involves the inclusion of six dummy atoms around the metal center , 
 as depicted in figure 3 , was shown to yield improved solvation energies .
md6 model complex used instead of mg in the 
 classical simulation of system ii , where the ds represent 
 the dummy atoms and the ls are the general ligands around mg .
subsequently , the mg ion and the surrounding ligands ( h2o , asp169 , asp141 , 
 and asn170 ) were included in the qm region ( system iii shown in figure 2 ) .
it should 
 be noted that only a part of the residues as shown in figure 2 was included in the qm region . 

inclusion of the complete residue did not alter the result significantly , 
 and hence , in all further calculations , we used chopped residues as 
 the ligands coordinated to the mg ion .
the change in 
 the barrier upon moving to system iii reflects qm charge 
 transfer effects , which are discussed below .
next , we considered 
 the effect of individual amino acid residues near the active site 
 by including them separately in the qm region , while still retaining 
 the residues in system iii ( figure 4 ) .
the different residues considered are 
 glu199 ( system iv ) , lys144 ( system v ) , tyr68 
 ( system vi ) , trp143 ( system vii ) , and met40 
 ( system viii ) . lys144 and
glu199 form h - bonds with the 
 catecholate ion and result in a charged qm system .
the free energy 
 barriers obtained for systems iv and v were 
 16.4 and 16.5 kcal / mol , respectively .
during the simulations , the 
 hydroxyl proton at the catecholate shuttles between the o atom of 
 the catecholate and the o atom of the glu residue .
as will be discussed 
 below , this corresponds to another reaction path and should be considered 
 accordingly .
the inclusion of tyr 68 ( which is considered to be an 
 important residue as its mutation decreases significantly the activity ) 
 yields a barrier of 15.8 kcal / mol .
the inclusion of met40 and trp143 
 yields activation free energies of 15.9 and 17.1 kcal / mol , respectively . 

we also included in the qm region tyr147 ( via ) and w38 
 ( viia ) , which are approximately 10  away from the 
 center of the substrates .
this yields barriers of 16.4 and 16.9 kcal / mol , 
 which are similar to the barriers in other systems where the residues 
 are much closer to the active site .
the activation free energies in kcal / mol are provided 
 in parentheses . after studying the effect 
 of individual residues
, we included both lys144 and glu199 ( system ix ) together with the mg metal and its ligands 
 ( figure 4 ) .
these two 
 residues lie closest to the substrates forming weak h - bonding interactions 
 and are therefore included in most of the systems described hereafter . 

further , we 
 included three residues together keeping the mg and its surrounding 
 ligands , lys144 and glu199 , constant . in system xiii
, 
 we studied the effect of another residue , trp43 , on the activation 
 barrier .
the different systems created ( figure 5 ) are ( a ) lys144 , glu199 , trp143 ( system x ) ; ( b ) lys144 , glu199 , tyr68 ( system xi ) ; ( c ) 
 met40 , tyr68 , trp143 ( system xii ) ; and ( d ) trp 143 , trp43 , 
 tyr68 ( system xiii ) . in the final and largest qm systems , 
 we included four residues together , namely , lys144 , glu199 ,
trp143 , 
 and met40 , resulting in system xiv and residues lys144 , 
 glu199 , trp143 , and tyr68 that generated system xv .
systems 
 that include in the qm region three residues along with the mg and its surrounding ligands .
the activation free energies 
 in kcal / mol are provided in parentheses . in total , 15 different qm systems have been considered and 
 it has been found that the free energy barrier lies in the range 14.519.4 
 kcal / mol . the largest system including seven residues does not show 
 much improvement over the system containing the mg and 
 ligands alone .
on the basis of our results , we propose that mg and its surrounding ligands ( system iii ) are 
 sufficient to be included in the qm region .
furthermore , we also carried 
 out qm / mm studies for the y68a mutant , where the experimental free 
 energy barrier for the mutant is 20.0 kcal / mol whereas that of the 
 wild type is 18.4 kcal / mol . in this mutant case , we only investigated 
 systems iii and i. it was found that , while 
 the calculated free energy for wild type with system i is 22.5 kcal / mol , the corresponding barrier for y68a is 26.8 kcal / mol 
 ( leading to an increase of 4.3 kcal / mol in the barrier ) . in the case 
 of system iii ,
the barriers for the wt and the mutant 
 are 14.5 and 16.6 kcal / mol , respectively ( leading to a 2.1 kcal / mol 
 increase in barrier ) .
it can be seen that the trend on increase in 
 barrier upon mutation is reproduced by both systems ( where system i gives a significant overestimate ) .
the results of 
 the above study are summarized in figures 6 and 7 as well as table 1 .
the analysis of 
 the results can start with the obvious finding that cases that do 
 not include the ion and its ligands in the qm systems ( cases i and ii ) present a major approximation , since 
 the corresponding qm treatment does involve major charge transfer 
 ( ct ) .
as we pointed out , the inclusion of the mg ion 
 in the qm region requires inclusion of the ligands of this ion , and 
 such a treatment can be very expensive when one uses a reasonable 
 level of the qm treatment . fortunately , treating the mg ion classically and using the proper parametrization can give reasonable 
 results .
the pmfs obtained for different systems for the pm6/mm potential , 
 which have been calculated using the wham procedure .
all systems include the substrates : sam and catechol ( shown in wires ) . 

another case , which requires special attention , is the inclusion 
 of glu199 in the qm system ( in addition to system iii ) without including additional residues . in this case
, we find a 
 partial proton transfer ( pt ) from the catecholate to glu199 .
however , 
 we actually have here a trivial case of another reaction mechanism 
 ( namely , a pt to glu199 ) . of course
, one should examine whether we 
 have a real pt or an artifact of the use of the pm6 method .
in fact , 
 the problem disappears already when we add lys144 or additional residues 
 to the qm region .
we also find that adding lys144 alone to system iii also leads to a deviation in the reaction free energy . 
 now
again we have a residue that is involved in the reaction , where 
 a pt from the catechol to the lys can be considered as an initial 
 step of the reaction ( see also below ) . 
 here a partial pt for the back reaction can appear in small qm regions , 
 but again adding more residues leads to a disappearance of the deviations .
we like to note that , although 
 lys144 is the most likely candidate to accept the proton from catechol 
 leading to the formation of catecholate ion , it is still unclear whether 
 the proton stays on lys144 or is transferred to some other residue 
 or water molecules . the close proximity of this residue to mg ( 3.2  ) in the crystal structure results in this ambiguity , 
 as it may lead to the low pka of these 
 residues .
in fact , the finding that the lys144ala mutant leads to 
 a minor reduction in kcat indicates that 
 the lys does not play a major role in the rate acceleration .
thus , although our previous evb calculations 
 gave a somewhat larger catalytic effect with neutral lys144 , we kept this residue ionized here ( in the cases 
 when it was in the classical region ) .
regardless of the actual protonation 
 state of the lys residue during the nucleophilic attack , our analysis 
 of the size dependence is still valid for the model system used .
overall , we note that , while the convergence of the activation free 
 energies is encouraging , the convergence for the reaction free energies 
 is not as good as that of the activation free energies ( table 1 ) , for single residues that 
 can be involved in side reactions like pt .
nevertheless , once other 
 qm residues were included with the problematic residues , we obtained 
 converging results .
this feature was found for glu199 , lys144 , and , 
 for less clear reason , trp143 , where in all cases adding more residues 
 leads to the disappearance of this problem .
qm / mm approaches provide a major tool for evaluating the energetics 
 of chemical reactions in complex molecular systems .
however , some 
 workers can assume that the results depend drastically on the size 
 of the qm system and thus presumably the qm / mm idea is inherently 
 a major approximation .
apparently , what is meant by  approximation  
 is not uniformly agreed upon .
for example , some may assume that an 
 accurate result can be obtained only when the protein and the solvent 
 are represented by a high level ab initio approach .
however , obtaining 
 reliable results also requires a major sampling and satisfying the 
 requirement of both a proper sampling and a reliable qm approach is 
 an extremely challenging task that would be hard to accomplish in 
 the near future .
there are ways to use large qm systems , including 
 our constraint dft ( cdft ) approach , or the use 
 of a fast processor or specialized computers that allow running with 
 large qm systems
. however , the requirement 
 of extensive sampling makes it crucial to exploit qm / mm approaches 
 with a limited size of the qm region , while exploring systematically 
 what are the limitations of such strategies .
this work explored the 
 approximations associated with the use of a relatively small size 
 for the qm region , by a systematic analysis of this issue while modeling 
 the free energy profile of the reaction of comt .
as seen from figure 8 , the convergence 
 of the activation free energy is quite reasonable ( except for models i and ii , as discussed below , the reasons for 
 those deviations are obvious ) .
convergence of the activation free energy 
 as a function of the number of qm residues .
the numbers in brackets 
 are the number of qm atoms for the given system .
the deviations between the results with different sizes of 
 the qm region are more pronounced when we compare the calculated exothermicity 
 ( table 1 ) . in analyzing 
 the corresponding results
first , we note that the cases that do not include the ion and its 
 ligands in the qm systems present a major approximation , since the 
 corresponding qm treatment does involve charge transfer ( ct ) .
thus , 
 we do have here an obvious case where one should try to include the 
 mg in the qm treatment .
fortunately , the ct effect can 
 be captured by changing the charge on the mg and the 
 ligands and by proper refinement of the van der waals parameter without 
 the use of computationally expensive methods .
the ct can also be represented 
 by including the mg in an evb treatment that can reproduce 
 this effect .
furthermore , even without including the mg in the evb region , we can easily reproduce the trend obtained with 
 the mg in the qm region by a standard evb treatment that 
 includes parametrization of the reaction exothermicity . as stated 
 in section iii ,
the convergence 
 for the reaction free energy is not as good as the activation free 
 energies .
this is in particular true for single residues ( i.e. , glu199 
 and lys144 ) that can be involved in side reactions such as pt reactions . 

however , once such residues are surrounded by other qm residues , we 
 obtain converging results .
such a behavior seems to reflect  instability  
 of the qm calculations , which do not follow a simple additive trend . 

this is reflected in the observation that the problem does not occur 
 when the problematic residues are part of the mm system .
it is also 
 likely that the size dependence will be reduced if we use a polarizable 
 force field .
it should also be noted 
 that our reaction is not a conventional sn2 reaction , as 
 it involves two charged species that react to give a neutral product , 
 thereby posing difficulties that arise in an sn1 reaction . 

thus , the solvation free energy of the reactive ionic species might 
 result in a sensitivity of the reaction free energy .
some workers 
 tend to attribute the main source of the size dependence problem to 
 the connection between the qm and mm link atom .
however , a part of this problem is drastically reduced by using 
 hybrid orbitals .
furthermore , the problem is drastically 
 reduced if one uses the same boundaries for reactions in water and 
 reactions in enzyme as is done with the evb method .
it is important 
 to point out that even with some dependence on the qm size it is possible 
 to explore quite reliably the catalytic and mutational effect .
this 
 issue has been demonstrated here with the wt to y68a mutations in 
 two different qm descriptions .
such a finding , for a mutant that is 
 the subject of significant controversy , is clearly important in view of the possible presumption that one 
 must use large qm regions in exploring mutational effects .
reference ( 12 ) obtained about 3 kcal / mol 
 deviation from the final result for about 10 residues but then suggested 
 that the convergence requires at least 30 residues ( where the deviation 
 from the converging result is about 23 kcal / mol ) .
ignoring 
 the issue that the calculations are based on a single point energy 
 evaluation , we like to address the implication that one needs to get 
 around 2 kcal / mol convergence before exploring catalytic problems . 

first , obtaining 
 2 kcal / mol errors due to the size of the system still allows exploring 
 enzyme catalysis , when the difference between the barrier in enzyme 
 and in water is frequently on the order of 10 kcal / mol .
of course , 
 it is also allowed ( as shown here ) to explore mutational effects . 

second , no first principle approach can at present give errors of 
 less than 2 kcal / mol .
that is , changing the level of the quantum method 
 can lead to several kcal / mol errors and adding a polarizable force 
 field can change the results by a few kcal / mol . obviously not performing 
 sampling can lead to enormous errors .
furthermore , very large errors can also be obtained from evaluation 
 of electrostatic contributions without proper relaxation of the protein . 

such calculations may produce an unphysical low dielectric that would 
 not provide sufficient screening for charge  charge interaction .
apparently , even proper free energy calculations 
 still involve several kcal / mol errors due to incomplete sampling , 
 which , for example , can miss water penetration effects .
the 
 above - mentioned errors are inherent to current calculations , and once 
 we are aware of such errors , we can better assess the effectiveness 
 of different approximations .
for example , the error can frequently 
 be reduced by using a well - defined reference ( e.g. , the wt barrier 
 in mutational studies and the reaction in water in exploring catalysis ) . 
 as usual ,
the best way of exploring different approximations is to 
 change the length of the run , the system size , and other factors and 
 to explore the stability of the calculations as well as the relationship 
 to the observed quantity .
it might also be pointed out here 
 that there are various implications that the catalytic effect deduced 
 from calculations with limited qm size ( mainly the conclusion that 
 the catalysis is associated with electrostatic effects ) might be questionable 
 and can be changed with larger qm sizes .
now , not only have evb calculations 
 seemed to give extremely stable conclusions , but also equally importantly , 
 other proposed catalytic factors ( e.g. , dynamical effects , quantum 
 effects , and more ) have never been reproduced by consistent calculations 
 regardless of the size used .
overall , the present study indicates that the 
 size dependence is not as crucial as one may suspect .
in particular , 
 the size dependence is significantly smaller than the catalytic effect . 
 of course
, the residues that are in direct contact with the reacting 
 atoms can be involved in a major charge transfer or even in another 
 reaction path .
however , the qm description of residues that are not 
 in direct contact do not change the activation barriers of chemical 
 transformation in a major way .
 although qm / mm calculations 
 are the primary current tool for modeling enzymatic reactions , the 
 reliability of such calculations can be limited by the size of the 
 qm region . 
 thus , we examine in this work the dependence of qm / mm calculations 
 on the size of the qm region , using the reaction of catechol - o - methyl transferase ( comt ) as a test case . 
 our study focuses 
 on the effect of adding residues to the qm region on the activation 
 free energy , obtained with extensive qm / mm sampling . 
 it is found that 
 the sensitivity of the activation barrier to the size of the qm is 
 rather limited , while the dependence of the reaction free energy is 
 somewhat larger . 
 of course , the results depend on the inclusion of 
 the first solvation shell in the qm regions . for example 
 , the inclusion 
 of the mg2 + ion can change the activation barrier due to 
 charge transfer effects . 
 however , such effects can easily be included 
 in semiempirical approaches by proper parametrization . 
 overall , we 
 establish that qm / mm calculations of activation barriers of enzymatic 
 reactions are not highly sensitive to the size of the qm region , beyond 
 the immediate region that describes the reacting atoms . 
Introduction
Computational Methods
Results
Discussion
Conclusion