text:
# Clues on the evolution of the Carina dwarf spheroidal galaxy from the color distribution of its red giant stars Based on data collected at E.S.O. La Silla, Chile, Proposals No. 63.N-0017 and 64.N–0512

## 1 Introduction

Since its discovery in 1977 (Cannon, Hawarder, & Tritton can+77 (1977)), the Carina dwarf spheroidal (dSph) galaxy has played an important role in our understanding of the nature of dwarf galaxies. Carina was in fact the first dSph galaxy in which a prominent intermediate age population was recognized (Mould & Aaronson mou+aar83 (1983)). The presence of a non-negligible component of old stars was subsequently revealed by the detection of a large number of RR Lyrae variable stars (Saha, Monet, & Seitzer sah+86 (1986); Dall’Ora et al. dall+03 (2003)). Mighell (migh90 (1990)) estimated the old population to constitute $`(17\pm 4)\%`$ of the galaxy population, a number confirmed by Smecker-Hane et al. (smec+94 (1994)) and Mighell (migh97 (1997)). A detailed description of the star-formation history (SFH) of Carina was given by Smecker-Hane et al. (smec+94 (1994), smec+96 (1996)). A deep $`BR`$ color-magnitude diagram (CMD) revealed the presence of three distinct turn-offs, pointing at three separate episodes of star formation approximately $`2`$, $`36`$ and $`1113`$ Gyr ago (see also Monelli et al. mone+03 (2003)). A quantitative assessment of this scenario was then provided by CMD simulation techniques applied to CTIO 4m and HST data (Hurley-Keller, Mateo, & Nemec hurl+98 (1998); Hernandez, Gilmore, & Valls-Gabaud her+00 (2000); Dolphin dol02 (2002)).

Quite surprisingly given this complex star-formation history, Carina shows a very thin red giant branch with all the different turn-offs connecting at the RGB base (Smecker-Hane et al. smec+94 (1994), smec+96 (1996)). This seems to imply that either no significant metal enrichment took place in Carina, or that age and chemical enrichment of the stellar populations played in opposite directions to yield a narrow RGB. While these alternative views have often been discussed, they have never been quantitatively addressed. The chemical evolution of dwarf spheroidal galaxies, in particular the role of galaxy-wide gas outflows and/or infall, remains a controversial issue and many basic questions are still open. Recent efforts to model the chemical evolution of dwarf spheroidal galaxies have been presented by Ikuta & Arimoto (ikut+arim02 (2002)), Carigi, Hernandez, & Gilmore (cari+02 (2002)), and Tolstoy et al. (tols+03 (2003)). All groups use known SFHs to infer the chemical-enrichment histories of dSph galaxies, but the scenarios differ in the role played by galactic winds and gas infall.

This Letter shows that realistic modeling of the color distribution of red giant stars by population synthesis methods can help constraining the chemical evolution of dwarf spheroidals galaxies. The intrinsic color width of the upper RGB in Carina is quantified here by exploiting the high statistics of a new $`BVI`$ database from wide-field observations. The star-formation history of Carina is then derived for a set of different metal-enrichment histories (the age-metallicity relation of this galaxy is still uncertain) by comparing observed and synthetic color-magnitude diagrams. The narrow RGB of Carina is shown to be consistent with the star-formation history derived from the observed CMDs, under reasonable assumptions on the chemical-enrichment history.

## 2 The data

Observations of Carina were carried out using the Wide Field Imager at the ESO/MPG 2.2m telescope. The camera consists of a mosaic of eight 2k $`\times `$ 4k EEV CCDs, in a square array of 2 $`\times `$ 4 detectors. In this configuration, the total number of pixels is 8k $`\times `$ 8k, with a scale of 0.238 arcsec pixel<sup>-1</sup> and a field of view $`32\mathrm{}\times 32\mathrm{}`$. Three different fields were observed: a central field and two outer fields with $`20\mathrm{}`$ offsets. Sky conditions were photometric, and standard star fields taken from the list of Landolt (lan92 (1992)) were observed for calibration purpose. A full description of the data and reduction methods will be presented in a forthcoming paper (E. V. Held et al. 2003, in prep.) while just a brief outline is given here. Pre-reduction of CCD images was performed within the IRAF environment. Each image was bias-subtracted and flat-fielded using twilight sky flats. After these steps, all images were astrometrically calibrated using the IRAF mscred package (Valdes vald98 (1998)) and the script package wfpred developed at the Padova Observatory.

Photometry was obtained using daophot ii/allstar (Stetson ste87 (1987)) on the co-added images, and the results were calibrated using calibration relations computed from standard stars observed in all CCDs. Similarly, aperture corrections were computed using a growth curve analysis performed on bright isolated stars in each CCD. Artificial star tests were done to evaluate the uncertainty of photometry and the degree of completeness of the data. The photometric errors in the $`V`$ band range from $`0.01`$ mag near the tip of the RGB, to $`0.3`$ mag near the detection limit.

## 3 Modeling the red giant branch of Carina

### 3.1 RGB width measurement

These wide field data allowed us to obtain a sound quantitative estimate of the RGB width in Carina. In order to measure the color spread of the RGB, a fiducial ridge line was fitted to the mean colors of red giants in $`0.5`$ mag bins. The distribution of color residuals about the fiducial line was then measured down to 1.5 mag below the tip. After quadratic subtraction of the contribution of photometric errors, the estimated intrinsic color spread of the upper RGB is:

$`\sigma _{(VI)}=0.021\pm 0.005`$

where the statistical error was estimated by bootstrap resampling of the RGB sample in 500 trials.

Using a direct comparison of the red giant branch of Carina with the ridge lines of Galactic globular clusters (Da Costa & Armandroff dac+arm90 (1990)), we estimated a mean metallicity $`[\mathrm{Fe}/\mathrm{H}]=2.08\pm 0.12`$ on the Zinn & West (zinn+west84 (1984)) scale (E. V. Held et al. 2003, in prep.). This value agrees well with previous photometric determinations (Mould & Aaronson mou+aar83 (1983); Mighell migh90 (1990); Smecker-Hane et al. smec+94 (1994); Mateo, Hurley-Keller, & Nemec mat+hur98 (1998)), as well as the low dispersion spectroscopic measurements of the Ca ii lines of Da Costa (dac94 (1994)) and Smecker-Hane et al. (smec+99 (1999)). A somewhat higher metallicity, $`[\mathrm{Fe}/\mathrm{H}]=1.6`$ with a spread of 0.5 dex, has recently been derived by Tolstoy et al. (tols+03 (2003)) from high-resolution spectroscopy of 5 red giant stars.

The difference between this new spectroscopic estimate and previous photometric \[Fe/H\] value can be only partly accounted for by age effects on the RGB color. The stars in Carina are, in fact, $`7`$ Gyr younger than the stellar populations of Galactic globular clusters (see below), and we know that younger isochrones are shifted blueward (see, e.g., Da Costa daco98 (1998)). Using Girardi’s et al. (gira+00 (2000)) models in the range $`Z=0.0001`$ to 0.0004, we find the effects of a younger age on the $`(VI)`$ color to be $`0.035`$ mag at the RGB level used for metallicity estimation. This bluer color mimics a lower metallicity, by $`0.17`$ dex (see also Smecker-Hane et al. smec+94 (1994)). By taking this effect into account, the age-corrected metallicity of Carina is then \[Fe/H\]$`=1.91\pm 0.12`$.

The color distribution of RGB stars formally corresponds (for an old population) to a metallicity dispersion $`\sigma _{[\mathrm{Fe}/\mathrm{H}]}=0.10\pm 0.04`$, which is significantly smaller than the quoted dispersion of spectroscopic measurements.

### 3.2 Population synthesis

Figure 3.1 (left panel) shows the CMD of the the central part of Carina, for which simulations were performed. An elliptical region with a 12 arcmin semi-major axis, $`65^{}`$ position angle, and $`0.33`$ ellipticity (Irwin & Hatzidimitriou irw+hat95 (1995)) was selected to this purpose. The intrinsically poor statistics of the red giant stars in Carina prevents reliable measurements of the RGB width in the outer regions. The observed CMD was simulated using a parametric technique similar to that employed by Aparicio, Gallart, & Bertelli (apa+97a (1997)), and a description of our methods is presented by Rizzi et al. (riz+01 (2002)). In brief, our technique assumes a number of known input parameters, such as the initial mass function (IMF), the distance to the galaxy, and the reddening. A Salpeter (sal55 (1955)) IMF is adopted, since recent work suggests a general validity of this IMF for masses larger than $`0.5M_{}`$ (Scalo sca98 (1998); Kroupa kro00 (2000)). A distance modulus $`(mM)_0=20.00\pm 0.06`$ and a reddening $`E(BV)=0.045\pm 0.01`$ were adopted (E. V. Held et al. 2003, in prep.). Using these input parameters, a set of synthetic stellar populations were generated using models from the library of Girardi et al. (gira+00 (2000)) and the simulation code ZVAR that linearly interpolates the evolutionary tracks both in age and metallicity, including the effects of binary stars. The observational effects were applied to the theoretical random populations using the photometric errors and completeness obtained from the results of artificial star experiments.

A basic ingredient of the simulation is the chemical-enrichment history (CEH), for which we chose to test a number of reasonable assumptions (Fig. 3.2). The first tested law is characterized by a prompt initial enrichment at ages of $`12`$ Gyr, which is responsible for producing much of the observed metals, followed by a slow metallicity increase leading to a final value of $`[\mathrm{Fe}/\mathrm{H}]1.9`$ (hereafter “PT law”, after the age-metallicity relation modeled by Pagel & Tautvaisiene pag+tau98 (1998) for the Large Magellanic Cloud). This trend is similar to the CEH recently derived from spectroscopic observations of dwarf spheroidal galaxies (Gallart et al. gal+02 (2002); Tolstoy et al. tols+03 (2003)). The final metallicity was constrained to yield a mean metallicity over the galaxy lifetime equal to the value derived from the photometry of Carina. This was done iteratively since the galaxy’s star-formation history is affected by the CEH normalization, and vice versa.

Other tested laws, illustrated in Fig. 3.2, include a constant metallicity law at $`[\mathrm{Fe}/\mathrm{H}]=1.91`$ and a linear law starting at $`[\mathrm{Fe}/\mathrm{H}]=2.3`$ and reaching $`[\mathrm{Fe}/\mathrm{H}]=1.76`$ at the present epoch. These laws sample a wide range of possibilities in the rate of metal enrichment, still yielding the correct mean \[Fe/H\] over the SFH. Other chemical evolution laws with a more metal-rich normalization would be inconsistent with the photometric data.

A total of 7 populations were generated using these input laws, with age bins centered at 1, 2.25, 4, 6, 8.5, 11 and 13 Gyr, each bin containing 30.000 stars. Tests were made with both a 30% fraction of binary stars and no binaries at all. The SFH that best fits the observed color-magnitude diagram was found by comparing the number of real and simulated stars inside a number of regions in the CMD (“boxes”). These were chosen to map (i) the luminosity function along the main sequence, (ii) the relative contribution of blue and red horizontal branch stars, (iii) the number of “red clump” stars, and (iv) the total number of stars along the RGB (see Fig. 3.1, left panel). The width of the RGB itself is not an input constraint for the best-fit CMD.

The best fitting synthetic diagram, obtained by minimizing the differences in star counts between the simulated and observed CMD with a downhill $`\chi ^2`$ minimization algorithm, is shown in the right panel of Fig. 3.1 for the PT law.

The derived star-formation histories are shown in Fig. 3.2, for the different assumptions on the metal-enrichment history. Our analysis of wide-field data agrees with previous results in that $`75`$% of the stars formed between 3 and 7 Gyr ago, with a small yet significant old stellar component. Note that, although the details of the SFH are different, the main episodes of star formation are relatively independent of the adopted CEH. In fact, our suitable choice of the CMD regions over which stars are counted makes the derived SFH more sensitive to the lifetimes of evolutionary phases than to the metal enrichment law (and uncertainties on the isochrone colors).

### 3.3 A thin RGB in Carina

We can now return to the problem of the origin of the extremely thin RGB in Carina, having in hand the necessary tools to quantitatively evaluate the effects of the chemical-enrichment and star-formation history. The color distribution of simulated RGB stars about a fiducial line was measured for each choice of the CEH, using the same methods as for the real data. We did not make use of the absolute colors of the synthetic RGBs, because of known uncertainties on the isochrone colors. The mean widths measured for the simulated diagrams are $`\sigma _{(VI)}=0.023\pm 0.003`$ for the linear and constant laws, and $`\sigma _{(VI)}=0.025\pm 0.003`$ for the PT law. The mean values and statistical errors result from 300 trials. These results are essentially independent of the adopted fraction of binary stars. The simulations thus indicate that the Carina’s narrow RGB naturally arises from the SFH derived from the data. The results for all the tested enrichment laws agree with the observed color spread within one standard deviation. Indeed, these values probably mark the lower limit to the measurable RGB width set by stochastic effects and photometric errors in the CMD simulations. As expected, no conclusions can be reached on the details of chemical enrichment from photometry alone. A large sample of spectroscopic measurements will help, provided that it is of the highest precision and accuracy.

## 4 Discussion and conclusions

In this Letter, we have reported new measurements of the mean metallicity (\[Fe/H\]=$`1.91\pm 0.12`$, after correction for the effects of age) and the color spread ($`\sigma _{VI}=0.021\pm 0.005`$) of red giant stars in Carina, based on new wide-field observations. Using synthetic color-magnitude diagrams to model the color distribution of RGB stars, we have quantitatively shown that the thinness of the red giant branch is a direct consequence of the derived star-formation history, for most chemical evolution laws that are plausible and consistent with the observed mean metallicity of the galaxy.

The most fundamental reason for the RGB being so thin is the fact that the predominant intermediate-age stellar population that makes up the RGB, was formed in a limited time interval ($`5\pm 2`$ Gyr ago), with only a minor contribution by an old population. In theory, both the age range and the metallicity evolution of the red giant stars contribute in some way to the color range of the RGB, but the effects are subtle and hardly measurable. In fact, (i) the age interval is too small for the isochrone separation to be measurable, and (ii) the sensitivity of the RGB color to changes in \[Fe/H\] is low for metal-poor systems, so the effect of modest metal enrichment during the main star-formation episode is negligible too.

Our photometric analysis bears some interesting implications on the chemical evolution history, by providing evidence that Carina could not experience a significant metal enrichment (i.e. above \[Fe/H\]$`1.9`$), certainly not before 3 Gyr ago. In fact, the lifetime-averaged metallicity inferred from photometric observations definitely rules out a metal-rich final composition. Little information is available on the most recent epoch, because stars younger than $`3`$ Gyr are only a trace population, but for the same reason no significant metal enrichment is expected in the last few Gyr. The most likely scenario remains that indicated by chemical evolution models, in which substantial metal enrichment is associated with the earliest star-formation epoch ($`12`$ Gyr).

In summary, both the episodic star-formation activity (with long-lasting episodes and quiescent phases) and the lack of a significant metal enrichment suggest that metal-enriched gas outflows play a role in the evolution of Carina and other low-mass galaxies (Vader vade86 (1986); Mac Low & Ferrara macl+ferr99 (1999); see also Gallagher & Wyse gal+wys94 (1994); Kunth & Östlin kun+ost00 (2000); Carigi et al. cari+02 (2002)). Closed-box models (e.g., Ikuta & Arimoto ikut+arim02 (2002)) appear unable to explain the very low metal abundances of dwarf galaxies like Carina. In a similar way, the oxygen abundances measured in the interstellar gas of dwarf irregular galaxies appear hardly consistent with a closed-box model (see, e.g., Saviane et al. savi+02 (2002), and references therein). Observational evidence of galactic winds provides support to this scenario for active starbursts, a recent example being the Chandra observation of a metal-enriched galactic wind in the dwarf galaxy NGC 1569 (Martin, Kobulnicky, & Heckman mart+02 (2002)).

We thank the referee for suggesting some significant improvements to the paper, Drs. L. Carigi and X. Hernadez for helpful remarks, and M. Tosi for useful conversations on dwarf galaxy evolution. Support for this work was provided by the Italian Ministry for University and Research through grants COFIN2001028897 and COFIN2002028935.


format: text

text:
# Ghosts of the Milky Way: a search for topology in new quasar catalogues

## 1 Introduction

The idea that we inhabit a compact multi-connected Universe has been presented by many authors (see review by Lachièze–Rey & Luminet (1995)). The simplest such Universe is one in which space is flat and may be visualised as tiled by a ‘fundamental parallelepiped’ with lengths $`\alpha _x`$, $`\alpha _y`$, $`\alpha _z`$. If the sides of the fundamental domain are significantly smaller than the horizon length, RH, observable effects would be present. One of the most interesting, albeit surreal, effects would be the appearance of ‘ghost images’. This is where radiation from an object, such as a galaxy, has traversed the fundamental domain, resulting in a repeat or ghost image of the object. The simplest ghost images to search for are ones of our own galaxy, since they appear at the nodes of a lattice of which the Milky Way is the origin. Therefore ghost images might be most simply identified by searching for antipodal pairs. In this paper we merge three new large quasar catalogues to search for ghost images of the Milky Way.

The best limits on the size of the fundamental domain come from analyses of the COBE measurement of the cosmic microwave background. The simplest space is the hypertorus. This is a spatially flat topology where opposite faces of the parallelepiped are identified. There are three other compact spaces, involving twists of $`\pi /2`$ or $`\pi `$ in pairing opposite faces of a parallelepiped (e.g. Levin 1998). There are also two compact models in which the fundamental polyhedron is a hexagonal prism with twists of $`\pi /3`$ and $`2\pi /3`$ (Lachièze–Rey & Luminet 1995). For equal sides of the parallelepiped, Levin (1998) obtained a minimum size of the side of the fundamental domain of 0.8RH, and a similar result for the hexagonal prism. Bond, Pogosyan & Souradeep (2000) obtained a somewhat larger limit for the hypertorus. For unequal sides the constraints are less strong, and there remains the possibility that there are many copies of the Universe within the horizon. The existence of antipodal pairs is a common feature of most topologies whether the geometry is hyperbolic, flat or spherical. To this extent, therefore, although we assume a flat geometry our analysis is also relevant to nearly-flat geometries.

On the basis of these COBE limits, it is clear that any search for ghosts needs to use objects at large redshifts, and so catalogues of quasars are an obvious source. Many authors have used quasars to probe topological structure (Fagundes & Wichoski 1987; Fagundes 1989; Roukema 1996), but have yet to constrain topology using these methods. Following the suggestion that many large galaxies go through a quasar phase, Fagundes & Wichoski (1987) searched a quasar catalogue for antipodal quasar pairs in the hope of finding ghost images of the Milky Way. In fact even with complete catalogues of quasars over the whole sky, good fortune would be required for success with such a search. This is because quasar lifetimes must be short relative to the age of the Universe, so the method relies on the quasar being alive at the particular time corresponding to the look-back time of any copies. Nevertheless quasar catalogues have grown enormously in the last few years, due in particular to the contributions of the 2-degree-field QSO survey (hereafter 2dF) and the Sloan Digital Sky Survey (hereafter SDSS). Therefore we considered it timely to repeat the search of Fagundes & Wichoski for antipodal quasar pairs.

The outline of the paper is as follows. In §2 we describe the quasar catalogues used. In §3 we explain the search method. The search is complicated by the peculiar velocity of the Milky Way $`v_{}`$, which displaces the ghost images from opposition by order $`v_{}/c`$, and so must be accounted for. A suitable search radius, then, depends on the uncertainty of the historical motion of the Milky Way, which depends, in turn, on the uncertainty of the cosmological density parameter. The results of the search are presented in §4, together with a brief discussion. We assume flat Universal geometry $`(\mathrm{\Omega }_\mathrm{m}+\mathrm{\Omega }_\mathrm{\Lambda }=1)`$ throughout and consider only topologies without axial twists.

## 2 Quasar Catalogues

The quasar catalogue used for this investigation was compiled from three sources:

1. The complete candidate list for the 2dF QSO survey (Smith et al., 2003, in prep.). This catalogue contains 47,768 objects, from which we extracted the 23,340 quasars with spectroscopic redshifts. A subset of 10,000 of the quasars has been published in Croom et al. (2001).
2. The early data release from the Sloan Digital Sky Survey Quasar Catalogue (Schneider et al., 2001) which contains 3,814 quasars.
3. The quasar catalogue of Véron-Cetty & Véron (2002)111available at: http://www.obs-hp.fr/www/catalogues/veron2\_10/veron2\_10.html which contains 23,760 quasars.

Duplicate quasars, primarily the 10,000 quasars in Croom et al. (2001) which are included in the Véron catalogue, need to be eliminated. We searched the merged list for pairs of objects within 2 arcsec and $`\mathrm{\Delta }z<0.02`$, and in each case eliminated one of the pair. This could, in fact, have included real pairs, but this would not affect our search for antipodal pairs, since it uses a much larger search radius. This left a final catalogue of 43,271 objects reaching a maximum redshift of $`z=5.8`$. The sky distribution of the quasars in this merged catalogue is illustrated in Figure 1.

## 3 Method

### 3.1 Peculiar Velocity

Because of the past motion of the local group, ghost images of the Milky Way will appear displaced from the lattice nodes. This angular displacement is of order $`v_{}/c7`$ arcmin, and must be corrected for. The effect is greatest for objects lying perpendicular to the direction of local group motion, and falls to zero for objects that lie on the path of local group motion. Therefore the idea is to apply an appropriate angular correction for each object in our catalogue to bring true antipodal ghost images into opposition. The uncertainty of the computed correction sets the search radius, since it is much greater than the uncertainty of the quasar coordinates, or any other relevant effect. Previous searches for antipodal quasar pairs have failed to correct for the motion of the local group.

To find this angular correction we turn to linear theory for an expression for peculiar velocity. From this we can find the comoving distance moved by the local group as a function of redshift, and from this compute the appropriate angular correction for each quasar.

The peculiar velocity, $`v_{pec}`$, is related to the peculiar acceleration, $`𝒈`$, by

$$𝒗_{pec}=\frac{a(z)}{a(0)}\dot{𝒙}=\frac{2f}{3\mathrm{\Omega }H}𝒈$$

(Peebles, 1993) where $`a(z)`$ is the scale factor at epoch $`z`$ and $`\dot{𝒙}`$ is a comoving velocity. Here

$$\mathrm{\Omega }=\mathrm{\Omega }_\mathrm{m}\left(\frac{H}{H_{}}\right)^2(1+z)^3.$$

The term $`f`$ is the dimensionless velocity factor, given by

$$f=\frac{d\mathrm{log}\delta }{d\mathrm{log}a}$$

where $`\delta `$ is the overdensity or density contrast. For a flat universe, a good approximation to $`f`$ is given by

$$f=\mathrm{\Omega }^{0.6}+\frac{1}{70}\left[1\frac{1}{2}\mathrm{\Omega }(1+\mathrm{\Omega })\right]$$

(Lahav et al., 1991).

In our linear approximation the peculiar acceleration is given by

$$𝒈=\frac{GM_p(𝒙_p𝒙_{lg})}{a^2|𝒙_p𝒙_{lg}|^3}$$

where $`𝒙`$ is a comoving distance and the subscripts refer to a single point-attractor, p, and the local group, lg. This equation provides the scaling of $`𝒈`$ with $`a`$, but the mass of the attractor, and its distance, are not needed explicitly, since we have a measurement of the peculiar velocity today.

The expression for peculiar velocity is numerically integrated to find the comoving distance moved by the local group as a function of redshift. This is then converted to an angular correction for a ghost image at redshift $`z`$. The only input variables to this correction algorithm are $`\mathrm{\Omega }_\mathrm{m}`$ and $`𝒗_{}`$. The values adopted for this analysis are:

1. $`\mathrm{\Omega }_\mathrm{m}=0.3\pm 0.1`$ (e.g. Verde et al. 2002).
2. $`𝒗_{}=627\pm 22`$km s-1 in direction $`(l,b)=(276\mathrm{°}\pm 3\mathrm{°},30\mathrm{°}\pm 3\mathrm{°})`$ from COBE data (Kogut et al., 1993).

The angular correction for each object depends on redshift and the angle between the quasar, the observer, and the direction of local group motion. For an object located perpendicular to the direction of local group motion (maximum correction) this correction is approximately 7 arcmin at $`z=0`$, and increases to approximately 12 arcmin at $`z=5`$. A plot of this angular correction is shown in Figure 2.

### 3.2 Search Radius

An appropriate search radius must be related to the uncertainties in the input data, and also involves a trade-off between contamination and incompleteness. Too large a tolerance and we recover too many false positives, swamping any real pairs, while with too small a tolerance we risk losing the few putative real pairs in the data.

We required that the redshifts of the pair agree within a tolerance $`\mathrm{\Delta }z=0.04`$. The large majority of the redshifts in the catalogue have been measured from intermediate resolution spectra, many by cross correlation against a template spectrum. For these the accuracy of the redshift is better than $`\mathrm{\Delta }z=0.01`$ (e.g. Croom et al. 2001). Therefore the chosen tolerance will include nearly all real pairs in the catalogue.

The position of any quasar in the catalogue, after correction for the peculiar velocity of the local group, is uncertain, because of the uncertainty of the parameters $`\mathrm{\Omega }_\mathrm{m}`$ and $`𝒗_{}`$. Gravitational lensing could also contribute, since the light paths from antipodal quasars traverse different regions of the fundamental domain. This is a much smaller source of uncertainty, however, and may be neglected. It would be very unlikely for a ghost pair to be lensed out of opposition by more than a few arcsec.

In principle we could compute a probability distribution function (pdf) for the coordinates of each quasar, and select a contour level at, say, the $`95`$ per cent confidence level. Then, in searching for an antipodal quasar, the size of the contour should be doubled. This is because the error in the peculiar velocity correction (due to poor choice of $`\mathrm{\Omega }_\mathrm{m}`$, $`𝒗_{}`$) is identical for each quasar of an antipodal pair. This procedure would be computationally expensive, and we followed a simpler Monte Carlo approach. The algorithm to correct the coordinates was run on the catalogue 200 times, each time using values of the input parameters $`\mathrm{\Omega }_\mathrm{m}`$, $`𝒗_{}`$, chosen from a Gaussian distribution about their mean. We computed the angular displacement from the mean position for each quasar in every simulation, and selected the angle which included $`95`$ per cent of all the points. This angle was $`67`$ arcsec. Therefore the procedure for the search was to correct the quasar coordinates using the best estimates for $`\mathrm{\Omega }_\mathrm{m}`$, and $`𝒗_{}`$, and then for each quasar to search for antipodal quasars in the catalogue within a search angular radius of $`134`$ arcsec. Because the coordinate pdfs for each quasar differ in size, shape, and orientation, in some cases this search radius will be unnecessarily large, and in some cases too small. But, by construction, it will include $`95`$ per cent of any real antipodal pairs.

We found that the largest source of error is from the uncertainties in our current knowledge of the speed and direction of the local group relative to the CMB. Our motion within the local group is a source of error unaccounted for in our search radius determination. As we are searching for a ghost of the Galaxy rather than a ghost of the local group, our motion within it should be considered. However, this motion is small compared with the large scale streaming motion of the local group itself.

## 4 Results and discussion

With a search radius of $`134`$ arcsec and $`\mathrm{\Delta }z=0.04`$ we recover seven antipodal quasars pairs. Details of the pairs are provided in Table 1. The table includes the coordinates of each quasar in the pair both before and after correction for peculiar velocity as well as the redshifts, brightnesses, the catalogue containing the quasar, the redshift difference, and the angular displacement from opposition of each pair.

The relative positions of the two quasars in each pair are plotted in Figure 3. For each pair the plot is centred on the antipode of the first quasar in the pair, and the points illustrate the offset of the second quasar from perfect anti-alignment, for the 200 Monte Carlo realisations of the parameters $`\mathrm{\Omega }_\mathrm{m}`$, and $`𝒗_{}`$. In each case the distribution of points overlaps perfect anti-alignment, indicating that, given the uncertainties in the parameters, the pairs could be real ghost images. As a further investigation, the angles between the axes of each quasar pair were computed to search for any specific angles that could indicate points on a lattice (such as $`\pi /2`$), but nothing striking was discovered. There are magnitude differences, up to $`\mathrm{\Delta }m=2`$, between the two quasars of a pair, but these might be due to variability because of the epoch differences of the observations, or time delay differences because the light paths to us traverse different regions of space, or the differences could be because each quasar is being viewed from a different direction. For the same reasons it would not be unreasonable for the spectra to differ also. A more important consideration is the number of pairs that might be expected by chance.

We investigated how the expected number of antipodal pairs varied with respect to the search tolerance. We varied the search angular radius from $`\mathrm{\Theta }=0`$ to $`\mathrm{\Theta }=10`$ arcmin, and the redshift tolerance from $`\mathrm{\Delta }z=0`$ to $`\mathrm{\Delta }z=0.4`$. Although the catalogue is inhomogeneous we found that the number of antipodal quasar pairs recovered followed closely the relations $`\mathrm{\Theta }^2`$ and $`\mathrm{\Delta }z`$ expected for unrelated pairs. Therefore we used a fit to the relations at large $`\mathrm{\Theta }`$, $`\mathrm{\Delta }z`$ to establish an estimate of the number of pairs expected by chance. For our search tolerance of $`\mathrm{\Theta }=134`$ arcsec we would have expected $`12\pm 3.5`$ quasar pairs by chance. Our result of seven pairs, therefore, is consistent with chance. While it remains possible that one or more of the pairs could be a real ghost pair, this cannot be proven from the existing data.

### 4.1 Discussion

We have presented a search for ghost images of our galaxy using the method of antipodal quasar pair detection, correcting for the angular distortion due to the history of motion of the local group. We have found no significant difference in the number of pairs found over the number we would expect purely by chance using our search tolerance on this catalogue. We have examined the angles between the axis of each quasar pair but failed to discover any indication of a lattice structure.

The uncertainties both with our present knowledge of cosmological parameters and the peculiar velocity of the local group make a search such as this increasingly difficult as quasar catalogues increase in size. To prove significance in a search such as this the search tolerance needs to be reduced to a value where we would be expecting only one or two detections purely by chance. In the future we can expect improvements in the accuracy of the parameters $`\mathrm{\Omega }_\mathrm{m}`$, $`𝒗_{}`$, that will allow a smaller search radius, that could compensate for the much larger future quasar catalogues. However, soon the uncertain correction for the peculiar motion of the Milky Way within the local group will become a significant factor. Given the COBE constraints on the size of the fundamental domain (and so the small possible number of ghost images), and the relatively short lifetime of quasars, it would seem unlikely that it will be possible to identify putative ghost images, unambiguously, by this method. Instead, extension of the search to configurations of three or more quasars on a lattice could yield a convincing detection.

## Acknowledgments

We are grateful to Andrew Jaffe, Thomas Babbedge, Matthew Pieri and Matthew Fox for their helpful comments. The 2dF QSO Redshift Survey (2QZ) was compiled by the 2QZ survey team from observations made with the 2-degree Field on the Anglo-Australian Telescope. Funding for the creation and distribution of the SDSS Archive has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, NASA, NSF, the U.S. Department of Energy, the Japanese Monbukagakusho, and the Max Planck Society. The SDSS website is http://www.sdss.org.


text:
# Shortcut to a Fermi-Degenerate Gas of Molecules via Cooperative Association

## Abstract

We theoretically examine the creation of a Fermi-degenerate gas of molecules by considering a photoassociation or Feshbach resonance applied to a degenerate Bose-Fermi mixture of atoms. This problem raises an interest because, unlike bosons, fermions in general do not behave cooperatively, so that the collective conversion of, say, two million atoms into one million molecules is not to be expected. Nevertheless, we find that the coupled Fermi system displays collective Rabi-like oscillations and adiabatic passage between atoms and molecules, thereby mimicking Bose-Einstein statistics. Cooperative association of a degenerate mixture of Bose and Fermi gases could therefore serve as a shortcut to a degenerate gas of Fermi molecules.

Pacs number(s): 03.75.Ss, 05.30.Fk, 34.10.+x

Laser cooling of molecules is difficult, if not impossible, to apply in practice. This bottleneck occurs because the removal of energy from the translational degrees of freedom is most often accompanied by excitation of ro-vibrational modes which, upon de-excitation, couple back into the translational motion in a manner that heats the system. A laser-cooling-based approach to a Fermi-degenerate gas of molecules is therefore unlikely to succeed. Moreover, buffer-gas cooling  and Stark-deceleration  techniques, while no doubt of practical interest in their own right, currently view quantum degeneracy from the horizon.

One possible route to ultracold molecules is photoassociation, and the formation of nondegenerate Fermi molecules was indeed recently observed. A shortcut to quantum degenerate molecules is in the early stages of development for Bose-Einstein condensates (BECs), so far delivering strongly enhanced molecule formation, precise measurements of the light-induced shift of the binding energy, and searches for a fundamental limit to the atom-molecule conversion rate. The magnetic field counterpart to photoassociation is the Feshbach resonance. Often seen as a tool for adjusting the scattering lengths, thereby enabling condensation in otherwise-incondensible systems, the latest studies of a BEC tuned near a Feshbach resonance have culminated in evidence for atom-molecule coherence.

Motivated by these exciting experiments, we report that a degenerate gas of Fermi molecules could be formed by applying either a photoassociation  or Feshbach resonance  to a degenerate mixture of Bose and Fermi atoms . The collective Fermi atom-molecule system is predicted to undergo Rabi-like oscillations, as well as adiabatic passage, from atoms to molecules, thereby mimicking cooperative behavior that was previously attributed to Bose statistics .

Photoassociation occurs when two atoms absorb a laser photon, thereby jumping from the free two-atom continuum to a bound molecular state. If the initial atoms form a Bose-Einstein condensate, then a photoassociation laser could conceivably be used to convert a BEC of atoms into a BEC of molecules . Collective free-bound photoassociation is theoretically identical to magnetoassociation, whereby a molecular condensate can be created when one atom from a condensed pair spin flips in the presence of a magnetic field tuned near a Feshbach resonance . Intuition developed in one instance is therefore generally applicable to the other, and we will often refer simply to collective association.

Analogous with coherent optical transients in few level atomic systems , photoassociation of a BEC has been predicted to induce Rabi-like oscillations between atomic and molecular condensates , whereby an entire gas of, say, two million Bose-condensed atoms are collectively converted into a million molecules that are, in turn, collectively converted back into (roughly) two million atoms, ad infinitum. Another interesting possibility arises because the ground state of the system is all atoms for large positive detunings (far below threshold) and all molecules for large negative detunings (far above threshold), so that a slow sweep of the laser detuning from one extreme to the other will collectively convert a BEC of atoms into a BEC of molecules . Incidentally, it was a particular combination of these two concepts, applied instead to magnetoassociation, that led to the observation of collective Ramsey fringes between an atomic condensate and a small fraction of molecular condensate dressed by dissociated atom pairs .

The statistics of neutral atoms is determined by the number of neutrons in the nucleus, which must be odd for fermionic atoms and even for bosonic atoms. The sum of the total number of neutrons in the nucleus of the constituent atoms similarly determines the statistics of neutral molecules. Accordingly, molecules formed by free-bound association of two fermions will necessarily result in a boson, whereas fermionic molecules are born from the union of a boson and a fermion. Given a degenerate mixture of Bose and Fermi gases , is it possible that collective free-bound association could serve as a source of degenerate Fermi molecules?

To address this question, we model a degenerate Bose-Fermi mixture of atoms  coupled by either a Feshbach or photoassociation resonance to a Fermi-degenerate gas of molecules. The initial bosonic \[fermionic\] atoms are denoted by the field $`\varphi (𝐫,t)`$ \[$`\psi _{}(𝐫,t)`$\], and the fermionic molecules by the field $`\psi _+(𝐫,t)`$. Neglecting particle-particle interactions, the Hamiltonian density for such an untrapped system can be written $`=_0+_I`$, where

$`{\displaystyle \frac{_0}{\mathrm{}}}`$ $`=`$ $`\varphi ^{}\left[{\displaystyle \frac{\mathrm{}^2}{2m_B}}\mu _B\right]\varphi +\psi _{}^{}\left[{\displaystyle \frac{\mathrm{}^2}{2m_{}}}\mu _{}\right]\psi _{}`$ (3)
$`+\psi _+^{}\left[{\displaystyle \frac{\mathrm{}^2}{2m_+}}+\delta _0\mu _+\right]\psi _+,`$
$`{\displaystyle \frac{_I}{\mathrm{}}}`$ $`=`$ $`\frac{1}{2}𝒦\left(\psi _+^{}\psi _{}\varphi +\varphi ^{}\psi _{}^{}\psi _+\right).`$ (4)
Here $`m_B`$ ($`m_{}`$) is the mass of a bosonic (fermionic) atom, $`m_+=m_B+m_{}`$ is the mass of a molecule, $`\mu _{B(\pm )}`$ is the so-called chemical potential that implicitly accounts for particle trapping, $`\delta _0`$ is the binding energy (detuning) of the molecular state relative to the dissociation threshold, and $`𝒦`$ is the (real) atom-molecule coupling.

We now make some simplifications to allow for ease of modeling. First, we scale the Fermi (Bose) atom fields as $`\psi _{}e^{i\mu _{}t}\psi _{}`$ ($`\varphi e^{i\mu _Bt}\varphi `$), and the molecule field as $`\psi _+e^{i(\mu _{}+\mu _B)t}\psi _+`$, which yields $`\mu _+\mu _+^{}=\mu _+(\mu _{}+\mu _B)`$; in turn, $`\mu _+^{}`$ is absorbed into the detuning. Second is that atom-molecule conversion occurs on a timescale much faster than the motion of the atoms in the trap, allowing us to neglect the kinetic energies, and justifying our omission of an explicit trapping potential for the particles. Third is that, whatever the Fermi energy associated with the fermionic components, it lies within the Wigner threshold regime, so that the coupling $`𝒦`$ can be taken as the same for all modes. Finally, we switch to momentum space, but retain only the $`𝐤=0`$ atomic condensate modes since, due to Bose stimulation, these transitions are favored over $`𝐤0`$ modes. The simplified Hamiltonian reads

$$H=\underset{𝐤}{}\left[\delta c_𝐤^{}c_𝐤\frac{1}{2}\kappa (c_𝐤^{}b_𝐤a+a^{}b_𝐤^{}c_𝐤)\right],$$
(5)
where $`a`$ annihilates a bosonic atom with wavevector $`𝐤=0`$, $`b_𝐤`$ ($`c_𝐤`$) annihilates an atom (molecule) with wavevector $`𝐤`$, and $`\kappa =𝒦/\sqrt{V}`$.

The system described by the Hamiltonian (5) is a particular nonlinear version of a two-level atom driven by an external field. Of course, the two-level atom is well known to undergo Rabi flopping , and similar behavior observed for the present system would in principle answer question of whether collective association will provide a shortcut to the formation of a Fermi-degenerate gas of molecules. We therefore consider an on-resonance system, for which $`\delta =0`$ and the Heisenberg equations of motion are

$`i\dot{a}=\kappa {\displaystyle \underset{𝐤}{}}b_𝐤^{}c_𝐤,`$ (7)
$`i\dot{b}_𝐤=\kappa a^{}c_𝐤,`$ (8)
$`i\dot{c}_𝐤=\kappa b_𝐤a.`$ (9)
If the system starts out with $`N_B`$ bosonic atoms, the operator $`a`$ will have a characteristic size $`\sqrt{N_B}`$. On the other hand, while the fermionic operators $`b_𝐤`$ and $`c_𝐤`$ are characteristically of order unity, the sum over momentum states will ultimately introduce the initial number of fermionic atoms. Hence, for equal numbers of bosons and fermions, $`N_B=N_F=N`$, it is intuitively obvious that the equations of motion (Shortcut to a Fermi-Degenerate Gas of Molecules via Cooperative Association) yield a system that evolves with the characteristic frequency $`\mathrm{\Omega }=\sqrt{N}\kappa `$.

Using Fock states, said intuition is confirmed by solving the time-dependent Schrödinger equation numerically, the results of which are shown in Fig. 1. For the simplest case of $`N=1`$, the system undergoes complete conversion from a doubly-degenerate Fermi-Bose system of atoms to a Fermi-degenerate gas of molecules in a time $`\mathrm{\Omega }^1`$. However, considering a larger number of initial atoms, $`N=5`$, we see that quantum many-body fluctuations not only frustrate complete conversion, but also adjust the oscillation frequency and lead to collapse and revival. Increasing the initial particle number to $`N=10`$ brings the amplitude of the first half-oscillation closer to the ideal case. This behavior is exactly analogous to the single-component bosonic case . Although limited computational resources preclude the explicit investigation of large particle number, based on the bosonic analogy  we fully expect the first half oscillation to be complete for large $`N`$, i.e., $`N1`$ initial Fermi atoms should be converted to $`N`$ Fermi molecules over a timescale $`\tau \mathrm{\Omega }^1`$.

Nevertheless, collisions between particles will shift the energy of a state with a given number of atoms and molecules in a manner that depends nonlinearly on the numbers, making it difficult to maintain exact resonance. Furthermore, once a the system has Rabi flopped to molecules, the depopulated atomic condensate makes Bose stimulation, and thus neglect of the $`𝐤0`$ BEC modes, questionable. And so we investigate the more robust possibility of rapid adiabatic passage . From the Hamiltonian (5), it should be clear that the system will favor all atoms for large positive detunings, while favoring all molecules for large negative detunings. With $`\mathrm{\Omega }`$ established as the characteristic frequency for collective atom-molecule conversion, changes in the detuning that are of the order of $`\mathrm{\Omega }`$, and occur over a time of order $`\mathrm{\Omega }^1`$, should qualify as adiabatic. Hence, if the detuning is swept in a linear fashion according to $`\delta =\xi \mathrm{\Omega }^2t`$, then dimensionless sweep rates $`\xi \stackrel{<}{}1`$ should enable rapid adiabatic passage to a degenerate gas of Fermi molecules.

Our suspicions are again corroborated by a Fock-state-based numerical solution to the Schrödinger equation, shown in Fig. 2. While many-body effects appear to rather seriously affect the efficiency of a marginally adiabatic sweep ($`\xi =1`$) compared to the $`N=1`$ case, the difference between $`N=5`$ and $`N=10`$ (not shown) is in fact small. Overall, many-body effects are expected to be weak for near-adiabatic sweeps ($`\xi 1`$), and vanishing for sweeps that are deep-adiabatic ($`\xi \stackrel{<}{}0.1`$), again in agreement with our BEC experience . On the matter of photodissociation to noncondensate modes–and the related pair correlations–we note that, in the all-boson case , such transitions can be neglected for a sweep directed as in Fig. 2, i.e., for $`\dot{\delta }<0`$. Given the success of the analogy so far, similar results are expected for a Fermi system.

Before closing, we estimate some explicit numbers. We eschew photoassociation because of the losses associated with the electronically-excited state , and focus on the atom-molecule coupling provided by the $`{}_{}{}^{87}\mathrm{Rb}`$-40K Feshbach resonance located at $`B_0=534`$, which has a width $`\mathrm{\Delta }_R=4`$G and an associated zero-field Fermi-Bose-atom scattering length $`a_{aa}^{FB}=17.8`$nm. Accordingly, the atom-molecule coupling is $`𝒦(4\pi |a_{aa}^{FB}|\mu _{\mathrm{Bohr}}\mathrm{\Delta }_R/m_r)^{1/2}=0.14\mathrm{cm}^{3/2}\mathrm{s}^1`$, where we have estimated the difference between the Fermi-Bose atom pair and molecular magnetic moments to be equal to the Bohr magneton $`\mu _{\mathrm{Bohr}}`$. Assuming $`N_B=10^5`$ condensate atoms in a trap with respective radial and axial frequencies $`\omega _r/2\pi =215`$Hz and $`\omega _a/2\pi =16.3`$Hz , the density of bosons is $`\rho _B=8.1\times 10^{13}\mathrm{cm}^3`$. As for the fermions, we assume a modest number, say, $`N_F=10^3`$, which has three consequences: (i) the atomic BEC will act as a reservoir, thus absorbing any heat created by holes in the Fermi sea ; (ii) barring an unfortunately large scattering length for Bose-atom and Fermi-molecule collisions, we can neglect the possibility of any Fermi-Bose collapse instabilities ; (iii) the size of the Fermi cloud ($`R_F=8.3\mu \mathrm{m}`$) is smaller than the BEC ($`R_B=10\mu \mathrm{m}`$), so that overlap is not an issue. Moreover, for $`N_BN_F`$, the timescale for atom-molecule conversion is $`\tau _{a2m}(\sqrt{\rho _B}𝒦)^1=8.2\times 10^7\mathrm{s}`$. This timescale is safely below the fastest timescale for trapped-atom motion $`\tau _t=(\omega _r/2\pi )^1=4.7\times 10^3\mathrm{s}`$, justifying our neglect of trap dynamics and the kinetic energy; physically put, this means that the Fermi energy is negligible compared to the atom-molecule coupling strength.

Lastly we discuss the neglected role of particle-particle interactions, i.e., collisions. These are described in terms of a single parameter, the $`s`$-wave scattering length $`a`$, where the value of $`a`$ differs depending on whether one has in mind Bose-Bose atom, Fermi-Bose atom, Fermi-Bose atom-molecule, or Fermi-Fermi atom-molecule collisions ($`s`$-wave Fermi-Fermi atom and molecule collisions are of course Pauli blocked). Whatever the species involved, the real part of $`a`$ describes elastic collisions between particles, and leads to so-called mean-field shifts. The imaginary part need only be considered for atom-molecule collisions terms, and arises because the molecules are vibrationally very hot, and therefore undergo vibrational relaxation induced by collisions with unconverted atoms. The collisional coupling strength is then $`\lambda =2\pi \mathrm{}a/m^{}`$, where $`m^{}`$ is the mass of the atom or the reduced mass of the atom-atom (atom-molecule) pair, and it is included in the theory as a density dependent detuning. The units of the collisional coupling are $`[\lambda ]=\mathrm{cm}^3/\mathrm{s}`$, so that an appropriate density factor is needed to arrive at a collisional rate. Estimating the real part of $`a`$ from the $`{}_{}{}^{87}\mathrm{Rb}`$-40K collisions, and the imaginary part from collisions , a typical density ($`\rho 10^{14}`$) leads to elastic and inelastic collisions that occur on a timescale $`\tau _c(\rho \lambda )^1\stackrel{>}{}1\mathrm{ms}`$. For rapid adiabatic passage, the contribution from the collisional interaction is negligible compared to the detuning, except near-resonance; however, the system only spends about $`\tau _{a2m}1\mu \mathrm{s}`$ in this region, which is short enough to expect collisions to be negligible. If, by chance, collisions do turn out to be a problem, the trap could be expanded to further reduce the rate of collisions with unconverted atoms.

In conclusion, we highlight that the term $`\mathrm{\Omega }=\sqrt{N}\kappa `$ was previously referred to as the Bose-enhanced free-bound coupling, and the detuning sweep, for example, referred to as Bose-stimulated rapid adiabatic passage from atoms to molecules. However, this behavior is now clearly independent of statistics, so that Bose stimulation of free-bound association has nothing whatsoever to do with Bose statistics, but is instead a many-body cooperative effect that applies equally well to Fermi-degenerate systems. We expect this situation to arise whenever the system is addressed as a unit, i.e., when the atom-molecule coupling strength is larger than the Fermi energy. To four-wave mixing in a Fermi gas of atoms , we therefore add a further example of fermions mimicking an effect that was previously attributed to Bose-Einstein statistics. Such behavior should allow for a shortcut to a degenerate gas of Fermi molecules via a collective photoassociation or Feshbach resonance.

The authors acknowledge Juha Javanainen for helpful discussions. This work supported by the Academy of Finland (MM and KAS, project 50314), and the Magnus Ehrnrooth foundation (OD).


text:
# Bibliography

April 8, 2003

Contact Order Dependent Protein Folding Rates:

Kinetic Consequences of a Cooperative Interplay

Between Favorable Nonlocal Interactions and

Local Conformational Preferences

Hüseyin KAYA and Hue Sun CHAN

Protein Engineering Network of Centres of Excellence (PENCE),

Department of Biochemistry, and Department of Medical Genetics & Microbiology, Faculty of Medicine, University of Toronto, Toronto, Ontario M5S 1A8, Canada

Running title: Physics of Contact-Order Dependent Protein Folding

Key words: calorimetry / chevron plot / Gō models / simple two-state kinetics /
single-domain proteins / nonadditivity

Corresponding author.
E-mail address of Hue Sun CHAN: chan@arrhenius.med.toronto.edu
Tel: (416)978-2697; Fax: (416)978-8548
Mailing address: Department of Biochemistry, University of Toronto, Medical Sciences Building – 5th Fl., 1 King’s College Circle, Toronto, Ontario M5S 1A8, Canada.

Abstract

Physical mechanisms underlying the empirical correlation between relative contact order (CO) and folding rate among naturally-occurring small single-domain proteins are investigated by evaluating postulated interaction schemes for a set of three-dimensional 27mer lattice protein models with 97 different CO values. Many-body interactions are constructed such that contact energies become more favorable when short chain segments sequentially adjacent to the contacting residues adopt native-like conformations. At a given interaction strength, this scheme leads to folding rates that are logarithmically well correlated with CO (correlation coefficient $`r=0.914`$) and span more than 2.5 orders of magnitude, whereas folding rates of the corresponding Gō models with additive contact energies have much less logarithmic correlation with CO and span only approximately one order of magnitude. The present protein chain models also exhibit calorimetric cooperativity and linear chevron plots similar to that observed experimentally for proteins with apparent simple two-state folding/unfolding kinetics. Thus, our findings suggest that CO-dependent folding rates of real proteins may arise partly from a significant positive coupling between nonlocal contact favorabilities and local conformational preferences.

INTRODUCTION

Generic protein properties as energetic constraints

The folding of many small single-domain proteins is well approximated by simple two-state thermodynamics and kinetics.1,2 In the past several years, we have shown that fundamental insights into protein energetics can be gained by using these general, apparently mundane properties as experimental constraints on protein chain models.3-10 This approach is based on the recognition that model interaction schemes capable of producing these commonly observed experimental properties are, somewhat surprisingly, not entirely straightforward to come up with. To date, much advance has been made by coarse-grained modeling of protein folding.7,11-15 Nonetheless, the interactions postulated by many existing models are insufficient for calorimetric two-state cooperativity.3,4 Furthermore, even common Gō models are not cooperative enough for simple two-state kinetics, their explicit native biases notwithstanding. Specifically, we recently found that several lattice6,9,10 and continuum (off-lattice)8 Gō-like formulations with essentially additive interaction schemes all led to chevron rollovers — a hallmark of folding kinetics that are often operationally referred to as non-two-state.9 Apparently, many-body interactions are needed to produce chevron plots with linear folding and unfolding arms consistent with a two-state description of equilibrium thermodynamics.10

Small single-domain proteins are characterized as well by a significant correlation between relative contact order (CO) and folding rate.16 Therefore, it is only logical to require a model protein interaction scheme to produce a similar correlation.17,18 Ising-like19,20 and other21,22 constructs without explicit chain representations have had successes in this regard. However, as for thermodynamic and kinetic cooperativities, achieving the CO dependence requirement in models with explicit chain representations appears to be a nontrivial task. Notably, an early lattice model study using a 20-letter alphabet suggested that proteins with higher CO should fold faster,23 thus predicting a trend opposite17 to that for real single-domain proteins.16,18 A more recent 20-letter lattice model investigation, on the other hand, found modest correlations between increasing CO and longer logarithmic folding time (correlation coefficient $`r0.70`$$`0.79`$ for chain lengths $`54`$).24 An earlier continuum Gō model studies of 18 proteins also found a modest correlation between increasing CO and slower logarithmic folding rates ($`r=0.69`$).25 But the corresponding dispersion in simulated folding rates covers only $`1.5`$ orders of magnitude, which is much narrower than the $`5`$ orders of magnitude covered by the real folding rates of the proteins in the given dataset. When a different potential function was used in a more recent continuum Gō model analysis, however, no correlation between CO and simulated folding rates was discerned.26

Recently, based on lattice 27mer simulations, Jewett et al.27 have proposed that enhanced thermodynamic cooperativity and many-body interactions — which are basic properties of individual two-state proteins to begin with1-10 — may also be a key to understand the correlation between CO and folding rate across different proteins. This is an attractive and insightful idea. However, the particular way in which thermodynamic cooperativity was enhanced by these authors led only to modest increases in folding rate dispersion relative to that for the corresponding lattice Gō models with pairwise additive contact energies. Both the dispersion in folding rates and the correlation of logarithmic folding rate with CO ($`r=0.75`$) for the most cooperative interaction scheme they reported were similar to that obtained from an earlier continuum Gō model study,25 as well as that from a recent simulation of 20-letter lattice models24 with only pairwise additive contact energies (see above). In our view, these results suggest that while CO-dependent folding may well derive from certain intraprotein interactions that are also responsible for high thermodynamic cooperativity, CO-dependent folding does not arise from thermodynamic cooperativity per se. In other words, how cooperativity is achieved can be critically important. Many a priori many-body mechanisms are consistent with high thermodynamic cooperativity. An example is the two rather different interaction schemes we considered in ref. 10 — one involves local-nonlocal coupling while the other assigns an extra favorable energy to the ground-state structure as a whole. But perhaps not all such mechanisms can mimic experimentally observed CO dependencies to the same degree. Therefore, to shed light on the physical mechanisms of CO-dependent folding, we endeavor to construct an interaction scheme that would provide larger dispersions in folding rates and better correlations with CO.

MODELS AND METHODS

The present study focuses on the idea of a cooperative interplay between local conformational preferences and the contact-like interactions that drive the packing of the protein core.3,5,6,10 We have shown that chain models embodying this idea can lead to calorimetric cooperativity and simple two-state kinetics,10 although our exploration thus far has been limited to model proteins that are mostly helical.3,5,6,10 Here we consider a general formulation of this idea, the basic ingredients of which are described by Fig. 1A. This hypothesis may be viewed as a synthesis of the local-dominant and the nonlocal-dominant perspectives.28 We were motivated by the recognition that both local29,30 and nonlocal31,32 intraprotein interactions are important determinants of protein structure and stability. Yet local conformational preferences alone are often insufficient for stable secondary structures under physiological conditions. Secondary structure formation is known to be context dependent;33 they are stable when packed in the core of a protein but are usually not stable in isolation (ref. 31 and references therein). Furthermore, conformational space grows exponentially with chain length, even when preferences arising from local excluded volume effects are taken into account.34 It follows that a large part of the stability and uniqueness of protein native structures cannot be explained by local interactions alone.35 On the other hand, our recent Gō-model studies have shown that nonlocal contact-like interactions by themselves are not cooperative enough for simple two-state kinetics6,8-10 if they are not coupled to local conformational propensities.

A simple model of local-nonlocal coupling

Here we explore the hypothesis in Fig. 1A by incorporating its form of local-nonlocal coupling into a new interaction scheme in Fig. 1B for explicit-chain models configured on three-dimensional simple cubic lattices. This allows the idea to be tested quantitatively. Fig. 1B may be viewed as a generalization of similar constructs we have used previously in the context of helical proteins.3,5,6,10 As a first step in our inquiry, we make the simplifying assumption that the interactions are native-centric,25-27,31,36-38 in that only native interactions are favored, while nonnative interactions are neutral (have zero energy). The local-nonlocal coupling in Fig. 1B involves nonadditive many-body interactions. A chain segment which is locally nativelike (with native bond and torsion angles) but make no native contact is not stabilized (contributing zero energy). On the other hand, nonlocal contact interactions between monomers far apart along the chain sequence are more favorable when the chain segments around the contacting residues are in their native conformations than when they are not. As such, the present model differs from models that additively combine contact energies and local favorabilities.39 The importance of nonadditive many-body effects in protein folding has been recognized,3,5,6,10,40-44 but they have not been used extensively to model calorimetric two-state cooperativity and linear chevron plots.3-10 Our aim here is to utilize extremely coarse-grained representations as a computationally efficient means to explore the general principles linking CO-dependent folding and proteinlike cooperativities. Many structural and energetic details of real proteins are beyond the scope of this work. In particular, the present work does not deal with the microscopic physical origins of local-nonlocal coupling. Instead we just presume that its presence in naturally occurring proteins could arise from evolutionary design. Because of these, the simple interaction scheme in Fig. 1B should be viewed only as a tentative model in this regard.

In order to examine the folding rates of a set of model proteins whose native structures cover a diverse range of CO values, we now consider chains of length $`n=27`$ configured on simple cubic lattices. For these 27mers, there are 103,346 distinct maximally compact conformations (not related by rotations or inversions)45,46 confined to a $`3\times 3\times 3`$ cube. The distribution of CO among these maximally compact conformations covers 97 different values27 from CO $`=`$ $`208/756=0.275`$ to $`402/756=0.532`$ (inset of Fig. 2A, where CO is computed using equation 1 of ref. 16). For each CO value, we randomly choose a maximally compact 27mer conformation as the native structure of a model protein (Table I).***Since the present choices of structures are independent of that by Jewett et al.,27 the structures listed in Table I do not necessarily coincide with those used in their study.

Folding and unfolding kinetics are modeled by standard Monte Carlo simulations using the Metropolis criterion and the elementary chain moves of end flips, corner flips, crankshafts, and rigid rotations. The relative frequencies of attempting these moves are 4.7%, 58.3%, 27%, and 10% respectively (c.f. ref. 6)The following typographical error in ref. 6 should be corrected. The relative attempt frequencies of corner flips and crankshafts used in this prior study of ours were, respectively, 60.6% and 27%, not the 27% and 60.6% stated on p. 901 of ref. 6. Time is measured by the number of attempted Monte Carlo moves for a given process. The set of elementary chain moves is chosen to mimic physically plausible processes. Lattice model kinetics are dependent on the choice of move set.12 Nonetheless, we expect the general trend predicted by the model is less sensitive to move set when kinetics are not dominated by trapping events,12 as is the case here and has been verified by Jewett et al.27 Progress towards the native state is tracked by the fractional number of native contacts $`Q`$ (ref. 3–6). To ascertain the implications of the local-nonlocal coupling we proposed, results from a highly cooperative interaction scheme with $`a=0.1`$ are compared with that from the additive scheme ($`a=1`$) of common Gō models (c.f. Fig. 1B). Folding trajectories are initiated at a randomly generated conformation; folding first passage time is defined by the formation of the $`Q=1`$ ground-state conformation. Unfolding trajectories are initiated at the ground-state conformation; unfolding first passage time is the time it takes for the chain to be left with three or fewer native contacts ($`Q3/28`$); $`Q=3/28`$ is chosen to define unfolding because it coorresponds approximately to the free energy minimum for the denatured state.

RESULTS

Sensitivity of folding rate on CO enhanced by local-nonlocal coupling

Fig. 2 provides the correlation between CO and folding rate among our 27mer models. It shows clearly that the local-nonlocal coupling mechanism postulated in Fig. 1 can lead to a significant enhancement of correlation as well as much increased sensitivity of folding rate to CO. Whereas the dispersion in folding rates among the common additive Gō models in Fig. 2A covers only approximately one order of magnitude (a factor of ten) and the logarithmic folding rates exhibit only a relatively weak correlation with CO (correlation coefficient $`r=0.63`$), the corresponding dispersion among the $`a=0.1`$ cooperative models in Fig. 2B covers approximately 2.5 to 3 orders of magnitude, with a strong correlation between CO and logarithmic folding rate ($`r=0.914`$) comparable to that observed among a selection of real, small single-domain proteins.18 Similar to the corresponding experimental situations,16,18 the comparisons in Fig. 2 were performed under conditions for which folding relaxation is essentially single-exponential, as is evident from the good agreements in Fig. 2 between median first passage time divided by $`\mathrm{ln}2`$ and the corresponding mean first passage time.6,47 To better delineate the effects of having weakened contact interactions when the chain segments locally adjacent to the contacting residues are nonnative, several $`a`$ values other than the $`a=0.1`$ used for the main plot are compared in the inset of Fig. 2B. It shows CO-dependent folding at different levels of local-nonlocal coupling (different $`a`$ values) for several 27mers with representative CO’s. The $`a=0`$ case here corresponds to complete interdependence between nonlocal contact and local structure. This inset indicates that sensitivity of folding rate to CO increases (the fitted line has a more negative slope) with decreasing $`a`$, and that the behavior of the $`a=0.1`$ models is very similar to that of the $`a=0`$ models. These results further affirm that local-nonlocal coupling is a key ingredient for the good correlation between CO and fold rate in these models. Nevertheless, as for real proteins,16,18 despite the good correlation, CO by itself cannot predict folding rates of the present models with high accuracy. Folding rates here can vary significantly for different structures with the same CO as well. For example, for the particular 27mer with CO $`=346/756=0.458`$ in Fig. 2B, the datapoint $`\mathrm{log}_{10}(\mathrm{folding}\mathrm{rate})=5.75`$ may be viewed as an “outlier” vis-à-vis the fitted line. However, for two other 27mers with the same CO but do not belong to the randomly chosen set in Table I (and therefore not plotted and not used in the correlation analysis of Fig. 2B), we found $`\mathrm{log}_{10}(\mathrm{folding}\mathrm{rate})=7.26`$ and $`7.60`$, which happen to be much closer to the fitted line in Fig. 2B. The reasons behind variations in folding rates among structures with same CO remain to be elucidated.

A consistent model of thermodynamic and kinetic cooperativity

Fig. 3 provides further analyses of the folding/unfolding kinetics of one example 27mer structure we choose to study in more detail. Consistent with our previous results,6,8-10 it shows that the model chevron plot48 predicted by the common additive Gō potential (upper plot) deviates significantly from simple two-state kinetics in that it exhibits a severe rollover under only moderately native conditions. More specifically, for this case rollover becomes significant at $`/k_\mathrm{B}T`$ values that are only slightly more negative (more favorable to folding) than that of the transition midpoint ($`/k_\mathrm{B}T1.43`$). In contrast, the chevron plot predicted by the model with a substantial local-nonlocal coupling (lower plot) is qualitatively similar to that of real, small single-domain proteins that fold and unfold with simple two-state kinetics.10 In particular, it has essentially linear folding and unfolding arms over an extended range of $`/k_\mathrm{B}T`$ values. We have also obtained for this model the equilibrium free energy of unfolding $`\mathrm{\Delta }G_\mathrm{u}`$ as a function of $`/k_\mathrm{B}T`$, where $`\mathrm{\Delta }G_\mathrm{u}`$ here is taken to be that between the unique $`Q=1`$ conformation and those with $`Q3/28`$. (The same definition is used for unfolding kinetics as stated above.) Because $`\mathrm{\Delta }G_\mathrm{u}`$ is essentially linear in $`/k_\mathrm{B}T`$, the linearity of the chevron arms over an extended $`/k_\mathrm{B}T`$ range implies an essentially linear relationship between folding/unfolding rates and $`\mathrm{\Delta }G_\mathrm{u}`$ within the corresponding regime (i.e., the model parameter $``$ may be eliminated in favor of the lower horizontal scale in Fig. 3). Furthermore, comparing the mean first passage times in Fig. 3 versus the corresponding median first passage times divided by $`\mathrm{ln}2`$ shows that folding or unfolding relaxation for this model is essentially single exponential6,47 for $`\mathrm{\Delta }G_\mathrm{u}<`$ $`10k_\mathrm{B}T`$. Essentially single-exponential folding under moderately folding conditions is further demonstrated by an approximately linear logarithmic distribution of first passage time8,9,49 shown in the inset. Similar to the cooperative models we recently investigated,10 for the model with local-nonlocal coupling in Fig. 3, the thermodynamic $`\mathrm{\Delta }G_\mathrm{u}`$ values matches well with the kinetically obtained quantity $`k_\mathrm{B}T\mathrm{ln}[(\mathrm{folding}\mathrm{rate})/(\mathrm{unfolding}\mathrm{rate})]`$ for $`\mathrm{\Delta }G_\mathrm{u}`$ ranging from $`10k_\mathrm{B}T`$ to $`6k_\mathrm{B}T`$ (lower V-shape). In other words, the folding/unfolding kinetics of this model is simple two-state6,8-10 within a $`\mathrm{\Delta }G_\mathrm{u}`$ range quite similar to that experimentally accessible to small single-domain proteins.10 Finally, the cooperative model in Fig. 3 is also calorimetrically two-state. Assuming that the interactions are temperature independent, the model’s van’t Hoff to calorimetric enthalpy ratio $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}`$ ($`\kappa _2`$ without baseline subtraction4) is determined to be $`0.992`$ (detailed calculation not shown), satisfying the requirement of $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}1`$ for two-state thermodynamics.3-5 Taken together, the above considerations imply that the local-nonlocal coupling mechanism for enhanced CO-dependent folding in Fig. 2B also provides — as it should — a consistent account of thermodynamic and kinetic cooperativities6,8-10 in simple two-state proteins (Fig. 3).

As it stands, the transition midpoints of all 27mers considered here with the local-nonlocal coupling parametrized by $`a=0.1`$ are very close to one another. This is because the interaction scheme in Fig. 1B assigns the same energy ($`=28`$) to every ground-state conformation. This is a simplifying assumption in the present modeling setup. Since the thermodynamic stabilities of real, small single-domain proteins are quite diverse,16,18 it is important to note that, in a broader perspective, our hypothesis that significant CO-dependent folding can emerge from local-nonlocal coupling is not contingent upon the different proteins in question having very similar thermodynamic stabilities. In more sophisticated models, for example, an extra favorable energy that differs from one 27mer to another may be assigned to the ground-state conformation (i.e., a different $`E_{\mathrm{gs}}`$ term as defined in ref. 10 for each 27mer). In that case, the thermodynamic stabilities of different 27mers can be very different, but their folding rates would not be affected by this extra feature of the model. In other words, the correlation between CO and folding rate in Fig. 2B would remain unchanged. As we have recently argued,10 such extra stabilizing energies for the ground state as a whole are physical plausible because experimental evidence50 indicates that in real proteins there is a partial separation between the driving forces for folding kinetics and the interactions responsible for thermodynamic stability.

DISCUSSION

Energy landscapes of the present models are further characterized in Fig. 4 for three representative structures with low, intermediate, and high CO values. In this figure, the low- and high-CO structures are, respectively, the fastest and slowest folding among the 97 structures in Table I, whereas the intermediate-CO structure is the one analyzed in Fig. 3. For the common additive Gō potential, energy $`E`$ is directly proportional to $`Q`$ ($`E=Q`$). However, for the cooperative models with local-nonlocal coupling, there are multiple energy levels for each $`Q`$, with $`E=Q`$ as the lower bound (left panels of Fig. 4). This means that, on average, the energetic separations between non-ground-state and ground-state conformations in the cooperative models with local-nonlocal coupling are larger than that in the additive Gō models. This feature is demonstrated directly in the right panels of Fig. 4, which show that the number of non-ground-state conformations within a given energy range is smaller for the cooperative models than for the additive Gō models except for the highest energies ($`E0`$). It follows that the overall thermodynamic cooperativities of the models with local-nonlocal coupling are substantially higher than that of the corresponding additive Gō models. This behavior is expected as well from our recent finding that simple two-state folding/unfolding kinetics (Fig. 3 above) requires “near-Levinthal” thermodynamic cooperativity.10 Indeed, for the three models in Fig. 4 with local-nonlocal coupling, the van’t Hoff to calorimetric enthalpy ratios $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}`$ are, from top to bottom, $`\kappa _2=`$ $`0.972`$, $`0.992`$, and $`0.998`$. These values are extremely high for model enthalpy ratios without baseline subtractions.4 In contrast, the corresponding additive Gō models are less cooperative, with $`\kappa _2=`$ $`0.751`$, $`0.861`$, and $`0.878`$. Here it is noteworthy that the additive Gō models’ $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}`$ ratios even after empirical baseline subtractions,4 $`\kappa _2^{(\mathrm{s})}=`$ $`0.885`$, $`0.961`$, and $`0.962`$, are lower than the $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}`$ ratios of the cooperative models in the absence of baseline subtractions.

Contact-order dependence indicative of special mechanisms of cooperativity

Obviously, thermodynamic cooperativity is a necessary ingredient for any protein chain model that purports to rationalize the generic properties of small single-domain proteins.3-10 For the particular interaction scheme we consider, the above analysis shows that features that give rise to significant CO-dependent folding also lead to high thermodynamic cooperativity. However, the converse is not necessarily true. More in-depth considerations and a comparison of the present results with that of Jewett et al.27 indicate that higher thermodynamic cooperativity per se does not necessarily give rise to more enhanced dependence of folding rate on CO. Our reasoning is as follows. First, for the present set of 27mer structures we have chosen randomly, the correlation between logarithmic folding rate and CO is quantified by $`r=0.63`$ ($`r^2=0.39`$) for the additive Gō interaction scheme. Despite that this correlation happens to be weaker than that of Jewett et al.’s collection of additive Gō models (their $`r^2=0.51`$), after cooperativity is enhanced by local-nonlocal coupling, the correlation between logarithmic folding rate and CO for our $`a=0.1`$ models is much higher ($`r^2=0.84`$, see Fig. 2 above, an improvement in $`r^2`$ value of $`0.33`$) Because all the model chains in the present study have the same length and the same number of native contacts, their correlation coefficient between folding rate and CO is the same as that between folding rate and the total contact distance (TCD) defined in ref. 51.  than the best case reported by Jewett et al.27 ($`r^2=0.57`$ for their $`s=3`$, an improvement in $`r^2`$ value of $`0.06`$ over that for their additive Gō models).§§§ If the $`s=3`$ interaction scheme of Jewett et al. is applied to the present set of structures and kinetic models, we found $`r^2=0.65`$ for the correlation between CO and folding rate. In this case, the folding rates span $`1.8`$ orders of magnitude; see ref. 52 for details.  Second, the folding rates of our cooperative models are much more sensitive to CO, covering 2.5 to 3 orders of magnitude, whereas those of Jewett et al. cover only approximately 1.3 orders of magnitude. This means that the present local-nonlocal coupling mechanism is significantly more effective in enhancing CO dependence than the nonlinear $`E`$$`Q`$ relationship postulated by Jewett et al. (equation 1 of ref. 27). Their interaction scheme does not make direct reference to chain conformations as such. Thermodynamic cooperativity is enhanced in their models by stipulating that the total contact energy $`E`$ (for a given conformation as a whole) does not decrease (does not become more favorable) linearly with increasing $`Q`$ as in common Gō models; but rather decreases at progressively faster and faster rates when $`Q`$ is closer to unity. Jewett et al. suggested that the “extraordinary cooperativity in protein folding” may originate from “three-body interactions.” But how three-body interactions might lead to their $`E`$$`Q`$ relationship remains to be elucidated. Third, in fact, if thermodynamic cooperativity is further increased in the interaction scheme of Jewett et al. by increasing their $`s`$ parameter, the energy landscape will eventually become a Levinthal golf course in the $`s\mathrm{}`$ limit. In that case, folding would be rate-limited by random conformational search and CO-dependence would be all but eliminated. Fourth, in this connection, we have recently considered three 27mer models with CO $`=0.28`$, $`0.40`$ and $`0.51`$ in a separate study. The thermodynamic cooperativity of these models are enhanced by assigning an extra stabilizing energy to the ground state but without local-nonlocal coupling.10 For the energetic parameters we considered, the folding rates of these models cover less than an order of magnitude.10 The same set of results also indicated that dispersion in folding rates under moderately folding conditions would decrease if thermodynamic cooperativity is increased by assigning an even stronger stabilizing energy to the ground state, in a manner similar to greatly increasing $`s`$ in Jewett et al.’s formulation. Taken together, these observations lead us to the conclusion that while thermodynamic cooperativity is certainly necessary, by itself it is not sufficient to guarantee CO-dependent folding rates similar to that observed experimentally16,18 if the underlying mechanism for thermodynamic cooperativity is not specified.

CO-dependent folding highlights the important role of local interactions in determining folding rates.16-18 It suggests that the mechanism of folding may involve relatively fast formation of local structure. In this regard, we note that under the general lattice scheme in Fig. 1B, formation of strong (unattenuated) native contacts with contact order $`|ji|=3`$ is relatively easier than formation of strong native contacts with higher contact orders. This is because in the $`|ji|=3`$ case there is an overlap between parts of the two local segments that have to be nativelike in order for the contact to be strong. Physically, how a general mechanism similar to that in Fig. 1 may arise in real proteins from solvent-mediated atomic interactions such as sidechain packing and hydrogen bonding remains to be elucidated. Many basic issues will have to be tackled to address this question. For example, correlations between backbone and sidechain rotamer conformations53 may contribute to such a mechanism. Another possibility is that aspects of anti-cooperativity of certain types of hydrophobic interactions54 may help disfavor premature nonspecific hydrophobic collapse (which would lead to kinetic trapping14) when the sidechains are locally less well packed than that in the native state. If this is the case, it could give rise to local-nonlocal coupling mechanisms similar to that postulated in Fig. 1.

In summary, while the models used in the present study are rudimentary, they provide strong evidence that a cooperative interplay between local conformational preferences and nonlocal favorable contact-like interactions is an important mechanism in accounting for experimentally observed CO-dependent folding of small single-domain proteins. We are optimistic that more rigorous applications of the CO-dependence constraint as well as the thermodynamic and kinetic cooperativity requirements would help further narrow down theoretical possibilities and thus contribute to a more realistic understanding of protein energetics.

Acknowledgments. We thank Robert L. Baldwin, Alan Davidson, Teresa Head-Gordon, Michael Levitt, Vijay Pande, Kevin Plaxco, Steve Plotkin, Wes Stites and Yaoqi Zhou for helpful discussions, and Vijay Pande and Kevin Plaxco for kindly sharing their work (ref. 27) before publication. The research reported here was partially supported by the Canadian Institutes of Health Research (CIHR grant no. MOP-15323), a Premier’s Research Excellence Award from the Province of Ontario, and the Ontario Centre for Genomic Computing at the Hospital for Sick Children in Toronto. H. S. C. is a Canada Research Chair in Biochemistry.

References

Table I

| $`\mathrm{\Delta }S_{ij}`$ | conformation | $`\mathrm{\Delta }S_{ij}`$ | conformation |
| --- | --- | --- | --- |
| 208 | uufddfuurddbuubddruufddfuu | 306 | uufrrbldrfdflurullddburdbr |
| 210 | uufddfuurbbdffdbbrffuubdbu | 308 | uufdfrbrbulddrfllfrruublfl |
| 212 | ufdfuubbrddffubufrddbuubdd | 310 | ufrulblfddrrbbllfuburdrfub |
| 214 | uuffdbdfrbufubbddrffuubdbu | 312 | uufrbbllffdrrdllbubdrurfdb |
| 216 | ufdfuubbrddfuufddruubbddfu | 314 | uufdrubbdfdfllbbuffubbrddr |
| 218 | ufdfuubbrdfufddbbruuffddbu | 316 | uffrddblbruufdllbdffrulubb |
| 220 | uuffddburfdbbuuffrddbuubdd | 318 | uufrbddbuullffdrrdllbubdru |
| 222 | uuffddburdfuubbddruuffdbdf | 320 | uufddfrruubbdfdbluuffldrbd |
| 224 | uufddrbufubrfdbdfflurulldd | 322 | uffdbrbrfufullbbrdrufldfdr |
| 226 | uufddfuurddbubdrffubbulfrf | 324 | ufrbddlfrflurullbbddffubrr |
| 228 | uffdbrbuffdrbbuffubbllfrfl | 326 | ufrfddluulddbbuufdrdbrfubu |
| 230 | ufdrbufublffddruurddbuubdd | 328 | ufrubdbuldldrrffllbufurblb |
| 232 | uufrrblddrufdluldfurdruull | 330 | uufdrubblddlfubuffddrrbbuf |
| 234 | uufddfuurddrbluurfdbbulddr | 332 | uuffrddruubbdfllfdbrrbluuf |
| 236 | ufddbbuurrflfrdlbdfrbubldr | 334 | ufrrdbdfluldbbrruuflbldrfd |
| 238 | uuffdbdfrrblbrulffrulbbrfd | 336 | uffrbrbuflblffrrddllbrrblu |
| 240 | uffdbrfurdbblufrbuffllbbrf | 338 | uffurrdldrbblurullfrrdldlf |
| 242 | uufrbddfflburflurrbbdffdbb | 340 | ufrrbbdffdlbrbllffurbubldr |
| 244 | ufdrrbluulfrdrbuffllddrrul | 342 | ufrullbrrblldrrdllfufdrrbu |
| 246 | ufdfurbdfruullbrrddblurull | 344 | uffdrdllbrbluuffdbrrdbuuff |
| 248 | ufdfrbuflurblbrrdldrffuubd | 346 | uuffddrrbbuufdfuldbubddflu |
| 250 | ufddbrfruublfdbrdblluurdru | 348 | uufrbbdlulddrrffllbuufdrrb |
| 252 | uffdbrrflurbbdlufufrbbllff | 350 | ufubrrdfdfuldblfuurrbldbdr |
| 254 | ufddbrblurrdfflubrfulbrbll | 352 | uffddrbllurrfubbddlluuffdd |
| 256 | ufdfurdruullbbrddrfluurdbu | 354 | ufrfdrbufubbllfrflddbrburd |
| 258 | ufrbdflfrrbbuullfrrdfulldr | 356 | uffrrbdbuullffrrbldbdflfrr |
| 260 | uufddfurbbrdlffrbufubblffl | 358 | uufdrrubddffuulldrdlbrbuuf |
| 262 | uuffdbrbufrfldrdllbrbrfubu | 360 | ufrddllfrruulldrblubddrruu |
| 264 | ufdrurddbuuldblurrddllffrb | 362 | ufubrrdffldrbblflfuurrbldb |
| 266 | ufdrurddllbrbluurrfldbrdfu | 364 | ufrfddlbblffubbuffrdbrdbuu |
| 268 | uuffdbrubrfddbuldflfrrulur | 366 | ufrfdbdfllbbuufdfurdbdbruu |
| 270 | uuffrrdllbdrbufrulbrddffll | 368 | uffurrbbddffllbrbuulfrdrfl |
| 272 | uufrdfuldbdfrruubblddfrubd | 370 | uffurrddbbuufllbrddflfrubr |
| 274 | ufdrubrfddllbbuurrdldrfuld | 372 | uufdrfdruubbddluufflddbrru |
| 276 | ufdrbdlfrrubdblluurffrbbdl | 374 | uffrdrbbuullffrrdbuldbdflf |

$`\mathrm{}`$ Table I to be cont’d

Table I $`\mathrm{}`$ (cont’d from last page)

| 278 | ufddrrbllbrrullurrfflbdfrb | 376 | uffrddllbuubddrfrbuufdlflu |
| --- | --- | --- | --- |
| 280 | uffrddlubdruubddllfubuffdd | 378 | ufdrfdlluubbdfdbrfrbuuffld |
| 282 | ufrbdffuldlubbddrrfflbuldf | 380 | ufrfddlbrbllffubbuffrdbrbu |
| 284 | uufrfldrrubbldrfdblfuldfrr | 382 | uufdrfdrbbuufdfullddbrbuuf |
| 286 | uffubbrddrffuldlbrurbufflb | 384 | ufrbbullddfuurrfllddrrbblu |
| 288 | ufdfrrubufldlubbrfdbdfrbuu | 386 | uffddrbllfuubbddrruuffdbll |
| 290 | uufrbbldrfdbllfubuffddrurd | 388 | uffrburbddffllbrbuulffrrdb |
| 292 | uffrddbbuufdldblffrulubbdf | 390 | ufrufddrbbuffubbllffddbrbu |
| 294 | ufdrdfulurbbddlluufddrfluu | 392 | ufrufddrbbuffubbllffddbrbu |
| 296 | ufdfrbdflbbruuffllddbbuufd | 394 | ufrrddlbburuflblddffurbrdb |
| 298 | ufrrdblblurfrbddffluldbrbl | 396 | uffrddblflbufubbddrruufdlf |
| 300 | ufdfrullddbuubddrffrbuubdd | 400 | ufrfddllubdrrblluuffrdbrbu |
| 302 | ufdfurddlluubbdfdbrfrbuufd | 402 | ufrufrbbddffllbrbuulffdrrb |
| 304 | ufrdlluurrbbdfdbllfuubdruf |  |  |

Table I. The ground-state 27mer conformations ($`n=27`$) used in this investigation are given by sequences of 26 bond directions, where r = right ($`+x`$), l = left ($`x`$), f = forward ($`+y`$), b = backward ($`y`$), u = up ($`+z`$), d = down ($`z`$). A structure is randomly selected for each of the 97 possible CO values amongst the compact 27mer structures with $`t_{\mathrm{max}}=28`$ contacts. Each integer $`\mathrm{\Delta }S_{ij}`$ is the sum of $`|ji|`$ over the $`(i,j)`$ nearest-neighbor contacts in the given conformation ($`ji3`$). Here CO $`=\mathrm{\Delta }S_{ij}/(nt_{\mathrm{max}})`$ $`=\mathrm{\Delta }S_{ij}/756`$.

Figure Captions

Figure 1. (A) Schematics of local-nonlocal cooperative energetics in protein folding. The conformation in the solid box represents the native (N) structure; the two filled circles depict a pair of nonlocal residues interacting favorably in the native state. The interaction strength between a residue pair is strong and essentially the same as that in the native structure if the chain segments sequentially local to both residues are nativelike, as in (i). \[Dotted boxes in (A) are used to mark nativelike chain segments.\] However, the interaction strength is weakened if one or two chain segments sequentially local to the interacting residues are not nativelike, as in examples (ii)–(iv). (B) A lattice implementation of this protein folding scenario. Here the favorable energy for every contact (between residues $`i`$ and $`j`$, $`|ji|3`$) in the ground-state native (N) structure is $``$ ($`<0`$) when the relative positions of the five residues centered at $`i`$ (residues $`i2`$, $`i1`$, $`i`$, $`i+1`$, and $`i+2`$) as well as the relative positions of five residues centered at $`j`$ (residues $`j2`$, $`j1`$, $`j`$, $`j+1`$, and $`j+2`$) are the same as that in N \[solid lines in (i)\], irrespective of the relative orientations of the two five-residue chain segments. However, if the local conformation of one or both sets of five contiguous residues is nonnative, the contact energy is weakened by an attentuation factor $`a`$ ($`0a<1`$). Examples of the latter situation is given by (ii)–(iv), where nonnative local chain segments are drawn as broken lines.

Figure 2. Correlation between the common (base 10) logarithm of folding rate and CO for the 97 structures in Table I under moderately folding conditions at $`/k_\mathrm{B}T=1.47`$, using (A) the common additive Gō potential and (B) the local-nonlocal cooperative interaction scheme with $`a=0.1`$. Solid lines are least-square fits. Here folding rate is the reciprocal of mean folding first passage time (folding rate $`=`$ 1/MFPT). Each MFPT is averaged from 500 trajectories. Associated with each value of $`\mathrm{log}_{10}(1/\mathrm{MFPT})`$ (filled circle) is an open circle marking the common logarithm of the median folding first passage time (FPT) divided by $`\mathrm{ln}2`$. If the kinetics is single-exponential, MFPT $`=`$ (median FPT)/$`\mathrm{ln}2`$. The inset in (A) is the distribution of CO among the 103,346 maximally compact 27mer conformations, wherein the number of conformations (vertical scale) is shown as a function of CO (horizontal scale). The inset in (B) uses six representative structures with different CO values ($`\mathrm{\Delta }S_{ij}=`$ 208, 224, 268, 310, 348, and 386 entries in Table I) to illustrate that $`\mathrm{log}_{10}(\mathrm{folding}\mathrm{rate})`$ (vertical scale) is more sensitive to CO (horizontal scale) when the local-nonlocal coupling is stronger. In this inset, different symbols denote different $`a`$ values; the lines fitted through the symbol are, from top to bottom, for $`a=1`$, $`0.75`$, $`0.5`$, $`0.25`$, $`0.1`$, and $`0.0`$.

Figure 3. Model chevron plots for a CO $`=0.410`$ structure ($`\mathrm{\Delta }S_{ij}=310`$ entry in Table I) are given by negative natural logarithm of MFPT as a function of $`/k_\mathrm{B}T`$ (filled symbols). Values of (median FPT)/$`\mathrm{ln}2`$ are shown by the open symbols. Squares (folding) and triangles (unfolding) are for the additive Gō potential ($`a=1`$, upper plot), whereas circles (folding) and diamonds (unfolding) are for the $`a=0.1`$ local-nonlocal cooperative interaction scheme (lower plot). Each MFPT is averaged from 500 trajectories, except for the model with local-nonlocal coupling at $`/k_\mathrm{B}T=1.47`$ (arrow). For this particular case, 7,500 folding trajectories were simulated to provide enriched statistics for the FPT distribution in the inset, wherein $`P(t)\mathrm{\Delta }t`$ is the fraction of trajectories with $`t\mathrm{\Delta }t/2<`$ FPT $`t+\mathrm{\Delta }t/2`$, and the bin size $`\mathrm{\Delta }t`$ for FPT is equal to $`5\times 10^6`$. The free energy of unfolding $`\mathrm{\Delta }G_\mathrm{u}`$ for the $`a=0.1`$ cooperative model is computed using Monte Carlo histogram techniques based on sampling at the transition midpoint $`/k_\mathrm{B}T=1.33`$. $`\mathrm{\Delta }G_\mathrm{u}`$ is essentially linear in $``$ (lower horizontal scale). The dotted V-shape, which fits well to the kinetic datapoints of the $`a=0.1`$ cooperative model over an extended regime, is an hypothetical simple two-state chevron plot consistent with the dependence of $`\mathrm{\Delta }G_\mathrm{u}`$ on $``$.

Figure 4. Energy landscapes of three representative models with local-nonlocal coupling ($`a=0.1`$, $`\mathrm{\Delta }S_{ij}=224`$, $`310`$, and $`386`$ entries in Table I; $`=1`$). The left panels show the correlation between $`E`$ and $`Q`$; each dot indicates that at least one conformation with the given $`(E,Q)`$ was encountered in our sampling. The right panels show these structures’ logarithmic densities of states, where $`g(E)`$ is the number of conformations with energy $`E`$ for the cooperative models ($`a=0.1`$, dots). Included for comparison are the $`\mathrm{ln}g(E)`$ values of the corresponding additive Gō models ($`a=1`$, open circles; $`=1`$). The densities of states here are estimated by Monte Carlo sampling at the models’ transition midpoints $`/k_\mathrm{B}T=1.33`$ ($`a=0.1`$) and $`/k_\mathrm{B}T=1.43`$ ($`a=1`$). Note that the cooperative models have more energy levels than the additive models. Therefore, to compare their densities of states on an equal footing, the open squares provide the natural logarithm of the number of conformations in the $`a=0.1`$ cooperative models with energies in the range $`m0.5E<m+0.5`$, where $`m=1,0,1,2,\mathrm{}`$ is an integer. Now the densities of states represented by the open squares ($`a=0.1`$) are directly comparable to that represented by the open circles ($`a=1`$) because their values are based upon the same unity bin size for $`E`$.