David Case

Frederick R Manby

The Ornstein-Zernike equation in molecular electronic structure theory

Keywords: Electron correlation, IEPA, Coupled cluster, Fluids, Ornstein-Zernike

The Ornstein-Zernike equation expresses a simple relationship between direct and full correlation effects of a many-body system. Although it was introduced in the context of fluids, it can be derived through simple physical arguments that are equally applicable to the electrons in a molecule. Direct correlation effects account for the correlation between two particles at a time. If a simple model for direct correlation can be found, the Ornstein-Zernike equation can be used to convert it to a treatment of the full correlation effect. We show that the independent electron pair approximation (IEPA) is a reasonable model for the description of direct correlation effects. IEPA followed by our Ornstein-Zernike treatment is closer to CCSD(T) than is CCSD for valence correlation energies of small molecules.

Introduction

In the Bakerian lecture of 1869, which incidentally appears never to have taken place [1], Thomas Andrews reported his observations on the liquefaction of gases under high pressure [2]. He had prepared carbon dioxide at or near its critical point and when the temperature or pressure was 'suddenly diminished' he observed 'flickering striae'.

His flickering striae were a manifestation of critical opalescence, and the phenomenon was later described by Smoluchowski [3] and Einstein [4] in terms of spontaneous density fluctuations in the fluid. In this early work 'mutual independence of the elements of volume' was assumed [START_REF] Ornstein | Proc. Akad. Sci[END_REF]: the theoretical apparatus needed to treat correlations between the density fluctuations was devised by Ornstein and Zernike [START_REF] Ornstein | Proc. Akad. Sci[END_REF].

Consider a simple (atomic), homogeneous fluid with density ρ. The radial distribution function g(r) expresses the extent to which the density is modified by correlations; so if there is a particle at position r 1 , the density at r 2 will be ρg(|r 1r 2 |) ≡ ρg 12 . In the long range the correlations die out, so g(r) → 1 as r → ∞. It is convenient to consider the shifted quantity h = g -1, called the total correlation function [START_REF] Hansen | Theory of Simple Liquids[END_REF].

In the Ornstein-Zernike (OZ) equation the total correlation function is related to a new quantity c 12 , the direct correlation function:

h 12 = c 12 + ρ dr 3 c 13 h 32 . (1)
The equation expresses the idea that the total correlation between particles arises † Corresponding author: fred.manby@bris.ac.uk primarily from the direct correlation effect, but with a correction arising from higher-order effects that bring in more and more additional particles. This leads to the interpretation of direct correlation as the correlation effect between just two particles in a homogenous background density. With this in mind, the OZ equation can be 'derived' based on the following physical arguments.

(1) Assume there is a particle at position r 1 (2) The direct correlation effect modifies the density at r 2 by an amount ρc 12 (3) The density at r 3 is modified as well, by the amount ρc 13 (4) This fluctuation at r 3 has an additional effect on the density at r 2 , with magnitude ρc 13 ρc 32 (5) Particle 3 can be anywhere, and after integrating over r 3 and cancelling one power of ρ, one obtains h 12 ≈ c 12 + ρ dr 3 c 13 c 32 [START_REF] Hansen | Theory of Simple Liquids[END_REF] The entire argument can be repeated with more intermediate particles, giving the infinite series h 12 = c 12 + ρ dr 3 c 13 c 32 + ρ 2 dr 3 c 13 dr 4 c 34 c 42 + • • • [START_REF] Evans | Density functionals in the theory of nonuniform fluids[END_REF] And finally, this series can be summed by a simple substitution to give the OZ equation.

The Ornstein-Zernike equation can be viewed as nothing more than a definition of c 12 in terms of h 12 . But if it turns out to be simpler to make approximate descriptions of the direct correlation function, the OZ equation can be used to convert c 12 to the full correlation function, h 12 . In the context of simple liquids it was found that the structure of c 12 was much simpler than that of h 12 . The full correlation function has characteristic long-range wiggles for a fluid close to its critical point, but the direct correlation function has a simple, short-range form at all T, P [START_REF] Hansen | Theory of Simple Liquids[END_REF].

Nothing in the physical argument above is special to simple fluids. In fact the OZ equation can just as well be applied to the electrons in a molecule. In electronic structure theory the familiar objects are two-particle densities (rather than correlation functions) and the system of interest is non-uniform. This suggests a recasting of the OZ equation in the form

ρ h 12 = ρ c 12 + dr 3 1 ρ(r 3 ) ρ c 13 ρ h 32 (2)
which is achieved by introducing spatial dependence in the one-particle densities and substituting the definitions ρ h 12 = ρ(r 1 )ρ(r 2 )h 12 and ρ c 12 = ρ(r 1 )ρ(r 2 )c 12 into Equation (1).

In Equation (2) the two-particle densities have the following interpretation: ρ h 12 is the correlation component of the exact reduced two-particle density; and ρ c 12 is the correlation component from a calculation that only includes the effects of direct correlation. If a simple model for ρ c 12 can be found, Equation ( 2) can be used to convert it into the total correlation component of the two-particle density.

The physically motivated 'derivation' above is quite compelling, but in any case there are various connections to more formal theory. The OZ equation emerges quite naturally in density functional theory (DFT) of fluids [START_REF] Evans | Density functionals in the theory of nonuniform fluids[END_REF] or indeed of electrons.

In fluids the direct correlation function is defined as a second functional derivative of the excess energy functional (that is the total energy minus the ideal or noninteracting contribution) [START_REF] Evans | Density functionals in the theory of nonuniform fluids[END_REF]:

c fluid 12 = - δ 2 β(F[ρ] -F id [ρ]) δρ(r 1 )δρ(r 2 ) . (3)
The quantum analogue of the OZ equation has been investigated by Chihara in the context of ionic liquids [START_REF] Chihara | [END_REF] and liquid metals [9]. In zero-temperature electronic density functional theory an Ornstein-Zernike equation can still be derived, but in this case the h 12 is related to a density response function, not a 2-particle density, and the c 12 is proportional to r -1 12 + f xc (r 1 , r 2 ). This is a potentially interesting connection, since this c 12 appears in the Dyson equation that connects the noninteracting response χ 0 to the full response χ. It is the quantity approximated by r -1

12 in the adiabatic-connection fluctuation-dissipation formulations of the random phase approximation [10]. However this connection will not be further pursued here. Instead we focus on the possible role of the OZ equation as an efficient means to build in inter-pair correlations to simple pair-based descriptions of molecular electronic structure theory.

2. The OZ equation in electronic structure theory

MP2 and direct correlation

The direct correlation function describes correlations in the system between just two particles, with the others remaining in their mean-field distribution. It should be possible to approximate the direct correlation function using what in electronic structure theory would be called a pair theory.

The simplest pair theory is MP2, and this choice also has the advantage that, to first order, the correlation contribution to the pair density is simply provided by the amplitudes:

ρ MP2 12 = T ia,jb ψ * i (r 1 )ψ a (r 1 )ψ * j (r 2 )ψ b (r 2 ) . (4)
Moreover, the correlation energy is in turn simple to extract from this this pair density:

E corr MP2 = dr 1 dr 2 ρ MP2 12 r 12 ≡ tr TK . (5)
(Here we assume that the amplitudes are properly antisymmetrized and the exchange integrals are not: K ia,jb = (ia|jb)).

The OZ equation in the basis of ia-pairs can be found by inserting the definition in Equation (4) into Equation (1), multiplying by ψ * i (r 1 )ψ a (r 1 )ψ * j (r 2 )ψ b (r 2 ), integrating over the coordinates of both electrons, and cancelling the overlap matrices that appear on the left and right of each term in the equation. This gives where

T h ia,jb = T c ia,jb + T c ia,kc X kc,ℓd T h ℓd,jb , (6)
X kc,ℓd = dr ρ(r) ψ k (r)ψ c (r)ψ ℓ (r)ψ d (r) . (7)
It should be noted that the derivation relies on the inversion of the overlap matrix

S ia,jb = dr ψ i (r)ψ a (r)ψ j (r)ψ b (r) (8)
which becomes ill-conditioned in the basis-set limit; this issue has been ignored in the present work, and in the triple-zeta basis sets used here, the overlap matrix is invertible in double precision. A further assumpinot is that it is acceptable to use the uncorrelated (Hartree-Fock) density in the definition of X.

It is perfectly possible to proceed directly with this equation, but the implementation can be made much more efficient using density fitting. Orbital pairs are defined in the usual way

|ia) = D ia µ |µ), (9)
with |µ)≡ Ξ µ a set of atom-centered fitting functions (we use the MP2FIT sets of Weigend et al. [11]). If density fitted MP2 is used, the exchange integrals are expressed as

K ia,jb = D ia µ Kµν D jb ν , (10)
where Kµν = (µ|ν), and the trace in Equation ( 5) can be written

E corr MP2 = T ia,jb D ia µ Kµν D jb ν . (11)
From this equation it is obvious that the amplitudes can be projected into the auxiliary basis set

Tµν = T ia,jb D ia µ D jb ν . (12)
Then the correlation energy can equally well be computed in the auxiliary basis as

E corr MP2 = tr T K . ( 13
)
The OZ equation can also be written in the fitting set,

Th = Th + Tc X Th (14)
where

Xµν = dr ρ(r) Ξ µ (r)Ξ ν (r) . ( 15
)
The X integrals (in either representation) are computed by numerical quadrature on DFT integration grids. We have checked that neither grid tolerances nor the threshold used to avoid divergences in 1/ρ have any effect on the conclusions presented.

Since both Equations ( 6) or ( 14) are formally equivalent we will now drop the tildes and consider the fitted and conventional versions together. The OZ equation can be simply rearranged for T h :

T h = T c (1 -T c X) -1 (16)
and a new correlation energy can be extracted as E h = tr T h K. For the small first-row molecules we used MP2 to obtain T c , computed T h by Equation ( 16), and then evaluated the correlation energy E h : the results are terrible. This, on proper reflection, is not altogether surprising. MP2 theory was selected as a candidate for producing the direct correlation density ρ c 12 on the grounds of being the simplest pair theory. But MP2 really is an exceptional theory, and its surprising accuracy arises from partial cancellation of two rather drastic approximations.

First, it is an independent-pair theory. In such theories the two electrons in a given pair can avoid each other without 'worrying' about the fact that there are other electrons around. This leads to an over-estimation of correlation effects. The exact independent electron pair theory, IEPA [12,13], overestimates the correlation energy by around 15% (see below), implying that the lack of treatment of inter-pair correlations raises the MP2 correlation energy by a comparable amount.

But in MP2 theory the virtual Hartree-Fock eigenvalues that appear in the energy denominator contain a spurious Coulomb and exchange interaction with all N (instead of N -1) electrons [14]. This results in a spurious elevation of the virtuals, an increase in the energy denominators, and, hence, an under-estimation of correlation effects [15]. The spurious interactions present in the Møller-Plesset denominator are explicitly subtracted in Epstein-Nesbet (EN) perturbation theory [13,[16][17][18]. EN2 performs very well for systems of two (correlated) electrons such as H 2 or LiH [15]. But for systems with more than two electrons (i.e. for systems where inter-pair correlations become relevant) EN2 overestimates the full configuration-interaction correlation energy by 15-20% [15,18].

Thus MP2 typically gives good correlation energies as a direct result of the cancellation of errors caused by neglecting inter-pair correlations and by using one-particle energies with spurious interactions.

IEPA and direct correlation

It should be clear from the previous section that a reasonable pair theory to use as a model for direct correlation is IEPA, that is, the pair theory in which the electron correlation within each pair is treated exactly.

To review the argument: direct correlation, as it appears in the physical motivation leading to Equation (1), is the correlation between pairs of electrons, and the OZ equation offers a way to introduce all of the inter-pair couplings. Therefore it seems reasonable to test the ansatz ρ c 12 = ρ IEPA

12

. There are two barriers to doing this directly, one technical and one theoretical. Technically it is non-trivial to obtain the 2-particle density in IEPA, although it certainly possible and we do intend to investigate this in the future. The theoretical problem is that even if one could obtain the exact 2-particle density, the correlation energy only depends on the part of it linear in the amplitudes, and it is not straightforward to extract this part from the whole 2-particle density. Here we propose to test the ansatz of Equation ( 6 For infinite systems it is reasonable for the entire density ρ to appear in the integral of Equation ( 1). For a finite system it becomes important to resolve the fact that the density in the integral should be the total density minus the density distribution of the two electrons under consideration. In the present case, with the ansatz T c = T IEPA , we are in a position to identify the densities that should be subtracted, and this suggests a modification of Equation ( 6)

T h ia,jb = T c ia,jb + T c ia,kc X (ij) kc,ℓd T h ℓd,jb , (17)
in which the X integrals contain a modified density:

X (ij) kc,ℓd = ρ -ψ 2 i -ψ 2 j ρ 2 ψ k ψ c ψ ℓ ψ d . ( 18
)
This theory has the pleasing feature of being exact for a single pair, or for any number of non-interacting pairs: in this case IEPA is exact, and the modified OZ equation results in T h = T c . On the negative side it has a high computational cost (O(o 4 v 3 ) for the multiplication in Equation ( 17)) and although this can be reduced by density fitting we will consider a much simpler approach. Given that in general one expects correlation effects to be dominated by 'diagonal' excitations (double excitations from the same doubly-occupied orbital), one can approximate Equation ( 17) by keeping the original form of Equation ( 6) but with the definition

X kc,ℓd = ρ -ψ 2 k -ψ 2 l ρ 2 ψ k ψ c ψ ℓ ψ d . (19)
This theory, which we will call OZ ′ , is also exact for isolated single pairs, but has much lower computational cost. It should be noted that both of these forms of the OZ equation break orbital variance; however, we will see later that this has little practical impact.

Numerical tests

To examine graphically the various quantities involved it is convenient to compute intracule quantities. In particular we compute

I(u) = T ia,jb (ia|δ(u -r 12 )|jb) (20)
for the various amplitudes involved, so that the integral

E = ∞ 0 du 4π u 2 I(u) u (21)
is the corresponding correlation energy. These quantities are shown in Figure 1 for water in a cc-pVTZ basis; all calculations were performed with frozen (1s 2 ) core electrons on non-hydrogen atoms.

The basic structure of the plots is as follows. The quantities vanish at the origin (because of the Jacobian) and integrate to zero (because i|a = 0). The intracules have a negative region at around 1 bohr, reflecting the avoidance of electron pairs as a result of correlation effects, and a balancing positive region at larger r. More particularly the data reveal several of the features discussed above. In panel (a) of Figure 1 it can be seen that the MP2 results lie close to CCD, although not, as discussed above, for entirely the right reasons. It is also clear that IEPA overestimates the correlation effect, producing a correlation intracule that is around 15-20% too large. In panel (b) the data CCD OZ -1 were obtained by setting T h = T CCD and solving the OZ equation (Equation ( 6)) for T c . This has the effect of removing inter-pair correlations. Comparing with CCD OZ -1 one can immediately see that MP2 is completely unsuitable as a model for direct correlation, but IEPA is rather close.

Finally, in panel (c), CCD is compared with results arising from processing MP2 and IEPA amplitudes with the OZ equation. The MP2 OZ result underestimates the correlation effects (here we have used the OZ equation to fold in inter-pair correlations, so the resulting underestimation of correlation is a pure symptom of the use of a one-particle energy spectrum with spuriously high virtual states). On the other hand IEPA OZ agrees quite closely with CCD, suggesting that if one starts of with an accurate description of the independent pair correlations, the simple OZ equation of Equation ( 6) can effectively incorporate inter-pair correlation effects.

The general features are reproduced for many molecules, and further examples are (more briefly) shown in Figure 2 for aldehyde, ketene, F 2 and H 2 . For hydrogen IEPA is already exactly equivalent to CCD, and processing the amplitudes with the (unmodified) Ornstein-Zernike equation introduces spurious many-electron effects to this two-electron system. This problem does not arise in the two modified OZ equations (Equation (17) or Equations ( 6) and ( 19)).

Correlation energies

To assess the accuracy of correlation energies obtained using the OZ equation we have performed calculations on selected small molecules. The valence correlation energies are shown in Table 1. The broad features are that IEPA overestimates CCD by an average of about 40 millihartree; the uncorrected OZ equation (using X as defined in Equation ( 7)) compensates for this, but often somewhat overshoots; and the modified OZ equation (using X as defined in Equation ( 19)) corrects for this over-shooting, deviating from CCD by an average of only 2 millihartree. For comparison MP2 differs from CCD by an average of around 8 millihartree.

The systematic over-estimation of correlation energies by IEPA is corrected by the OZ equation, but as always the spread of the data is also important. Here the performance is less good, with the standard deviation in the OZ ′ method relative to CCD being on the order of 20 millihartree. In considering this result, it should be noted that the OZ equation certainly has no direct relation to CCD, and the standard deviation of CCD relative to the more accurate method of CCSD(T) for these molecules is itself 11 millihartree. Incidentally, it is instructive to look at the correlation energies per correlated pair: the average error for MP2 valence correlation energies per pair relative to CCD over all molecules except H 2 is around 1 millihartree. For H 2 the MP2 correlation energy is 7.5 millihartree too small, and here the true problem with MP2 is revealed: since H 2 has only one pair, there is no compensation between the use of the 'wrong' denominators and the neglect of inter-pair correlations.

IEPA is not invariant to unitary rotations amongst the occupied orbitals. It is therefore obvious to ask how the results vary with different representations of the occupied space, and in particular whether the results are improved using localized orbitals. All calculations have therefore been repeated with Boys [START_REF] Boys | Localized Orbitals and Localized Adjustment Functions[END_REF] and Pipek-Mezey [START_REF] Pipek | [END_REF] localized orbitals, and the data are presented in has been previously observed, we note that with local orbitals IEPA correlation energies are even more overestimated than with canonical orbitals [21]. With Boys localized orbitals, the mean signed error in IEPA relative to CCSD(T) is -74.2 millihartree; the OZ equation reduces this to just -9.6 millihartree, with a standard deviation of around 9 millihartree. With Pipek-Mezey orbitals, the agreement is even better, with a mean signed error of -4.8 millihartree and a standard deviation around 6 millihartree. IEPA followed by the OZ equation agrees more closely with CCSD(T) than does CCSD, which has a mean error and standard deviation of +20 millihartree and 9 millihartree for the same set of molecules.

The modified OZ ′ method does not perform so well, with correlation energies typically falling about halfway between IEPA and unmodified OZ results. The spread of errors is also increased with OZ ′ , indicating that the approximate treatment suggested Equation [START_REF] Boys | Localized Orbitals and Localized Adjustment Functions[END_REF] does not perform particularly well. Implementation and testing of the full theory corrected for spurious inter-pair correlations (Equa- 1 and2 relative to CCSD(T)/cc-pVTZ. The horizontal axis is marked in millihartree, and the data presented are for CCSD (solid), MP2 (dashed), IEPA (dot-dash) and OZ (dot-dot-dash). The IEPA and subsequent OZ calculations used Pipek-Mezey local orbitals.

200

100 100 tion (17)) will be investigated in the future.

In any case, the OZ equation does appear to do a reasonable job of correcting IEPA for inter-pair correlations. For further striking evidence one can view the statistics graphically (Figure 3). IEPA using Pipek-Mezey orbitals followed by the OZ equation is closer to CCSD(T) than is CCSD.

Conclusions

The Ornstein-Zernike equation connects the direct correlation function to the full correlation function of a fluid. A simple physical argument leads to the equation if 'direct correlation' is taken to mean correlation within a single pair of particles, neglecting all others. The electronic structure theory that corresponds to that de- Table 1. Valence correlation energies (in millihartree) in the cc-pVTZ basis for small first-row molecules using a variety of correlation treatments. OZ and OZ ′ refer to the OZ equation applied to IEPA amplitudes using X defined in Equations ( 7) and [START_REF] Boys | Localized Orbitals and Localized Adjustment Functions[END_REF], respectively. The mean error and standard deviations are relative to CCD. One can interpret the OZ energy expression diagrammatically. Using dots to denote the X integrals, a solid line for the IEPA amplitudes and dashes for exchange integrals, the OZ (or OZ ′ ) energy expression can be written as

Molecule

E h = + + + • • • . ( 22
)
From this expansion it is clear that the energy expression is topologically a sum of ring diagrams. In other contexts this kind of structure has the effect of incorporating some part of the inter-pair correlation to all orders. The method described in this paper has numerous theoretical deficiencies, some of which are listed below:

(1) The OZ equation for which there is a clear physical argument contains 2-particle densities, not amplitudes. The use of amplitudes is justified only for MP2 (because the amplitudes provide directly the first-order contribution the 2-particle density). Even if two-particle densities were used, there is no straightforward procedure for extracting the correlation energy.

(2) If amplitudes are used in the OZ equation, then the density ρ that appears in the integral term should probably be more properly considered as a density of excitations, and the whole theory be considered as a theory of a fluid of (bosonic) excitons i → a.

(3) We have made no consideration of spin. A proper treatment of spin would presumably further increase accuracy, and will be investigated in the future.

Despite these shortcomings we believe that the evidence provided shows three conclusions very clearly. First the OZ equation can be applied to the electrons in a molecule, because there is nothing particular to classical systems or to systems at finite temperature in the physical argument presented in the Introduction. Second, IEPA is a reasonable model for direct correlation, as defined by the OZ equation. This conclusion is equivalent to saying that the OZ equation appears to be a reasonable approach for incorporating inter-pair correlations. Lastly, the accuracy of MP2 arises from a compensation between the underestimation of pair-correlation energies and the overestimation arising from the neglect of inter-pair correlations. that trip and numerous insights on this project are very gratefully acknowledged. We are very grateful to Prof Hans-Joachim Werner for converting the CCSD code in Molpro into an IEPA code in a period of about 20 minutes. On a personal note, it has been my great pleasure and honour to collaborate with Prof Werner on local, density fitted and explicitly correlated electronic structure theory over the past years, and I thank him for his support, insight and friendship -FRM.

2 . 3 .

 23 ) with the definition T c = T IEPA . Page 5 of 23 URL: http://mc.manuscriptcentral.com/tandf/Correction for finite number of particles

Figure 1 .

 1 Figure 1. Comparison of correlation intracules 4πr 2 I(r) for water in cc-pVTZ with several correlation methods. Panel (a): comparison of CCD, MP2 and IEPA. The MP2 correlation intracule closely agrees with that from CCD, but correlation effects in IEPA are clearly overestimated. Panel (b): comparison of MP2 and IEPA with the direct correlation intracule obtained by processing the CCD amplitudes using the OZ equation in reverse. Panel (c): comparison of the CCD correlation intracule with intracules obtained by processing MP2 and IEPA through the OZ equation.

9 Figure 2 .Figure 3 .

 923 Figure 2. Correlation intracules 4πr 2 I(r) for four small molecules using a cc-pVTZ basis and the correlation methods CCD (solid), IEPA (dot-dash) and IEPA-OZ (dash). Panels show intracules for (a) acetaldehyde, (b) fluorine molecule (c) ketene and (d) hydrogen molecule. Note that for H 2 CCD and IEPA are exactly equivalent.

Page 2 of 23 URL: http://mc.manuscriptcentral.com/tandf/tmph

	October 19, 2009 12:23 Molecular Physics oz-iepa
	Molecular Physics
	Ornstein-Zernike equation	3
	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y

Page 3 of 23 URL: http://mc.manuscriptcentral.com/tandf/tmph

	October 19, 2009 12:23 Molecular Physics oz-iepa
		Molecular Physics
	4	D Case and F R Manby
	F o
	r
	P
	e e r
		R e v i e w
		O n
		l y

Table 2 .

 2 First, as

Page 7 of 23 URL: http://mc.manuscriptcentral.com/tandf/tmph

	October 19, 2009 12:23 Molecular Physics oz-iepa
		Molecular Physics
	8	D Case and F R Manby
	F o
	r
	P
	e e r
		R e v i e w
		O n
		l y

Page 10 of 23 URL: http://mc.manuscriptcentral.com/tandf/tmph

	October 19, 2009 12:23 Molecular Physics oz-iepa
	Molecular Physics
	Ornstein-Zernike equation	11
	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y

Table 2 .

 2 Valence correlation energies (in millihartree) in the cc-pVTZ basis for small first-row molecules using a variety of correlation treatments. IEPA and IEPA plus OZ equation were performed using Boys and using Pipek-Mezey localized orbitals, and errors are relative to CCSD(T).

			Boys		Pipek-Mezey
	Molecule	CCSD(T)	IEPA	OZ	IEPA	OZ
	C 2 H 2	-339.3	-395.3 -341.8	-383.4 -333.8
	C 2 H 4	-374.9	-422.3 -377.3	-416.7 -373.2
	C 2 H 4 O	-618.2	-710.5 -635.8	-705.6 -632.0
	C 2 H 5 OH	-653.2	-738.0 -660.1	-735.6 -658.2
	C 2 H 6	-414.5	-458.1 -413.4	-458.1 -413.4
	CH 3 CHO	-612.6	-702.7 -625.6	-693.9 -619.0
	CH 3 OH	-461.6	-522.5 -465.2	-520.0 -463.4
	CH 4	-224.7	-244.8 -221.0	-244.8 -221.0
	CO	-376.7	-443.4 -387.3	-434.5 -380.9
	CO 2	-622.5	-732.5 -642.4	-720.6 -633.3
	H 2	-39.3	-39.2	-36.8	-39.2	-36.8
	H 2 CCO	-574.9	-667.7 -591.1	-655.6 -581.8
	H 2 O	-275.4	-314.6 -274.5	-312.4 -272.9
	H 2 O 2	-523.2	-613.3 -539.1	-604.0 -532.1
	HCHO	-422.9	-490.6 -433.5	-481.5 -426.6
	HCONH 2	-640.7	-738.3 -655.3	-729.8 -648.7
	HCOOCH 3	-850.8	-980.2 -874.1	-970.5 -866.6
	HCOOH	-661.6	-768.1 -679.2	-758.5 -671.8
	HNCO	-605.3	-708.6 -623.6	-698.5 -615.8
	NH 2 CONH 2	-857.3	-986.0 -876.1	-978.9 -870.6
	NH 3	-255.3	-286.9 -253.2	-286.8 -253.2
	Mean error	-	-74.2	-9.6	-67.8	-4.8
	Std. deviation	-	35.2	8.8	32.3	6.4

Acknowledgements

FRM is grateful for discussions on the Ornstein-Zernike equation with distinguished experts on the statistical mechanics of fluids, namely Prof Jean-Pierre Hansen (Cambridge), Prof Bob Evans (Bristol) and Prof Ben Widom (Cornell). Some of the ideas presented here were developed while FRM was a visiting scholar at Cornell, hosted by Prof Garnet Chan, and Prof Chan's financial contribution to


Johannes Kästner
email: kaestner@theochem.uni-stuttgart.de

Paul Sherwood

The ribosome catalyses peptide bond formation by providing high ionic strength

Keywords: Ribosome, catalytic activity, QM/MM, ionic strength, umbrella sampling Ribosome, catalytic activity, QM/MM, ionic strength, umbrella sampling

The ribosome catalyses peptide bond formation by

providing high ionic strength

Introduction

The ribosomes provides the link in the translation between the main alphabets in biology, the DNA sequence and amino-acid sequence of proteins. It catalyzes the formation of peptide bonds between amino acids bound to tRNA.

Each ribosome consists of a large subunit and a small subunit (50S and 30S in procaryotes) assembled around mRNA. Two molecules of tRNA attach to the mRNA and reach into the large subunit. They carry the amino acids which are to be connected by the ribosome. The three-dimensional structure of the ribosome of different species is known from crystallography [1][START_REF] Hansen | Proc. Natl. Acad. Sci[END_REF][START_REF] Bashan | [END_REF][4][5][6] at resolutions between 2.4 Å and 3.5 Å.

The formation of the peptide bond happens in the so-called peptidyl transferase center (PTC). There, the tip of one tRNA molecules locate a new amino acid in the A-site (amino site) of the ribosome. This is spatially close to the P-site (peptidyl site) where another tRNA is connected to the evolving peptide string. The latter is attached to the one amino acid on the A-site by forming a peptide bond which elongates the peptide by one amino acid. Then, the now empty tRNA molecule on the P-site leaves the ribosome, the whole ribosome moves in a process called translocation, so that the remaining tRNA with the peptide chain occupies the P-site. Then, a new tRNA molecule arrives with another amino acid to be attached to the peptidyl chain. The correct sequence alignment is achieved through the interaction with mRNA.

An intriguing structural feature of the environment of the PTC is that it exclusively consists of chemically inert moieties: RNA, water, and counter-ions, most probably sodium or potassium. In this study we aim at explaining how such a chemically inert environment is able to catalyze the peptide bond formation in the ribosome. Figure 1 gives an overview of the atoms likely to be involved in the peptide bond formation. The new bond is formed between the N-atom of the one amino acid connected to the A-site and the C-atom of the peptide connected to the P-site. The bond from the latter atom to O3' of the P-site (which is part of the ribose of the last RNA residue, Ade 77 of tRNA) is broken. Additionally, a hydrogen atom leaves N of the P-site and the same or a different hydrogen is added to O3' of the P-site. The hydrogen atom may be directly transferred, in what will be referred to as the direct mechanism in this paper, or there may be a proton shuttle mechanism. The proton shuttle can proceed in a concerted manner with the change in the bonds between the heavy atoms or in a general acid/base mechanism.

Temperature-dependent kinetic measurements of peptide formation in a model system as well as the full ribosome [7][8][9][10] showed that the ribosome actually increases the enthalpic part of the reaction barrier compared to the reaction in water, but decreases the entropic part. Results are summarized in Table 1.

An influence of Ade 2451 (E. coli notation is used throughout this paper) on the catalytic activity as well as charge relay involving Gua 2447 were postulated [1]. This would be consistent with a general acid/base mechanism. However, replacement of Ade 2451 does not lead to a complete loss of function [7,[START_REF] Thompson | Proc. Natl. Acad. Sci[END_REF][START_REF] Polacek | [END_REF], and Gua 2447 can be replaced by Ade without changing the the pK a value of 7.5 [13]. A general acid-base mechanism involving proton abstraction from the amino group of the A-site to Ade 2451 was found unfavorable by theoretical investigations [14].

Density functional calculations of the termini of tRNA with two amino acids in vacuum [15] resulted in an activation barrier of 148.6 kJ mol -1 . The tRNA bases were found to rotate with respect to each other during the reaction. Two hydrogen bonds were found to stabilize the transition state: one between the inactive amino-H of the A-site and the peptidyl-O of the A-site (HB1), the other one between H2' of Ade 77 of the A-site and peptidyl-O of the P-site (HB2). Both are indicated in Figure 1. Including the ribosome environment in a force field description and treating the reactive part with the empirical valence bond method (EVB) [START_REF] Warshel | Computer Modeling of Chemical Reactions in Enzymes and Solutions[END_REF][START_REF] Åqvist | [END_REF] resulted in one barrier of 41.9 kJ mol -1 towards a

water [8] water [10] model system [8] ribosome [ tetrahedral intermediate and another barrier of 71.2 kJ mol -1 relative to the reactant for breaking the remaining C-O bond [14,18]. In the tetrahedral intermediate, the aminonitrogen of the A-site is already bound to the carboxyl-carbon of the P-site, but its bond to O3' of the P-site is not yet broken.

A proton shuttle mechanism via the -O2'-H group was found by computer simulations [18]. The same shuttle mechanism was found [19] not to catalyze the respective reaction in water.

Different theoretical work [20] addressed the decoding, the correct alignment of tRNA at mRNA in order to correctly transfer information from mRNA to the protein sequence. All-atom molecular dynamics (MD) simulations of the whole ribosome (2.6 Mio atoms) of Thermus thermophilus were used.

Here, we investigate the formation of the peptide bond and the acompanying proton transfer in the direct mechanism as well as in a shuttle mechanism via O2' of Ade 77 of the P-site. We used density functional theory (DFT) calculations, coupled to the ribosome environment described with a classical force field in a QM/MM approach. To cover the entropic part, excessive sampling of the system was necessary which is prohibitive with DFT. Thus, we used the semiempirical self-consistent-charge density-functionaltight-binding (SCC-DFTB) method [21] with parameters given in the same citation to describe the bond breaking and formation. SCC-DFTB was shown to provide good geometries in biochemical systems. The energies obtained in this way are less reliable than those from DFT, however, SCC-DFTB is about two orders of magnitude more efficient in terms of CPU time, enabling finite-temperature MD sampling.

Methods

Level of theory

Our model of the ribosome was described by a QM/MM approach [22][23][24][START_REF] Senn | Atomistic Approaches in Modern Biology -From Quantum Chemistry to Molecular Simulations[END_REF][START_REF] Senn | [END_REF][27][28]. Using electrostatic embedding, the MM charges of the force field polarize the QM part. Covalent bonds between the QM part and the MM part were truncated on the QM side by hydrogen link atoms. The charge-shift scheme [24] was employed in order to avoid over-polarization of the QM density near the links. The CHARMM22 force field [29][30][31][32] was used for the MM part. Two different approaches were used for handling the bond-breaking and formation in the QM-part: DFT and SCC-DFTB. The BP86 functional [33][34][35][36][37] with a DZP basis [38] as implemented in GAMESS-UK was used for the DFT calculations. The BP86 functional has been shown to agree well with higher-order ab initio calculations on the peptide bond formation in water [39]. The SCC-DFTB, on the other hand was used to study an enzymatic process of cleaving a peptide bond with good qualitative agreement to experiment, especially for ∆∆G values [40]. Thus, both levels seem appropriate for our purposes. Convergence of the results with respect to the basis set was tested by comparing energy minima obtained with larger basis sets, see the Supporting Information. Coulomb fitting [41] was used with the Ahlrichs [42,43] basis set. The DFT calculations were performed with GAMESS-UK [44,[START_REF]GAMESS-UK is a package of ab initio programs[END_REF], the QM/MM coupling with ChemShell [24,[START_REF] Chemshell | a Computational Chemistry Shell[END_REF], and the force field calculations with DL POLY as included in ChemShell.

If not explicitly stated otherwise, the QM part contained the atoms of the ribose rings of both adenolyl-termini of the tRNA and the whole amino acids connected to those. The truncation of the QM subsystem was done by cutting through the bonds from C4'-C5', and C1'-N9, the nitrogen atom connecting the purine ring system to the ribose. This resulted in a total of 58 QM atoms plus 4 hydrogen link atoms. This choice of the truncation is well justified because only single bonds have been cut, charge transfer between the regions designed as QM and MM regions is not to be expected, and the distance between the cuts and the chemically active atoms is 3 bonds at the minimum. Effects of the truncation, however, have been tested by enlarging (and shrinking) the QM-region and calculating the difference between two energy minima, see the Supporting Information.

System preparation

The ribosome of Thermus thermophilus was modeled based on the computer model of the whole ribosome [START_REF] Tung | [END_REF] (PDB entries 1TWT and 1TWV). This model includes both subunits, two tRNA molecules, and a short strand of mRNA. The amino acid phenylalanine was modeled on the A site based on the transition-state analogue RAP characterized by crystallography [48] (pdb entry 1VQP). Glycine with its N-terminus acetylated to mimic the continuing protein chain was modeled on the P-site, also based on 1VQP.

All residues with at least one atom within 30 Å of the central atom (N) were included in our model, all other residues were deleted. All atoms more than 30 Å away from the central atom were frozen. This resulted in 3 fragments of protein chains, located close to the surface of the sphere, 228 bases of ribosomal RNA of the longest chain of the large subunit, 5 bases of tRNA of the A-site, and 6 bases of the P-site, being included. The system was superimposed with a sphere of pre-equilibrated TIP3P [49] water molecules of a radius of 30 Å. All water molecules overlapping with any atoms of the model were deleted. The system is highly negatively charged due to the charge of -1 of each RNA unit. The fragments of the protein chains contain 13 positively charged amino acids and 1 negatively charged one. Charge neutrality was achieved by replacing 225 water molecules, each within 7 Å of a phosphorus ion, by sodium ions. The final setup contained a total of 13,310 atoms.

The system was truncated by a spherical boundary potential, i.e. a flat potential within a cutoff distance r c from a fixed point in space (the reaction center), and a potential of the form E(r) = k 2 r 2 for r > r c . The values r c = 30 Å and k = 10 -3 Hartree Bohr -2 (about 9.4 kJ mol -1 Å-2 ) were used.

The system was equilibrated extensively on an MM-only level in order to make sure all cavities within the RNA were filled with water and that the ions were distributed consistently with the charges on the RNA phosphate groups. These equilibration runs, as well as the umbrella sampling simulations (see below) were performed at a temperature of T = 300 K in a canonic ensemble achieved by a Nosé-Hoover chain thermostat [50][51][52][53] with a chain length of 4 and a characteristic period of 20 fs, corresponding to a thermostat wavenumber of 375 cm -1 . Newton's equations of motion were integrated with a reversible noniterative leapfrog-type integrator [54] with a time step of 1 fs. All atoms except the frozen atoms were allowed to move. Internal degrees of freedom of all water atoms were constrained, consistent with the TIP3P model. Towards the end of the equilibration runs, restraints to keep the N atom of the A-site in proximity to the C atom of the P-site were introduced. In order to provide an ensemble of starting geometries close to the assumed transition state, restraints of the form E(ξ) = k 2 (ξ -ξ ref ) 2 along the reaction coordinates:

ξ C = d(C-O) -d(C-N) F o r P e e r
ξ H = d(N-H) -d(O-H)
with reference values of ξ ref C = ξ ref H = 0 and force constants of k = 0.05 Hartree Bohr -2 (about 469 kJ mol -1 Å-2 ) were used. Here and in the following, N and H stand for the atoms of the amino group of the A-site, and C and O stand for the peptidyl-carbon and O3' of the P-site. Next, an equilibration run of 145 ps was performed on the QM/MM level with DFTB as the QM method. The same restraints as in the final MM equilibration described above were used (ξ ref C = ξ ref H = 0 and k = 0.05 Hartree Bohr -2 ). Statistical tests [55] applied as described previously [56], showed that the system was equilibrated with respect to ξ C and ξ H already after 25 ps. Ten snapshots were taken in intervals of 10 ps from 55 ps to 145 ps. These served as starting configurations for the following geometry optimisations and transition state searches.

Geometry Optimisation

In each of these snapshots, we aimed at finding a transition state and subsequently minimised the energy to find the connected reactant structure. However, a transition state search in such a large system, 28,812 degrees of freedom, is difficult to perform. So, first, we kept the restraints in place and minimised the total system on a DFTB/MM level of theory using an L-BFGS [57,58] (limited-memory version of a Broyden-Fletcher-Goldfarb-Shanno) minimiser in hybrid-delocalised internal coordinates (HDLC) [59]. This makes sure all vibrational elongations distributed over the whole system are damped out, which saves time and effort for the following transition-state search. We removed the restraints and used the super-linearly converging variant [60] of the dimer method [61][62][63] in HDLC coordinates to locate the transition states. This lead to transition states corresponding to the direct mechanism. We used weights of 1 for all atoms in the QM part and weights of 0 for all other atoms in an algorithm described previously [60] which effectively restricts the transition mode to the QM atoms. The energy is minimised with respect to the coordinates of all atoms with weight 0.

Starting from these transition states for the direct mechanism optimised at the DFTB/Charmm level, we optimised transition states for both the direct mechanism and the proton shuttling mechanism at the DFT(BP86)/Charmm level. The dimer method was used with the same settings as described above. We also distorted the system slightly into the direction of the transition mode and minimised the energy to find the structures of the tetrahedral intermediate (where stable) and further to the reactant at the DFTB/Charmm level. The geometries of the transition states at the DFT(BP86)/Charmm level were used equivalently to find the reactant state structures on that level.

To locate the transition states of the proton shuttle mechanism, we started out from the TS of the direct mechanism and used restraints to move the system in the vicinity of the TS of the proton shuttle mechanism. Then we removed all restraints and used the dimer method in HDLC coordinates and with weights as described above to find the TS. This was successful in 7 of the 10 snapshots. In the remaining snapshots we used the improvedtangent nudged-elastic band method [64][65][START_REF] Jónsson | Nudged Elastic Band Method for Finding Minimum Energy Paths of Transitions[END_REF] with a climbing image with the same weights as for the dimer method in the version implemented [START_REF] Goumans | [END_REF] in DL-FIND. This resulted in better starting configurations and led to convergence of the subsequent dimer searches in one case (snapshot 5). In the two remaining snapshots (3 and 9) we were unable to locate transition states for the proton-shuttle mechanism.

Reaction barriers were consistently calculated with respect to the energy of the reactant state obtained for the same snapshot. The optimiser DL-FIND [68] was used for all geometry optimisations and transition state searches.

Influence of individual residues on the reaction barrier

After having located the transition states and the reactant states connected to them, we estimated the electrostatic influence of individual residues (RNA units, water molecules, and sodium ions) on the activation barrier using the static structures. The full QM density was replaced by restricted electrostatically fitted charges (RESP), [69] fitted to reproduce the electrostatic potential of the full DFT (BP86) density, as these can be expected to result in a quite accurate electrostatic energy at hugely reduced cost compared to the full QM density. The change in the activation barrier ∆∆ ‡ E i due to the charge on residue i can be determined as

∆∆ ‡ E i = ∆ ‡ E 0 -∆ ‡ E i (1)
where ∆ ‡ E 0 is the electrostatic component of the activation energy calculated using ESP charges instead of the QM density, and ∆ ‡ E i is the electrostatic component of the activation energy with all charges on residue i set to zero. In these calculations, the geometries of both the RS and the TS are kept unchanged. Thus, the self-energy of the QM-part is constant and drops out of ∆∆ ‡ E i . Note that ∆∆ ‡ E i contains contributions from MM-MM interactions as well as from QM-MM interactions. Thus, if ∆∆ ‡ E i is positive, the atom charges of residue i increase the barrier (destabilise the TS). Otherwise these charges stabilise the TS.

∆∆ ‡ E i estimated in this way is a semiquantitative measure of the electrostatic influence of individual residues, helpful to determine which residues play a role in the catalytic activity. However, ∆∆ ‡ E i is certainly too crude an approximation to be compared with the effect of a mutation of the respective residue on the reactivity. Among the effects not covered by ∆∆ ‡ E i are the substitution of residue i by other moieties, changes in the geometries of RS and TS by the change, as well as changes in the polarisation of the QM part.

Umbrella sampling

Biased MD simulations [70,71] on the DFTB/Charmm level were used to investigate the free-energy changes and, thus, the role of entropy along the reaction. A bias of the form

E(ξ) = k 2 (ξ -ξ ref ) 2
was used in both reaction coordinates ξ C and ξ H with a force constant of k = 0.05 as described above. A total of 238 windows were sampled in steps of 0.5 Bohr in both of the reaction coordinates. Each window was sampled until 10 ps of the trajectory were trend-free in both reaction coordinates, often leading to significantly longer equilibration times in each window. The results were analysed with two-dimensional umbrella integration [72]. Despite rather long total production time of 2.4 ns, however, the statistical significance of the results was not satisfactory. Especially the reduced quality due to a lack of stability of the sampling around the transition state deteriorated the accuracy. From an approximate potential energy surface, whose contours are shown in Figure 2, however, it is clear that the minimum free-energy path proceeds roughly along the diagonal of the graph. Thus, the sum ξ CH = ξ C +ξ H serves as a good one-dimensional reaction coordinate for umbrella sampling simulations.

One-dimensional umbrella sampling simulations along ξ CH were performed in 74 windows from ξ ref CH = -11.3 Bohr to 10.6 in steps of 0.3. Each window was sampled until 10 ps of trend-free date were achieved. Umbrella integration analysis [73,74]

resulted in ∆ ‡ G = 228.1 ± 3.4 kJ mol -1 .
Visualisation was found to be essential for detecting errors, understand the results, and inspire interpretation. VMD [75] versions 1.8.6 and 1.8.7 were used for molecular visualisation tasks as well as to create Figures 3 and5.

ξ H ( Å) 4 3 2 1 0 -1 -2 -3 4 3 2 1 0 -1 -2 -3
Figure 2. Contour lines (in distances of 50 kJ mol -1 of the free energy with respect to the reaction coordinates ξ C and ξ H (solid lines) and the boundaries of the sampled area (thin dashed line). The free-energy surface clearly shows that ξ CH = ξ C + ξ H can be expected to be a good reaction coordinate. Table 2. Activation barriers in kJ mol -1 for the direct mechanism and the proton-shuttle mechanism. Possible reasons for outliers are discussed in section 3.5.

Snapshot ∆ ‡ E direct ∆ ‡ E shuttle ∆ ‡ E shuttle -∆ ‡ E direct 1

Results and Discussion

Reaction Mechanism

The available data, although obtained from a high level of theory and at least some sampling, preclude an answer as to weather the direct mechanism or the proton-shuttling mechanism occur in the real system. It may well be that both mechanisms play a role. The activation barriers (potential energies) are given in Table 2. For three snapshots, the direct mechanism shows a lower barrier than the proton-shuttle mechanism. For the remaining five snapshots for which comparison is possible, the proton-shuttle mechanism is (partially very slightly) preferable. The average barriers are with 104.2 kJ mol -1 and 106.6 kJ mol -1 practically identical.

An exponential average between the snapshots (which is close to the lowest energybarrier) may be more appropriate than the direct average to estimate a turnover rate and to compare to the experimental value for the barrier. This is justified because the relevant quantity determining the reaction rate is the free energy rather than the potential energy. Following Jarzynski's equation [76] the free energy ∆A of a process equals the exponential average, ∆A = -k B T ln exp(-∆W /k B T ) , of the work ∆W (which can be approximated by potential energy differences ∆E) drawn from a canonic ensemble. Here, k B refers to Boltzmann's constant and T to the absolute temperature. Of course, a sample of ten snapshots is nowhere near a complete set and, therefore, can not be expected to provide an accurate average. Furthermore, our barriers have not directly been obtained from a canonic ensemble. Instead, starting configurations have been drawn from a (true) canonic ensemble, but then the reactants and transition states were optimised. Still, it may be more relevant to take the exponential average (which is close to the lowest barrier) into consideration rather than the direct average. A different reasoning for the lowest barrier found being most important is transitionstate theory. If many paths lead from the reactant to the product, the one with the lowest barrier has the highest reaction rate. If the barriers differ significantly (by more than k B T ), as they do in our results, then the reaction will almost exclusively proceed via the fastest route. Starting from different MD-snapshots is only a means of finding different transition states.

The exponential average for the activation energy at T = 300 K is with 78.7 kJ mol -1 noticeably smaller for the direct mechanism than for the proton-shuttle mechanism with 90.1 kJ mol -1 . The lowest barriers contribute most to the exponential average. So, a similar difference is obtained when just comparing the lowest values for the two mechanisms: 73.1 kJ mol -1 for the direct mechanism and 85.9 kJ mol -1 for the proton-shuttle mechanism. The exponential average for the 18 values of both mechanisms taken together is 80.2 kJ mol -1 . While these data show some distinction between the two mechanisms under investigation, in the real system still both mechanisms might occur. There is a substantial spread in the activation energies between the different snapshots. The difference between the highest and lowest values is as much as 51.5 kJ mol -1 for the direct mechanism and even 71.6 kJ mol -1 for the proton-shuttle mechanism. Even if these values may be dominated by outliers, they show the importance of calculating the barrier at more than one snapshot. Differences between the snapshots in the geometry of the reacting atoms are discussed here, while differences in the coordination, which may explain the main outlier (snapshot 10) are discussed in section 3.5. The spread also shows that the ribosome is a very flexible system with lots of geometrical changes occurring even during limited sampling time.

Tables 3 and4 show the main geometric data of the transition states in the two mechanisms under investigation. In the direct mechanism, the C-N distance is almost constant between the different snapshots, with a standard deviation of only 0.021 Å. This distance is, with 1.562 Å on average, also already close to the product, with typically 1.365 Å. The low spread and short bond in the TS can be interpreted as the C-N bond being mostly formed, maybe only as a single bond, already in the TS. In the PS, the C-N bond has partially double-bond character. By contrast, the C-O distance, and especially the O-H distance vary substantially between the snapshots. These bonds also dominate the transi- tion mode.

The geometric data of the atoms directly involved in the transition are uncorrelated with the activation barrier. Thus, at least for the direct mechanism, the differences in the activation barriers between the snapshots must be caused by interactions with the environment, rather than by early or late transition states.

The transition modes vary rather substantially between the snapshots. While in the first snapshot, the transition is dominated by hydrogen transfer, the transition mode significantly extends to the sugar rings and a nearby water in snapshot 4, the one with the lowest barrier. Correspondingly, the imaginary frequency is rather high (1208 cm -1 ) in snapshot 1, while it is quite low (201 cm -1 ) in snapshot 4.

The situation is similar for the geometric data of the proton shuttle mechanism, see Table 4. Here, the C-N and the N-H distances show the smallest standard deviation, with 0.044 Å and 0.032 Å, respectively. These distances are also similar in values to fully established bonds. The main transition (and also main variation between the snapshots) is the breaking of the C-O bond and, concertedly, the formation of the bond O-H2'. As with the direct mechanism, a correlation between geometric features of the reactive centre and the activation energy could not be found.

Note that neither of the two hydrogen bounds previously found in a gas-phase model [15] are present in any of the transition states we found by QM/MM. There was just no local minimum of that kind found in the geometry optimisation. In section 3.6, we will report on tetrahedral intermediates that have been found in 6 of the 10 snapshots when using SCC-DFTB as QM method in static calculations. However, no minimum corresponding to a tetrahedral intermediate was found on the free-energy surface using that QM method. Thus, if a tetrahedral intermediate exists in the real system, it is a very shallow minimum with it energy similar to that of the reactant state. Others [14,18] have found the barrier from the RS to the tetrahedral intermediate to be significantly smaller than the one from the latter to PS. Thus, a tetrahedral intermediate is irrelevant for the reaction mechanism.

F

Decomposition of the QM/MM energies

In order to assess the influence of the ribosome environment on the PTC, we simulated the equivalent model reaction in the gas phase: rib-phe(gas) + rib-ala(gas) → rib(gas) + rib-ala-phe(gas) where rib stands for the ribose moiety (without the RNA base which would be present in adenosyl), and ala and phe for the amino acids used in our model, alanine (with acetylated N-terminus) and phenylalanine, respectively. This model corresponds to the QM part of our QM/MM simulations. The activation barrier for the gas-phase reaction is ∆ ‡ E(gas) = 118.0 kJ mol -1 . To obtain this number, the optimisations were started from the QM/MM-RS and TS of snapshot 4. While it is unlikely that other local minima would be reached in gas phase optimisations of the other snapshots, this possibility cannot be excluded. Thus, for this snapshot, the ribosome decreases the activation barrier by 44.9 kJ mol -1 compared to the gas-phase reaction.

In the gas phase, the ribose rings of Ade 77 connected to the peptides take over some of the tasks of the whole surrounding. Hydrogen bonds between the peptidyl-O atoms and the hydrogen at the O2' position keep the reactants in place, a finding that confirms previous gas-phase studies [15].

To investigate the geometric influence of the environment, an energy decomposition scheme, similar to a previously used one [77], is applied. As a first step, we consider the distortion of the gas-phase reactant complex by the ribosome. Thus, we compare the energy of the reactant complex in the gas-phase geometry to the reactant complex in the QM/MM geometry. This-strictly positive-energy contribution E RS (gas//env) -E RS (gas//gas) = 113.6 kJ mol -1 . It is mainly caused by breaking of intramolecular hydrogen bonds in the gas-phase complex, which are to be replaced by interactions with the environment. In this comparison, only the geometry of the QM/MM model is used. All interactions with the ribosome environment, also the polarisation of the electron density, are omitted.

The same decomposition can be performed for the transition state, which provides

the activation barrier calculated in the gas-phase at the QM/MM-optimised structures, ∆ ‡ E(gas//env), corresponding to the process RS(gas//env) → TS(gas//env)

We obtain ∆ ‡ E(gas//env) =156.4 kJ mol -1 , substantially (83.3 kJ mol -1 ) higher than the ribosome-catalyzed barrier of ∆ ‡ E(env//env) =73.1 kJ mol -1 for the same snapshot. This comparison clearly shows, that the transition state is significantly stabilised by the ribosome compared to the reactant state. This stabilisation is due to interactions with the environment of the reactive centre within the ribosome rather than due to just placing the reactants in a favourable position next to each other. Of course, here the comparison was done with respect to the gas phase. A fairer comparison may be to the reaction in water. However, we doubt the validity of common continuum solvation models for this system, as much of the interactions are hydrogen bonds and interaction with ions, see below, which cannot be covered by continuum solvation models. Taking solvation in water or salt solutions into account explicitly is computationally rather demanding and, therefore, outside of the scope of this work.

From comparison to the gas-phase reaction, we learned that the ribosome has to distort the reactants in order to reduce the barrier. Now, we will discuss which specific interactions to the environment lead to the significant reduction of the barrier by 83.3 kJ mol -1 . The electrostatic influence of individual residues on the catalytic activity is assessed by calculating ∆∆ ‡ E i for each residue for the direct mechanism obtained with snapshot 4 (DFT(BP86)/Charmm), the one with the lowest activation energy. Visualisation by means of showing only residues with |∆∆ ‡ E i | larger than a modifiable threshold and colouring them according to ∆∆ ‡ E i for both the RS and the TS geometries helped to understand the interplays between the residues involved.

Electrostatic influence of the environment of the PTC on the reaction

The largest electrostatic influence of a single residue (∆∆ ‡ E i = 39.4 kJ mol -1 ) is caused by a water molecule which forms a hydrogen bond to the RS, but not to the TS. The H-bond is formed to the peptidyl-O of the P-site. In the TS, instead, the same water molecule H-binds to another water, which causes a decrease of the barrier by ∆∆ ‡ E i = -14.4 kJ mol -1 . The latter water molecule, in turn, is pulled away from a sodium by the former one when the TS is formed, which causes the sodium to increase the barrier by ∆∆ ‡ E i = 29.6 kJ mol -1 . One sodium atom stabilises the transition state by ∆∆ ‡ E i = -39.2 kJ mol -1 by moving closer to O2' of the P-site (2.223 Å at the TS vs. 2.287 Å at the RS). The movement, however, is only possible, because a hydrogen bond between O2' of the P-site and O of the A-site which holds O2' away from the same sodium in the RS, is broken in the TS. The same sodium atom is also stabilised by O4' of Ade 2451 which is supposed to have catalytic activity. However, ∆∆ ‡ E i of Ade 2451 is only +3.4 kJ mol -1 , so it has hardly any net electrostatic influence.

Another water molecule and a sodium ion stabilise the TS by ∆∆ ‡ E i = -29.5 kJ mol -1 and -17.9 kJ mol -1 , respectively. The sodium is tighter bound between the peptidyl-O of the P-site, O of the (truncated) amino acid chain of the P-site, and O2' of the A-site's Ade 76 as well as by O2' of Ura 2584. The water molecule and Ura 2584 more efficiently stabilise the sodium in the TS and thus reduce the barrier. The energetic contribution of Ura 2584 is ∆∆ ‡ E i = -24.6 kJ mol -1 , mainly due to better stabilisation of the sodium ion (distance 2.38 Å in the TS and 2.49 Å in the RS).

Concerning the RNA environment, Cyt 75 of the P-site destabilises the TS by ∆∆ ‡ E i = 13.5 kJ mol -1 . One rational may be the slightly weaker hydrogen bond between its NH 2 group and a phosphate-oxygen of Ade 76 of the A-site.

It is striking that there are not less than 5 sodium ions and 5 water molecules with an individual electrostatic influence on the barrier of more than 12 kJ mol -1 , compared to only two RNA residues.

The electrostatic influence of these individual residues seems quite pronounced compared to the overall barrier. However, one has to keep in mind, that any residue, would it be removed artificially (by mutagenesis or any other method) would be replaced by other moieties, which might have a similar influence on the reaction.

It should also be noted that only the electrostatic interaction at static structures of the RS and the TS enter the analysis based on ∆∆ ‡ E i . The rigidity of the static structures leads to rather large numbers for ∆∆ ‡ E i . In a dynamic situation at finite temperature, the system would accommodate changes in the environment by altering details in the reaction path which likely reduce the numerical values. So, rather than the absolute numbers, the most important feature of the analysis based on ∆∆ ‡ E i is, that water and, especially, sodium ions play a major role in reducing the barrier. The apparent importance of ions for the reaction warrants to review the ionic strength, the concentration of ions in our model. In a sphere with a radius of 30 Å we have 238 singly positively charged ions (and as many negatively charges ones), resulting in an ionic strength of I = 3.5 M (M standing for mol dm -3 ). The ionic strength found in the ribosome is, thus, about 25 times higher than the physiological ionic strength of about 134 mM [78,79]. Note that this can not stem from a deficiency in our model. We cut out a specified volume of the experimental structure of the ribosome. Thus, the concentration of RNA-phosphate, which causes the high ionic strength, is taken from experiment. We just added the counter-ions. These, however have to be present in the real system as well because huge accumulations of negative charge would be unstable. Since the negatively charge phosphate groups are distributed rather evenly within the volume of the ribosome (each RNA residue contains one phosphate group) one can also expect a more or less uniform distribution of the cations. We could have used doubly-charged calcium instead of sodium, which would even increase the ionic strength. We, thus, conclude that the ribosome, and likely other ribozymes as well, catalyse reactions by placing the reactants in an environment with very high salt concentration.

Variability of the stabilisation of the TS by the ribosome

We studied the environment of the PTC in the different snapshots, searching for interactions of the active atoms with their surrounding. Notably, these were rather different between the snapshots. This causes the large spread in activation energies. Certainly, for systems as flexible as the ribosome, it is essential to calculate the activation energy in independent geometries rather than only using one state.

We presented the analysis based on ∆∆ ‡ E i only for snapshot 4. However, for the other snapshots we observed similar changes in the interaction of the active region with sodium ions between the RS and the TS. These will be discussed on more qualitative grounds in the following.

In snapshot 4, the peptidyl-O of the P-site is H-bound to a water molecule in the RS, but not in the TS, see above. The same is true for snapshot 6. In snapshot 10, the same atom is H-bound to two water molecules in the TS, but to only one in the RS. The RS, but not the TS, is stabilised by a H-bond from the peptidyl-O of the A-site to the amine group of Ade 2602 in snapshot 5.

O3' of the P-site, the bond of which to the peptidyl-C is broken in the main reaction, is stabilised in the TS by a H-bond/interaction to a water in snapshot 5, to H-O2' in snapshots 4 and 5, as well as to a sodium ion in snapshot 3. These interactions are absent in the RSs of the respective snapshots and thus lower the barrier.

An interaction of the hydrogen atom to be transferred in the main reaction with O2' of the P-site is present in the TS but not the RS of snapshots 2 and 5, and in the RS but not the TS of snapshot 8. This H-bond is the key feature of the proton shuttle mechanism.

The peptidyl-O of the A-site is stabilised by a water in the TS, but not RS, of snapshot 4, and in the RS but not the TS of snapshot 8.

In all snapshots but 9 and 10, the peptidyl-O of the P-site is stabilised by a sodium ion in the reactant state as well as in the transition state. In the last two snapshots (9 and 10) the sodium ion had moved away and is absent in both RS and TS (both in the TSs of the direct mechanism and the proton-shuttle mechanism). This may be one of the reasons for the high activation barrier of snapshots 9 and 10. Also it was observed that the lonepair of the amino group of the A-site, which is to attack the carbon atom of the P-site, already points towards this carbon in many of the reactant states. Notably, this is not the case in snapshots 2, 5, and 8-10. This may also contribute to the high activation energy of This incomplete list already hints the relatively low activation energy of snapshot 4: the most selectively catalytically active interactions of the reactive centre with the environment are found in snapshot 4. It, moreover, again shows the importance of sodium ions in the catalytic activity.

DFT compared to SCC-DFTB

In order to allow for extensive sampling of the free-energy profile, we augmented our DFT calculations by SCC-DFTB results. The latter method is a couple of orders of magnitude faster than DFT and, thus, allows for molecular dynamics sampling with proper equilibration. Before results of the free-energy simulations can be discussed, however, we show that the geometries of the stationary points on the PES are very similar to those found by DFT. We performed the SCC-DFTB calculations for the direct mechanism only.

Bond distances at the transition states of the 10 snapshots as well as the activation energies obtained at the SCC-DFTB/Charmm level are given in Table 5. The activation energies are, as expected, with a direct average of 195.0 kJ mol -1 and an exponential average of 161.9 kJ mol -1 noticeably higher than those obtained with DFT(BP86) (104.2 kJ mol -1 and 78.7 kJ mol -1 , respectively). A tetrahedral intermediate was found in 6 of the 10 snapshots. The energies of these tetrahedral intermediates relative to the reactant states are given in Table 5.

The geometries are, however, quite similar to the ones obtained with DFT(BP86)/Charmm. The C-N bond at the transition state is predicted about 0.36 Å longer on average, the N-H bond 0.11 Å longer and the O-H bond 0.21 Å shorter than by DFT(BP86)/Charmm, see Table 5. This means, concerning the carbon-"transfer", the SCC-DFTB/Charmm transition state is slightly earlier, concerning the hydrogen transfer, it is slightly later than the DFT(BP86)/Charmm transition state.

The similarity in the location of the transition states between DFT and SCC-DFTB is also illustrated in Figure 6. Using the reaction coordinates ξ C and ξ H , all transition states are found in a quite narrow region of about ±0.5 Å in both directions. While the regions between the energy expressions largely overlap, it is noticeable that SCC-DFTB/Charmm is shifted towards slightly lower values in ξ C and slightly higher values in ξ H . Still, the structural similarity is surprisingly high. Since the entropy is strongly related to the area in configuration space available to the reaction, thus on geometric grounds, free-energy sampling based on SCC-DFTB/Charmm can be expected to give valuable information about the entropy of the reaction.

Free energy of activation

We calculated the free energy of activation by umbrella sampling. The entropy can then be estimated by comparison with activation energy obtained from static structures. To obtain reliable values for the free energy in systems as flexible as the ribosome, extensive sampling is unavoidable. This forces us to use the semiempirical SCC-DFTB method rather than more accurate DFT calculations. We used umbrella sampling simulations along ξ CH = ξ C + ξ H with umbrella integration analysis. We obtained a free energy of activation ∆ ‡ G = 228.1 ± 3.4 kJ mol -1 based on the SCC-DFTB method as QM method. As in the static SCC-DFTB calculations, this value is significantly higher than the experimental value of ∆ ‡ G = 63 kJ mol -1 . The reactant state was found at ξ CH = -5.349 Å and the transition state at ξ CH = 0.135 Å. A tetrahedral intermediate was not found on the free-energy surface. Possible values of ξ C and ξ H consistent with the result of the one-dimensional free-energy simulations along ξ CH are indicated in Figure 6 by a red line. The location of the transition state on the free-energy surface is, thus, fully consistent with the locations of the transition states of the static calculations.

Note that the activation-free-energy obtained from umbrella sampling simulations is not related to a particular snapshot. It was obtained from 740 ps of MD-simulation; significantly longer than the time interval of 90 ps along which the snapshots were taken.

These value for ∆ ‡ G has to be compared to ∆ ‡ E obtained with the same Hamiltonian as the umbrella sampling simulations. Using the snapshots described above, we obtained activation energies from 156.2 to 217.6 (and in one case even 280.4) kJ mol -1 . Assuming again that the exponential average is the relevant one, this results in an entropic contribution of T ∆ ‡ S = -66.2 kJ mol -1 . Even if we consider the direct average of ∆ ‡ E for the enthalpy, this leaves an entropic contribution of T ∆ ‡ S = -33.1 kJ mol -1 , while experimentally it is roughly zero.

The apparent discrepancy with the experimental data can be explained by one or more of the following reasons. (1) There may well be finite-temperature effects additionally to the entropy which are not covered by ∆ ‡ E obtained from static calculations. (2) The zero-point energy (ZPE) may lower the barrier. An estimate for this contribution has been obtained by the harmonic approximation for the vibrational zero-point contributions of only the atoms of the QM part for snapshot 1. In this case, the ZPE lowers the barrier by 13.7 kJ mol -1 . (3) The SCC-DFTB/Charmm Hamiltonian may be inappropriate for the system. (4) Here, we only sampled the direct mechanism. While the activation energies for the direct mechanism and the proton-shuttle mechanism were found to be similar, the free energies of activation may be different. Even if both mechanism occur competingly in the real system, this fact would increase the activation entropy. To address these questions, further simulations will be necessary.

Summary

Our QM/MM simulations of the formation of the peptide bond in the ribosome have shown the importance of a high concentration of ions in this ribozyme. By averaging over 10 snapshots taken from molecular dynamics simulations, we found reaction barriers in good agreement with the experimental data. The variation between the snapshots was found to be substantial, indicating a high degree of flexibility of the environment of the PTC. Comparison to the gas-phase reaction shows that the ribosome reduces the barrier by electrostatic influence of the environment rather then by just favourably positioning of the reactants. The high concentration of mobile ions (counter-ions to the negatively charged rRNA and tRNA) in the ribosome, about 25 times the physiological concentration, was found to be key to the catalytic activity of the ribosome. This finding can probably be generalised to other ribozymes as well.

We studied two reaction mechanisms, a direct proton transfer including a 4-membered transition state, and a proton shuttle via the O2'-H group of tRNA in a 6-membered transition state. Both were found to have similar activation energies. They may compete in the real system.

A tetrahedral intermediate, an energy-minimum in which the C-N bond of the peptide bond to be formed already exists, but the C-O bond is not yet broken, was found to be absent on the potential energy surfaces of most snapshots, as well as on the free-energy surface.

The free-energy simulations for the direct proton transfer mechanism resulted in a significantly higher free energy of activation than the potential energy barrier. This may be due to additional thermal effects which are not covered by the static potential energy calculations, or by a reduced free-energy barrier for the proton-shuttle mechanism. Further investigations, which are out of the scope of this work, are required to address this question.

In summary, we have found two possible mechanisms of peptide bond formation in the ribosome. The contribution of the ribosome to lowering of the energy barrier was found to be its increased concentration of mobile ions rather than just spatial alignment of the reactants. This also puts the catalytic mechanism of other ribozymes into a new light.

Effect of the basis set

To test the effect of the basis set, the energy differences between the reactant state and the transition state of the direct mechanism were calculated with different basis sets. These calculations were done during feasibility studies preceding the current work. Thus, they were done with a smaller model with only those residues free to move that had at least one atom within 3 Å of the QM-region. Thus, they do not directly correspond to one of the snapshots in the paper. Nevertheless, the results given in Table 1 show the suitability of the chosen basis set. One can clearly see from Table 1 that polarisation functions are essential for obtaining reliable energetics. However, a DZP basis set is sufficient.

Effect of the size of the QM region

The energy differences between the same states as for the basis set convergence tests above were calculated with different sets of atoms included in the QM part. In these calculations, only residues with at least one atom within 3 Å of the reactive centre (the two amino acids) were optimised. In all energy differences used in the main paper, more atoms were optimised. The standard QM region as used in all results given in the main manuscript contained 57 atoms. The energetic effect is depicted in Figure 1. Explanations of the different QM regions are:

• 26 atoms: in order to push the limit, only the backbone of the two amino acids and C3' with its hydrogen atom were included in the QM part. In contrast to the 35-atom model, the rest of the P-site ribose ring was treated by MM.

• 35 atoms: only the P-site ribose ring, the P-site amino acid (glycine), the backbone of the A-site amino acid, and C3' (and the hydrogen bound to it) of the A-site ribose ring were included in the QM part. In contrast to the 44-atom model, the rest of the A-site ribose ring was treated by MM.

• 44 atoms: the phenyl ring of the phenylalanine residue was treated by MM. Somewhat surprisingly, this reduces the energy difference by about 6 kJ mol -1 .

• 68 and 79 atoms: the A-site phosphate, and both phosphates (A and P sites) were treated by QM, respectively. This only has a minor effect.

• 83 atoms: the purine bases of both Ade 76 were treated by QM. This slightly reduces the energy difference by 3 kJ mol -1 .

• 105 atoms: both phosphates and both purine bases were treated by QM. Also only a minor effect of 3 kJ mol -1 compared to the standard QM region.

Overall, one can clearly see that for sensibly chosen sizes of the QM regions, the result does not significantly depend on the size.

Figure 1 .

 1 Figure 1. Atoms in the active center of the ribosome, the peptidyl transferase center (PTC). The blue dashed line separates the A-site and the P-site. The two short red arrows indicate attacking lone pairs in the hydrogen exchange mechanism, the longer green arrow indicates the direct mechanism.

Figure 3 .

 3 Figure 3. The transition-state structure of the lowest-energy snapshot (4). Indicated are atoms and residues which interact with the reaction core. The colour code is: H (white), C (gray), N (blue), O (red), Na (yellow). The structure is truncated at numerous positions for better visibility.

Figure 4 .

 4 Figure 4. Electrostatic influence of individual residues on the activation energy (∆∆ ‡ E i ) obtained at snapshot 4. Positive values indicate an increase of the activation energy by the respective residues, negative values indicate a decrease.

Figure 5 .

 5 Figure 5. All residues with |∆∆ ‡ E i | > 12 kJ mol -1 are shown as thick sticks. Numbers correspond to ∆∆ ‡ E i in kJ mol -1 . Additionally, both termini of the tRNA, which contain the active centre, are shown (ball-and-stick). The 4 atoms mainly included in the reaction (NCHO) are shown in larger spheres for better orientation. Residues shown in green stabilise the transition state, residues shown in red destabilise the transition state with respect to the resting state. Higher colour saturation indicates a stronger effect.

Figure 6 .

 6 Figure 6. Location of transition states in DFT and SCC-DFTB (QM/MM each) compared to umbrella sampling (US) with SCC-DFTB. The circles indicate the transition states with the lowest activation energies for the respective method. The diagonal line indicates possible location of the transition state in the one-dimension umbrella sampling simulation.

  Supporting information to:Ribosome catalyses peptide bond formation by providing high ionic strength by Johannes Kästner and Paul Sherwood

1 Page

 1

Figure 1 : 2 Page 19 of 19 URL

 1219 Figure 1: Relative energy depending on the size of the QM region.

Table 1 .

 1 Energetics (activation enthalpies and entropies) in kJ mol -1 of peptide-bond formation in water, in a model system, and in the full ribosome at T = 298 K.

	9]

Table 3 .

 3 The transition states (geometric data in Å and activation energy in kJ mol -1 ) of the ten snapshots optimised with BP86/Charmm.

	Snapshot d(C-O) d(C-N) d(N-H) d(H-O2') d(O2'-H2') d(H2'-O)
	1	1.992	1.608	1.081	1.710	1.335	1.130
	2	2.119	1.610	1.067	1.808	1.198	1.251
	3
	4	1.951	1.613	1.066	1.818	1.216	1.236
	5	2.186	1.585	1.071	1.772	1.120	1.376
	6	2.229	1.552	1.111	1.584	1.139	1.335
	7	1.992	1.587	1.114	1.555	1.202	1.256
	8	1.917	1.704	1.083	1.647	1.231	1.213
	9
	10	1.956	1.586	1.158	1.507	1.568	1.040
	avg.	2.043	1.606	1.094	1.675	1.251	1.230
	RS	1.357	-	1.031	-	0.999	-
	PS	-	1.365	-	0.987	-	0.980

Table 4 .

 4 The geometric data in Å of the transition states for the proton shuttle mechanism optimised with BP86/Charmm.

  The low spread in th C-N distance of the transition states indicates the existence of a tetrahedral intermediate, a geometry where the C-N bond is (partially) formed, but the N-H bond (and the C-O) bond not yet broken. However, when optimising downwards from the transition states, in only one snapshot (number 4) was it possible to locate a tetrahedral intermediate. In all other cases, the system converged towards the reactant state. The tetrahedral intermediate found for snapshot 4 is only 3.3 kJ mol -1 higher in energy than the RS. The C-N distance is with 1.668 Å still slightly longer than in the transition state,

	October 23, 2009 15:2 Molecular Physics rib
	3.2. Tetrahedral intermediate
	d(C-O)=1.452 Å, d(N-H)=1.041 Å, and d(O-H)=2.086 Å. Of course, this does not mean
	that no tetrahedral intermediate exists on the potential energy surfaces of the other snap-
	shots.
	o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

Table 5 .

 5 Geometric data of the transition states for the direct mechanism in Å, activation energy, and energy of the tetrahedral intermediate E TI (where existing) in kJ mol -1 optimised with SCC-DFTB/Charmm. MAE and MSE refer to the mean absolute error and the mean signed error in Å, respectively, compared to the DFT(BP86)/Charmm geometries reported in Table3. snapshots 2, 8, and 10.

	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

Table 1 :

 1 Relative energy with different basis sets. Single point refers to the geometry optimised at the DZP level. Relaxed refers to a structure optimised in the given basis set. Energies in kJ mol -1 . The names of the basis sets correspond to the keywords in GAMESS-UK.

	Basis set	single point relaxed
	SV AHLRICHS	99.58	96.47
	DZ AHLRICHS	97.79	94.57
	DZP	-	154.6
	DZP AHLRICHS	147.69	147.96
	TZVP AHLRICHS	162.94	162.94

Taylor & Francis and I.T. Consultant

Acknowledgment. We acknowledge support by HECToR, and CSED of STFC for access to its IBM Blue Gene/P system. J.K. thanks the German Research Foundation (DFG) for financial support of the project within the Cluster of Excellence in Simulation Technology (EXC 310/1) at the University of Stuttgart.