18 Berendsen applied to single and multiple proton transfer reactions by Hammes-Schiffer and collaborators [119,120,121,122,103]. It seems that surface hopping (also called Molecular Dynamics with Quantum Transitions,MDQT)is a rather heavy tool to simulate proton dynamics.A recent and promising development is path integral centroid dynamics [123]that provides approximate dynamics of the centroid of the wavefunctions.Several improvements and applications have been published [123,124,125,126127,128]. We have applied the DME method to several proton transfer processes in aqueous solution (hydrogen malonate [111,114],HIV protease [129,130], and a slow neutral ester hydrolysis [131).The latter case is a slow reaction with a rate constant of about once per minute,where the rate-limiting step is the transfer of a proton from a water molecule (which is in close proximity to the ester carbon,see Fig.la)to another water molecule,concerted with electron redistribution that leads to the tetrahedral intermediate (Fig.1b). We used a combination of semi-empirical QM(to obtain the energy of the reactant state as a function of the distance between ester carbon and water oxygen and to obtain partial charges to be used in MD for intermediate states),thermodynamic integration(to obtain the free energy difference in solution between the equilibrium reactant state and the activated state from which proton transfer is probable),and DME (to obtain the initial rate of proton transfer).For the latter the proton potential along the transfer path was calculated during a biased MD simulation. This proton potential,which is rapidly and heavily fluctuating (Fig.2), drives the evolution of the proton wave function.The rate of transfer was then deduced from the coarse-grained average increase of proton population at the reactant side.Only the initial rate can be determined this way;this has the advantage that the protonic motion does not yet 'backreact'on the envi- ronment,but it has the disadvantage that no information is obtained about recrossing events.A good agreement with experimental rates was found,and a deuterium isotope effect of 3.9 was simulated,compared to the experimental value of 3.2.The picture is that proton transfer is driven by fluctuations due to solvent motion which transiently provide good tunneling conditions;the probability of transfer is however almost always negligibly low except when the carbon-water distance happens to be favourable for the reaction.The free energy of activation to reach this favourable state is one of the main reasons why the reaction is so slow.While the method is still crude,it can provide us with good estimates for reactions rates,including enzyme catalysis,when proton transfer is the rate-limiting step. 3.3 Reduction of Complexity:Averaging over Degrees of Freedom Substantial headway towards longer time scales and larger systems can only be expected from reduction of system complexity.It is here where future
MD:The Limits and Beyond 19 H CCHCl2 CH3O (a)Reactant state 日 O±-H H H CHC12 CH3O (b)Product state Fig.1.The rate-determining step in the neutral hydrolysis of paramethoxy-phenyl dichloroacetate.In the reactant state (a)a water molecule is in proximity of the carbonyl carbon;after concerted proton transfer to a second water molecule and elcctron redistribution,a tetrahedral intermediate (b)is formed
20 Berendsen 30.0 20.0 10.0 0.0 -10.0 -2000.0 20.0 40.0 60.0 80.0 t(ps) Fig.2.The fluctuating difference between the proton potential at the product side relative to that at the reactant side (the difference between the two wells in a double-well proton potential).Whenever this difference is close to zero,tunneling conditions are favourable. innovations are expected to be most fruitful.All such reductions concern in one way or another the omission of the explicit description of a (large) fraction of the degrees of freedom in the system. The first requirement is the definition of a low-dimensional space of 're- action coordinates'that still captures the essential dynamics of the processes we consider.Motions in the perpendicular null space'should have irrele- vant detail and equilibrate fast,preferably on a time scale that is separated from the time scale of the essential'motions.Motions in the two spaces are separated much like is done in the Born-Oppenheimer approximation.The average influence of the fast motions on the 'essential'degrees of freedom must be taken into account:this concerns (i)correlations with positions ex- pressed in a potential of mean force,(i)correlations with velocities expressed in frictional terms,and (i)an uncorrelated remainder that can be modeled by stochastic terms.Of course,this scheme is the general idea behind the well-known Langevin and Brownian dynamics. In special cases(as in colloidal solutions)some particles can be consid- ered as 'essential'and other particles as 'irrelevant',but in most cases the essential space will itself consist of collective degrees of freedom.A reaction coordinate for a chemical reaction is an example where not a particle,but some function of the distance between atoms is considered.In a simulation of the permeability of a lipid bilayer membrane for water [132]the reaction coordinate'was taken as the distance,in the direction perpendicular to the bilayer,between the center of mass of a water molecule and the center of mass of the rest of the system.In proteins (see below)a few collective de- grees of freedom involving all atoms of the molecule,describe almost all the
MD:The Limits and Beyond 21 macromolecular motion.One can also consider local densities of molecular species or molecular groups as essential degrees of freedom;in that case local fux densities are considered that respond to gradients of the thermodynamic potential (see below). 30.0 25.0 20.0 15.0 10.0 5.0 0.0 0.0 0.9 1.8 2.7 3.6 4.5 z(nm) Fig.3.The potential of mean force for a water molecule(open symbols)and an NH3 molecule (filled symbols)penetrating into a lipid bilayer membrane.The aqueous phases are on both sides (regions 1),the headgroups are predominantly in regions 2,regions 3 contain the most ordered part of the hydrocarbon phase,and region 4 represents the middle of the hydrocarbon phase containing the end groups.The large circles with error bars were determined by integrating average constraint forces;the smaller points in regions 1 and 2 were determined directly from the observed water density and the smaller points in the middle were obtained by particle insertion (Marrink et al.,J.Phys.Chem.98 (1994)4155 and ibid.100 (1996)16729). Potentials of Mean Force In principle the potential of mean force in which the essential'coordinates move can be determined from detailed simulations including all degrees of freedom,in which the essential degrees of freedom are either treated as constraints or are restrained with an umbrella potential. In simple,low-dimensional cases this actually works:for example,in the simulation of the transport of water and other small molecules through a lipid bilayer [132,133]the potential of mean force for a small molecule as function of the depth z in the membrane could be derived from integration of the mean force acting on a molecule constrained at a given z,using several constrained MD runs (see Fig.3).The results could be checked in this case
22 Berendsen by particle insertion in the middle of the bilayer.From analysis of the force fuctuations in the same simulations it was possible to derive the diffusion coefficient in the z-direction as a function of z,thus yielding a complete one- dimensional picture of the transfer path and allowing the computation of the permeation coefficient.A similar procedure was applied to a single file of hydrogen-bonded water molecules [134. However,in many applications the essential space cannot be reduced to only one degree of freedom,and the statistics of the force fluctuation or of the spatial distribution may appear to be too poor to allow for an accurate determination of a multidimensional potential of mean force.An example is the potential of mean force between two ions in aqueous solution:the momentaneous forces are two orders of magnitude larger than their average which means that an error of 1%in the average requires a simulation length of 108 times the correlation time of the fluctuating force.This is in practice prohibitive.The errors do not result from incorrect force fields,but they are of a statistical nature;even an exact force field would not suffice. Thus one must rely on macroscopic theories and empirical adjustments for the determination of potentials of mean force.Such empirical adjustments use free energy data as solubilities,partition coefficients,virial coefficients, phase diagrams,etc.,while the frictional terms are derived from diffusion coefficients and macroscopic theories for hydrodynamic interactions.In this whole field of enquiry progress is slow and much work(and thought!)will be needed in the future. Essential Dynamics and Rigid Bodies in Proteins Proteins happen to be a very special kind of polymer,constructed such that they fold into specific structures,necessary for specific functions.But they are also flexible in a functional way,because functions as substrate binding,catalysis,allosteric effects and regulation require specific dynamical properties.It is therefore not too surprising that proteins are found to exhibit a collective motions with large amplitude in a low-dimensional subspace,described by 10 to 20 collective degrees of freedom,while the other 10,000 degrees of freedom seem to perform insignificant,low-amplitude and gaussian-distributed fluctuations [135,137. The 'essential degrees of freedom'are found by a principal component analysis of the position correlation matrix Cij of the cartesian coordinate displacements xi with respect to their averages (xi),as gathered during a long MD run: C=(ci-()(xj-(cj）)》 Each eigenvalue of this matrix indicates the contribution to the total fluctuation of the collective degree of freedom given by the corresponding eigenvector.If eigenvectors are ordered in descending order,it is generally observed [135,136,138]that the first,20 eigenvalues contribute some 80%of