text:
# Topological and energetic factors: what determines the structural details of the transition state ensemble and “on–route” intermediates for protein folding? An investigation for small globular proteins.

## I Introduction

Our understanding of the protein folding problem has been thoroughly changed by the new view that has emerged in the last decade. This new view, based on the energy landscape theory and funnel concept , describes folding as the progressive evolution of an ensemble of partially folded structures through which the protein moves on its way to the native structure. The existence of a deep energy funnel in natural proteins and the relatively simple connectivity between most conformational states which are structurally close makes this description possible even when only a few simple reaction coordinates that measure similarity to the native structure are used. The folding mechanism is controlled by both the shape of this free energy landscape and the roughness on it, which arises from the conflicts among interactions that stabilize the folded state and therefore can create non–native conformational traps .

The energetic roughness, however, is not the only limiting factor in determining a sequence’s foldability. Even if the energetic roughness could be completely removed, the folding landscape would not be completely smooth. Theoretical  and experimental  advances indicate that the final structure of the protein also plays a major role in determining a protein’s foldability. Some particular folding motifs may be intrinsically more designable than others. To address this difference in foldability which is not dependent on energetic frustration, we have introduced the concept of “topological frustration” .

Let us imagine an ideal situation for which the order of native contact formation during folding is not biased. In this “ideal” situation, there are an enormously large number of equivalent folding pathways, and an analysis of the transition state ensemble would show that for this ensemble nearly all parts of the protein have a similar probability of participation. The structure in the transition ensemble has been estimated by analogy with minimalist lattice models made to reproduce the global landscape features of small, fast folding proteins: similar Levinthal entropies, stabilities and energetic roughness as gauged by the glass transition temperatures. These models show a transition state ensemble about half way through the unfolded and folded states . In this ideal case, all the contacts in this transition ensemble would exist with the same probability.

Although the average amount of native formation in the transition ensemble is about 50%, the lattice simulations show that, even when the sequence is designed to have substantially reduced energetic frustration, there are variations in the amount of nativeness of specific contacts in the transition state ensemble . Real proteins display similar heterogeneity in contact formation. In systems with no energetic frustration and equal native interactions, these variations in the transition state ensemble are created solely by the folding motif and polymeric constraints that make certain contacts more geometrically accessible and stable than others. This variation in frequency that some contacts are made in the transition state ensemble generally reduces the entropy of the transition state and, when determined by the native motif, is a gauge of the amount of “topological frustration” in the system. Although this type of frustration can be modified by some design tricks , it cannot be completely eliminated: it reflects an intrinsic difficulty in folding to a particularly chosen shape. Minimalist models have shown how this heterogeneity leads to a transition ensemble that is a collection of diffuse nuclei which have various levels of native contact participation . The minimalist models calibrated to real proteins show similar overall levels of contact heterogeneity as in real proteins . This picture of a transition state composed of several diffuse nuclei has been confirmed by other lattice and off–lattice studies . In addition to selecting sequences which have low levels of energetic frustration, evolution appears to have selected for a particular set of folding motifs which have reduced levels of “topological frustration”, discarding other structures to which it is too difficult to fold .

Guided by theoretical folding studies on lattice, off–lattice, and all–atom simulations (see for instance , , , , , ) as well as recent experimental evidence , we suggest that real proteins, and especially small, fast folding (sub–millisecond), two–state like proteins, have sequences with a sufficiently reduced level of energetic frustration that the experimentally observed “structural polarization” of the transition state ensemble (viz. the variation in the amount of local native structure) is primarily determined by the topological constraints. That is, in well designed sequences, the variations are more determined by the type of native fold than by differences in sequence which leave the native fold relatively unchanged and the energetic frustration small.

The amount of native local structure in the transition state can be experimentally measured by using single and double point mutants as probes in the $`\mathrm{\Phi }`$ value technique . If the topology is a dominant source of heterogeneity in transition state structure, then the majority of evolved sequences which fold to the same motif would exhibit similar local structure in the transition state ensemble. We provide evidence in this paper that not only is this the case, that much of the transition state ensemble is determined by the final folded form, but, also for larger proteins that are not two–state folders, some “on–route” intermediates are determined by topological effects as well. Thus it appears that the dominance of topology in folding extends even into some larger, slower folding proteins with intermediates. This fact is consistent with some recent observations by Plaxco and collaborators that reveal a substantial correlation between the average sequence separation between contacting residues in the native structure and the folding rates for single domain proteins .

To ascertain the extent of topological control of the folding behavior, we create several simplified energetic models of small, globular proteins using potentials created to minimize energetic frustration. We show that these energetically unfrustrated models reproduce nearly all the known global features of the transition states of the real proteins on whose native structures they are based, including the structure of folding intermediates. We directly compare the structure of the transition state ensemble experimentally determined by $`\mathrm{\Phi }`$ value measurements with the numerically determined one. The simulated transition state ensemble is inferred from structures sampled in equilibrium around the free energy barrier between the folded and unfolded states. This free energy is computed as a function of a single reaction coordinate that measures the fraction of formed native contacts. The validity of this method has been demonstrated in references .

The organization of the paper is as follows: in section II we present in some detail the physical concepts underlying this work in the light of recent experimental results. In section III we present results for a sample of five small, globular proteins, and compared these results against the available experimental data. The off–lattice model used in our study is presented in the Appendix. In order to investigate the relevance of the topology, we chose a model which reproduces the topological features of a given real protein and eliminates most of the energetic frustration and variations in the strength of native residue–residue contacts. The predicted transition state for these proteins are in good agreement with experimental evidences, supporting our hypothesis of the major role played by topology.

## II Checking the Folding Mechanism by Analyzing the Transition State Ensemble

How do we know what the folding transition state ensemble looks like? Experimental analysis of folding transition state ensembles has been largely performed using the $`\mathrm{\Phi }`$–value analysis technique introduced by Fersht and co–workers . $`\mathrm{\Phi }`$ values measure the effects that a mutation at a given position along the chain has on the folding rate and stability:

$$\mathrm{\Phi }\frac{RT\mathrm{ln}(k_{mut}/k_{wt})}{\mathrm{\Delta }\mathrm{\Delta }G^0}$$
(1)
where $`k_{mut}`$ and $`k_{wt}`$ are the mutant and wild–type folding rate respectively, $`R`$ is the ideal gas constant, $`T`$ is the absolute temperature, and $`\mathrm{\Delta }\mathrm{\Delta }G^0`$ is the difference in the total stability between the mutant and wild–type proteins in $`kcal/mol`$.

Because the folding event of small fast folding proteins is well described as a diffusive process over a barrier determined by the free energy profile, the folding rate can be written as a Kramer’s–like equation

$$k=k_0\mathrm{exp}[\mathrm{\Delta }G^{}/RT]$$
(2)
where $`k_0`$ is a factor depending on the barrier shape and the configurational diffusion coefficient of the system. If $`k_0`$ is insensitive to small sequence changes, what appears to be true for reasonably unfrustrated sequences  the $`\mathrm{\Phi }`$ value is then seen to be a ratio of free energy changes of the folding barrier to stability:

$$\mathrm{\Phi }=\frac{\mathrm{\Delta }\mathrm{\Delta }G^{}}{\mathrm{\Delta }\mathrm{\Delta }G^0}$$
(3)
where $`\mathrm{\Delta }\mathrm{\Delta }G^{}`$ is given by

$$\mathrm{\Delta }\mathrm{\Delta }G^{}=\mathrm{\Delta }G_{mut}^{}\mathrm{\Delta }G_{wt}^{}=RT\mathrm{ln}\frac{k_{mut}}{k_{wt}}.$$
(4)
When this relationship is valid and the mutation can be considered a small perturbation, the $`\mathrm{\Phi }`$ value is a convenient measure of the fraction of native structure which is formed in the transition state ensemble around the site of the mutation. A $`\mathrm{\Phi }`$ value close to 1 means that the free energy change between the mutant and the wild type is almost the same in the transition state and native state, indicating that native contacts involving the mutated residue are already formed at the transition state. Inversely, a $`\mathrm{\Phi }`$ close to $`0`$ means that the free energy change is the same in the transition state and unfolded states, so the local environment of the residue is probably unfolded–like. A detailed analysis of the mutation is needed to determine exactly what contacts are disrupted under mutation. Ideally, mutations are made which eliminate small hydrophobic side–groups. Studies using $`\mathrm{\Phi }`$ values with multiple same–site mutations generally support the accuracy of $`\mathrm{\Phi }`$ value as a structural measurement of the transition ensemble , although sizable changes in the transition state structure have been induced in at least one protein through a single point mutation . In interpreting $`\mathrm{\Phi }`$ values, it is also important to remember that they only measure the relative change in structure, not the absolute amount of structure. This leads to the possibility that some mutants with low $`\mathrm{\Phi }`$ values may have nearly native local environments in the transition state, a possibility seen clearly in the experimental studies of Procarboxypeptidase A2 .

The validity of $`\mathrm{\Phi }`$ values as structural measurements clearly supports the Kramer’s–like description of the folding rate and the fact that the $`\mathrm{\Phi }`$ can be properly understood as a ratio of the free energy change of the transition ensemble over the change of the native ensemble (equation 3). This latter equation is very convenient as a starting point for computing $`\mathrm{\Phi }`$ values. In several recent simulation papers for lattice and off–lattice protein models, we have investigated this issue at length . All these studies concluded that as long as the systems present a weak or moderate level of energetic frustration (such as the Gō–like models in this work), $`\mathrm{\Phi }`$ values determined from changes in the free energy barrier,determined using a single simple reaction coordinate, yield quantitatively correct $`\mathrm{\Phi }`$ values. Therefore, all the calculations performed in this work were done utilizing eq. 3 — no actual kinetics was performed but only the appropriate sampling of the protein configurational space (see Appendix and refs. ,, ,, for example, for details). Technically, as long as the folding barriers are of a few $`k_BT`$ or more and the displacement of the barrier position along this reaction coordinate under mutation is sufficiently small, the $`\mathrm{\Phi }`$ values can be computed using free energy perturbation:

$$\mathrm{\Phi }=\frac{\mathrm{\Delta }\mathrm{\Delta }G^{TS}\mathrm{\Delta }\mathrm{\Delta }G^U}{\mathrm{\Delta }\mathrm{\Delta }G^F\mathrm{\Delta }\mathrm{\Delta }G^U}=\frac{\mathrm{ln}e^{\mathrm{\Delta }E/RT}_{TS}\mathrm{ln}e^{\mathrm{\Delta }E/RT}_U}{\mathrm{ln}e^{\mathrm{\Delta }E/RT}_F\mathrm{ln}e^{\mathrm{\Delta }E/RT}_U}.$$
(5)
We use equation 5 to compute $`\mathrm{\Phi }`$ values for our protein models using fixed transition, unfolded, and folded regions identified by the free energy profile viewed using a single order parameter: $`Q`$, the fraction of native contacts formed in a given conformation.

What experimental evidence exists as to the role of topology in determining the average structure in the folding transition state ensemble? The clearest evidence to date of the role of topology comes from comparisons of the transition state structure of two homologues of the SH3 domain (src SH3 and $`\alpha `$–spectrin SH3). These two homologues have only weak identity ($`30`$% identity with gaps), but $`\mathrm{\Phi }`$ values at corresponding sequence positions are highly correlated , supporting the degeneracy in the folding behavior for these two sequences. Furthermore, one of these sequences has a strained $`\mathrm{\Phi }`$$`\mathrm{\Psi }`$ conformation in the high $`\mathrm{\Phi }`$ region of the distal turn. The fact that this strain does not detectably lower the $`\mathrm{\Phi }`$ values in the local neighborhood , suggests that the sequence details and local stability are less important for determining how structured a region is in the transition state ensemble than its location in the final folded conformation. Other evidence indicates that these results may be more generally applicable than simply for SH3 or $`\beta `$–sheet proteins. Sequence conservation has been shown not to correlate with $`\mathrm{\Phi }`$ values , indicating that in general sequence changes at a given position in a protein weakly affect the $`\mathrm{\Phi }`$ value at that position.

Results for some small fast folding proteins (such as CI2 and the $`\lambda `$repressor) suggest that the transition state is an expanded version of the native state, with a certain degree of additional inhomogeneity over the structure  (similar to the theoretical predications for small $`\alpha `$–helical proteins ), while results for other proteins (as the $`\beta `$–sheet SH3 domain) show apparently larger structural heterogeneity in the transition state . This difference in the degree of “structural polarization” that is emerging between small $`\alpha `$–helix and $`\beta `$–sheet proteins suggests that the folding mechanism of a given protein is fundamentally tied to the type of secondary structural elements and their native arrangement. Current studies using $`\mathrm{\Phi }`$ value technique have been made of src SH3 , $`\alpha `$spectrin SH3 , CI2 , Barnase , Barstar , $`\lambda `$repressor , CheY , protein L , Procarboxypeptidase A2 , RNase H  and the tetrameric protein domain from tumor suppressor p53 .

In this paper, we analyze five proteins (SH3, CI2, Barnase, RNase H and CheY) that have been extensively studied experimentally and for which, therefore, details of their transition state ensemble are quite well known. We generate sequences (and potentials) for simulating these different globular proteins. These sequences have the native backbone folds of real experimentally studied globular proteins but sequence and potential interactions designed to drastically reduce the energetic frustration and heterogeneity in residue–residue interactions. By comparing the transition state structures of these unfrustrated models with the experimental studies of their real protein cousins, we quantify the effects of the native topology. If topology completely determines how folding occurs, then the model and real proteins should have identical folding behavior and $`\mathrm{\Phi }`$ values. If energetic frustration and heterogeneity are critical for determining the folding mechanism, then the variations in $`\mathrm{\Phi }`$ values with position should bear far reduced similarity to those in the real proteins on which the computer homologues are based.

Two of the five studied proteins are simple two–state like fast folding proteins (SH3 and CI2), while the other three (Barnase, RNase H and CheY) are known to fold through the formation of an intermediate state. We show not only that our simple models can reproduce most of the $`\mathrm{\Phi }`$ value structure, but also that models for Barnase, RNase H and CheY correctly reproduce the folding intermediates of these proteins, suggesting that many of the “on–route” intermediates are also largely determined by the type of native fold.

We represent the five globular proteins using a simplified Cα model with a Gō–like  Hamiltonian as detailed in the Appendix. This potential is in its details unlike that of real proteins, which have residue–residue interactions with many components (Coulomb interactions, hydrogen bonding, solvent mediated interactions, etc., etc.). The crucial features of this potential are its low level of energetic frustration, that characterizes good folders and a native conformation equal to the real protein. The ability of this model to reproduce features of the real transition state ensemble and real folding intermediates is a strong indication that the retention of the topology is enough to determine the global features of their folding mechanism. Using these models, we simulate the dynamics of a protein starting from its native structure, for several temperatures. To monitor the thermodynamics of the system, we group the configurations obtained during a simulation as a function of the reaction coordinate, $`Q`$, defined as the fraction of the native contacts formed in a conformation ($`Q=0`$ at the fully unfolded state and $`Q=1`$ at the folded state). The choice of $`Q`$ as order parameter for the folding is motivated by the fact that in a funnel–like energy landscape, a well designed sequence has the energy of its conformations reasonably correlated to degree of nativeness, and the parameter $`Q`$ is a good measure of the degree of similarity with the native structure. Our Gō–like potential is minimally frustrated for the chosen native structure, and the prediction of transition state ensemble structures and folding rates for these Gō–like systems has been shown to be quite accurate . From the free energy profile as a function of $`Q`$, it is easy to locate the unfolded, folded and transition state ensembles, as it is shown in next section. Since these models consider totally unfrustrated sequences, they may not reproduce the precise energetics of the real proteins, such as the value of the barrier heights and the stability of the intermediates, nonetheless they are able to determine the general structure of these ensembles.

In order to compare the folding process simulated using our model to the actual process for a given protein (as obtained from experimental $`\mathrm{\Phi }`$–values analysis), we need to choose a “mutation” protocol to compute $`\mathrm{\Phi }`$ values. Experimentally, the ideal mutation is typically one that removes a small hydrophobic side–group such as a methyl group that makes well–defined and identifiable residue–residue contacts in the native state. The $`\mathrm{\Phi }`$ value is then sensitive to this known contact. Our computational mutation is the removal of a single native bond, so our computer $`\mathrm{\Phi }`$ values are sensitive to the fractional formation of this bond $`Q_{ij}`$ between residues $`i`$ and $`j`$. We make these mutations because, as in most real mutations, they are sensitive to the formation of specific contacts, rather than being averages over interactions with many parts of the native structure. They mostly resemble the interaction $`\mathrm{\Phi }_{int}`$ value made by making double cycle mutants . $`\mathrm{\Phi }`$ values are computed from equation 5. In an ideal, perfectly smooth funnel–like energy landscape, all the $`\mathrm{\Phi }`$ values should be equal; in an energetically unfrustrated situation, $`\mathrm{\Phi }`$ values variations are due to the structure of the native conformation.

## III Determining the Transition State Ensemble of Small Globular Proteins

We have discussed the idea of “topological frustration” and its role in determining the structural heterogeneity of the transition state ensemble. We explore its role directly by creating protein models which drastically reduce the energetic frustration and energetic heterogeneity among residue–residue native interactions leaving the topology as the primary source of the residual frustration. Results obtained with these models, constructed using a Cα level of resolution with a Gō–like potential designed to fold to the native trace of chosen proteins, are then compared against the experimental data of those proteins. Five proteins with different folding motifs and different amounts of transition state heterogeneity (variation in $`\mathrm{\Phi }`$ values) and/or intermediates have been investigated.

We first analyze Chymotrypsin Inhibitor II (CI2), a mixed $`\alpha `$$`\beta `$ protein with a broad distribution of $`\mathrm{\Phi }`$ values (nearly uniform from 0 to 1). Then we present an analysis for the src SH3 domain, a largely $`\beta `$–sheet protein with a more polarized transition state structure (a substantial number of large $`\mathrm{\Phi }`$ values). We then apply the same technique to Barnase, Ribonuclease H (RNase H) and CheY, three other mixed $`\alpha `$$`\beta `$ proteins which fold via a folding intermediate. Although these proteins are not two–state folding proteins, we demonstrate that topology is also the dominant determinant of their folding behavior. We show that the topology plays a major role not only in the transition state ensemble, but it is also largely responsible for the existence and general structure of the folding intermediate. This result may be quite common for “on–route” folding intermediates and could provide a computational method for distinguishing between “on–pathway” and “off–pathway” structures which are inferred from experiments. To check the applicability of this method, the same approach presented in this paper has been extended elsewhere  to a pair of larger proteins (Dihydrofolate Reductase and Interleukin–1$`\beta `$). Even for these very large proteins we found that the overall structure of the transition state and intermediate ensembles experimentally observed can be obtained utilizing similar simplified models.

### A Analysis of two–state folders: CI2 and SH3

#### 1 CI2

The Chymotrypsin Inhibitor 2 (CI2) protein is a 64 residue protein, consisting of six $`\beta `$–sheets packed against an $`\alpha `$–helix to form a hydrophobic core. Experimental studies  have established that CI2 folding and unfolding can be modeled by simple two–state kinetics. The structure of the transition state for this protein has been extensively characterized by protein engineering , by free energy functional approaches , by a geometrical variational principle , and by all–atom molecular dynamics simulations . These studies have shown the transition state has roughly half of the native interactions formed in the transition state ensemble and a broad distribution of $`\mathrm{\Phi }`$ values in agreement with the general predictions of the energy landscape theory used with a law of corresponding states for small proteins . The broad distribution of $`\mathrm{\Phi }`$ values suggests that most hydrophobic contacts are represented at a level of about 50% in the transition state ensemble.

We constructed a Gō–like Cα model of CI2 as described in the Appendix. Several fixed temperature simulations were made and combined using the WHAM algorithm  to generate a specific heat versus temperature profile and a plot of the potential of mean force as a function of the folding order parameter $`Q`$ (see figure 1). From the free energy profile, we identified the dominant barrier, and used the thermal ensemble of states at its location to generate $`\mathrm{\Phi }`$ values from equation 5. The ranges of values of $`Q`$ used to determine each of these ensembles are shaded in figure 1. The mutations have been implemented by the removal of single attractive interactions (they are replaced with the same short ranged repulsive interactions used between residues without native interactions). The values computed via this method are shown in figure 2. Also shown in this figure is the fractional formation of individual native contacts in the transition state. The small difference between these two figures is primarily due to the fact that in the $`\mathrm{\Phi }`$ calculations the native contact formation at the folded and unfolded states are also taken into account. Because of the higher concentration of contacts between residues near–by in sequence and the local conformational preferences, the unfolded state shows a high level of local structure. The inaccurate representation of local contacts in the unfolded state makes the short range $`\mathrm{\Phi }`$ values less reliable as transition structure estimates than long range $`\mathrm{\Phi }`$ values.

From the calculations, we detect three significant regions of large $`\mathrm{\Phi }`$ values: the $`\alpha `$–helix, the mini–core defined by strands 3 and 4 and their connecting loop, and between the C–terminus of strand 4 and the N–terminus of strand 5. These regions generally have $`\mathrm{\Phi }`$ values in excess of $`0.6`$. Slightly smaller values of about $`0.5`$ exist for the short range contacts between the N–terminal of strand 3 and the C–terminal of the $`\alpha `$–helix and for contacts between strand 3 and strand 4. All other regions lack a consistent set of large $`\mathrm{\Phi }`$ values. Despite the large number of native contacts between strands 1 and 2 and the $`\alpha `$–helix and between strands 5 and 6 and the $`\alpha `$–helix, only low $`\mathrm{\Phi }`$ values are observed in this region (nearly all below $`0.2`$ in value). A comparison between these data and the exhaustive analysis of Fersht and colleagues  shows excellent overall agreement. They have found that “$`\beta `$–strands 1, 5 and 6 … are not structured in the transition state….”. Strand 2 also shows a highly reduced amount of structure. Furthermore, “the central residues of $`\beta `$–strands 3 and 4 interact with the $`\alpha `$–helix to form the major hydrophobic core of CI2.” The hydrophobic mini–core in this region (defined as the cluster formed by side–chains of residues 32, 38, and 50) is detected by single mutant and double mutant $`\mathrm{\Phi }`$–values  to be at least 30% formed in transition ensemble. Similarly, they found the $`\alpha `$–helix, particularly the N–capping region, to be highly ordered.

In summary, we see a quite good overall agreement except for a discrepancy in the short range interactions in the loop region between strands 4 and 5. This protein shows generally higher $`\mathrm{\Phi }`$ values between interactions which are more local in sequence and lower $`\mathrm{\Phi }`$ values between interactions which are distant in sequence. The results are thus consistent with the picture of the transition state as a collection of non–specific and somewhat diffuse nuclei . This overall low level of frustration suggests a low level of “topological frustration” in this model as well and a particularly designable motif.

#### 2 src SH3 domain

Src SH3 is the 57 residue fragment of Tyrosine–Protein Kinase that stretches from T84 to S140. It has five $`\beta `$ strands (and a short 3–10 helix) in an anti–parallel arrangement, forming a partial $`\beta `$ sandwich. Experimental measurements have shown that the SH3 domain folds using a rapid, apparently two–state mechanism. A $`\mathrm{\Phi }`$ value analysis  reveals that the distal loop hairpin and diverging turn regions are both highly structured and docked together at the transition state; the hydrophobic interactions between the base of the hairpin and the strand following the diverging turn are partially formed, while other regions of src SH3 appear only weakly ordered in the transition state ensemble. The overall representation of the transition state structure of src SH3 —having the distal loop and diverging turn largely formed and other regions weakly formed— agrees with studies of $`\alpha `$spectrin SH3 which has a similar backbone structure but a dissimilar sequence (approx 30% identity with gaps). This observed similarity along with evidence of a strained backbone conformation in the distal loop of the $`\alpha `$spectrin SH3  supports the concept of “topological” dominance in folding .

Fig. 3 shows the folding behavior as obtained from our dynamics simulations of the Gō–like analogous of the src SH3. The free energy barrier defining the transition state location is evident in the figure. As before, we have computed $`\mathrm{\Phi }`$ values from equation 5 by mutating (removing) every native residue–residue attractive contact. The results of this calculation are shown in figure 4. In addition to $`\mathrm{\Phi }`$ values, the contact formation probability at the transition state ensemble have been calculated. Our previous caveats concerning $`\mathrm{\Phi }`$ values for local interactions still apply. We observe the highest collection of off–diagonal (long range) $`\mathrm{\Phi }`$ values is in the diverging turn —distal loop interaction exactly as seen from the experimental $`\mathrm{\Phi }`$ value measurements. We see very low values in the RT loop region, in accord with the two mutants in this loop. We also see medium to high values between the two $`\beta `$ strands which are connected by the distal loop. The transition state structure of the SH3 presents a substantially larger degree of structural polarization than CI2, where the $`\mathrm{\Phi }`$ values are much more uniform. This suggests that SH3 has a backbone conformation which is intrinsically more difficult to fold, i.e., there is a greater level of “topological frustration” in this structure. Nevertheless the transition state composition is well reproduced for both the two proteins.

### B Analysis of three proteins which fold throughout the formation of an intermediate state: Barnase, RNase H and CheY

Barnase, RNase H and CheY are three small $`\alpha `$$`\beta `$ proteins (although larger than the previous two proteins): Barnase is a 110 residue protein, composed by three $`\alpha `$–helices (located in the first 45 residues) followed by five $`\beta `$–strands; RNase H consists of 155 residues which arrange themselves in five $`\alpha `$–helices and five $`\beta `$–strands; CheY is a 129 residues, classic $`\alpha \beta `$ parallel fold in which five $`\beta `$–strands are surrounded by five $`\alpha `$–helices. Experimental results show that these three proteins do not fold by following a simple two–state kinetics directly from the unfolded state to the native structure, but fold through the formation of a metastable intermediate which interconverts into the native state. This brings up an interesting question: is topology alone able to determine the presence of an intermediate in the folding process? In Figs. 5, 7 and 9 we show evidence for the first time that such intermediates can be created solely from a Gō–like minimalist model which preserves the native topology. The presence of this intermediate during these protein’s folding events is a requirement of the native protein motifs. The free energy changes upon mutations of a wild–type three–state protein are experimentally measured both for the intermediate and the transition state, to define two different sets of $`\mathrm{\Phi }`$–values for the protein:

$$\begin{array}{c}\mathrm{\Phi }_I=\frac{\mathrm{\Delta }\mathrm{\Delta }G_I\mathrm{\Delta }\mathrm{\Delta }G_U}{\mathrm{\Delta }\mathrm{\Delta }G_F\mathrm{\Delta }\mathrm{\Delta }G_U}\\ \mathrm{\Phi }_{TS}=\frac{\mathrm{\Delta }\mathrm{\Delta }G_{TS}\mathrm{\Delta }\mathrm{\Delta }G_U}{\mathrm{\Delta }\mathrm{\Delta }G_F\mathrm{\Delta }\mathrm{\Delta }G_U}\end{array}$$
(6)
where $`\mathrm{\Phi }_I`$ provides information about the structural composition of the intermediate state (I), and $`\mathrm{\Phi }_{TS}`$ of the transition state (TS). Following we discuss in some details the results for the three proteins. Since, as for the first two proteins, the $`\mathrm{\Phi }`$–values and the native contact probabilities provide somewhat similar information, for simplicity, we show only the results obtained for the native contact probabilities (for safety we have checked the $`\mathrm{\Phi }`$–values and determined that similar information is recovered).

#### 1 Barnase

The analysis of experimentally obtained $`\mathrm{\Phi }`$ values  for the Barnase shows that some relevant regions of the structure are fully unfolded in the intermediate while other regions are fully folded.

Fig. 6 shows the intermediate and the transition state structure obtained from the Gō–like model. The intermediate shows substantial structural heterogeneity: there are very high probability values for interactions within the $`\beta `$–sheet region and its included loops, and very low values for interactions within the $`\alpha `$–helices and their loops and between the $`\alpha `$–helical and $`\beta `$–sheet regions. Some local short range helical interactions are formed. The transition state ensemble structure shows the same structure as the intermediate with the addition of strong interactions within helices 2 and 3; between helix 2, helix 3, the first $`\beta `$–strand, and the intervening loops; and between the second $`\beta `$–strand and the second helix.

Comparing these simulation results with extensive mutagenesis studies of reference , we observe a good qualitative agreement. The $`\beta `$–sheet region is highly structured in the intermediate as it is the core region 3 (consisting of the packing of loop 3, that joins strands 1 and 2, and of loop 5, that joins strands 4 and 5, with the other side of the $`\beta `$–sheet). In agreement with experiments, the earliest formed part of the protein appears to be the $`\beta `$–sheet region. Also the core region 2 (formed by the hydrophobic residues from helix 2, helix 3, the first strand, and the first two loops) is found to be only weakly formed in the intermediate and the transition state.

There are two minor discrepancies between the Barnase model and the experimental data. First, we slightly overestimate the formation of core region 2 in the transition state ensemble. Second, we underestimate the amount of structure in core region 1 (formed by the packing of the first helix against a side of the $`\beta `$–sheet) in both intermediate and transition ensemble. In particular, we under–represent the interaction between helix 1 and the $`\beta `$–sheet region. The experimentally observed early packing of helix 1 against the rest of the structure is not reproduced by our model. Clearly there are some important energetic factors which have been neglected by the simple model. These may be inferred from the Barnase crystal structure. For example, one can see that helix 1 is largely solvent exposed, with interactions between it and the remainder of the protein formed by only five of the eleven helix residues. 83 % of the interactions reside on the hydrophobic residues PHE7, ALA11, LEU14 and GLN15, and the 17 % of the interactions are formed by the charged residues ASP8 and ASP12, while the solvent exposed part of the helix is composed of polar residues. Large stabilizing interactions other than tertiary (most hydrophobic) interactions are neglected in the model, being probably responsible for the failure in predicting the formation of the structural parts involving helix 1. In this structural detail, it appears that the topological factors are not the leading determinant of the folding behavior.

#### 2 Ribonuclease H

Kinetic studies of the wild–type RNase H have shown that an intermediate state is populated in the folding process, and the structure of this intermediate has been extensively investigated by circular dichroism, fluorescence and hydrogen exchange methods  and by protein engineering . Fig. 7 shows that, consistently with the experimental evidences, we find an intermediate state in the folding process of the RNase H model. Experimental results indicate that the most stable region of the protein intermediate involves the $`\alpha `$–helix 1, the strand 4, the $`\alpha `$–helix 4 and the $`\alpha `$–helix 2. Hydrogen exchange experiments have shown that the $`\alpha `$–helix 1 is the region of the protein most protected from exchange, suggesting that most of the interactions involving the $`\alpha `$–helix 1 are already significantly formed at the intermediate state of the folding process. The helix 4 and the $`\beta `$–strand 4 are the next most protected regions, while the $`\alpha `$–helix 5 has low to moderate level of protection. After the completion of the this intermediate structure, the rate–limiting transition state involves the ordering of the $`\beta `$–sheet and the $`\alpha `$–helix 5. The packing of helix 5 across the sheet is found to be the latest folding event.

The results of the model for RNase H show a good agreement with the experimental evidences. As shown in Fig. 8, we find that the formation of contacts involving the helix 1 is the earliest event in the folding process. Contacts arising from the $`\alpha `$–helix 4 and the $`\beta `$–strand 4 are then formed at the intermediate state and consolidated at the transition state. In agreement with the experimental results, we find that, at the transition state, interactions between the $`\alpha `$–helix 1, the strand 4 and the rest of the protein are mostly formed; the $`\alpha `$–helix 4 is also well structured and interactions between the helix 4 and the other parts of the protein are partly formed. Interactions among the strands are almost all formed, but the sheet is not yet docked to the helix 5.

#### 3 CheY

Utilizing protein engineering , the transition state of CheY has been characterized and it can be described as a combination of two subdomains: the first half of the protein (subdomain 1), comprising the $`\alpha `$–helices 1 and 2 and the $`\beta `$–strands 1–3, is substantially folded whereas the second half (subdomain 2) is completely disorganized. The helix 1 seems to play the role of a nucleation site around which subdomain 1 begins to form. Moreover, an intermediate has been detected at the early stage of the folding process where all the five $`\alpha `$–helices are rather structured. The last two helices, however, are very unstructured in the later occurring transition state. From this result it has been suggested that a misfolded species is visited at the beginning of the folding process.

Our simple model detects two possible intermediates for this protein, one of them is an “on–route” intermediate that is short–living and occurs just before the transition state ensemble ($`Q`$ around 0.6 in Fig. 9). Surprisingly, the unfrustrated model is also able to detected a “misfolded” trap in the folding of CheY. Since non–native interactions are not allowed in the model, this trap is a long–living partially folded state created by the topological constrains. There is no direct connection between this trap state and the fully folded state. The structure of this trap is shown in Fig. 10 and it agrees with the experimental observation of all helices well structured. Differently from the previously discussed proteins, the model of CheY seems to have a tendency to first form a “wrong” part of the protein and, when this happens, a partial unfolding must occur before the folding can be completed.

Finally, analyzing the transition state structure, we find a good agreement with the experimental data. As shown in Fig. 10, the first part of the protein (subdomain 1) is almost fully folded at the transition state ensemble, while subdomain 2 is completely unfolded.

## IV Conclusions

Recent theoretical studies and experimental results suggest that the folding mechanism for small fast folding proteins is strongly determined by the native state topology. The amount of energetic frustration, arising from the residual conflict among the amino–acid interactions, appears largely reduced for these proteins so that topological constraints are important factors in governing the folding process. Towards exploring this topological influence in real proteins, we analyzed the folding process of the Gō–like analogous of five real proteins. Since we have used Gō–like potentials, the energetic frustration is effectively removed from the system, while the native fold topology is taken into account. It is important to highlight that the results from such studies exhibit the overall topological features of the folding mechanism, although we do not expect the precise energetic values for barrier heights and intermediate state stabilities. For example, real proteins are not necessarily totally unfrustrated and they have only to minimize energetic frustration to a sufficiently reduced level in order to be good folders. Also, as long as energetic frustration is small enough, creating some heterogeneity at the native interactions may help to reduce topological frustration , and that will energetically favor some contacts over others.

The effective use of a small number of global order parameters as reaction coordinates, in interpreting real data or studying more detailed protein folding model, depends critically on the degree of frustration present in real proteins . Since our results show that general structural features of the transition state ensemble in real proteins, at least for this class of fast folding proteins, is reproducible by using a substantially unfrustrated potential, several different global order parameters should work to explain the folding mechanism. For this reason, it should not be a surprise the fact that, utilizing energy landscape ideas and the funnel concept, some very simple models with approximate order parameters determined by single or few sequence approximation  have been successful in predicting qualitative features of the transition state ensemble.

Again, we have compared in details the structure of the transition state ensemble of the five proteins resulting from our simulations with experimental data. The agreement between our results and the experimental data supports the idea that energetic frustration is indeed sufficiently reduced and the protein folding mechanism, at least for small globular proteins, is strongly dependent on topological effects. The structure of the transition state ensemble of the CI2 presents a broad distribution of $`\mathrm{\Phi }`$ values —i.e. a reduced degree of structural polarization— in agreement with predictions based on the energy landscape theory (see , ). On the other hand, the structure of the SH3 transition state ensemble shows a higher degree of polarization. Nevertheless, by using our simplified Gō–like model, we have reproduced the transition state composition for both proteins, demonstrating that topology is largely responsible for the observed experimental differences. The last three proteins we have analyzed, (Barnase, RNase H and CheY) are known to fold through a three–state kinetics, involving the formation of an intermediate structure. Our Gō–like model of these proteins also fold with a three–state kinetics with intermediates that are analogous to the ones detected experimentally. This fact suggests that topology is also a dominant factor in determining the “on–route” intermediates.

## V acknowledgments

We thank Viara Grantcharova and David Baker for informations about the SH3 structure. We also thank Vladimir Sobolev for the CSU software. We are indebted to Angel García, Peter Wolynes, Steve Plotkin, Jorge Chahine, Joan Shea, Margaret Cheung, Charlie Brooks, Amos Maritan and Jayanth Banavar for helpful discussions. One of us (C.C.) expresses her gratitude to Giovanni Fossati for his suggestions and for carefully reading the manuscript, and to the Center for Astrophysics & Space Sciences of UCSD for the usage of graphics facilities and computer time. This work has been supported by the NSF (Grant #96–03839), by the La Jolla Interfaces in Science program (sponsored by the Burroughs Wellcome Fund) and by the Molecular Biophysics training grant program (NIH T32 GN08326).

## Model and Method

In order to investigate how the native state topology affects the folding of a given protein we follow the dynamics of the protein by using a Gō–like Hamiltonian  to describe the energy of the protein in a given configuration. A Gō–like Hamiltonian takes into account only native interactions, and each of these interactions enters in the energy balance with the same weight. It means that the system gains energy as much as any amino acid pair involved in a native contact is close to its native configuration, no matter how strong the actual interaction is in the real protein. Residues in a given protein are represented as single beads centered in their C–$`\alpha `$ positions. Adjacent beads are strung together into a polymer chain by mean of bond and angle interactions, while the geometry of the native state is encoded in the dihedral angle potential and a non–local potential. The energy of a configuration $`\mathrm{\Gamma }`$ of a protein having the configuration $`\mathrm{\Gamma }_0`$ as its native state is thus given by the expression:

$$\begin{array}{c}E(\mathrm{\Gamma },\mathrm{\Gamma }_0)=_{bonds}K_r\left(rr_0\right)^2+_{angles}K_\theta \left(\theta \theta _0\right)^2+\\ _{dihedral}K_\varphi ^{(n)}\left[1+\mathrm{cos}\left(n\times (\varphi \varphi _0)\right)\right]+\\ _{i<j3}\{ϵ(i,j)[5\left(\frac{\sigma _{ij}}{r_{ij}}\right)^{12}6\left(\frac{\sigma _{ij}}{r_{ij}}\right)^{10}]+ϵ_2(i,j)\left(\frac{\sigma _{ij}}{r_{ij}}\right)^{12}\}.\end{array}$$
(7)
In the previous expression $`r`$ and $`r_0`$ represent the distances between two subsequent residues at, respectively, the configuration $`\mathrm{\Gamma }`$ and the native state $`\mathrm{\Gamma }_0`$. Analogously, $`\theta `$ ($`\theta _0`$) and $`\varphi `$ ($`\varphi _0`$) represent the angles formed by three subsequent residues and the dihedral angle defined by four subsequent residues along the chain at the configuration $`\mathrm{\Gamma }`$ ($`\mathrm{\Gamma }_0`$). The dihedral potential consists of a sum of two terms for every four adjacent $`C_\alpha `$ atoms, one with period $`n=1`$ and one with $`n=3`$. The last term in Eq. (7) contains the non–local native interactions and a short range repulsive term for non–native pairs (i.e. $`ϵ(i,j)=constant>0`$ and $`ϵ_2(i,j)=0`$ if $`i`$$`j`$ is a native pair while $`ϵ(i,j)=0`$ and $`ϵ_2(i,j)=constant>0`$ if $`i`$$`j`$ is a non–native pair). The parameter $`\sigma _{ij}`$ is taken equal to $`i`$$`j`$ distance at the native state for native interactions, while $`\sigma _{ij}=4`$ Å for non–native (i.e. repulsive) interactions. Parameters $`K_r`$, $`K_\theta `$, $`K_\varphi `$, $`ϵ`$ weight the relative strength of each kind of interaction entering in the energy and they are taken to be $`K_r=100ϵ`$, $`K_\theta =20ϵ`$, $`K_\varphi ^{(1)}=ϵ`$ and $`K_\varphi ^{(3)}=0.5ϵ`$. With this choice of the parameters we found that the stabilizing energy residing in the tertiary contacts is approximately twice the stabilizing energy residing in the torsional degrees of freedom. This balance among the energy terms is optimal to study the folding of our Gō–like protein models. The native contact map of a protein is derived with the CSU software based upon the approach developed in ref. . Native contacts between pairs of residues $`(i,j)`$ with $`ji+3`$ are discarded from the native map as any three and four subsequent residues are already interacting in the angle and dihedral terms. A contact between two residues $`(i,j)`$ is considered formed if the distance between the $`C_\alpha `$’s is shorter than $`\gamma `$ times their native distance $`\sigma _{ij}`$. It has been shown  that the results are not strongly dependent on the choice made for the cut–off distance $`\gamma `$. In this work we used $`\gamma =1.2`$. We have used Molecular Dynamics (entailing the numerical integration of Newton’s laws of motion) for simulating the kinetics of the protein models. We employed the simulation package AMBER (Version 4.1)  at constant temperature, i.e. using Berendsen algorithm for coupling the system to an external bath . Both temperature and energy are measured in units of the folding temperature $`T_f`$ in the simulations.

For each protein model, several constant temperature simulations were made and combined using the WHAM algorithm  to generate a specific heat profile versus temperature and a free energy $`F(Q)`$ as a function of the folding reaction coordinate Q. This algorithm is based on the fact that the logarithm of probability distribution $`P(Q)`$ of the values taken by a certain variable Q (e.g. the order parameter) at fixed temperature T may serve as an estimate for the the free energy profile $`F(Q)`$ at that temperature. In fact, the probability to have a certain value $`Q_1`$ for the variable Q, at temperature $`T=1/\beta `$, in the canonical ensemble is given by:

$$P_\beta (Q_1)=\frac{W(Q_1)e^{\beta E(Q_1)}}{Z_\beta }$$
(8)
where $`W(Q)`$ is the density of configurations at a point $`Q`$ in the configurational space, $`Z_\beta `$ is the canonical partition function at temperature $`T=1/\beta `$ and $`E(Q)`$ is the energy of the system at the value Q of the reaction coordinate***Since our model is almost energetically unfrustrated, the energy fluctuations for a set of configurations with fixed Q are strongly reduced such that the energy in a given configuration could be considered as a function of $`Q`$.. Since the free energy $`F`$ is

$$F(Q)=E(Q)TS(Q)$$
(9)
and the entropy S(Q) is related to the configurational density $`W(Q)`$

$$W(Q)e^{S(Q)/k}$$
(10)
where $`k`$ is the Boltzmann constant, it follows that

$$\frac{P_\beta (Q_1)}{P_\beta (Q_2)}=\frac{e^{\beta F(Q_1)}}{e^{\beta F(Q_2)}}$$
(11)
and free energy differences can be computed by

$$\beta (F(Q_1)F(Q_2))=\mathrm{log}\frac{P_\beta (Q_1)}{P_\beta (Q_2)}.$$
(12)
By using the procedure of refs. , data from a finite set of simulations can be used to obtain complete thermodynamic information over a large parameter region.

Probability distributions are obtained by sampling the configurational space during Molecular Dynamics simulations.

For the smaller proteins (CI2 and SH3) we have determined the errors on the estimates of the transition temperature and contact probabilities (or $`\mathrm{\Phi }`$ values). This has been accomplished by computing these quantities from several (more than 10) uncorrelated sets of simulations. We found that the standard deviation for each single contact probability is 0.06 for CI2 and 0.05 for SH3, while the transition temperature is determined in both cases with an uncertainty smaller than 0.5%. These errors are obtained using about 200 uncorrelated conformations in the transition state ensemble. Since Barnase, RNase H and CheY have twice to three times the number of tertiary contacts of SH3 and CI2, in order to have appropriate statistics, we have sampled about 500 uncorrelated conformations (thermally weighted) for every transition state ensemble or intermediate.

### Captions to the figures

Fig. 1. (a) Free energy $`F(Q)`$ as a function of the reaction coordinate $`Q`$ around the folding temperature for the model of CI2. Free energies are measured in units of $`k_BT_f`$, where $`T_f`$ is the folding temperature. The unfolded, folded and transition state regions are shown in the light blue shaded areas. (b) A typical sample simulation at a temperature around the folding temperature. The reaction coordinate Q as a function of time (measured in arbitrary unit of molecular dynamics steps) is shown. The two–state behaviour is apparent from the data. The unfolded and folded states are equally populated at the folding temperature. (c) Heat capacity as a function of the temperature (units of folding temperature).

Fig. 2. The results for the transition state structure from the simulations for CI2. The probability of native contact formation at the transition state (left panel), and bond $`\mathrm{\Phi }`$–values (right panel) are shown. Different colors indicate different values from 0 to 1, as quantified by the color scale. The $`\alpha `$–helix, the interactions between the strands 4 and 5, and the minicore (i.e. interactions between residues 32,38 and 50) are the parts of the structure formed with the highest probability, although they are not fully formed. Overall, the transition state ensemble appears as an expanded version of the native state where most contacts have a similar probability of participation, but some interactions are less like to occur. These results agree with the transition state structure experimentally obtained.

Fig. 3. (a) Free energy $`F(Q)`$ as a function of the reaction coordinate $`Q`$ for a set of temperatures around the folding temperature. Free energies are measured in units of $`k_BT_f`$. The choices for the unfolded, folded and transition state regions are marked as shaded regions. (b) The reaction coordinate Q as a function of time (unit of molecular dynamics steps), from a typical sample simulation around the folding temperature. As in Fig. 1, the two–state behaviour is apparent. At the transition temperature the model protein has equal probability to be found in the unfolded or in the folded state. (c) Heat capacity as a function of the temperature, in units of folding temperature.

Fig. 4. The transition state structure as obtained from the simulations for SH3. Panel in the left represents the probability for a native contact to be formed at the transition state, while the panel in the right shows the results for bond $`\mathrm{\Phi }`$–values. Different colors indicate different values from 0 to 1, as quantified by the color scale. Diverging turn and distal loop are marked on the contact map. The interactions within and between these two parts of the protein chain appear to be formed with high probability. The interactions between the two strands joined by the distal loop are partially formed, while the contacts involving the first 20 residues do not contribute to the transition state structure. This description of the transition state is in agreement with experimental results.

Fig. 5. (a) Free energy $`F(Q)`$ of Barnase protein as a function of the reaction coordinate $`Q`$ around the folding temperature. Free energies are measured in units of $`k_BT_f`$. The unfolded, folded and intermediate state regions are marked in green, while the top of the two barriers are marked in light blue. The local minimum in the free energy profile between the unfolded and folded minima locates the folding intermediate state. The presence of a folding intermediate state is also evident from panel (b), where the order parameter $`Q`$ is plotted as a function of time for a typical molecular dynamics simulation around the folding temperature. In the interval $`Q(0.40.5)`$, the same state (i.e. with the same average structure) is visited both from the unfolded and folded structures.

Fig. 6. The probability of native contact formation for the intermediate (left panel) and transition state (right panel) structures as obtained from our simulations of Barnase. Different colors indicate different values from 0 to 1, as quantified by the color scale. The earliest formed part of the protein appears to be the $`\beta `$–sheet region, in agreement with experimental results. The core 3 (formed by loops 3 and 5 to the $`\beta `$–sheet) is formed at the intermediate and transition state, while core 1 (the packing of the helix 1 against the $`\beta `$–sheet) and the core 2 (the interactions between the hydrophobic residues from the helices 2 and 3, the strand 1, and the first two loops) start to form only after the transition state. The formation of the $`\alpha `$–helix 1 occurs as a late event of the folding from our simulations, while from experimental results it seems to be already formed at the intermediate and transition state. The early formation of the $`\alpha `$–helix is most probably due to energetic factors rather then from topology requirements (and then beyond the prediction possibility of this model), as detailed in the text.

Fig. 7. (a) Free energy $`F(Q)`$ of the model of RNase H as a function of the reaction coordinate $`Q`$ around the folding temperature. Free energies are measured in units of $`k_BT_f`$. The regions corresponding to the unfolded, folded and intermediate state are marked in green, while the top of the two barriers are marked in light blue. A folding intermediate is detected as a local minimum in the free energy between the unfolded and folded minima. In panel (b) the fraction of native contacts formed, $`Q`$, is plotted versus the simulation time for a sample of our simulations (at a temperature $`T=0.99T_f`$) where the transition from unfolded to folded state is observed. The local minimum of panel (a) corresponds to a transiently populated intermediate (located at $`Q`$ around 0.4) that later evolves to the fully folded state.

Fig. 8. The probability of native contact formation at the intermediate (left panel) and transition state (right panel) structure, as observed for the RNase H model. Different colors indicate different values from 0 to 1, as quantified by the color scale. In agreement with experimental results, we found that interactions involving the $`\alpha `$–helix 1 are the first formed in the folding process. Contacts between the $`\alpha `$–helix 1 and the strand 4 are highly probably formed at the intermediate. Also the $`\alpha `$–helix 4 is well structured and the $`\beta `$–sheet is partly formed. These interactions strengthen at the transition state where also the $`\beta `$–sheet is almost completely formed, while the packing of helix 5 across the sheet is not yet accomplished.

Fig. 9. (a) Free energy $`F(Q)`$ profile for the model of CheY plotted as a function of the reaction coordinate $`Q`$ for a set of temperatures around the folding temperature. Free energies are measured in units of $`k_BT_f`$. Differently from the corresponding figures of Barnase (Fig. 5) and RNase H (Fig. 7), two different structures are populated between the folded and unfolded states. In addition to the “on–route” intermediate state (marked in green as the regions corresponding to the folded and unfolded states), a “misfolded” intermediate structure (marked in brown at Q around 0.4) is transiently visited from the unfolded state. The top of the two barriers are marked in light blue. In agreement with experimental results, we found that in this “misfolded” structure, all the five $`\alpha `$–helices are rather structured while, in the later occurring “on–route” intermediate and transition state ensemble, the helices 4–5 are completely unstructured (see fig. 10). Panel (b) shows a typical sample of the simulation around the folding temperature, in a region where the folding occurs. The first transiently populated intermediate state corresponds to a structure where all the helices are formed. Before to proceed to the folded state, a partial unfolding occurs.

Fig. 10. The probability of the native CheY contacts to be formed in the “misfolded” intermediate (left panel) and transition state (right panel) for the model protein. Different colors indicate different values from 0 to 1, as quantified by the color scale. In agreement with experimental data, all the helices are mostly formed in the transiently populated “misfolded” structure, while helices 4 and 5 are rather unstructured at the transition state. The two subdomains experimentally detected in the CheY transition state  are evident in the figure: the first part of the protein (all interactions arising from the $`\alpha `$–helices 1–2 and the $`\beta `$–strands 1–3) is folded, while the second part (interactions among the $`\alpha `$–helices 4–5 and the $`\beta `$–strands 4–5) is completely unfolded. The helix 3 is structured but the interactions between the helix 3 and the rest of the protein are not completely formed.


format: text

text:
# Measurements of Light Nuclei Production in 11.5 A GeV/c Au + Pb Heavy-Ion Collisions

## I Introduction

Relativistic heavy ion collisions are believed to reach energy densities an order of magnitude greater than that of normal nuclear matter. These collisions allow the examination of the strong interaction in a novel environment as well as providing a possible doorway to new states of matter. In order to understand the dynamics of the collision system, one must use the only available tools-the species and momenta of the particles which exit the collision region. The use of emitted hadrons to probe the collision system is complicated because these hadrons rescatter many times as they traverse the collision region, and consequently lose some of their direct information about the earlier stages of the time evolution. However, the final space-time extent of the system at freeze-out (the time when strong interactions cease) and position-momentum correlations of the emitted particles contain much information about the entire time evolution. In principle, these carry information about the equation-of-state of the early collision region.

In order to extract information about both the momentum and position distributions of the source at freeze-out, it is necessary to measure multi-particle correlations. One widely used technique is Hanbury-Brown-Twiss interferometry  for which the correlations between particles are due to quantum statistics. Another method is through the measured yields of light nuclei, which are formed by the coalescence of individual nucleons.

Because of the violence of heavy ion collisions, it is highly improbable for a nuclear cluster near center-of-mass rapidity $`y_{cm}`$=1.6 in a collision at these energies to be a fragment of the beam or target nucleus  . This would involve a cluster suffering a momentum loss of several GeV/c per nucleon that does not destroy the cluster, which is typically bound by only a few MeV per nucleon. These nuclei then are formed by coalescence and so represent correlations of several nucleons. As the mass of measured nuclei increases, of course, so does the number of particles involved in the correlation, and so does the sensitivity to features of the freeze-out distribution.

In part due to the fragility of these states, the observed light nuclei are believed to be formed only near freeze-out of the collision system, at which time the mean free path of a bound cluster is long enough for it to escape without further collision. It is this notion that gives rise to a class of models of light nuclei production, the coalescence models (for example, ). In general, these models assume a phase space distribution of nucleons at freeze-out and impose some coalescence conditions on the freeze-out positions and momenta of the nucleons in order to calculate the yields of nuclei. These models differ both in their assumptions about the phase space profiles and coalescence conditions. Their differences are often characterized by their predictions of the invariant coalescence, or $`B_A`$, parameters which are defined as

$$B_A\frac{(E\frac{d^3N_A}{dP^3})}{(E\frac{d^3N_{neutron}}{dp^3})^N(E\frac{d^3N_{proton}}{dp^3})^Z}$$
(1)
where a nucleus with baryon number $`A`$ and momentum $`P=pA`$ is formed out of $`Z`$ protons and $`N`$ neutrons.

In the simplest momentum space coalescence models, coalescence is assumed to take place between any nucleons with a small enough momentum difference. Early experimental results at the Bevalac with beam energies of $`500A`$MeV and high energy proton induced reactions revealed $`B_A`$ parameters which were approximately constant for these different collision systems ( contain useful summaries). The assumption of such simple models is that if the collision spatial volume is similar to size of the cluster, all nucleons whose momentum difference is less than a fixed value will fuse to form a nuclear cluster. Thus the experimentally observed constant $`B_A`$ values may indicate collision volumes in these systems that are not substantially larger than the RMS size of a deuteron.

In more advanced models, assuming a quantum mechanical sudden approximation  and using a density matrix formalism, accounting for both the positions and momenta of the nucleons , the $`B_A`$ parameters take on a relationship with the source volume $`V`$ of $`B_A(1/V)^{(A1)}`$. Heavy-ion experimental results at higher energies at the AGS ($`10AGeV`$) and CERN ($`160AGeV`$) revealed $`B_A`$ values that decreased with beam energy. This observation was understood to be a sign of significant expansion in the collision volume before freeze-out. This larger source volume creates a situation where some nucleons with small relative momentum will have too large a spatial separation to coalesce, thus reducing $`B_A`$.

The density matrix formalism  assumes that although the collision volume can expand significantly, there is no correlation between the momentum and the position of a given particle. This assumption then leads to a prediction of no kinematic dependence of the $`B_A`$ parameter. However, there is a great deal of evidence that collective motion is present leading to expansion of the collision volume and significant position-momentum correlations  . Although the overall expansion of the system tends to decrease $`B_A`$ values by spatially isolating nucleons from each other, collective motion makes it more likely that nucleons that are spatially close together also have similar momenta, which to some extent works in the opposite direction by increasing $`B_A`$. Other coalescence models have made an effort to include the effects of both larger source volumes and collective motion, by including coalescence as a ’afterburner’ in collision cascade models such as RQMD  and by analytical calculations  .

Light nuclei production can also be calculated from thermal models   which assume at least local thermal equilibrium of the system and thus that particle production (including, in this case, composite particles) is governed by a single temperature and chemical potentials. This gives rise to expressions for the $`B_A`$ with the same $`(1/V)^{(A1)}`$ dependence found in some coalescence models. Collective expansion can be included in these models ; it affects only the amount of energy available in local rest frames for particle production.

There has been much experimental effort in the measurement of light nuclei in heavy ion reactions. Previous measurements at AGS energies, including clusters of $`A4`$, show values for $`B_A`$ that are considerably lower than at Bevalac energies , indicating a much greater expansion of the system. This is also consistent with AGS results showing that the $`B_A`$ become smaller with more central collisions and larger target nuclei  . At CERN-SPS energies, production of secondary particles (chiefly pions) is several times as large as at the AGS, leading to a larger expansion before hadronic freeze-out. Coalescence of high mass clusters is thus much less probable, as indicated by yet lower values for the $`B_A`$ as measured for example by experiment NA52 .

In this paper we describe and report the results of measurements by E864 of the yields of light nuclei in collisions of $`{}_{}{}^{197}Au`$ with beam momentum of 11.5 A GeV/c on targets of $`{}_{}{}^{208}Pb`$ and $`{}_{}{}^{197}Pt`$. The yields are measured for nuclei from baryon number A=1-7. In Section II we briefly describe the experimental apparatus and the analyses used to produce our final invariant multiplicities. In Section III we report the results and compare them with measurements of other experiments where such measurements overlap. Finally in Section IV we examine trends in the data and discuss interpretations in the context of several different models of light nuclei production.

## II Experiment 864

### A Apparatus

Brookhaven AGS Experiment 864 is an open geometry, high data rate spectrometer which was chiefly designed to search for rarely produced objects in Au+Pb collisions. Figure 1 shows a schematic view of the experimental apparatus, a thorough description of which is given in . Event centrality (impact parameter) is characterized by the charged particle multiplicity measured in an annular scintillator counter  located approximately 10 cm downstream of the target, which subtends an angular region from 16.6o to 45o in azimuth when viewed from the target. The products of an interaction travel downstream through two dipole spectrometer magnets, M1 and M2. Charged particle identification is performed using information from the scintillator hodoscope walls (H1, H2, and H3) and the straw tube stations (S2 and S3). The hodoscopes walls each consist of 206 1cm thick scintillator slats placed vertically. They provide information about the charge, time-of-flight, and position of each charged particle hit, and this information is used to identify candidate charged particle tracks. These tracks are then rejected or confirmed and further refined by spatial hit information provided by the straw tube stations. Under the assumption that the track originates in the target, a rigidity is assigned to the track by a look-up table generated from a GEANT simulation of the experimental apparatus (using a technique described in  as applied for the PHENIX experiment at RHIC). With information on rigidity, time-of-flight, and charge, a mass can then be assigned to the track, providing particle identification. A typical charge one mass distribution with a field of 0.45 Tesla in our spectrometer magnets is shown in Figure 2. Mass resolutions of 3 to 7% RMS are typical for particles with velocity $`\beta .985`$.

At the downstream end of the apparatus is our hadronic calorimeter ; an array of 754 towers, each measuring 10cm x 10cm on the front face, each of which provides energy and time information. This is the essential piece of the apparatus in our analyses of yields of neutral particles, for which the tracking detectors serve only to provide a veto. In addition, the energy and timing information from each tower is used to provide the input for a level II high mass trigger (the Late-Energy Trigger or LET ), which provides an enhancement of approximately a factor of 50 in our searches for rare high mass states.

Measurements reported in this paper are from a variety of experimental conditions, including different trigger conditions and different magnetic field settings in M1 and M2. The different data sets and their experimental conditions are listed in Table I. Because of the large acceptance open geometry design of the experiment, the different data sets often have significant regions of overlap with one another, allowing a consistency check on the measurements-see for example the measurements of alpha particles introduced in Section III.

### B Data Analysis

In order to measure the yield of a given species, a mass plot analogous to Figure 2 is made for each kinematic bin. The number of tracks which lie within its mass peak are determined. Background under each peak is then estimated, generally with fits to signal plus background, and subtracted away from this count. The invariant multiplicity in a given kinematic bin is then determined by correcting this number of raw counts for the geometric acceptance, trigger efficiency, charge cut efficiency, track quality cut efficiency, and detector efficiency.

Geometric acceptances are generally 25% or lower as shown in Figure 3 for protons. Charge cut and track quality cut efficiencies are typically 90 to 95%. The three redundant charge measurements for each track allow easy calculation of the charge cut efficiency in each hodoscope simply by examining charge measurements made in each hodoscope against results from the other two; charge misidentification in each of the three hodoscopes is less than one percent, so fewer than one track in a million is assigned an incorrect charge. In general, the track quality cuts are determined by comparison with monte-carlo simulations. In cases where the efficiencies are particularly high, they are determined directly from the data. The total track detector efficiency is 85 to 90%; this is determined by excluding each detector in turn from the track-finding process.

The trigger efficiency is significant only for those measurements which were made using the LET. This mass and momentum dependent efficiency ranges from approximately 40% to 90% for the measurements reported here. The LET efficiency is determined by two methods. The first, used mainly for slow, high mass states for which the efficiency is very high, is done simply using a monte-carlo simulation of the shower generated by the object and knowledge of the LET look-up table in each tower. The second method is from the equation $`ϵ_{LET}=N_{LET}/(N_{LET}+R\times N_{nonLET})`$ where R is the rejection factor provided by the LET (i.e. number of events which fire the LET divided by the total number of events). $`N_{LET}`$ and $`N_{nonLET}`$ are the numbers of particles of interest which do and do not produce LET triggers in LET triggered events, respectively.

Sources of systematic error that we have quantified include possible error in the determination of the efficiencies listed above as well as error in background subtraction (particularly relevant for the deuterons and alpha particles at their highest rapidities). Also examined were the effects of changing the assumed input distribution for each particle species in the determination of geometrical acceptances and efficiencies and effects of possible differences in the magnetic field with the field maps that were used in reconstruction of tracks.

Overall errors are generally dominated by systematics, particularly for the lighter states. Statistical errors can be significant for the heavier states, particularly in the determination of LET efficiencies in which the number of particles of interest which do not fire the LET generally has the largest statistical error.

## III results

Measurements of invariant multiplicities for protons, neutrons, deuterons, $`{}_{}{}^{3}He`$, and $`{}_{}{}^{4}He`$ are shown in Figures 4 through 9, and the values of the data shown in these figures are listed in tables in the Appendix. Figure 4 displays proton invariant multiplicities for three different bins of collision centrality. Figure 5 shows only proton yields from 10% most central Au+Pb collisions along with neutron multiplicities measured by E864 (see Reference ) for comparison. Figures 6 and 7 display deuteron and $`{}_{}{}^{3}He`$ invariant multiplicities for the same centrality bins used for the proton measurements. Tritons, $`{}_{}{}^{4}He`$, $`{}_{}{}^{6}He`$, $`{}_{}{}^{6}Li`$, $`{}_{}{}^{7}Li`$, and $`{}_{}{}^{7}Be`$ are measured by E864 only in 10% most central collisions; yields for tritons and alpha particles are shown in Figure 8 and Figure 9, respectively, while yields for the heavier nuclei are listed in tables in the Appendix. Added detail concerning most of these measurements may be found in the Ph.D. theses listed in Reference .

### A Contributions from Hyperon Decays

For comparison to other experimental results and calculations of coalescence parameters, it is important to quantify the contribution to proton yields that is made by protons which come from decays of hyperon states, which to E864 are indistinguishable from primordial protons. There are three dominant hyperon decays which produce protons: $`\mathrm{\Lambda }p+\pi ^{}`$, $`\mathrm{\Sigma }^0\mathrm{\Lambda }+\gamma p+\pi ^{}+\gamma `$, and $`\mathrm{\Sigma }^+p+\pi ^0`$. The contributions from these decays were evaluated using a GEANT simulation of the experiment with an input distribution taken from measurements of E891 for the Lambda  and an input distribution from RQMDv2.3  for the Sigma. From this simulation it was determined that protons from hyperon decays account for approximately 12% of the measured yields of protons with only a slight kinematic dependence. The proton and neutron yields from E864 shown in Figures 4 and 5 have not been corrected for hyperon feed-down (nor have the values listed in the appendix tables).

### B Comparisons with Other Experimental Results

Figure 10 shows a comparison of the light nuclei measurements from AGS experiments E864, E877 , and E878 . Because of the different beam momenta of the experiments (10.8 A GeV/c in E878), the yields are plotted versus beam normalized rapidity. The yields shown for E864 and E877 are average yields at approximately $`p_T/A`$ = 150 MeV/c, while the E878 yields are measurements at $`p_T0`$. Other caveats to the comparisons of the yields shown in Figure 10 are noted in the figure caption.

Proton yields measured by E864 are clearly higher than measurements by E878. Some of this difference can be attributed to protons which are feed-down from hyperon decay. The acceptance of E878 for these feed-down protons is only about 10% of what it is for primordial protons, while in E864 the two acceptances are nearly the same. When this difference is taken into account (see Section III A) and the E864 yields of primordial protons are lowered by approximately 12%, the results of the two experiments are different by approximately 25% at midrapidity; this lies within the range of systematic errors for the two experiments. At higher rapidities, the agreement is better.

Comparison of the three experiments’ measurements of deuterons and tritons shows close agreement and the measurements of E864 and E878 of alpha particle yields also agree within errors.

## IV Discussion

### A General Trends of the Spectra

#### 1 Transverse Dependence of Yields

E864 has sufficient coverage in transverse momentum for us to extract measurements of inverse slope parameters of the different light nuclei yields in the rapidity range $`2.2y2.4`$. Shown in Figure 11 are yields of protons, deuterons, $`{}_{}{}^{3}He`$, and $`{}_{}{}^{4}He`$ as a function of transverse mass $`m_Tm_0=\sqrt{p_T^2+m_0^2}m_0`$ in this rapidity slice. Overlayed on the measurements are fits of each species to a Boltzmann distribution in transverse mass, from each of which we extract an inverse slope parameter, $`T`$, as noted in Figure 11. The fits from which we extract $`T`$ are linear fits to the log of the invariant multiplicities divided by $`m_T`$; this is the same fitting method used for determining the slope parameters for neutrons in . The $`\chi ^2`$ values are less than one per degree of freedom for all these fits.

These slope parameters, as well as those for neutrons in this same rapidity range, are displayed in Figure 12 as a function of mass number. Polleri et. al.  have demonstrated the sensitivity of these trends in the inverse slopes to the density and velocity profiles of the nucleons at the time when coalescence occurs. To make this point, they have performed calculations of the behaviour of these trends for different assumptions about the source distributions. Two of these assumptions give rise to the two curves shown overlayed on the data in Figure 12. The first has a ’box’ spatial profile and a linear velocity profile, and the second has a Gaussian spatial distribution and a velocity profile $`v(r_T)r_T^{(1/2)}`$. In determining the curves shown in Figure 12 we have fit the data to these two different functional forms with the numerical constraint of T=100 MeV for zero mass. Neither of these sets of model parameters provides an adequate description of the data. These curves are meant only to illustrate the sensitivity of these measurements: clearly, the sensitivity to the differences in these assumptions increases with increasing mass of the measured nuclei.

Shown in Figure 13 are these same trends in light nuclei slope parameters including centralities other than the 10% most central collisions. For the more peripheral events, the trends are consistent with a linear dependence of slope parameter on mass (matching the calculation in Reference  including box density profiles). Only for the most central events is there a clear rollover in the slope as a function of mass.

#### 2 Longitudinal Dependence of Yields

In order to examine the rapidity dependence of the yields of light nuclei, we observe the trends in multiplicities in the $`p_T/A`$ range from 100 to 200 MeV/c. (Our transverse coverage is not sufficient to integrate in transverse momentum for a full measurement of dN/dy over this entire rapidity range.) In Figure 14 we show the invariant multiplicities of deuterons for this low $`p_T`$ range for three different centralities as a function of rapidity. The yields are concave as a function of rapidity (i.e. they are lowest at center-of-mass and increase toward beam or target rapidity) and become more concave for the more peripheral events. We can parameterize this concavity by fitting each set of yields to a quadratic $`a+b(yy_{cm})^2`$. These fits are shown overlayed on the data in Figure 14. The ratio of coefficients, $`b/a`$, from this parameterization serves as a measure of the relative concavity of a species’ spectrum as a function of rapidity at low $`p_T`$ .

Values of this ratio $`b/a`$ are plotted in Figure 15 as a function of mass number $`A`$ for protons, deuterons and $`{}_{}{}^{3}He`$ in three different collision centralities. We observe that the concavity of a spectrum increases with the mass of a species and with collision centrality. This can be understood as some form of collective motion in the longitudinal direction, either expansion or incomplete stopping, which increases in peripheral events. It also is consistent with a rapidity dependent transverse expansion which pushes the nuclei out of the transverse momentum range measured here.

### B $`B_A`$ parameters

The trends noted in Section IV A will also be observable in a study of the behaviour of the $`B_A`$ parameters. For simplicity, we can from our neutron measurements  characterize the neutron spectrum as a factor of 1.19$`\pm `$.08 greater than the proton spectrum (feed-down from hyperons is taken into account in determining this ratio), and so we evaluate $`B_A`$ as

$$B_A\frac{(E\frac{d^3N_A}{dP^3})}{(1.19)^N(E\frac{d^3N_{proton}}{dp^3})^A}.$$
(2)
for a nucleus of baryon number $`A`$ with $`N`$ neutrons and $`Z`$ protons. Again, all invariant yields are evaluated at a common velocity.

In Figures 16 and 17 we show measurements of $`B_2`$ and $`B_3`$ as a function of rapidity and transverse momentum. Because nucleons which are feed-down from weak decays are not available as nucleons for coalescence, we have subtracted the contribution to proton and neutron invariant yields from hyperon feed down for the calculation of $`B_A`$.

The values of the $`B_A`$ in Figures 16 and 17 are clearly not momentum independent as is assumed in many early coalescence models ( and , for example). This is expected given the large amount of evidence for collective flow in heavy-ion collisions . The values of the $`B_A`$ parameters seem to increase slightly with increasing transverse momentum and there is a clear increase away from center-of-mass rapidity. Both of these increases in $`B_A`$ away from center-of-mass momentum are consistent with expectations from an expanding source, although the longitudinal dependence can also be interpreted as a sign of incomplete stopping. Indeed the fact that the invariant yields of antiprotons near $`p_T`$=0 are strongly peaked near center-of-mass rapidity as measured by E864  while the proton yields are essentially flat may be taken as further evidence of incomplete stopping.

### C Source Size Calculations

The $`B_A`$ parameters can be related through coalescence models to source sizes. Following the model of Sato and Yazaki  which uses a density matrix representation of the source distribution and projects it onto a representation of the deuteron wave function, we can relate the $`B_A`$ parameters to the RMS source radius through

$$B_A=(\frac{2J_A+1}{2^A})\frac{A^{5/2}}{m^{(A1)}}(4\pi \frac{\nu _A\nu }{\nu _A+\nu })^{(3/2)(A1)}$$
(3)
with $`\nu _A`$ the size parameter for a cluster with baryon number $`A`$ and spin $`J_A`$ ($`m`$ represents the nucleon mass). This model assumes the absence of collective motion of the nucleons, and therefore a radius (given by $`R_{RMS}=\sqrt{3/2\nu }`$) which is independent of momentum. We evaluate the source size using $`B_2`$, $`B_3`$, $`B_4`$, $`B_6`$, and $`B_7`$ values from our measurements nearest the collision center-of-mass (again assuming a neutron to proton ratio of 1.19) and list the results in Table II. We have used values for the $`\nu _A`$ from  for $`A`$ = 2,3 and 4 and following  have done a polynomial extrapolation to determine $`\nu _6`$ and $`\nu _7`$. This extrapolation for $`\nu _A`$ may be suspect particularly for the halo nucleus $`{}_{}{}^{6}He`$, but the final value for $`R_{RMS}`$ is quite insensitive to the value for $`\nu _6`$. The extracted radius for deuterons $`A`$=2 is considerably larger than the initial size of the colliding nuclei, and the radius parameters decrease for clusters of increasing mass.

A model by Scheibl and Heinz  which includes the effect of collective flow in a density matrix prescription for coalescence leads to an expression for source dimensions:

$$B_2=\frac{3\pi ^{3/2}<C_d>}{2m_TR_{}^2(m_T)R_{}(m_T)}e^{2(m_Tm)(\frac{1}{T_p^{}}\frac{1}{T_d^{}})}$$
(4)
where $`R_{}`$ and $`R_{}`$ are the transverse and longitudinal dimensions of the fraction of the source which contributes to deuteron emission (comparable to the radius parameters extracted in the YKP parameterization of HBT interferometry), $`<C_d>`$ is a quantum mechanical correction factor for the finite size of the deuteron which is evaluated by the authors under various assumptions about the source, and $`T_p^{}`$ and $`T_d^{}`$ are the inverse slope parameters for protons and deuterons. Equation 4 as written assumes a box density profile for the source; a Gaussian profile would result in the absence of the final exponential factor. Plugging in our results for $`B_2`$, we can extract values for $`(R_{}^2(m_T)R_{}(m_T))^{1/3}`$ which are shown in Figure 18 as the solid circles. For the calculations shown here we have used 0.75 for $`<C_d>`$ and the values for $`T_d^{}`$ and $`T_p^{}`$ as measured at rapidity 2.3.

Shown also in Figure 18 as hollow circles are the results of source size calculations with the similar fragment coalescence model of Llope et. al.  via the equation

$$R_c^3=\pi ^{3/2}\frac{(2J_c+1)}{(2J_a+1)(2J_b+1)}\frac{m_c}{m_am_b}\frac{E_a\frac{d^3N_a}{dp_a^3}E_b\frac{d^3N_b}{dp_b^3}}{E_c\frac{d^3N_c}{dp_c^3}}$$
(5)
which relates the effective source radius $`R`$ in the frame of a cluster $`c`$ which may be formed through the coalescence of smaller clusters $`a`$ and $`b`$. (Note that the radius parameters $`R`$ shown in Figure 18 are meant to describe sources of the form $`\rho (r)exp(r^2/2R^2)`$ and so correspond to RMS radii of $`R_{RMS}=\sqrt{3}R`$.)

We observe in Figure 18 that there is a decrease in source size with increasing distance away from the center-of-mass (again, as expected for a source with radial expansion), but our measurements do not give a clear picture concerning scaling of the source size with transverse mass  such as has been noted in results for sizes from HBT two-particle correlations and in measurements of $`B_2`$ at the CERN SPS .

### D Comparison with RQMD

We can also compare our results with predictions from the cascade model RQMD  version 2.3. The complex many-body processes by which light nuclei are formed are not included in RQMD, rather an afterburner  is used to calculate the coalescence of these states based upon an input of the positions and momenta of the nucleons at freeze-out. This model then explicitly includes the position-momentum correlations due to expansion, etc. that are present in RQMD.

Figure 19 displays the E864 measurements along with predictions from RQMD with the coalescence afterburner for comparison over the transverse momentum range $`0.1p_T/A0.2`$ GeV/c as a function of rapidity. For comparison, RQMDv2.3 was run under two different conditions, one including the effect of repulsive mean-field potentials (potentials mode) and one not (cascade mode). We note that there is an increasing disagreement with increasing mass as noted previously in Reference . Calculations of $`{}_{}{}^{4}He`$ production in cascade mode are at a level of approximately 100 lower than our measurements; with potentials on, the discrepancy is still larger by about a factor of two. The level of disagreement varies somewhat in with transverse momentum (the slopes are in fact generally better predicted with potentials on than off) but particularly for the heavier states the change is slight compared to the overall level of disagreement.

### E Scaling of Yields Versus Mass

Shown in Figure 20 are the invariant yields in a small kinematic region at or near $`y=1.9,p_T300MeV`$. Over ten orders of magnitude, the yields in this kinematic bin fit very closely to an exponential dependence with a penalty factor of approximately 48 for each nucleon added (see  and references therein for a discussion of such exponential behavior of cluster yields at lower collision energies). Of course, as we have seen above, the yields for different species have different kinematic dependences due to collective motion, and so this should not be taken as a penalty factor which governs the integrated yields of the various species, which from the kinematic dependences of the $`B_A`$ parameters discussed in Section IV B can clearly be different.

We can make a crude estimate of the difference in these two penalty factors by parameterizing the rapidity and mass dependences of the inverse slope parameters as follows: from our neutron measurements  which extend up to beam rapidity, we can roughly parameterize the distribution of T in rapidity as a Gaussian with a width of 1.1 units. From Figure 13, we also can make a rough parameterization of the mass dependence of the inverse slope as $`T(A+1.0)`$ With these two parameterizations and our penalty factor of 48 at $`y=1.9,p_T300MeV`$, we calculate a penalty factor of 25 in the integrated yields for the addition of an extra baryon to a coalesced state. This should be considered as a lower limit on such a penalty factor since as previously noted the inverse slope parameter rolls over as a function of mass rather than following the linear dependence we have used in this estimate .

We can also estimate this penalty factor in overall yields by following section VI of Reference . Equation 6.10 in this reference allows us effectively to make an estimate of the difference between the penalty factor near $`p_T`$ = 0 and the overall penalty factor in two limiting cases: a static, homogeneous fireball (which would translate our penalty factor of 48 at $`p_T`$ =0 to an overall penalty factor of 72) and a rapidly expanding system (which would result in an overall penalty factor of 39).

### F Spin and Isospin Dependences of Yields

In Figure 21, we display three ratios as a function of transverse momentum in the rapidity ranges y=1.8-2.0 and y=2.0-2.2. The three ratios are the ratios of invariant yield of neutrons over invariant yield of protons, the ratio of yield of tritons over yield of $`{}_{}{}^{3}He`$ and the ratio of $`{}_{}{}^{6}He`$ to $`{}_{}{}^{6}Li`$. The $`n/p`$ and $`t/^3He`$ ratios are consistent with a value of approximately 1.2 (in Reference , we extract values of $`1.19\pm .08`$ and $`1.23\pm .04`$, respectively, for the two ratios in the range $`1.6y2.4`$.). In contrast, the $`{}_{}{}^{6}He/^6Li`$ is much nearer to a value of 0.3. $`p`$, $`n`$, $`{}_{}{}^{3}He`$ and $`t`$ are all spin J=1/2 states while $`{}_{}{}^{6}He`$ is spin 0 and $`{}_{}{}^{6}Li`$ is spin 1. We take this as evidence that the yields scale as the degeneracy factor 2J+1 which is commonly predicted in thermal and coalescence models.

With the dependences upon mass number, isospin and spin divided away one can examine the yields for other dependences. This topic, including a possible dependence on binding energy per nucleon with an inverse slope parameter dependence of a few MeV, is discussed in Reference .

## V Summary

We have shown results of measurements of light nuclei from $`A`$=1 to $`A`$=7. The increase with mass of light nuclei inverse slope parameters appears to roll over at approximately $`A`$=3 in central events but not necessarily in more peripheral events. Also we have shown that the yields near $`p_T`$=0 are concave as a function of rapidity and that this relative concavity increases in more peripheral events and in higher mass nuclei, consistent with both radial expansion and incomplete stopping. These trends are also evident from our observations of the kinematic dependences of the $`B_A`$ parameters. From these parameters we have extracted source dimensions from various models. Efforts to extract more quantitative information about the source from these measurements using the cascade model RQMD with a coalescence afterburner were unsuccessful as predictions of the model differ from our results by an amount that increases with mass and reaches a level of 100 or more by $`A`$=4.

We have also examined the overall scaling of the yields up to $`A`$=7, extracting a penalty factor of about 48 to add a nucleon to a coalesced state near midrapidity at low transverse momentum. This likely translates into a somewhat smaller penalty factor in overall yields for the addition of a nucleon, but we have argued that this is unlikely to differ by as much as a factor of two from our measured penalty of 48.

## VI Acknowledgements

We gratefully acknowledge the efforts of the AGS staff in providing the beam. This work was supported in part by grants from the Department of Energy (DOE) High Energy Physics Division, the DOE Nuclear Physics Division, and the National Science Foundation.

## VII APPENDIX

Shown in Table III through Table XVI are results of measurements by E864 of invariant multiplicities of light nuclei from $`A`$=1 through $`A`$=7 in 10% most central Au+Pb collisions. Results for less central data are also listed for protons, deuterons and $`{}_{}{}^{3}He`$.