text:
# Crucial stages of protein folding through a solvable model: predicting target sites for enzyme-inhibiting drugs.

## I Introduction

One of the open fundamental questions in molecular biology is how to predict the folded state of a protein from the knowledge of its sequence. Despite a large increase in available computing power in the past years, it has been impossible to answer this question by means of computer simulations of various degrees of complexity and detail. However an increasing amount of experimental Fersht95 ; sh3A ; sh3B ; Chiti99 ; Plaxco  and theoretical results Miche99b ; M00 ; baker ; Clem ; cieplak  supports the view that the folding of natural proteins into their native state is largely influenced by the native-state topology (for a brief review see bakernature ). Accordingly, the folding process is regarded as a well defined sequence of obligatory steps to be taken in order to reach the native state. Even if protein sequences have evolved to fold efficiently, the kinetics en-route to the native state might be hindered by the realization of particularly difficult (rate-limiting) steps, such as the formation of non-local amino acid interactions (contacts) that usually requires the overcoming of large entropy barriers. Some non-local native contacts are rather crucial for the folding process, because their formation helps establishing further native interactions and leads to a rapid progress along the folding pathway until another barrier is met. Their formation is associated to bottlenecks for the entire folding process. Strikingly, the amino acids involved in such crucial contacts are those for which the largest changes in the folding kinetics are observed in site-directed mutagenesis experiments Fersht95 , as first proven for CI2 and Barnase Miche99b . This suggests that protein sequences have been carefully optimised so to exploit the conformational entropy reduction accompanying the folding process Funnel2  through the selection of the key amino acids. The number and importance of bottlenecks depends significantly on several factors. Among the most important are the contact order of the proteinbaker  and whether it folds in two or more stages Jackson .

In previous studies hiv ; Settanni , we have shown how the most delicate folding stages can be identified within a molecular dynamics approach, by monitoring the formation probability of native and non-native contacts from the unfolded to the native state. This can either be done as a function of time at a fixed temperature around the folding temperature or working at thermal equilibrium for a succession of decreasing temperatures (annealing). In principle, the two approaches need not to be equivalent but, for the quantities we have investigated, they give consistent results. Then, concerning the indentification of crucial contacts, one can safely concentrate on studying thermodynamic equilibrium at various temperatures. The main limitation of Molecular Dynamics (MD) and MonteCarlo (MC) simulations, especially for long protein chains, is that they are extremely time-demanding and plagued with statistical errors which can affect the predictions based on the study of the relative sensitivity of contact formation. Therefore it would be highly desirable to develop a suitable theoretical model, amenable to a deterministic (and computationally fast) treatment, thus resulting in a deeper understanding of the problem. Ideally, such a model should encompass all the “necessary ingredients” that usually are included in computer simulations: peptide chain constraints, effective interactions between residues, favourable monomeric positions, etc. In the following we describe a recently developed theoretical scheme gaussian , that, while being very simplified and approximate compared to other schemes based on MD or MC simulations, can be treated analytically, leading to expressions that can be evaluated exactly. The calculated quantities rival those obtained through more sophisticated but computationally demanding MC and MD techniques. The purpose of the present paper is to show how the model can be employed to yield helpful observables to identify the folding bottlenecks. In particular we apply the method to the HIV-1 Protease (HIV-1 PR), an enzyme which is crucially involved in the HIV infection condra . In general, the accurate knowledge of bottlenecks has important pharmaceutical ramifications because their knowledge may be exploited in a rational drug design. Due to the large amount of available clinical data, HIV-1 PR is a natural choice for a stringent test of our automated predictive scheme.

## II Theory

The model we adopt builds on the importance of the native state topology in steering the folding process, that is in bringing into contact pairs of amino acids that are found in interaction in the native state. A primary quantity of interest that we shall calculate is the probability that a given native contact is established at a definite stage of the folding process. Probably, the oldest attempt to calculate such quantity dates back to Flory who tried to estimate the probability $`p_{ij}`$ that two sites $`i`$ and $`j`$ in a long harmonic chain (the peptide) are in contact flory . The approximation introduced by Flory was to neglect correlations between residues, which amounts to considering the chain embedded in a highly-dimensional space. As a result, the $`p_{ij}`$’s are a decreasing function of the sequence separation $`|ij|`$. Clearly, this approximation is not apt to pinpoint the key folding sites, since it exploits the native topology at the simplest level; in fact it takes into account only the contact order of native interactions. The Flory approach, however, can be refined by incorporating correlations between the formation of pairs, triplets etc. of contacts chan90 ; goddard ; CT2 . Here we use a recently introduced energy function that allows to calculate the $`p_{ij}`$’s within a self-consistent analytic scheme. The strategy is similar in spirit to that of Go and Scheraga Go  where only the formation of native interactions is energetically rewarded and is common to all recent approaches which exploits the native state topology Miche99b ; M00 ; baker ; Clem ; cieplak .

We describe the proteins by the coordinates $`𝐫_i`$ of the Cα atom of the $`i`$-th amino acids. The simplified energy functional for the chain of $`N`$ residues is

$$H=\frac{KT}{2}\underset{i=1}{\overset{N1}{}}(𝐫_{i,i+1}𝐫_{i,i+1}^0)^2+\frac{1}{2}\underset{ij}{}\mathrm{\Delta }_{ij}[(𝐫_{ij}𝐫_{ij}^0)^2R^2]\theta _{ij}$$
(1)
where $`K`$ is the strenght of the peptide bonds, assumed to be harmonic, and $`T`$ is the absolute temperature in units of the Boltzmann constant.

The relative position between amino acid centroids is denoted by $`𝐫_{ij}=𝐫_i𝐫_j`$ and the corresponding native positions are indicated with the superscript $`0`$. $`\mathrm{\Delta }`$ is the contact matrix, whose element $`\mathrm{\Delta }_{ij}`$ is 1 if residues $`i`$ and $`j`$ are in contact in the native state (i.e. their Cα separation is below the cutoff $`c=6.5`$ Å) and 0 otherwise. The matrix $`\mathrm{\Delta }_{ij}`$ along with the set $`𝐫_{ij}^0`$ encodes the topology of the protein. The factor $`\theta _{ij}`$ has the form

$$\theta _{ij}=\mathrm{\Theta }(R^2(𝐫_{ij}𝐫_{ij}^0)^2)$$
(2)
where $`\mathrm{\Theta }(x)`$ is the unitary step function and $`R`$ is a distance cutoff defining the range of the interaction between non consecutive amino acids. In standard off-lattice approaches, the interaction $`V(d)`$ between non-bonded amino acids at a distance $`d`$, is taken to be a square well potential, or some type of Lennard-Jones interaction. Our choice in Eq. (1) is a sort of “harmonic well” which, while being physically sound and viable, is suitable for a self-consistent treatment, as explained below. The location of the outer rim of the well is controlled by $`R`$, which can be set to a few Angstroms ($`R=3`$ Å in the present study) to penalise conformations where the the separation of two residues differs too much from the native one. In the native state each $`\theta _{ij}`$ is close to $`1`$ while in the denaturated state case $`\theta _{ij}`$ are usually negligible.

While the present form of the model does not accurately describe the effects of self-avoidance this does not lead to a qualitatively wrong behaviour in the highly-denatured ensemble (large $`T`$). The treatment of steric effects becomes progressively more accurate as temperature is lowered. In fact, the model guarantees that the native state is the true ground state and therefore protein conformations found at low temperature inherit the native self-avoidance. The connectedness of the chain, as well as its entropy, are captured in a simple but non-trivial manner. The most significant advantage of the model is that it can be used to explore the equilibrium thermodynamics without being hampered by inaccurate or sluggish dynamics.

Two limit cases of model 1 are worthy of notice. In the absence of any bias towards the target structure (i.e. when both $`\mathrm{\Delta }_{ij}`$ and the $`\{r^0\}`$’s are removed) the model reduces to the standard gaussian polymer model whose behaviour is exactly known flory ; Klo99 . Furthermore, the limit when $`T0`$ (when all native contacts are established and the bonded-energy term fluctuations are negligible) the model reduces to the gaussian network model that has been introduced and used to study the near-native vibrational properties of several proteins bah97 ; hiv99 ; kes00 ; ani01 .

The thermodynamics of the model is fully determined by the partition function

$$𝒵(T)=\underset{i=1}{\overset{N}{}}d^3r_i\mathrm{exp}(H/T).$$
(3)
In the integral of equation (3) and in the following, it is always meant that translational invariance is always explicitly broken by fixing, for example, the centre of mass of the system (see Appendix).

The integral (3) is still hard to treat analytically, due to the presence of non-quadratic interactions in the last term of Hamiltonian (1). We thus perform a further, but non-trivial, simplification by replacing $`H`$ with the variational hamiltonian $`H_0`$

$$H_0=\frac{KT}{2}\underset{i=1}{\overset{N1}{}}(𝐫_{i,i+1}𝐫_{i,i+1}^0)^2+\frac{1}{2}\underset{ij}{}\mathrm{\Delta }_{ij}[(𝐫_{ij}𝐫_{ij}^0)^2R^2]p_{ij}$$
(4)
where the factors $`\theta _{ij}`$ are now subtituted by parameters $`p_{ij}`$ independent of the coordinates. Due to its quadratic form, the model described by Eq. (4) can be solved with the standard techniques for Gaussian integrals. Such parameters have to be optimally determined so to ensure self-consistency:

$$p_{ij}=\mathrm{\Theta }((𝐫_{ij}𝐫_{ij}^0)^2R^2)_0.$$
(5)
The symbol $`\mathrm{}_0`$ indicates that the thermal averages are performed through the Hamiltonian $`H_0`$. Now in such self consistent approach the problem is fully solved and we can compute the resulting partition function from which we extract all the thermal properties and averages. In particular the logarithm of the partition funzion, $`\mathrm{ln}(𝒵)`$, has the following explicit expression:

$$\mathrm{ln}(𝒵)=\frac{3(N1)}{2}\mathrm{ln}(2\pi )\frac{3}{2}\mathrm{ln}N\frac{3}{2}\mathrm{ln}(det{}_{}{}^{}M)+\frac{R^2}{2T}\underset{lm}{}\mathrm{\Delta }_{lm}p_{lm}$$
(6)
where the matrix $`M`$ is defined as:

$$M_{i,j}=\{\begin{array}{cc}\hfill K(2\delta _{i,1}\delta _{i,N})+2_l\mathrm{\Delta }_{i,l}p_{i,l}/T& \mathrm{for}i=j\hfill \\ & \\ \hfill 2p_{i,j}\mathrm{\Delta }_{i,j}/TK\left[\delta _{i,j+1}+\delta _{i,j1}\right]& \mathrm{for}ij.\hfill \end{array}$$
(7)
and the prime in (6) denotes that the zero eigenvalue of $`M`$ has to be omitted (see Appendix).

The quantities $`p_{ij}`$ in Eq. (5) represent precisely the occurrence probability of a contact between residues $`i`$ and $`j`$ and indicate the frequency with which that native contact is established. At thermal equilibrium their dependence on temperature reflect the status of compactness of the protein molecule. For instance, well below the folding temperature, $`T_F`$, each $`p_{ij}(T)`$ is expected to assume a value close to unity, as all native contacts are already formed. Instead, for temperatures much larger than $`T_F`$, all $`p_{ij}(T)`$ tend to be very small, reflecting the low propensity of the protein to establish contacts.

Thermodynamics quantities can be easily derived from the $`p_{ij}`$’s. Another quantity necessary to characterize the folding transition is the specific heat, which exhibits one or more peaks in correspondence of significant structural rearrangements of the protein conformation. Since every energy change is mainly associated to the formation of native interactions, we address the question of which native contacts contribute mainly to the peak(s) of the specific heat. A clear answer to this question is readily found in the temperature behaviour of frequencies $`p_{ij}`$. Indeed, each $`p_{ij}(T)`$ exhibits a sigmoidal shape, and the modulus of its derivative develops a sharp maximum in correspondence of the point of inflection (crossover temperature). The importance of every native contact $`ij`$ turns out to be characterized by the crossover temperature and the maximum slope of its $`p_{ij}`$, which can be regarded as an indicator of its degree of cooperativity. In fact, the most important contacts are those with high crossover temperature and associated high cooperativity.

This fact allows a complete identification and classification of the bottlenecks, because we are now able to indentify those contacts that are termodinamically relevant to peaks and shoulders of the specific heat.

## III Application to HIV-1 Protease

The human immunodeficiency virus (HIV) encodes a protease, HIV-1 PR, whose inhibition is crucial to prevent the maturation of infectious HIV particles condra . The role of the Protease in the infection spreading is to act as ”molecular scissor” cleaving inactive viral polyproteins into smaller, functional proteins. In the presence of protease inhibitors, viral particles are unable to mature and are rapidly cleared. Extensive clinical trials have lead to the development of five HIV-1 PR inhibitors approved by the Food and Drug Administration (FDA): Saquinavir mesylate (SAQ), Ritonavir (RIT), Indinavir sulfate (IND), Nelfinavir mesylate (NLF) and Amprenavir (APR) BIOCH88 . Such drugs are particularly effective in short-term treatments, while their long-term efficacy is limited by resistance. Indeed mutants resistant to protease-inhibitors can emerge in vivo already after less than one year condra . Table 1 summarises the list of HIV-1 PR known mutating sites causing drug resistance.

In an earlier work, the study of the near-native harmonic vibrations of the HIV-1 PR has shown that a number of sites that are paramount to the stability of the native enzyme are close to some of the residue of Table 1 hiv99 . The self-consistent scheme of eqn. 4 allows to extend this result by modelling the partially-folded ensemble at finite temperature.

In particular, we will be concerned in the characterization of such ensemble near the folding transition temperature. The motivation to do so stems from a recent study hiv  where we have shown that such mutating amino acids correspond, with high statistical significance, to sites involved in the folding kinetic bottlenecks. The rationale for this finding is that the most effective drugs can be eluded only by mutations occurring in correspondence of the key sites. Due to the sensitivity of the folded native conformation to these sites, only fine-tuned mutations are allowed in correspondence of these sites. Such mutation have to result in a native-like enzymatic activity and in the avoidance of the drug action. These constraint act as a severe selective pressure on the mutated proteases that the HIV virus is able to express. As a result, the mutations that will ultimately cause drug-resistence are expected to occur in correspondence of the crucial sites. These residues are heavily influenced by the native topology, and hence should display little dependence on the particular (effective) drug to be eluded.

It is therefore our purpose to apply the scheme introduced in the previous section and identify the key residues within our topology-based scheme. The method, being completely analytic, is free from statistical uncertainty, common to all MC and MD simulation methods, or from difficulty (due to spatial restraints) to reach the target native state below the folding temperature.

## IV Results and discussion

The structural model at the basis of our analysis is the free enzyme condra . It is a homodimer with C2 symmetry, each subunit being composed by 99 residues (Fig. 1). Previous studies hiv  have shown that geometrically important residue positions can be obtained considering a single monomer. Indeed the specific heat of the whole homodimer on decreasing the temperature shows a peak in correspondence of the folding of each sub-unit and then at lower temperature another peak signals the aggregation of the two sub-units. Thus, in the following, we will be concerned only with a single monomer. The specific heat is obtained through numeric differentiation of the average internal energy, which has the following explicit analytic expression in terms of the $`p_{ij}(T)`$’s and the quantities introduced before:

$$E=\frac{3(N1)T}{2}\frac{R^2}{2}\underset{ij}{}\mathrm{\Delta }_{ij}p_{ij}(T).$$
(8)
The study of Go and Scheraga Go  showed that systems described by energy-scoring-functions that reward the formation of native contacts display cooperative (all-or-none) folding transitions with an associated peak(s) in the specific heat. Consistently with these expectations, the specific heat calculated by differentiating Eq. (8) with respect to $`T`$ shows a single peak, see Fig. 2, thus providing an unambiguous criterion for identifying the folding transition temperature, $`T_F`$. The width of the specific heat peak at the folding transition in Figure 2 is larger than the typical one found in experimental Jackson  and theoretical studies Kaya1 ; kaya2 . It is possible to enhance the cooperativity of the transition by intervening on the actual value of $`K`$ in Eq. (1); in fact, a decrease of $`K`$ leads to sharper transitions. An alternative criterion for fixing the value of $`K`$ is provided by its influence on the average amount of native structure that is formed at the native state. Since we are particularly interested in monitoring the progressive establishment of native contacts, we adopted this second possibility to set the value of $`K`$. In fact, by choosing $`K=1/15`$ in (1), we ensure that, at $`T_F`$, the average fractional occupation of native contacts, $`q`$:

$$q=\frac{\underset{i,j}{\overset{}{}}\mathrm{\Delta }_{i,j}p_{i,j}}{_{i,j}^{}\mathrm{\Delta }_{i,j}}$$
(9)
is about 50 % (see Fig. 2), as established in several experiments and numerical studies. The primed summation symbol indicates that the sum is not carried out over consecutive pairs. The degree of native similarity, $`q`$ is a useful overall indicator to monitor the progress towards the native state in a folding process CT ; Funnel5 . While the ultimate quantities of interest are the $`p_{i,j}`$’s, it is useful to consider an intermediate level of description and focus on the whole network of contacts that a given site takes part to. A natural order parameter is provided by the “average environment formation” Finkel ; karp  which, for a generic site $`i`$ is defined as:

$$P_i=\frac{\underset{j}{\overset{}{}}\mathrm{\Delta }_{i,j}p_{i,j}}{_j^{}\mathrm{\Delta }_{i,j}}.$$
(10)
$`P_i`$ is a measure of the fraction of established native contacts the $`i`$-th residue partecipate to (clearly, $`P_i`$ is defined only when the denominator of eqn. 10 is non-zero). The environment profiles for three different temperatures are shown in Fig. 3. The irregular behaviour of the profiles results from a complex interplay of the burial of the sites and the locality of their contacts. The hierarchical formation of secondary structures at high temperature is clearly visible. It is instructive to correlate the location of the sites known to cause resistance to drug treatments (see Table 1) with the features of the profiles. In particular, several mutating sites responsible for drug resistance (see Table 1) can be found in correspondence of the peaks of the environments (see in particular sites 20,63,71,77,84). The most precise way to identify the key residues is, however, through the analysis of the fractional occupation of native contacts and not through the environments, since they only carry averaged information. Typical $`p_{i,j}`$ curves as a function of temperatures are shown in Fig. 4.

As anticipated in section Theory, all $`p_{i,j}`$’s have monotonic sigmoidal shapes which mainly reflect the sequence separation, $`|ij|`$ and the native burial of each of the residues. In general, each contact is established at a different crossover temperature and with different intensity hiv . The data relative to the frequencies of native-contact formation is conveniently summarised in the color-coded contact maps of Fig. 5. A bright red color is used to highlight those contacts with the largest crossover temperatures above $`T_F`$, see Fig. 5a, or highest intensity in Fig. 5b. Both these intuitive notions can be used to identify the key folding contacts. The inspection of Fig. 5 reveals that several kinetic bottlenecks (red regions) are located three-four contacts downstream the three $`\beta `$-turns in HIV-1 PR. In addition, the formation of contacts around residues $`84`$ and $`30`$, despite being so far away along the sequence, appears to be a crucial folding stage since it allows the collapse of the individual secondary structure motifs. It is striking that these results make an excellent parallel with those of Ref. hiv , where long and delicate MD simulations of the unfolding/refolding of HIV-1 PR were carried out using a much more sophisticated energy-scoring function. This provides a cross validation for the robustness of the results obtained both in the stochastic and the present, analytic, scheme. The emphasis is on the exactness of the present approach that allows to determine easily the $`p_{i,j}`$’s with an arbitrary accuracy. The absence of stochastic noise allows to compile Table 2 which shows the top contacts ranked according to crossover temperature and intensity. Sites that are known to cause drug resistance through mutations are highlighted in boldface. It is apparent that a high fraction of the top key folding contacts do, indeed, contain key mutating sites. To test the significance of such matches we compare the number of marked mutating sites contained in each column of Table II with the number of those contained in a randomly compiled table. We expect a random list of $`t`$ elements extracted among $`N`$, $`m`$ of which are marked, to contain an average of $`tm/N`$ marked elements with a square deviation of $`tm(Nm)(Nt)/(N^2(N1))`$. For the case of HIV-1 PR the total number of contacts (excluding consecutive residues) within a cutoff radius of 6.5 Å is $`N=180`$ and the number of those which include at least one known mutating site is $`m=60`$. An analysis of the contacts of Table 2 (selected according to crossover temperature or cooperativity of formation) shows that the number of matches observed among the top sites tipically exceeds that expected from a random choice by a standard deviation (the precise amount depend of how many top ranking contacts are considered. An alternative and apparently more stringent approach is to identify independent groups of highly correlated contacts, and then search for the key residues in each group. To a first approximation, the correlated sets of interacting pairs may be identified with the clusters in the contact map. This leads to define six main groups, the three $`\beta `$-sheets, the helix and the two sets of long-range contacts, around contacts 14-60 and 23-84, respectively (see Fig. 5). The four contacts in each group with the highest intensity of formation above $`T_F`$ are summarised in Table 3. Out of the 24 contacts, 12 of them involve a key site, which is two standard deviations away from the number of matches expected on a random basis ($`7.9\pm 2.1`$). Again, this testifies both the reliability of the general scheme followed here and also its robustness in the different possible implementations.

Interestingly, the results of Table 3 account better than those of Table 2 for the heterogeneous location of the key folding sites. The emerging conclusion is that a complete description of the crucial contacts can be obtained only by monitoring all the key stages of the folding process. In standard MC and MD simulations of protein unfolding/refolding, it is the simulated dynamics that reveals which, and how many, delicate stages exists. In the present approach, the folding process is characterised analytically, thus the complete set of folding bottlenecks follows from the study of distinct groups of interrelated contacts.

Finally, we remark that the determination of the key contacts does not uniquely provide the key folding sites, since two sites are involved in each pairwise contact. This ambiguity can, in several cases, be resolved either by selecting those sites that take part in several crucial contacts, or by examining their distribution on the three-dimensional native structure for clues that may help breaking the ambiguity.

## V Conclusions

We have used an analytical technique to study and characterize the folding process of globular proteins. This deterministic method allows the automated identification of contacts involved in folding rate-limiting steps. As a result, the whole folding process is particularly sensitive to mutations occuring at sites involved in such crucial contacts. We test our scheme and its usefulness in pinpointing the crucial sites by applying it to HIV-1 protease. For this enzyme, extensive clinical trials have allowed the identification of several sites involved in drug-resistance mutations. Such sites have a meaningful overlap with the key folding sites predicted by our scheme with a modest computational effort compared to more sophisticated stochastic simulations techniques. This indicates that the available inhibiting drugs are quite effective since they can be eluded only by mutations of the (sensitive) key sites of the protease.

The proposed approach to identify the crucial residues is quite general and ought to be useful to identify the kinetic bottlenecks of other viral enzymes of pharmaceutical interest, thus aiding the development of novel effective inhibitors.

We expect to focus our future efforts on improving the present approach by taking into account the propensities of different amino acids to form contacting pairs. This limitation can be overcome by introducing physically viable (attractive) pairwise interactions Seno:98:prl ; Sippl95 ; Miyazawa96 ; Maiorov:92 ; stabloc . In the present approach this possibility was deliberately avoided to highlight the influence of the native state topology alone on the kinetic bottlenecks, irrespective of the different chemical nature and strength of the effective amino acid interactions. We expect that the inclusion of such effects, while not distorting the overall picture presented here, may change the relative strength of spatially-close contacts. This may improve the agreement between Table 1 and tables 2-3 by resolving those cases were a site adjacent to a mutating one is selected.

We are indebted to Paolo Carloni for several illuminating discussions and for having stimulated the present work. This work was supported by INFM, Murst Cofin2001.

## Appendix A Appendix

In this appendix we discuss how the translational invariance of a quadratic energy scoring function can be explicitly broken by fixing the center of mass of the system in the origin. The constrained partition function is written as:

$$𝒵=\underset{i=1}{\overset{N}{}}d^3x_ie^{1/2_{i,j}𝐱_iA_{ij}𝐱_j}\delta ^3(\underset{i}{}𝐱_i)$$
(11)
where the matrix $`A`$ incorporates the quadratic dependence of $`H_0`$ in eqn (4) from the space co-ordinates (and also includes the $`1/T`$ factor to yield the usual Boltzmann weight). The translational invariance of $`H_0`$ implies that $`A`$ satisfies the property: $`_jA_{ij}=0`$, which amounts to say that the uniform vector, $`𝐯_1N^{\frac{1}{2}}(1,1,1,1\mathrm{},1)`$ is an eigenvector of $`A`$ with eigenvalue $`\lambda _1=0`$. We assume that $`H_0`$ is invariant only for the simultaneous translation of all the coordinates, $`\{𝐱_i\}`$. In this case all other eigenvalues, $`\{\lambda _{i>1}\}`$ are strictly positive and the corresponding eigenvectors $`𝐯_{i>1}`$ are all orthogonal to the zero mode $`𝐯_1`$.

By rewriting the Dirac-delta constraint as

$$\delta ^3(𝐳)=\underset{c\mathrm{}}{lim}\left(\frac{c}{2\pi }\right)^{\frac{3}{2}}e^{c𝐳𝐳/2}$$
(12)
the partition function takes on the form $`𝒵=lim_c\mathrm{}𝒵_c`$, where

$$𝒵_c=\left(\frac{c}{2\pi }\right)^{\frac{3}{2}}\underset{i=1}{\overset{N}{}}d^3x_ie^{1/2_{i,j}𝐱_iA_{ij}^{}𝐱_j},$$
(13)
where $`A_{ij}^{}=A_{ij}+c`$. It is straightforward to see that $`A^{}`$ admits the same eigenvectors of $`A`$. Only the zero mode eigenvalue will change from zero to $`cN`$, while all others will be unmodified. Upon performing the Gaussian integrations in $`𝒵_c`$ we obtain:

$`𝒵_c`$ $`=`$ $`\left({\displaystyle \frac{c}{2\pi }}\right)^{\frac{3}{2}}\left({\displaystyle \frac{2\pi }{cN}}\right)^{\frac{3}{2}}{\displaystyle \underset{i=2}{\overset{N}{}}}\left({\displaystyle \frac{2\pi }{\lambda _i}}\right)^{\frac{3}{2}}`$ (14)
$`=`$ $`{\displaystyle \frac{1}{N^{3/2}}}\left(2\pi \right)^{\frac{3(N1)}{2}}{\displaystyle \underset{i=2}{\overset{N}{}}}\lambda _{i}^{}{}_{}{}^{\frac{3}{2}}.`$
This shows that $`𝒵_c`$ is effectively independent of $`c`$ and, therefore, the partition function $`𝒵`$ simplifies to

$$𝒵=N^{\frac{3}{2}}\left(2\pi \right)^{\frac{3(N1)}{2}}\left(det{}_{}{}^{}A\right)^{\frac{3}{2}},$$
(15)
where the prime denotes that the determinant is calculated omitting the zero mode eigenvalue.


format: text

text:
# How Well Do We Know the Electromagnetic Form Factors of the Proton?

## 1 Introduction

Elastic form factors measurements probe the charge and magnetization distributions of the nucleon, and provide strong constraints on models of nucleon structure. Prior to the year 2000, all of the high-$`Q^2`$ proton form factor data came from cross section measurements, utilizing the Rosenbluth technique to separate the electric and magnetic form factors, $`G_E`$ and $`G_M`$. A global analysis of the large body of data on elastic electron-proton scattering indicated that $`G_M`$ follows the dipole form, $`G_M=1/(Q^2+0.71)^2`$, with $``$5% deviations (Fig. 1). While the measurements of $`G_E`$ at high $`Q^2`$ are significantly less precise, the extracted ratio of $`G_E`$ to $`G_M`$ is roughly constant.

More recent $`G_E/G_M`$ results, from measurements of the polarization of the recoil protons, show that $`G_E`$ falls more rapidly with $`Q^2`$ halla  . There are significant deviations from the global Rosenbluth analysis above $`Q^2`$=1 GeV2, as shown in Fig. 1. This discrepancy indicates that something is wrong with one of these two techniques, or one or more of the experiments. I will briefly review the two techniques, focusing on potential sources of systematic uncertainties. Next, I will give an overview of a recently completed JLab experiment designed to test the compatibility of the two techniques. Finally, I will present an analysis of the previous data designed to look for possible sources of the discrepancy between the two techniques. The goal is to determine what kind of errors would have to exist in the cross section measurements to explain the discrepancies in the form factors ($`>`$100% error in $`G_E`$, few % in $`G_M`$) and what impact these errors might have on our knowledge of the proton form factors, as well as other measurements.

## 2 Extractions of the Elastic Form Factors

Rosenbluth extractions (L/T-separations) of the form factors are performed by measuring elastic electron-proton scattering at fixed virtual photon energy and momentum ($`\nu ,\stackrel{}{q}`$), while varying the electron energy and scattering angle to vary the virtual photon polarization, $`ϵ`$. The reduced cross section, $`\sigma _R`$, can then be expressed in terms of the form factors, which depend only on $`Q^2`$ ($`=q^2\nu ^2`$):

$$\sigma _R\frac{d\sigma }{d\mathrm{\Omega }}\frac{ϵ(1+\tau )}{\sigma _{Mott}}=\tau G_M^2(Q^2)+ϵG_E^2(Q^2),$$
(1)
where $`\tau =Q^2/(4M_p^2)`$. $`G_M`$ is then extracted from the reduced cross section at $`ϵ=0`$, while $`G_E`$ is extracted from the $`ϵ`$-dependence. Due to the $`ϵ/\tau `$ weighting of the electric term relative to the magnetic term, the contribution from $`G_E`$ decreases as $`1/Q^2`$ for a fixed ratio of $`G_E/G_M`$, and isolating the contribution of $`G_E`$ becomes increasingly difficult as $`Q^2`$ increases. Because of this, it is important, and increasingly difficult as $`Q^2`$ increases, to ensure that the $`ϵ`$-dependent systematics do not overwhelm the uncertainties in the extraction. Because $`ϵ`$ is correlated with beam energy, scattering angle, and scattered electron energy for a fixed value of $`Q^2`$, and because the Mott cross section varies rapidly with angle at fixed $`Q^2`$, there are several potential sources of $`ϵ`$-dependent errors which might affect the extracted form factors.

An alternative technique involves using polarized electrons and measuring the polarization of the recoiling proton. The ratio of the transverse to longitudinal components of this polarization is directly related to $`G_E/G_M`$. Measuring a ratio of two polarization components means that uncertainties in the cross section, beam polarization, and detector analyzing power all cancel out, significantly reducing the dominant sources of systematic uncertainty. While this method is clearly superior at large $`Q^2`$ values, the discrepancy between the Rosenbluth and recoil polarization measurements occur at $`Q^2`$ values as low as $``$1 GeV2, where both techniques give precise measurements.

## 3 Jefferson Lab Experiment E01-001

Jefferson Lab experiment E01-001 was designed to test the consistency of two techniques. In the global analysis of the cross section measurements, a major concern is the relative normalization of the difference experiments. While a normalization factor for each experiment is determined from the best fit to the entire data sets, there is still room for the normalizations to vary which could lead to a change in the extracted form factors. Single experiment extractions eliminate the effect of normalization uncertainties, which can dominate the uncertainty in a global analysis. For the existing single experiment L/T extractions, the dominant sources of uncertainty come from possible errors that could be correlated with $`ϵ`$. The largest such uncertainties come from rate dependent corrections, as the Mott cross section varies rapidly with scattering angle for fixed $`Q^2`$, and kinematic-dependent corrections, which may be especially important for extremely large or small values of the scattered electron momentum, where effects such as multiple scattering or non-linearities in magnetic spectrometers may become important. The goal of E01-001 was to make an extremely precise Rosenbluth extraction of $`G_E/G_M`$ in a single measurement, with careful checks on systematic uncertainties. Data was taken at three $`Q^2`$ values from 2.6 to 4.1 GeV2 to see if the two techniques give consistent results for the ratio of $`G_E/G_M`$ in the region where the previous L/T separations disagree with the new polarization transfer measurements.

E01-001 differs from previous experiments in two main respects: first, we measured the elastic cross section by detection of the struck proton rather than the scattered electron, and second, we made simultaneous measurements at high and low $`Q^2`$ values for each beam energy. At fixed $`Q^2`$, the proton momentum stays fixed, and so there are no momentum-dependent corrections for the protons. In addition, the rate-dependence is dramatically reduced when detecting protons. For the kinematics of our experiment, the proton cross section variation is $`<`$50% over the full $`ϵ`$-range at each $`Q^2`$ value. If the electrons were detected at the corresponding kinematics, the cross section would vary by 1-2 orders of magnitude between high and low values of $`ϵ`$. Finally, the cross section is typically a factor of 2-4 less sensitive to uncertainties in beam energy and scattering angles, making the measurement less sensitive to small uncertainties in the scattering kinematics.

The main measurement is compared to a normalization point at low $`Q^2`$ (0.5 GeV2), where the ratio of $`G_E/G_M`$ is well known and, more importantly, where the $`ϵ`$-range of the measurement is very small ($`\mathrm{\Delta }ϵ0.05`$). Because $`ϵ`$ is nearly constant for the low $`Q^2`$ measurement, the reduced cross section has a very small, and well known, $`ϵ`$-dependence, and we can use this data as a luminosity monitor to correct for the beam current, target thickness, and target density fluctuations.

The experiment was run in May 2002, and data was taken at $`Q^2`$=2.64, 3.2, and 4.1 GeV2. The two lower $`Q^2`$ points will give the most precise results, and each should provide a better than 7$`\sigma `$ differentiation between form factor scaling, $`\mu _pG_E/G_M=1`$, and the decrease in $`G_E/G_M`$ measured by the polarization measurements. The $`Q^2=4.1`$ point has larger uncertainties, due to both reduced statistics and increased background, but should still give a 4$`\sigma `$ separation between the two results. The existing data are shown in Fig. 1 along with projected results for E01-001 under two different assumptions for the extracted ratio. Because the rate-dependence and kinematic sensitivities are greatly reduced compared to previous measurements, the errors are dominated by uncertainties which are uncorrelated at the different $`Q^2`$ values, making the three measurements largely independent. In a typical L/T-separation measurement, any rate-dependent or momentum-dependent errors would be likely to give similar effects for the extractions at all $`Q^2`$ values, and might change the overall trend of the data for all $`Q^2`$ points.

## 4 Global Reanalysis of Cross Section Data

Because of the difficulties in performing L/T-separations at high $`Q^2`$, where $`G_E`$ contributes only a few percent to the cross section, the recoil polarization technique is more reliable at large $`Q^2`$. However, the disagreement between the two techniques extends to $`Q^21`$ GeV2. In this range, the electric form factor contributes 20-30% to the cross section at $`ϵ=1`$, and Rosenbluth separations give a precise measurement of $`G_E/G_M`$. As mentioned above, the experimental normalization factors and $`ϵ`$-dependent systematics are very important for these extractions, and it is possible that the disagreement is due to problems with some subset of the cross section measurements in the global analysis. In fact, it has been noted that the results from different experiments that have extracted $`G_E`$ are inconsistent at large $`Q^2`$ values. We will first present the individual Rosenbluth separations, and show that the inconsistencies between different data sets appear to be a combination of the assumptions in the analyses along with an error in one of the data sets. We then will try to determine if the disagreement between the two techniques can be explained by a simple problem with one or more data sets in the analysis, or any problems in the analysis itself.

### 4.1 Analysis of Individual Rosenbluth Measurements

Figure 2 shows the measurements of $`\mu _pG_E/G_M`$ for several different Rosenbluth separation measurements walker ; single\_lt . The data has been binned into five $`Q^2`$ bins, and the solid and dashed lines show the weighted average for each bin along with the 1$`\sigma `$ uncertainties. While the combined data set shows approximate form factor scaling, with a decrease of $``$10% at moderate $`Q^2`$ values, the individual measurements have significant scatter about this average. Comparing each data point to the average value for its $`Q^2`$ bin, we get $`\chi _\nu ^2=1.26`$ for 40 degrees of freedom (i.e. a confidence level (CL) of 13%). The disagreement is more obvious if we focus on the high $`Q^2`$ data: $`\chi _\nu ^2=1.63`$ for 17 degrees of freedom for data above $`Q^2=1.5`$ GeV2 (5% CL). The scatter of these results has been used to argue that the experiments are inconsistent, and that these results should be discarded.

Before concluding that the Rosenbluth extractions are not reliable, we should examine these data more carefully. There are two problems in this comparison of ‘single-experiment’ extractions. First, the Walker data has a correction at small angles that was discovered by a later SLAC experiment, but was not taken into account in the original analysis. Second, the other extractions shown in Fig. 2 are not really single-experiment measurements. In three cases (Litt, Price, and Berger), the values of $`G_E`$ and $`G_M`$ come from combining a new set of cross section measurements with older data at different $`ϵ`$ values. Various procedures have been used to determine a relative normalization between the two experiments, but the uncertainties from this determination are either ignored altogether, or applied without the taking into account the fact that adjusting the normalization of one data set leads to uncertainties that are highly correlated between the different $`Q^2`$ values. For the Bartel and Andivahis data, the form factors are extracted using only the new data, but that data comes from multiple data sets (taken using different spectrometers, or detecting protons rather than electrons). For these experiments, a direct measurement of the relative normalizations at identical kinematics was possible, and so the normalization factors should be better determined than in the previous cases, where the normalization factors had to be determined from data sets that did not have any kinematical overlap. However, while the normalization factors should be better determined, the correlated nature of the uncertainties from this determination was not taken into account. Thus, the extracted form factor ratios shown in Fig. 2 are not a proper basis for determining the consistency of the different cross section measurements.

While the ratios shown in Fig. 2 are correct for the normalizations used in these analyses, the ratio will increase or decrease for all $`Q^2`$ values if these normalization factors are varied. Because the ratio is sensitive to the $`ϵ`$-dependence, and because this dependence decreases with $`Q^2`$, a shift in the normalization factors would have the largest effect at higher $`Q^2`$, and so each of these data sets could be caused to rise or fall with $`Q^2`$ by varying the normalization factors. While it may require a large change in the normalization correction, these analyses can all be adjusted to reproduce the falloff seen in the polarization transfer measurement (e.g. for the most precise measurement, the Andivahis data, a 4% change in the measured correction for the low-$`ϵ`$ data brings the ratio into agreement with the Hall A data, but at the cost of introducing a 6$`\sigma `$ disagreement between the two spectrometers for measurements at identical kinematics). In order to study the consistency of the data, we must avoid the large uncertainties related to the normalization factors. To do this, we must either fully take into account the correlated uncertainties arising from the choice of normalization procedure, or else extract the form factors from single experiments, where these normalization issues do not arise. Alternatively, in a global analysis the normalization factors can be better determined, as the data sets will have significantly more overlap than in the case where two data sets, one at high $`ϵ`$ and the other at low $`ϵ`$, are combined. This should lead to a more precise determination of the normalizations. We can then see if adjustments of the normalization factors within these uncertainties can significantly change the results. In the following sections, we will both analyze the single-experiment Rosenbluth extractions and perform a global analysis to examine the discrepancy with the polarization measurements.

We begin by repeating the extraction of $`G_E/G_M`$ for only those experiments where the $`ϵ`$-range was adequate to perform an L/T-separation using the data from a single detector. For the Walker data, cross sections taken below 15 were excluded to avoid the error from the missing correction. For the Andivahis data, we use only the data from the 8 GeV spectrometer, and exclude the 1.6 GeV spectrometer data (which always provided a single low-$`ϵ`$ point). Of the other experiments, only the Litt and Berger data had enough $`ϵ`$ range to perform a stand alone Rosenbluth extraction. The extractions from this limited data set are shown in Fig. 3. The average again yields a ratio that close to unity, but the data sets are now in better agreement: $`\chi _\nu ^2=1.18`$ for 10 degrees of freedom (30% CL), for data above $`Q^2`$=1.5 GeV2. The average is clearly well above the polarization data; in fact, all 20 data points lie above the Hall A fit.

### 4.2 Global Fit to Cross Section Data

While these few stand-alone extractions are self-consistent, we would like to examine the full body of data to determine if the disagreement between the global analysis and the new polarization results can be at least partially explained by some problem in the data or analysis. The fit may be affected by inclusion of bad data points or data sets. It may even be that while the best fit yields a roughly constant ratio of $`G_E`$ to $`G_M`$, adjustment of the normalization factors for the experiments may bring down this ratio at high $`Q^2`$ without significantly decreasing the overall quality of the fit. To test such explanations of the discrepancy between the two techniques, a new global analysis of the cross section measurements is presented which can be used to test the above possibilities.

The global analysis is largely a repeat of the analysis performed in Ref. walker . We use the same data sets, and perform a combined fit to the electric and magnetic form factors (with $`1/G_E`$ and $`1/G_M`$ parameterized with 6th order polynomials in $`q=\sqrt{Q^2}`$) along with a normalization factor for each of the data sets. However, there are some differences in the data sets and fitting procedures. For the Walker experiment, we exclude the data below 15 degrees, as discussed in the previous section. For the Andivahis measurement, we use the final published cross sections, which were not available at the time of the previous global analysis. For the Andivahis and Berger experiments, we break up the data into subsets, one for each detector configuration. Thus, data taken from 13 experiments is broken up into 16 subsets, each with it’s own normalization constant. Finally, we exclude data below $`Q^2=0.3`$ GeV2 and above $`Q^2=10`$ GeV2, as we are mainly interested in the Rosenbluth results in the region where we have polarization transfer measurements: $`0.5<Q^2<6`$ GeV2. This initial fit gives similar results to the Bosted parameterization bosted  of the Walker global analysis. The top two curves in Fig. 4 show the Bosted fit (dashed line) and the result of the new fit (uppermost solid line).

Removal of the low angle Walker data combined with breaking up the Andivahis and Bartel data sets slightly reduced the extracted ratio, but it is still just 5-10% below scaling, well above the parameterization of the Hall A data (dotted line). In examining the individual experiments, there are no data sets that have unacceptably large $`\chi ^2`$ values compared to the global fit, and no data point that lie beyond 3$`\sigma `$ from the fit. These simple tests do not indicate any clear problems with the data sets. However, if a data set has an $`ϵ`$-dependent systematic uncertainty that is not too large compared to the uncorrelated uncertainties, then this data set may bias the extracted ratio over a range in $`Q^2`$ without having an unusually large $`\chi ^2`$. In addition, the total $`\chi ^2`$ for the fit is quite low: $`\chi ^2`$ = 220.4 for 274 degrees of freedom, indicating that some of the experiments have overestimated the uncertainties. Thus, the simple statistical test above may not be sufficient, and we need additional checks for possible bad data sets that could have a large impact on the overall result.

The first test involved repeating the fit 16 times, with a different data set removed each time. These fits had only very small changes in the extracted ratio. The fit was repeated, this time with the removal of the three data sets whose exclusion caused the largest reduction in the extracted ratio. Even this ‘worst-case’ removal of three data sets leads to a small reduction ($`<`$10%) in the ratio (middle curve in Fig. 4). Thus, if the discrepancy is caused by errors in the cross section measurements, it must be a systematic, $`ϵ`$-dependent error that impacts several data sets rather than just a problem with one or two experiments.

Finally, we wish to see if small modifications to the normalization factors can remove the inconsistency between the two techniques without making the overall fit significantly worse. First, we perform a constrained fit to the data. We use the same 16 data sets as in our original fit, allowing $`G_M`$ and the normalization factors for each experiment to vary, but requiring $`G_E`$ to match the polarization data by constraining the ratio $`\mu _pG_E/G_M=10.13(Q^20.04)`$. In this way, the normalization factors will be optimized to reproduce this ratio, as well as maximizing the consistency between the different data sets. If only small adjustments to the normalization factors are required, then the $`\chi ^2`$ for this fit should be only slightly worse than the unconstrained fit. When the ratio is constrained, the total $`\chi ^2`$ of the fit increases by 60.5 while the number of degrees of freedom increases by 6 (from 274 to 280). While the total $`\chi ^2`$ is still close to one, due to the overly conservative error estimates in some of the data sets, the increase in $`\chi _\nu ^2`$ is 0.20, which is extremely large for a fit to more than 300 data points. Constraining $`G_E/G_M`$ to match the Hall A fit clearly gives too much weight to the polarization data, and so the unconstrained fit was repeated one more time, but this time fitting to both the cross section data and the recoil polarization data from Refs. halla  (bottom solid line in Fig. 4). Again, the $`\chi ^2`$ increase is significant: $`\chi ^2`$ increases by 65 while the number of degrees of freedom increases by 26 (the number of additional data points). In addition to the fact that the overall fit quality is worse ($`\mathrm{\Delta }\chi _\nu ^2=0.15`$ for $``$300 degrees of freedom), the recoil polarization data has larger deviations from the global fit than any of the other data sets.

From the above tests we conclude that it is not possible to explain the discrepancy between the two techniques without significant errors in several data sets, or modifications to the cross section normalization factors that lead to a significantly worse fit. This is not too surprising, as the extractions from single experiments, which do not suffer from uncertainties in the overall normalization, were in agreement with each other but did not agree with the recoil polarization measurements.

## 5 Conclusions

The discrepancy between electric form factors extracted from recent polarization transfer measurements and older cross section data implies either a fundamental flaw in one of the techniques, or a problem in one or more experiments. We have shown that the Rosenbluth extractions of $`G_E`$ from previous measurements are self-consistent, providing that one looks only at analyses that use a single set of data (and after removing the small angle Walker data). In our new global analysis of the cross section data, we find no simple explanation for difference between L/T and polarization measurements. A problem with the cross section measurement would have to be an $`ϵ`$-dependent error involving several different data sets.

A problem with either the polarization transfer or L/T data might have significant implications for other experiments using these techniques. If the problems are shown to be in the cross section data, then the implications are not limited to other L/T measurements. To explain the observed ($`>`$100%) inconsistencies in $`G_E`$, the cross section measurements must have $`ϵ`$-dependent systematic errors on the scale of a few % or more (which could be rate-dependent, angle-dependent, or energy-dependent errors). While recoil polarization measurements can determine the ratio $`G_E/G_M`$, the cross sections are still needed to extract the individual form factors. Errors of a few % or more in an unknown subset of the cross section measurements could lead to errors in the form factors at the few % level, and could even imply similar errors in the $`Q^2`$-dependence. In addition, elastic cross section measurements are often used as a benchmark to determine normalization for a variety of measurements.

Until we know which result is correct, we cannot be certain of our knowledge of $`G_E`$ at $`Q^2>1`$ GeV2. If E01-001 verifies the polarization transfer measurements, we must still understand the source of the discrepancy in the cross section data in order to know the form factors to better than the few % level. While it seems likely that such errors would come from cross section measurements at ‘extreme’ kinematical conditions, and thus would thus have the largest impact on measurements of the $`ϵ`$-dependence, the actual consequence of the discrepancy between the new and old form factor measurements will not be clear until we understand why they disagree.

This work is supported in part by the U. S. Department of Energy, Nuclear Physics Division, under contract W-31-109-ENG-38.


text:
# What can rapid terrestrial biogenesis tell us about life in the universe?

## References

Caldeira, K. and Kasting, J.F. (1992) The life span of the biosphere revisited. Nature 360, 721-723.

Carter, B. (1983) The Anthropic Principle and its Implications for Biological Evolution. In Proc. R. Soc Discussion Meeting on the Constants of Physics, edited by W.H. McRea and M.J. Rees, R. Soc., London and in Philos. Trans. R. Soc. London, A310, 347-355.

Hanson, B. (2002) working paper on line at http://hanson.gmu.edu/vita.html Must Early Life be Easy? The Rhythm of Major Evolutionary Transitions.

Lineweaver, C.H. and Davis, T.M. (2002) Does the rapid appearance of life on Earth suggest that life is common in the Universe? Astrobiology in press, astro-ph/0205014.

Livio, M. (1999) How rare are extraterrestrial civilizations and when did they emerge? Astrophys. J. 511, 429-431.