text:
# The far-infrared energy distributions of Seyfert and starburst galaxies in the Local Universe: ISO photometry of the 12𝜇m active galaxy sample1footnote 11footnote 1Based on observations with ISO, an ESA project with instruments funded by ESA Member States (especially the PI countries: France, Germany, the Netherlands and the United Kingdom) with the participation of ISAS and NASA

## 1 INTRODUCTION

One of the main goals of selecting a complete sample of galaxies using the IRAS 12$`\mu m`$ flux was the definition of a complete and largely unbiased sample of active galaxies in the local universe. This selection was done twice, the first one producing a list of 26 Seyfert type 1 galaxies (hereinafter Seyfert 1’s) and 32 Seyfert 2’s (Spinoglio & Malkan, 1989), out of a sample of 390 galaxies from the IRAS PSC (1988). The second was a selection of 53 Seyfert 1’s, 63 Seyfert 2’s, two blazars and 38 high-luminosity non-Seyferts (i.e. galaxies with $`L_{IR}10^{44}\mathrm{erg}\mathrm{s}^1`$ without an obvious Seyfert type optical spectrum) out of a sample of 893 galaxies (Rush, Malkan & Spinoglio, 1993) (RMS) from the IRAS Faint Source Catalog (Moshir et al., 1992). It was found that the 12$`\mu m`$ flux is approximately a constant fraction ($``$ 1/5) of the bolometric flux in active galaxies. Moreover, also for non-Seyferts (mostly spirals) 12$`\mu m`$ selection was found to be the closest available approximation to selection by a limiting bolometric flux, which is about 14 times $`\nu F_\nu `$ at 12$`\mu m`$ (Spinoglio et al., 1995) (hereinafter S95). It therefore follows that deep surveys at 12$`\mu m`$ should provide complete samples at different bolometric flux levels of normal and active galaxies, which will not suffer the strong selection effects present both in the optical-UV and far-infrared. We refer to the introduction of the paper presenting the original sample (Spinoglio & Malkan, 1989) for the arguments in favor of the mid-infrared selection. After the determination by RMS, two more works derived the local 12$`\mu m`$ luminosity function (Xu et al., 1998; Fang et al., 1998) and a new selection of galaxies from the IRAS database at 12$`\mu m`$ is in progress (Alexander & Aussel, 1999).

The study of the infrared energy distributions of active and starburst galaxies in the local universe is not only important by itself, but also because it is needed to understand galaxy evolution. Deep ISO surveys have shown that indeed a strong evolution is present at z $``$ 1, indicating that many galaxies at that epoch were experiencing phases of extremely high luminosities in the far-infrared, likely to represent violent events of star formation \[see e.g. the reviews of Malkan (2000), Genzel & Cesarsky (2000), and Franceschini et al. (2001)\]. Also, the peak of the energy output of galaxies lies in the far-infrared. This makes redshifted galaxies appear relatively brighter at the longest far-infrared and sub-millimeter wavelengths; surveys in these wavebands are therefore capable of discovering very distant galaxies at substantial cosmic look-back times (see e.g. Malkan, 2001). It is thus even more urgent to understand the local infrared emission, to provide a firm base for comparing the local universe “activity” with the properties and evolution observed in the recent and upcoming cosmological surveys (e.g. Franceschini et al., 2001; Malkan, 2001).

Several papers have reported ISO photometric data on galaxies in the unexplored range from 100-200$`\mu m`$, but most of these are on few individual sources, often the prototypes of the various classes of active, starburst and normal galaxies. The only exceptions are the study of 42 objects (20 Seyfert 1’s and 22 Seyfert 2’s) (Pérez García & Rodríguez Espinosa, 2000)(PGRE) of the CfA Seyfert galaxy sample (Huchra & Burg, 1992), the study of a random sample (Haas et al., 2000) of 17 Palomar Green quasars (Schmidt & Green, 1983), the study of 22 radio-loud and radio-quiet quasars (Polletta et al., 2000) and one of 10 galaxy/quasar pairs from the 3CR catalog (Meisenheimer et al., 2001). The CfA Seyfert galaxy sample is smaller (48 objects) than the 12$`\mu m`$ active galaxy sample (158 objects), and is selected by the optical light of the galaxy, which may bias the results significantly. For example, a recent work on high resolution radio observations of the Seyfert galaxies of the 12$`\mu m`$ sample (Thean et al., 2000, 2001), report a detailed statistical analysis of the redshift distributions and 12$`\mu m`$ flux distributions of Seyfert 1’s and 2’s and compared these with the CfA sample. In the 12$`\mu m`$ active galaxy sample the two types of Seyferts are equally distributed according to redshift, their 12$`\mu m`$ flux density distributions are well matched and therefore Seyfert 1’s and 2’s are equally luminous in the IRAS 12$`\mu m`$ band. On the contrary, the CfA Seyfert sample has different redshift distributions for the two Seyfert types, and is biased against including distant Seyfert 2’s.

The three latter works on quasars (Haas et al., 2000; Polletta et al., 2000; Meisenheimer et al., 2001) are complementary to ours since most of the sources have much higher redshifts ($`0.1<\mathrm{z}<2.0`$). However, these studies are not suited to studying the population of active galaxies in the local universe.

The present work is based on observations of 90 galaxies, randomly selected from the 12$`\mu m`$ galaxy sample, including 29 Seyfert 1’s, 35 Seyfert 2’s, 12 high luminosity non-Seyferts (hereafter called starburst galaxies) and 14 normal galaxies for comparison. It is the largest sample of active galaxies for which far-infrared photometry is available out to 200$`\mu m`$. Since it contains both active and starburst galaxies it is expected to be the basis for firm and statistically significant conclusions on the far-infrared behaviour of galaxies in the local universe.

The paper is organized as follows: section 2 describes the observations collected and reduced longward of 100$`\mu m`$; section 3 shows the results in terms of average energy distributions and presents a derivation of “pure” infrared emission components, for quiescent disks, starbursts, and Seyfert nuclei; section 4 shows the color-color diagrams of the various classes of galaxies that can be explained by thermal emission from dust and presents a correlation between the far infrared color \[200 - 100\] and the 60$`\mu m`$ excess; section 5 shows the correlations among the different luminosities and between colors and luminosities. Section 6 presents our conclusions. Finally, in the Appendix we present the results obtained shortward of 100$`\mu m`$.

## 2 OBSERVATIONS

The Infrared Space Observatory (ISO) (Kessler et al., 1996) data presented here were obtained in the open time proposal “IR energy distributions and imaging of the complete sample of 12$`\mu m`$ active galaxies”. Photometric data at several wavelengths between 60 and 200 $`\mu m`$ were collected with the two array cameras, C100 and C200 of the ISOPHOT imaging photo-polarimeter (Lemke et al., 1996). In particular, we present in this section the ISOPHOT photometric results obtained with the C200 detector array in the spectral range 120-200$`\mu m`$ on 39 objects belonging to the 12$`\mu m`$ active galaxies sample. The data obtained with the C100 array camera are presented in Appendix 1, for two reasons: first, these data are less precise than the C200 data, having larger uncertainties, due either to a difficult subtraction of the background and to a less accurate calibration (see the discussion in Appendix 1); second, these data are less important, because all the galaxies of the sample have good IRAS detections at 12, 25, 60 and 100$`\mu `$m. To increase the statistics, we have also analyzed ISOPHOT C200 archive data on other 19 galaxies and we have included in the discussion ISOPHOT literature data on another 32 galaxies, bringing the present sample to 90 objects. This is therefore currently the largest sample of active galaxies in the local universe studied to 200$`\mu `$m.

The journal of the C200 ISOPHOT observations is given in Table 1, which presents the details of each single observation that we analysed: for each source are given the coordinates of the observations, the redshift, the galaxy type (hereafter we call starburst nuclei sbn the high luminosity non Seyfert galaxies), the ISO observation identification number and the total integration time spent on-source.

All of the detections are consistent with that expected for unresolved point sources. That is, none of the galaxies were significantly resolved by the large ISOPHOT PSF, which is estimated to be $``$ 100 $`\mathrm{}`$ FWHM in diameter at 120 $`\mu `$m \[e.g. see Fig.4.8 of Laureijs et al. (2001)\].

The C200 ISOPHOT observations were collected in chopped photometric mode in the five filters centered at 120, 150, 170, 180 and 200$`\mu m`$. The data reduction has been performed using the PIA222PIA is a joint development by the ESA Astrophysics Division and the ISOPHOT consortium led by the MPI fur Astronomie, Heidelberg. Contributing institutes are DIAS, RAL, AIP, MPIK and MPIA. software, version 7.2 (Gabriel et al., 1997) and the results were verified with version 8.0. Data were deglitched, corrected for non-linearity and then for each ramp a straight line was fitted. In most cases the first signal per chopper plateau was discarded since the detector response was not stabilized. Then the data were corrected for dark current and for signal dependence on the ramp integration time and calibrated with the internal Fine Calibration Source (FCS1). The same reduction procedure was applied also to the calibration data. An independent reduction of several of these observations was made by one of us (MM) at IPAC333IPAC (Infrared Processing and Analysis Center) is part of the Division of Physics, Mathematics and Astronomy at Caltech, Pasadena, and the Space and Earth Sciences Program Directorate at JPL, Pasadena., those results are in good agreement with the values published here. We realize there are faint outer parts to some of these galaxies, and this was even measured and discussed in the appendix of RMS. However, we find that only a small fraction of the total flux is actually missed due to the central concentration of the light, and that our colors should be very reliable anyway, one of the main points of this study.

The ISOPHOT flux densities have been corrected using the color correction, from the PIA (Gabriel et al., 1997), corresponding to black-body functions with an emissivity law proportional to the inverse square of the wavelength. The chosen dust temperatures were the lowest color temperature that we derive from the \[200 - 100\] color of each class of galaxies, derived from the positions in the far-infrared color-color diagrams presented in section 4. They are 20K for Seyfert 1’s and normal galaxies, 25K for Seyfert 2’s and 30K for high luminosity non-Seyferts. At these temperatures, the color correction is, however, quite small (less than 15%).

Because all the galaxies have good detections in all the IRAS bands, we have constructed their spectral energy distributions using the new ISOPHOT data together with the IRAS data given in RMS. In Fig. 1 we present the spectral energy distributions of each individual object. Table 2 presents the measured flux densities and 1 $`\sigma `$ statistical errors in the five ISOPHOT C200 wavebands for each object. Absolute calibration uncertainty is in the range 10—30%.

For most of the sources a good general agreement is found between IRAS and ISO measurements and the C200 values place themselves on a smooth curve running from the IRAS to the ISO photometry. Data for Arp 220, NGC 6240 and IZW 1 were taken from the archive. The results reported here agree well with those published in the literature for the first two sources (Klaas et al., 1997) while the C200 data of IZW 1 of Table 2 agree with those published by Haas et al. (2000) but not with those by PGRE. The measurements were taken with different observational setups, the ones reported here were gathered in chopping mode while those by PGRE in staring mode and the resulting values are much lower in this latter case. We ascribe this difference to the different background subtraction procedure and it could be due to an overestimation of the background value by this latter authors. They indeed did not detect IZw 1 at 200$`\mu `$m while we obtained a signal-to-noise ratio larger than 18 at the same wavelength. Data taken towards NGC 6090 (MK 496) reported here are lower than those published by Calzetti et al. (2000) by a factor of 1.3-1.7, however they agree within the large errors quoted by these authors.

In the following analysis we also included the ISO results published by several authors on galaxies belonging to the 12$`\mu `$m galaxy sample: 17 CfA Seyfert galaxies (PGRE), 10 normal galaxies (Siebenmorgen, Krugel & Chini, 1999), 4 nearby spiral galaxies (Alton et al., 1998) and the Seyfert 2 NGC 7582 (Radovich et al., 1999). The CfA Seyfert 1’s included are: NGC 3227, NGC 3516, NGC 4051, NGC 4151, MK 766, MK 231, NGC 5033, NGC 5548, MK 817 and NGC 7469, while the CfA Seyfert 2’s are: NGC 1143/44, NGC 3079, NGC 3982, NGC 4388, NGC 5256, NGC 5929 and NGC 1068. The normal galaxies included are: MK 323, MK 332, MK 555, MK 538, MK 799, NGC 5719, NGC 6918, NGC 7083, UGC 2936 and UGC 2982 and the nearby spiral galaxies are: NGC 134, NGC 660, NGC 5194 and NGC 5236.

Table 3 gives the integrated mid-to-far infrared luminosity from 8 to 250$`\mu `$m for all 90 galaxies, including those taken from the literature. This luminosity has been computed using the formula:

$$L_{FIR}=4\pi (\frac{cz}{H_o})^2\times \underset{i}{}(S_{\nu _i}\times \mathrm{\Delta }\nu _i)\times 10^{11}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2$$
(1)
with the $`S_{\nu _i}`$ the flux densities in Jy and $`\mathrm{\Delta }\nu _i`$ = 13.48, 5.16, 2.58, 1.0, 1.04, 1.02, 0.95, 0.65, 0.48, for $`\lambda _i`$ = 12, 25, 60, 100, 120, 150, 170, 180 and 200 $`\mu m`$, respectively444The four frequency intervals ($`\mathrm{\Delta }\nu _i`$) relative to the IRAS wave-bands have been taken from Rowan-Robinson and Crawford (1989), while the others have been computed from Table 2.7 of Laureijs et al. (2001).. We adopted a value of $`\mathrm{H}_\mathrm{o}=75\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$.

To show how the present sub-sample is representative of the 12$`\mu `$m active galaxy sample, we have computed average redshift and completeness of the sub-sample and compared with that of the original sample (RMS). The average redshift of the Seyfert 1’s in the present sub-sample is $`z=.029\pm .034`$, compared to the value of $`z=.030\pm .038`$ of the original sample; that of the Seyfert 2’s is $`z=.021\pm .018`$, which is exactly the same of the value of the original sample; that of the starburst galaxies is $`z=.034\pm .014`$, compared to the value of $`z=.030\pm .010`$ of the original sample. The completeness test was applied only to the Seyfert galaxies, because the starburt nuclei of the 12$`\mu `$m sample are not supposed to be a complete sample (see RMS): the present sub-sample contains 63% of the Seyfert 1’s and 57% of the Seyfert 2’s of the original sample down to a flux limit of 0.30 Jy at 12$`\mu `$m and 61% of the Seyfert 1’s and 64% of the Seyfert 2’s of the original sample down to 0.40 Jy. We conclude from this analysis that the present sub-sample represents fairly well the 12$`\mu `$m active galaxy sample, and there is no apparent bias.

## 3 SPECTRAL ENERGY DISTRIBUTIONS

### 3.1 Average energy distributions

Fig. 2 shows the average observed power $`\nu \times F_\nu `$ from 4400 Å to 200$`\mu m`$ of each class of galaxies (except for the very few nearby spiral galaxies, for which the statistics is poor) normalized to the 12$`\mu m`$ power observed by IRAS. The optical and near-infrared data have been corrected for aperture sizes to represent the total fluxes (S95). The shape of the optical-to-far-infrared spectral energy distribution (hereinafter SED) changes from one class to another: Seyfert 1’s show a general decrease from the optical to the far-infrared with only a small bump around the L (3.6$`\mu m`$) band; Seyfert 2’s show two pronounced peaks with about the same power around the J-H bands (1.2-1.6$`\mu m`$) and at 60$`\mu m`$; starburst galaxies have again the same two peaks in their SEDs, but the far-infrared peak is much brighter than the near-infrared peak; finally normal galaxies show two peaks with about the same power.

To better show the far-infrared turnover, we show in Fig. 3 the mean far-infrared SED for each class of galaxies, except nearby spirals, normalized to the observed power at 60$`\mu m`$.

The main results here are: i) Seyfert 1’s have a flatter SED shortward of 60$`\mu m`$, compared to all other classes of galaxies; ii) the two Seyfert types have virtually identical spectra from 60 to 200$`\mu m`$, while they appear different only shortward of 60$`\mu m`$ ; iii) the SED of starburst galaxies show the steepest drop off beyond 60$`\mu m`$; iv) the normal galaxies in the sample observed show stronger emission at wavelengths longer than 150$`\mu m`$ than the other types of galaxies.

The various types of galaxies show a sequence in the slope of the short wavelength part of the SED (from 12 to 60$`\mu m`$) from the very flat Seyfert 1’s, through the intermediate Seyfert 2’s, to the steep starburst and normal galaxies. The flatter SED of Seyfert 1’s in the range 12-120$`\mu m`$ probably arises from the warmer dust heated by the active nucleus.

### 3.2 Testing unified models of Seyfert 1’s and 2’s

Unified models claim that the observational differences between Seyfert 1’s and 2’s can be attributed to the different orientation of a hypothetical dusty torus (e.g. Antonucci, 1993). To test this hypothesis, we have fitted in Figure 4 the average slope of the 12-200$`\mu `$m far-infrared SEDs of Seyfert 1’s and Seyfert 2’s with the sum of an optically thick dusty torus seen face and edge on respectively, from the models by Granato & Danese (1994) and grey-body thermal emission at 25K with inverse square wavelength dependence of the dust emissivity (meant to represent the extended dust emission from the galactic disk). As shown in Figure 4, the cool component produces about 40% of the observed fluxes at 60$`\mu `$m. Although any detailed model fitting is beyond the scope of the present paper, there is a rough qualitative consistency between the data and this simple model. However, in detail the model fails because it does not predict strong enough flux differences at 12 and 25$`\mu `$m as the torus orientation shifts from face-on to edge-on. With a more realistic range of torus orientations in Seyfert 1’s and 2’s, the disagreement with observations would be even more significant.

### 3.3 Spectral Decomposition of Active and Quiescent Components

In this section we investigate if the decomposition of the observed SEDs in physically distinct spectral components that was proposed by Rowan-Robinson and Crawford (1989) (RRC) for the IRAS data can be extended to the longer wavelength ISO results.

First, we discuss the SED as defined by the IRAS data only: if we simply take the average flux ratios among the four IRAS wavelengths ($`F_{12\mu m}`$: $`F_{25\mu m}`$: $`F_{60\mu m}`$: $`F_{100\mu m}`$) of the different types of galaxies, from our sample we obtain for Seyfert 1’s: $`1:2.3:10.3:21.7`$, for Seyfert 2’s: $`1:3.3:20:29.2`$, for starburst galaxies: $`1:4.5:27:35`$ and for normal spiral galaxies: $`1:2.4:16:30`$. Subtracting the normal galaxies ratios from those of starburst galaxies, we are able to identify a “warm” starburst component with ratios: $`1:2:11:5`$, peaking at 60$`\mu `$m. As an exercise, if we add together the ratios of Seyfert 1’s with those of this warm starburst component, we can roughly reproduce the ratios of Seyfert 2’s. This suggests that there might be a difference between the two types of Seyfert nuclei in the star forming activity of their host galaxies, confirming the finding that Seyfert 2 nuclei lie preferentially in galaxies experiencing more enhanced star forming activity compared to Seyfert 1’s (Maiolino et al., 1995). It is possible that this results from a selection effect: since Seyfert 2 nuclei are relatively weaker at 12$`\mu m`$ than Seyfert 1’s, they are more likely to fall into our flux-limited sample if their brightness at this wavelength is enhanced by strong star formation.

We want now to extend the analysis including the longer wavelength ISOPHOT data. Our broad-band SEDs are probably too crude a description of the complexities of galaxies to allow a full principal component decomposition into physically distinct emission components. Nonetheless, we have followed RRC in using the 12-25-60-100$`\mu `$m colors to identify those galaxies in our sample which closely resemble the SEDs of the “quiescent cirrus” disk, the “starburst” component, and the “pure Seyfert” nucleus. For each of the three types of galaxies, the normal spirals, the starburst galaxies and the Seyfert 1’s, we have selected those objects lying in the two IRAS color-color diagrams close (i.e., within 0.2 magnitudes) to the colors of “pure disc”, “starburst” and “Seyfert” components of RRC (labeled D, B and S, respectively in Fig.1 of RRC). We have then plotted in Fig.5 the complete 12-200$`\mu `$m SEDs of the selected galaxies individually and their average value for each class in Fig.6.

As shown in Fig.5, five normal spiral galaxies have nearly the RRC colors of the “pure cirrus disk” component: $`F_{12\mu m}`$: $`F_{25\mu m}`$: $`F_{60\mu m}`$: $`F_{100\mu m}`$ = 1:1:12.7:40. These are the only normal galaxies in our sample that have the color \[60 - 25\]$``$0.8 and the color \[100 - 60\]$``$0.3. Similarly seven starburst galaxies have IRAS colors close to the “pure starburst” spectrum with the four IRAS fluxes scaling as 1:5:24.5:24.5. These galaxies have been selected because they are the only starburst galaxies in our sample with 0.6$``$\[60 - 25\]$``$0.9 and -0.05$``$\[100 - 60\]$``$0.25. Seven Seyfert 1’s galaxies have IRAS colors near the “pure Seyfert nucleus” and are the only Seyfert 1’s in our sample with \[60 - 25\]$``$0.25 and \[25 - 12\]$``$0. As also confirmed from our discussion of the two colors diagrams (see subsection 4.1), the strong distinction between quiescent disks and starbursts remains clear out to 200$`\mu `$m. The cirrus and starburst spectra probably represent extremes of minimal and maximal recent star formation, that tend to be found in the least and most luminous galaxies, respectively. The pure Seyfert spectrum is rather similar to the pure starburst spectrum between 100 and 200$`\mu `$m. Both show a relative lack of cold dust, and the Seyfert’s tend to be weaker at 120$`\mu `$m.

The Seyfert 2’s are spread all around the IRAS color-color diagrams. As can be seen in Fig.5, some of them have IRAS spectra close to the pure starburst template (we have selected four in Fig.5). And indeed their ISOPHOT far-infrared spectra also match the pure starburst spectrum well, since their infrared continuum appears to be dominated from dust around star forming regions. Those Seyfert 2’s with IRAS colors like quiescent cirrus disks also resemble pure disks in the 100–200$`\mu `$m region. Again it appears that the Seyfert 2’s nucleus contributes a minor fraction of the observed far-infrared luminosity in those objects.

## 4 COLOR-COLOR ANALYSIS

### 4.1 Far-infrared color-color diagrams

We have constructed color-color diagrams with the ISO and IRAS photometry to see how the SEDs of the different classes of galaxies differ, and if these can be used to segregate them. Fig.7 and 8 show, respectively, the \[200 - 100\] versus \[60 - 25\] and the \[200 - 100\] versus \[100 - 60\] color-color diagrams555We define color
$$[\lambda _1\lambda _2]=Log(F_\nu (\lambda _1)/F_\nu (\lambda _2))$$
where $`F_\nu (\lambda _i)`$ is the flux density in Jy at $`\lambda _i`$. . The averaged colors for each type of galaxies are reported in Table 4.

From the two diagrams we conclude:

* For non-Seyfert galaxies, our interpretation of the different positions of the different types of galaxies in Fig.8 is that the \[100 - 60\] and especially the \[200 - 100\] colors define a correlated sequence running from galaxies with the weakest star formation up to the starbursts, in agreement with the SED behaviour as shown in Fig.3; Table 4 can be used to quantify this result: the \[100 - 60\] colors of starburst and nearby spiral galaxies do not overlap each other at the 1$`\sigma `$ level; even better, the \[200 - 100\] color is able to separate starburst from normal and normal from nearby spiral galaxies, while the colors of Seyfert’s are intermediate between starburst and normal galaxies.
* The star formation rate does not, however, influence the \[60 - 25\] colors; from Table 4 we note that the \[60 - 25\] color is virtually the same among starburst, normal and nearby spiral galaxies. The \[60 - 25\] color only shows a small displacement between Seyfert 1’s and 2’s, and Seyfert galaxies in general from normal and starburst galaxies; its average value decreases moving to Seyfert galaxies and becomes the lowest for Seyfert 1’s (Table 4). However, as can be seen from the Table 4, the scatter in individual galaxies is comparable to the amount of this weak shift, making it marginally significant.
* Seyfert 1’s and 2’s are not distinguishable from normal galaxies, or from each other in either \[100 - 60\] or \[200 - 100\] colors; as Table 4 shows, for both colors Seyfert and normal galaxies overlap at the 1$`\sigma `$ level.
* As can be seen from both Fig. 7 and 8 and from Table 4, all the infrared colors of the CfA Seyfert’s are more similar to those of normal inactive galaxies than are those of our 12$`\mu `$m Seyfert’s.

The first result indicates that both \[200 - 100\] and \[100 - 60\] colors decrease when star formation activity increases. We will confirm (see subsection 4.4) that the former color also correlates with the 60$`\mu `$m excess, that we consider a measure of the star formation activity in a galaxy. The second result was already known (e.g. S95). It comes about because recent star formation boosts both the 25 and 60$`\mu `$m luminosities of galaxies in the same proportions found in quiescent spiral galaxies. The ratio of these two luminosities remains constant. The best explanation for the third result is that none of the Seyfert nuclei in these galaxies make a substantial contribution to the total fluxes observed at wavelengths of 60$`\mu `$m and longward. The emission from Seyfert nuclei is relatively stronger at 25$`\mu `$m, compared with that of normal galactic disks. This indicates that the thermally emitting dust in Seyfert nuclei has a substantially hotter temperature distribution than in normal galaxies, even those with strong starbursts. Furthermore, the dust temperatures in Seyfert 1’s tend to be even higher than those in Seyfert 2’s, as was already appreciated by Edelson & Malkan (1986) and Edelson, Malkan & Rieke (1987). The fourth empirical result illustrates a general difference between the CfA and our 12$`\mu `$m Seyfert samples: in the former the active nuclei tend to be relatively less luminous than their host galaxies. This results in far-infrared colors in the CfA Seyfert’s which are more “normal”, i.e. similar to those of non-active spiral galaxies. One prediction is that the mid-infrared (10–25$`\mu `$m) “compactness” \[defined in Edelson & Malkan (1986) as the ratio between the flux in a small beam and the total integrated flux\] of the CfA Seyfert’s is less than that of the 12$`\mu `$m Seyfert’s, even at a given redshift.

Figures 9 and 10 show the far-infrared colors \[150 - 100\] and \[200 \- 150\] versus the IRAS colors \[60 - 25\] and \[100 - 60\], respectively. The \[200 - 150\] color is able to separate starburst from normal galaxies, as is the \[200 - 100\] color, while the \[150 \- 100\] color is not. This means that the separation in mainly due

to the longer wavelength slope. Finally the separation between CfA and non-CfA Seyfert galaxies is even more apparent in the \[200 - 150\] color than in the \[200 - 100\] color. The CfA Seyfert’s lie in the upper part of the diagrams at \[200 - 150\] $`>`$ -0.4, towards the position of normal galaxies, while most of the non-CfA Seyfert’s are located at \[200 - 150\] $`<`$ -0.4.

### 4.2 Spectral Curvature

Color-color diagrams can also be used to understand if the spectra between 100–150 and 150–200 $`\mu `$m can or cannot be described by simple power laws. In Fig. 11 we compare the 120 to 100$`\mu `$m flux ratio with the broader wavelength baseline of 100 to 150$`\mu `$m. The solid line shows the values that would be observed for pure power laws $`F_\nu \nu ^\alpha `$ with:

$$\alpha =\frac{log(F_{\nu _1}/F_{\nu _2})}{log(\nu _1/\nu _2)}=\frac{log(F_{\nu _1}/F_{\nu _3})}{log(\nu _1/\nu _3)}$$
(2)
where $`\nu _1`$, $`\nu _2`$ and $`\nu _3`$ are three generic frequencies and $`F_{\nu _i}`$ the correspondent flux densities.

However, most galaxy flux ratios are better fitted with curved 100–120–150$`\mu `$m spectra. Points lying to the right of this line show downward spectral curvature. Most of the sample — except for 4 normal galaxies and a few low-luminosity Seyfert’s — do in fact have stronger 120$`\mu `$m emission than would be predicted by a power-law interpolation between 100 and 150$`\mu `$m. Relative to their 100 and 150$`\mu `$m fluxes, only a few galaxies show deficits of 120$`\mu `$m emission. The dashed regression line (Fig.11),

$$[150100]=1.06\times [120100]0.09$$
(3)
shows relatively little scatter about it. The dispersion of 0.09 dex corresponds to random errors in each flux of 16%, not much larger than our typical uncertainties. The line correspondent to pure power law behaviour for these colors has a slope of 2.22 compared to the flatter slope of 1.06 fitting the data.

Similarly, the regression line fitting the \[170 - 150\] color versus the \[200 - 150\] color is given by:

$$[170150]=0.09\times [200150]0.03$$
(4)
The line correspondent to a pure power law for these two latter colors has a slope of 0.43 compared to the flatter slope of 0.09 fitting the data.

The situation appears even simpler for the 180$`\mu `$m flux, which generally does fall near a power law interpolation between 150 and 200$`\mu `$m (as shown by the solid line in Fig. 12). Thus we find that for most galaxies, measurements of 100, 150 and 200$`\mu `$m would be sufficient to predict the 120, 170 and 180$`\mu `$m fluxes accurately, but not with a simple power law interpolation.

### 4.3 Dust Color Temperatures

Another important outcome that can be derived from the color-color diagrams is the estimate of the temperature(s) that is (are) dominating the emission in the far-infrared. To perform this analysis we have chosen the \[200 - 100\] versus \[60 - 25\] color-color diagram and we adopted for fitting the combination of two grey-bodies, leaving to vary the dust emissivity law. We show this color-color diagram separately for the different types of galaxies in Fig. 13. The emission is broad-band: fitting it from 25 to 200$`\mu `$m requires a range of dust temperatures. To compare the model with the observed colors, we have overplotted the combination of the two grey-bodies in Fig. 13. This introduces the minimum necessary number of fitting parameters. To fit the long-wavelength color a steep dependence of the dust emissivity law is needed,

$$ϵ(\lambda )\times B(\lambda ,T)\mathrm{with}ϵ\lambda ^2$$
(5)
A much flatter wavelength dependence would not fit most of the 200$`\mu m`$ data for any type of galaxies.

The non-Seyfert galaxies show a clear sequence of increasing average dust temperature going from the nearest spirals up to the most luminous starbursts. The correlation of the two colors is:

$$[200100]=1.125\times [6025]1.252$$
(6)
(the linear regression coefficient is 0.33 for 25 data points) or alternately, for the \[100 - 60\] color:

$$[200100]=1.626\times [10060]0.700$$
(7)
(the linear regression coefficient is 0.73 for 26 data points)

These correlations are essentially re-statements of the correlation first presented by S95, in their Appendix B (see their Figure 13). The corresponding shift in grey-body temperature that we find is the same as what they found, after correction for the fact that their shallower wavelength dependence of dust emissivity made their fitted temperatures systematically higher than ours for a given observed spectrum. Our interpretation is that higher rates of star formation raise the average dust temperature, shifting the peak of the thermal dust emission to higher frequencies, and tilting the entire far-infrared spectrum toward the blue. Splitting the 100–200$`\mu `$m spectrum into two pieces, we find that most of this curvature occurs between 150 and 200$`\mu `$m.

In the normal spirals belonging to the 12$`\mu m`$ galaxy sample and observed by Siebenmorgen, Krugel & Chini (1999) and by Alton et al. (1998), dust at even lower temperature (15 K) is detected, confirming the previous finding of these latter authors.

Most of the Seyfert 1’s far-infrared data occupy a region that can be fitted by the mixture of two grey-bodies with $`\mathrm{T}_{\mathrm{min}}=22\mathrm{K}`$ and $`\mathrm{T}_{\mathrm{max}}=5570\mathrm{K}`$. The three nearby CfA Seyfert 1’s galaxies (NGC 3227, NGC 5033, NGC 4051, with data taken from PGRE) with an average redshift of z=0.0031, show the highest value of the \[200 - 100\] color, requiring a lower minimum dust temperature (the fit in Fig. 13a gives $`\mathrm{T}_{\mathrm{min}}=18\mathrm{K}`$).

Almost all Seyfert 2’s far-infrared data (Fig. 13b) can be fitted by the mixture of two grey-bodies with $`\mathrm{T}_{\mathrm{min}}=25\mathrm{K}`$ and $`\mathrm{T}_{\mathrm{max}}=5565\mathrm{K}`$. They show a somewhat narrower temperature distribution compared to Seyfert 1’s, and do not extend to the higher values. Again the four reddest CfA Seyfert 2’s — NGC 4388, NGC 3982, NGC 1143/44 and NGC 3079 — (PGRE) with the coolest dust ($`\mathrm{T}_{\mathrm{min}}=20\mathrm{K}`$), have an average redshift of 0.011, as opposed to 0.019 for all the other Seyfert 2’s plotted.

Comparing the three classes of galaxies, there is a trend with the most active objects (Seyfert 1’s) have the widest temperature range (22–70K); Seyfert 2’s are intermediate, while starbursts have the warmest low temperature component (30K), together with the coldest high temperature one (60K).

A similar result has been found by Klaas et al. (1997) analyzing the ISOPHOT data of the three interacting galaxies Arp 244, NGC 6240 and Arp 220: the increase of the highest dust temperature component is accompanied by the decrease of the lowest dust temperature component they detect. In four narrow-line Seyfert 1’s Polletta & Courvoisier (1999) found that the lowest dust temperature is in the range 20-40K.

We are aware that there could be an additional selection effect in nearby galaxies if there is very extended emission from the galaxy discs. Measurements not performed in mapping mode could have lost some of the long wavelength emission. The detection of very cold dust only in nearby objects could thus be explained with the very low surface brightness of its thermal emission. This dust component has been found by Alton et al. (1998) to be surprisingly extended. We cannot exclude therefore that also for our objects there might be emission from cold dust (T $``$ 15K) that escaped detection with our observations.

### 4.4 The far-infrared turnover versus the 60$`\mu `$m excess

In this subsection we want to relate the strength of the star formation activity in a galaxy with the steepness of the far-infrared turnover, measured by the \[200 - 100\] color, with the aim to see how the different classes of galaxies behave and if there is a separation between galaxies with and without active nuclei.

We have chosen as the indicator of enhanced recent star formation, which warms dust around HII regions, the 60$`\mu `$m “excess” as the ratio of the observed 60$`\mu `$m flux to the flux that a source would have at 60$`\mu `$m from power-law interpolation of the flux between 12 and 100$`\mu `$m. We have already seen (see subsection 4.1) that the starburst galaxies of our sample are characterized by a low value of their \[200 - 100\] color. In Fig.14 we plot the \[200 - 100\] color index as a function of the 60$`\mu `$m excess. We find a correlation between the steepness of the far-infrared (100–200$`\mu `$m) turnover and the strength of the 60$`\mu `$m “excess”. While this diagram does not perfectly separate the galaxies of different classes, it nevertheless shows that they cluster preferentially in different regions of the diagram. Seyfert 1’s (excluding six objects of the CfA sample) cluster in a no-60$`\mu `$m excess region with a color \[200 - 100\] $`<`$ 0. The starburst galaxies cluster in the central area of the diagram and have all 60$`\mu `$m excess. Normal galaxies and nearby spirals have no 60$`\mu `$m excess (except 2 objects) and a color $`0.5<[200100]<+0.5`$. Seyfert 2’s are widely spread all over the diagram, but with a 60$`\mu `$m excess generally higher than Seyfert 1’s. Four especially infrared-luminous galaxies — MK 273, Arp 220, MK 938 and FSC05189-2524 lie in the “starburst” zone. In fact detailed mid-infrared spectroscopy (in the range 3-30$`\mu `$m) suggests that much of the total infrared luminosity in these extreme objects is indeed powered by star formation (Genzel et al., 1998). As we already noted, many of the CfA Seyfert galaxies have higher excess and/or higher \[200 - 100\] color, with respect to the other Seyfert galaxies. This is presumably because their faint Seyfert nuclei are too feeble to influence the overall far-infrared continuum emitted by their bright host galactic disks.

As stated above, the 60$`\mu `$m excess is measured with respect to the underlying 12 and 100$`\mu `$m fluxes, which are dominated by the “cirrus” component powered by the integrated stellar luminosity in the galaxy. When the 60$`\mu `$m bump is especially strong, the warm dust associated with recent star formation even contributes to the observed 100$`\mu `$m flux, making a relative excess at that wavelength compared with the emission from cold dust at longer wavelengths. The result is that the more actively star-forming galaxies, such as infrared-luminous starbursts, have the sharpest far-infrared turnover. They show the steepest drop off from 100 to 200$`\mu `$m, with an average slope over that wavelength range of $`\alpha _{100200\mu m}`$ = + 2.2, compared to $`\alpha _{100200\mu m}`$ = + 1.4, $`\alpha _{100200\mu m}`$ = + 1.6 and $`\alpha _{100200\mu m}`$ = + 0.8 for Seyfert 1’s, 2’s and normal galaxies, respectively (see Table 4).

A simple least squares fit to the data show a correlation of the color \[200 - 100\] with the 60$`\mu `$m excess with a linear regression coefficient of R=-0.48 (for 83 objects), corresponding to a probability of 99.9996%. This correlation improves by excluding the Seyfert 1’s (R=-0.61 for 65 objects).

We suggest that that the diagram shown in Fig.14 can be used to separate the starburst dominated objects from the AGN dominated ones. Objects located in the upper right part of the diagram are “more” starburst dominated, while those at the left, having a fainter excess, are the AGN dominated objects. We suggest that starburst activity in galaxies, i.e. with high rates of current star formation, results in excess emission in the 60$`\mu `$m band accompanied by a general heating of the galactic ISM and thus a decrease of the \[200 - 100\] color.

We note that many of the Seyfert galaxies of the CfA sample are shifted upwards in Figure 14 with respect to the bulk of the 12$`\mu `$m Seyfert’s. Again we interpret this as an indication that their stronger long-wavelength emission arises from the dominance of quiescent disk dust.

## 5 LUMINOSITY ANALYSIS

### 5.1 Luminosity Correlations

In this subsection we analyze the correlations of the 200$`\mu `$m luminosity with the luminosity in each of the IRAS bands and the total mid- to far-infrared luminosity. Although diagrams which correlate one luminosity with another in the same object should be used with caution, we have a large sample with good selection criteria in which the luminosities span over 4 orders of magnitude. In Figures 15 and 16 we show such correlations and give the formulae for the best-fitting regression lines. As expected, the slopes in the correlations of the 200$`\mu `$m luminosity against the 25 and 60$`\mu `$m luminosities are flatter than those against the 12 and 100$`\mu `$m luminosities (see also S95). This is because the 25 and 60$`\mu `$m luminosities are preferentially elevated in those galaxies with more active recent star formation. The 100$`\mu `$m and 200$`\mu `$m emission, in contrast, are relatively stronger in the quiescent disk emission (see below), and correlate almost linearly with each other.

These figures are a logical extension of the findings of S95. They showed that galaxies with higher bolometric luminosity are relatively brighter in the 25 and 60$`\mu `$m bands, due to the elevated “starburst” component in them, while they appear relatively fainter at both shorter (12$`\mu `$m) and longer (100$`\mu `$m) wavelengths. The “pivot points” at which the relative luminosity neither increases nor decreases, occur around 12 and 100$`\mu `$m. This is our interpretation of why the 12 and 100$`\mu `$m luminosities have a correlation with the bolometric luminosity with a flatter slope (1.06 and 1.13, respectively) compared to the slope of the correlations of the 25 and 60$`\mu `$m luminosities with the bolometric luminosity ($``$ 1.2)(see Fig.7 of S95). For wavelengths much shorter than 12$`\mu `$m or longer than 100$`\mu `$m we therefore expect the luminosity to track less than linearly with the bolometric luminosity. That is what these figures are showing.

It is not particularly surprising that the 200$`\mu `$m luminosity correlates most closely with the 100$`\mu `$m luminosity, given that these are the two closest wavelengths. The fact that the Seyfert galaxies (of both types 1 and 2) lie on the same track as the non-Seyfert’s once again indicates that the Seyfert nucleus makes a negligible contribution at these wavelengths. Since the 200$`\mu `$m luminosity is dominated by the “quiescent cirrus” component, it does not increase as fast as the 25 or 60$`\mu `$m luminosities, which are strongly effected by starbursts. Thus the slopes of the correlation between the 200$`\mu `$m luminosity and those at these latter wavelengths are much flatter than 1 (0.83 and 0.86 respectively).

The 200$`\mu `$m luminosity even decreases less rapidly than the total luminosity. This is seen by the flat slope of only 0.79 in Figure 16, where the X-axis is the “total” mid- to far-infrared luminosity, defined in equation (1) and given in Table 3. We note that a similar slope (0.87) is relating the optical blue luminosity with the bolometric luminosity (see Fig.10 of SM95), confirming the above interpretation.

### 5.2 Color-luminosity diagrams

As discussed above, it is believed that the more luminous non-Seyfert galaxies have higher light/mass ratios. They are luminous not so much because they are bigger or more massive, but principally because they have higher proportions of recently formed stars, which are highly luminous. A higher proportion of their total luminosity emerges in the far-infrared, and in particular from the warmed dust grains which emit strongly at 25–60$`\mu `$m– what we describe as the “starburst” component (e.g. S95). One consequence of this is that the most actively star-forming galaxies should have the sharpest far-infrared turnovers, as shown in subsection 4.4, because their dust has a warmer temperature distribution, as seen is subsection 4.3. Thus we expect a systematic trend of infrared SEDs with total luminosity. Spinoglio et al. (1995) already showed that the higher the bolometric luminosity of a galaxy, the higher its ratio of starburst/cool dust components.

Indeed we do find that the far-infrared turnover becomes systematically sharper in the more luminous galaxies. This can be seen in the inverse correlation between the color index \[200 - 100\] and the 12$`\mu `$m luminosity, shown in Fig. 17. In this figure we see that this relation holds in general, irrespective of the type of the galaxy, with a regression coefficient of R=-0.48 for 81 objects, for all the observed galaxies in the 12$`\mu `$m sample, excluding the 3C and PG quasars. The fit to the non-Seyfert galaxies only is:

$$[200100]=13.5130.317\times Log(L_{12\mu m})$$
(8)
(the linear regression coefficient is -0.52 for 26 data points).

We also plot in Fig.18 the inverse correlation between the color index \[200 - 100\] and the 60$`\mu `$m luminosity, which gives a regression coefficient of R=-0.48 for 83 objects, for all the galaxies in our sample, excluding the 3C and PG quasars. The fit for the non-Seyfert’s only gives:

$$[200100]=14.9830.346\times Log(L_{60\mu m})$$
(9)
(the linear regression coefficient is -0.63 for 26 data points)

The four most noticeable outliers are the two 3C objects and the two PG quasars, which have unexpectedly red \[200 - 100\] colors . There is therefore an indication that even at low redshift (z $``$ 0.2) quasars (either the optically selected PG and the radio selected 3C) might have an excess in cold dust emission.

In Figures 19 and 20, we separate this color-luminosity trend into the slopes between 100 and 150$`\mu `$m and 150 to 200$`\mu `$m. Clearly the trend is due entirely to the steepening of the 200–150$`\mu `$m slope in the more luminous galaxies, which is the same whether the luminosity is measured at 12 or at 60$`\mu `$m. The 100–150$`\mu `$m slope in fact shows no significant luminosity dependence.

## 6 SUMMARY

New ISOPHOT photometry of a large sample of nearby active and normal galaxies shows that:

* The 60-200$`\mu m`$ SEDs of active and normal galaxies are similar. The 100 to 150 to 200$`\mu `$m spectra have significant curvature, generally downward. The trends are sufficiently systematic that measurements at these wavelengths can give accurate predictions of the fluxes at the intermediate wavelengths of 120, 170 and 180 $`\mu `$m.
* In the non-Seyfert galaxies, higher luminosities (presumably generated by higher rates of recent star formation) are correlated with hotter average dust temperatures at all wavelengths from 12 to 200$`\mu `$m. In particular, the 150 to 200$`\mu `$m slope steepens in more luminous galaxies.
* The decomposition of the observed SEDs into physically distinct spectral components has been extended from the IRAS data (RRC) to our ISO results: a strong distinction is apparent out to 200$`\mu m`$ between the the quiescent disk component and the starburst component. We suggest that these components represent the extremes of minimal and maximal recent star formation, found in the least and most luminous galaxies, respectively.
* The \[200 - 100\] versus \[100 - 60\] color-color diagram is able to separate the three types of galaxies: starburst galaxies, normal galaxies and nearby spirals, while Seyfert’s are at intermediate locations. The large-aperture 60–200 $`\mu `$m spectra of Seyfert galaxies we measured are dominated by the emission from the host galaxy. They therefore depend on the relative contributions of “starburst” and “quiescent disk” emission in the galaxy, and are independent of the presence of an active nucleus.
* The \[200 - 100\] vs \[60 - 25\] color-color diagram shows that a mixture of two black-bodies with a warm and a cold component is able to fit all the far-infrared data, only if we assume the dust emissivity law proportional to the inverse square of the wavelength. Comparing active galaxies, Seyfert 1’s show the larger range in temperature (22-70K), Seyfert 2’s are intermediate, while starburst galaxies have the narrower range in dust temperature (30-60K). Nearby normal spirals and some of the CfA Seyfert galaxies show very cold dust components (15-18K).
* A correlation is found between the slope of the far-infrared turnover (as measured from the \[200 - 100\] color) and both the 60$`\mu `$m excess and the 12 and 60$`\mu `$m luminosities. We have confirmed and extended the finding of S95 that increasing star formation rates cause the more luminous non-Seyfert galaxies to have relatively stronger 25–60$`\mu `$m emission, resulting in systematically bluer far-infrared colors. We suggest that the plane defined by the \[200 - 100\] color versus the 60$`\mu `$m excess can be used to separate starburst dominated galaxies from AGN dominated ones.

This research was funded in Italy from the Italian Space Agency (ASI). This work benefited from extensive use of the NASA Extragalactic Database (NED) and from the help given by the IPAC staff in planning the ISO observations. We thank Brian Rush and Nao Suzuki for help in preparing the observing command templates. We also thank the anonymous referee, whose thorough comments helped in improving this article.

## Appendix A DATA OBTAINED WITH THE ISOPHOT C100 ARRAY CAMERA

Some C100 observations, namely for objects FSC00521-7054, NGC 7674, MCG+1-33-36, NGC 262, FSC03362-1642, ESO253-G3, NGC 1365, NGC 4501, NGC 4922A/B were performed in staring mode, however only for the latter three an adjacent empty sky position was observed to measure the sky background. In the former cases sky background was estimated as follows: we assume that the sources are point-like and therefore only the central pixel (# 5) contains the source flux. The adjacent pixels were averaged to get a background value, which was then subtracted from the signal of the central pixel. This procedure is, however, affected by large uncertainties since it requires a good flat-fielding of the array because of its unstable behavior and in two cases, for NGC 7674 and NGC 262, the values inferred are much lower than the corresponding IRAS photometry. Note that the uncertainties quoted here are only the statistical ones. The results of the C100 observations are reported in Table 5.


format: text

text:
# Optimal prediction of folding rates and transition state placement from native state geometry

## I Introduction

In the past three decades, there has been a growing effort of the scientific community for studying and understanding the principles that govern the folding process of a sequence of amino acids in the corresponding native structure Anfinsen (1973); Creighton (1993); Branden and Tooze (1991). In recent years, several proteins, in particular those folding via a two-state mechanism Jackson and Fersht (1991) have provided an extraordinary benchmark for experimental and theoretical characterization of the folding pathways. The significant amount of experimental data available for several structurally unrelated proteins Jackson and Fersht (1991); Jackson (1998); Huang and Oas (1995); Burton et al. (1996); Kragelund et al. (1995, 1996); Mines et al. (1996); Schonbrunner et al. (1997); Reid et al. (1998); Grantcharova and Baker (1997); Villegas et al. (1995); Khorasanizadeh et al. (1993); Scalley et al. (1997); Silow and Oliveberg (1997); Van Nuland et al. (1998); Taddei et al. (1999); Ferguson et al. (1999); Smith et al. (1996); Kuhlman et al. (1998); Hamill et al. (1998); Tan et al. (1996), has opened the possibility to identify and isolate the factors that influence the folding rate. Besides considering detailed chemical interaction, such as those affecting free-energy barriers, an appealing and elegant line of investigation has focused on the effects of the native state structure on the folding process Plaxco et al. (1998a); Debe and Goddard III (1999); Dinner and Karplus (2001); Ivankov and Finkelstein (2001).

From a qualitative point of view, the influences of structural effects was traditionally summarised in the tenet that proteins with high helical content fold faster than proteins with mixed alpha/beta content, the slowest folding being for the all-beta ones. This useful and intuitive rule of thumb, fails to account for the very different rates observed between proteins in each of the alpha, alpha/beta or beta families Jackson (1998); Aurora et al. (1997); Maritan et al. (2000a); Capaldi and Radford (1998). A deep insight into this problem was provided by the work of Plaxco et al. Plaxco et al. (1998a) who introduced the concept of contact order, which captures, quantitatively, features beyond the mere secondary structure motifs. The highly significant correlation of contact order and experimental folding rates shows the extent to which the mere topology of native state can influence the folding process. However, the highly organised native structure of proteins is too rich to be captured by a single parameter such as the contact order. Indeed, the latter cannot account in the same satisfactory way for the transition state placement, three-state folding rates or the diversity of folding rates among structurally similar proteins Clarke et al. (1999).

In the present study we investigate how the topology of the native state can be further exploited to provide optimised predictions for protein folding rates and the transition state placement. To do so we consider, among others, one particular topological descriptor that is crucial for characterising the connection and interactions of native contacts: the clustering coefficient , or cliquishness. Such parameter, heavily studied in the context of graph theory Bollabas (1985); Watts and Strogatz (1998); Strogatz (2001) is shown to have highly significant correlation with folding rates. The advantage of using this topologic descriptor is that it allows to capture the cooperative formation of native interactions, as proved by its statistically relevant correlations with the transition state placement. Further, we discuss how the different topologic aspects captured by the cliquishness and contact order can be combined to yield optimal correlations higher than for the individual descriptors.

## II Theory and Results

Customarily, at the heart of theoretical or numerical studies of topology-based folding models is the contact matrix (or map) Go and Scheraga (1976) which will be used extensively also in the present context. The generic entry of the contact map, $`\mathrm{\Delta }_{ij}`$, takes on the value 1 if residues $`i`$ and $`j`$ are in contact and zero otherwise. Several criteria can be adopted to define a contact; in the present study we consider two amino acids in interaction if any pair of heavy atoms in the two amino acids are at a distance below a certain cutoff, $`d`$. All values of $`d`$ between 3.5 Å and 8 Å have been considered and reported. The contact map provides a representation for the spatial distribution of contacts in the native structures that is both concise and often reversible (since native structures can be recovered when appropriate values of $`d`$ are used). Plaxco and coworkersPlaxco et al. (1998a) have used the contact map to describe and characterize the presence and organization of secondary motifs in protein structures. The parameter that was introduced, the relative contact order, provides a measure of the average sequence separation of contacting residues and is defined as

$$\text{relative contact order}=\frac{1}{L}\frac{\underset{ij}{}\mathrm{\Delta }_{ij}|ij|w_{ij}}{_{ij}\mathrm{\Delta }_{ij}w_{ij}},$$
(1)
where $`i`$ and $`j`$ run over the sequence indeces, $`w_{ij}`$ is the contact degeneracy (i.e. the number of pairs of heavy atoms in interaction) and $`L`$ is the protein length. Remarkably, the contact order was shown have a highly significant linear correlation with experimental folding rates. The result of Plaxco and coworkers can be explained, a posteriori, with intuitive arguments: a high contact order corresponds to few local interactions. One may thus expect that the route from the unfolded ensemble to the native state is slow, being hindered by the overcoming of several barriers Bryngelson et al. (1995); Wolynes et al. (1995); Dill and Chan (1997); Dobson et al. (1998); Sali et al. (1994) due to spacial restraints, as recently analysed by Debe and Goddard Debe and Goddard III (1999) and previously by Chan and Dill Chan and Dill (1990) and also observed in topology-based numeric studies Cecconi et al. (2001). These considerations are based purely on geometric arguments and do not take into account the influence of specific interactions between the residues. In principle, the latter may well override the topological influence on the folding process, but surprisingly, as remarked in a recent review article Baker (2000) this is often not the case Cecconi et al. (2001); Micheletti et al. (1999); Clementi et al. (2000); Galzitskaya and Finkelstein (1999); Munoz et al. (1999); Alm and Baker (1999); Klimov and Thirumalai (2000).

Our aim is to exploit as much as possible the topologic information contained in the native state to improve both the accuracy of predictions for folding rates and gain more fundamental insight into the process. To this purpose we have considered additional topologic descriptors besides the contact order. The one that appeared most significant is a parameter termed cliquishness or clustering coefficient Bollabas (1985); Watts and Strogatz (1998); Strogatz (2001). For a given site, $`i`$, the cliquishness is defined as:

$$\text{cliquishness}(i)=\frac{\underset{jl}{}\mathrm{\Delta }_{ij}\mathrm{\Delta }_{il}\mathrm{\Delta }_{lj}}{_{jl}\mathrm{\Delta }_{ij}\mathrm{\Delta }_{il}}=\frac{\underset{jl}{}\mathrm{\Delta }_{ij}\mathrm{\Delta }_{il}\mathrm{\Delta }_{lj}}{N_c(N_c1)/2},$$
(2)
where $`N_c`$ is the number of contacts to which site $`i`$ takes part to. As for the contact order, also the cliquishness has an intuitive meaning; in fact it provides a measure of the extent to which different sites interacting with $`i`$ are also interacting with each other. Of course, the cliquishness is properly defined only if site $`i`$ is connected to, at least, two other sites. To ensure this, we included also the covalently bonded interactions $`[i,i\pm 1]`$ in (2). The importance of taking the cliquishness into account for discriminating fast/slow folders can be anticipated since a higher interdependency of contacts (large cliquishness) will likely result in a more cooperative folding process. In fact, the formation of a fraction of interactions will result in the establishment of a whole network of them. Consistently with this intuitive picture one should also expect that a large/small cliquishness will affect in different ways the amount of native-like content of the transition state.

We have tested and verified these expectations by calculating the average cliquishness for 40 proteins for which folding rates and transition state placement, $`\theta _m`$, have been measured. $`\theta _m`$ is deduced from the variation of folding/refolding rates upon change of denaturant concentration ($`m_F^{}`$ and $`m_U^{}`$) and provides an indirect indication of how much the solvent-exposed surface of the transition state is similar to that of the native one. $`\theta _m`$ ranges between 0 and 1; higher values denote stronger similarity with the native state. It is worth pointing out that, although the model underlying the calculation of $`\theta _m`$ relies on a two-state analysis, an effective $`\theta _m`$ can be inferred for three-state folders as well Jackson (1998). Since reliable $`\theta _m`$’s are not available for all proteins, the number of entries used to correlate the cliquishness and $`\theta _m`$ (see Tables 1 and 2) is slightly smaller than that used for tge logarigthm of refolding rates, $`\mathrm{ln}K_F`$. The set of proteins used, shown in Tables 1 and 2, was built up from experimental data collected in previous studies and predictions (often topology-based) of folding rates Jackson (1998); Plaxco et al. (1998a); Debe and Goddard III (1999); Dinner and Karplus (2001); Ivankov and Finkelstein (2001).

As indicated, the entries include both two-state and three-state folders, proteins belonging to the same structural family as well as proteins under different experimental conditions. This allows to examine to what extent predicted folding rates are consistent with the wide variations of folding velocities observed in structurally-related proteins and in different experimental conditions. As discussed in detail below, when the comprehensive set of Table 1 and 2 is used, the correlation found between cliquishness and folding rates is 0.71, with a statistical significance of $`t=10^5`$, more relevant than the one between a suitably defined contact order and folding velocities ($`r=0.66`$, $`t=510^5`$). As will be shown, the predicting power of the two quantities can be combined to achieve the optimal correlation of 0.74.

The prediction of the transition state placement, turns out to be more difficult when either of the two topologic parameters is used. While for the contact order it is equal to 0.23, the cliquishness yields the value of 0.48 which is not significantly improved by combining the two descriptors. Though the linear correlation of the clustering coefficient and the transition state placement is not as high as for the folding-rate case it is nevertheless statistically meaningful, having a probability of 0.004 to have arisen by chance.

### II.1 Two- and three-state folders

Before considering the more general case of all entries in Tables 1 and 2, we focus on two-state folders, i.e. proteins with a cooperative (all-or-none) transition between the unfolded and folded states. The neatness of this process, due to the absence of any significantly populated intermediate state, makes them ideal candidates for identifying and isolating the factors that influence the folding rate.

In the present context this separate test is important since it appears that the relative contact order is a much stronger descriptor for two-state folders, than for the general case. As a matter of fact, when both two- and three-state folders are considered, the influence of the average sequence separation of native contacts on folding properties is better captured by a different version of the contact order, which we shall term “absolute”, obtained when the r.h.s. of eq. (1) is not divided by $`L`$:

$$\text{absolute contact order}=\frac{\underset{ij}{}\mathrm{\Delta }_{ij}|ij|w_{ij}}{_{ij}\mathrm{\Delta }_{ij}w_{ij}}.$$
(3)
In the following we shall report and compare the performance of both parameters; furthermore we shall always consider the absolute value of the linear correlation coefficients, $`|r|`$, without regard to its sign, which can be easily inferred from the plots.

The original definition of contact order has an unrivaled performance in the prediction of folding rates for the two-state folders of Table 1. As visible in Fig. 1, it gives a stable correlation for cutoffs in the range 5 Å$`d`$ 7 Å, with the maximum value of $`r=0.80`$ for the cutoff $`d=4.5`$ Å. The statistical significance of such correlation can be quantified through a calculation of the probability, $`t`$, to observe by pure chance a correlation higher than the measured one (in modulus). The standard model underlying such estimates relies on the hypothesys of normal distribution of the deviates of the correlated quantities. As a rule of thumb, the upper value of $`t=0.05`$ is taken as a threshold for statistically meaningful correlations. For the value of $`r=0.80`$ reported above, this probability is $`t=310^5`$, which is, therefore, extremely significant.

Consistently with previous results, we found that the transition state placement is a much more elusive quantity to predict than folding rates. In fact, all topologic descriptors yield a poorer correlation compared to $`\mathrm{ln}K_F`$ (see Fig. 2). For the relative contact order, the best $`r`$ is 0.48 (for $`d=6.0`$) with an associated $`t=0.02`$. As anticipated, the performance of the absolute contact order in this particular context is significantly inferior then the relative one (see Figs. 1 and 2) and hence will not be further commented.

Concerning the performance of the novel parameter under scrutiny, the cliquishness, it can be seen from Figs. 1 and 2 that it is statistically meaningful for both folding rates and transition state placement. There are, however, significant differences with respect to the contact order. For folding rates the optimal $`r`$ is 0.67 ($`d=4.6`$ Å) and the associates value of $`t`$ is 5 $`10^4`$, one order of magnitude larger than for the relative contact order. For $`\theta _m`$ the situation is reversed since the optimal value of $`r=0.58`$ (for $`d=3.8`$ Å) has the statistical relevance of $`t=0.004`$, with a marked improvement over the previous case.

It is also interesting to note that cliquishness-based correlations have a non-trivial dependence on the cutoff $`d`$. In fact, due to the overall compactness and steric effects, the degree of dispersion of the cliquishness values for different sites in the same or different proteins is much more limited compared, e.g. to the average sequence separation of contacts. This leads to the observed decay of the correlations when the cutoff $`d`$ is increased.

The applicability of topology-based models are not limited to two-state folder, but can be extended to include three-state folders as well Micheletti et al. (1999); Clementi et al. (2000). Despite the addition of the 11 entries corresponding to three-state folders, the performance of cliquishness-based predictions for folding rates and $`\theta _m`$ improves from the values reported for two-state folders. As shown in Figs 3 and 4 the associated optimal correlations for $`\mathrm{ln}K_F`$ and $`\theta _m`$ are $`r=0.71`$ and 0.49, again observed for the same cutoff values ($`d`$) mentioned for the two-state case. The corresponding statistical significances are now, $`t=110^5`$ and $`t=0.004`$, which, despite the enlargement of the experimental set, show even an improvement over the two-state case.

From Figs. 3 and 4 it can be noticed that the performance of the relative contact order is noticeably poorer than the absolute contact order which, being a much better descriptor, becomes the focus of our analysis. The corresponding measured correlations are, in fact, $`r=0.66`$ for $`\mathrm{ln}K_F`$ and $`r=0.20`$ for $`\theta _m`$ with corresponding values of $`t=510^5`$ and $`t=0.23`$.

A direct comparison of how the clustering coefficient and the absolute contact order correlate with $`\mathrm{ln}K_F`$ and $`\theta _m`$ can be made by inspecting the plots of Figs. 5 and 6. It is worth pointing out that the analysis of the deviations from the linear trends of Figs. 5 and 6 reveals that a particular protein, 1urn, is among the top outliers for both cliquishness and contact order-based analysis, although no simple explanations is available for this singular behaviour. Although for both folding parameters the cliquishness gives a more significant correlation than contact order, the difference is particularly dramatic for the transition state placement which is notoriously difficult to capture with topology-based predictions Plaxco et al. (1998a).

An important conclusion stemming out of this observation is that the transition state structure (and hence $`\theta _m`$) is more influenced by the degree of interdependency of native contacts than their average sequence separation. This is in accord with the intuition that highly interdependent contacts may mutually enhance their probability of formation, thus facilitating the progress towards the native state during the folding process. This is, indeed, consistent with the negative correlation observed between cliquishness and native content, $`\theta _m`$, at the transition state. It is important to stress that the presence and effects of the cooperative formation of native interactions cannot be captured by parameters based on measures of contact locality. This highlights the importance of considering all viable topologic descriptors to characterize the folding process, since they do not impact in the same way on various folding properties.

### II.2 Optimal combined correlation

A natural question that arises is whether it is possible to combine the predicting power of cliquishness and contact order to achieve correlations with experimental folding rates and transition state placements that are better than the individual cases.

Indeed, as shown in Appendix A, it is straightforward to combine in an optimal linear way the two quantities to improve the prediction accuracy. The quantitative increment in the correlation is clearly related to the amount of independent information contained in the two topologic descriptors. Hence, an important issue is to what extent cliquishness and contact order are mutually correlated.

If, in place of a physical contact map, $`\mathrm{\Delta }_{ij}`$, one uses a random symmetric matrix, no meaningful correlation will be found. The contact maps of real proteins, however, display features that are highly non-random which reflect both (i) the physical constraints to which a compact three-dimensional chain is subject and (ii) the presence and organisation of secondary motifs Chothia (1984, 1992); Levitt and Chothia (1976); Rose and Seltzer (1977).

With the aid of numeric simulations it was possible to assess the degree of interdependency of clustering coefficient of native contacts and their average sequence separation resulting from the first of the mentioned effects. This was accomplished by considering, in place of the proteins of tables 1 and 2, 150 computer-generated compact structures respecting basic steric constraints found in real proteins (details can be found in the Methods section). As visible in the plot of Fig. 8 the level of mutual contact order-cliquishness correlation observed in these artificial structures is $`r=0.25`$ which is significantly smaller than the actual correlation of the two quantities found in real proteins. In fact, the typical correlation for cliquishness and contact order (either relative or absolute) is around 0.65. Such non trivial correlation can been ascribed to the special topologic properties of naturally occurring proteins whose ramifications have been investigated in a variety of contexts Maritan et al. (2000a, b); Denton and Marshall (2001); Micheletti et al. (1999); Hunt et al. (1994); Socci et al. (1994). Thus, the very presence and organization of secondary motifs in proteins makes it possible, on one hand, to exploit the native topology to predict e.g. folding rates, while on the other it limits the amount of independent information contained in different topologic descriptors.

Nevertheless, since the mutual correlation is not perfect, it is still possible to achieve, by definition, better predictions by combining cliquishness and contact order. The degree of enhancement depends also on the statistical significance of the individual starting correlations. For these reasons, the improvement is noticeable for folding rates, while it is not significant for transition state placement. For the case of two-state folders, the optimal combination yields correlations of $`r=0.86`$ while for the more general case of two and three-state folding rates one has $`r=0.74`$ which leads to a discernible improvement over previous cases, as visible in Fig. 7. To the best of our knowledge, this is the highest correlation recorder among similar studies involving a comparable number of entries (also including non-linear prediction schemes Dinner and Karplus (2001)). Due to the fact that the optimal combined correlations are found a posteriori, the associated values of $`t`$ are no more meaningful indicators of statistical significance.

Besides the cliquishness, we have investigated other parameters that are routinely used to characterise general networks (networks of contacts in our case). In particular, we considered the “diameter” of the contact map, defined as the largest degree of separation between any two residues, and also its average value. The diameter measures the maximum number of contact that need to be traversed to connect an arbitrary pair of distinct residues. Although the contact-map diameter is an abstract object, it conveys relevant topological information about protein structure, since it measures the long-range structural organisation. We found, a posteriori, that neither the maximum, nor the average diameter, correlate in a significant manner with the folding rate or transition state placement.

## III Conclusions

We have analysed important topological descriptors of organised networks (in our case the spatial network of native contacts) that could be used, individually, or in mutual conjunction, to describe and predict experimental parameters used to characterize the folding process. It is found that, besides the previously introduced contact order, a topologic parameter, termed cliquishness or clustering coefficient, is a powerful indicator of both the folding velocity and the transition state placement for two- and three-state folders. The predicting power of the cluquishness is that it takes into account the presence and organisation of clusters of interdependent contacts that are putatively responsible for the cooperative formation of native-like regions. This property appears well-suited to reproduce important features in the transition state that are otherwise elusive to other topologic analysis. The high statistical significance of the observed correlations testifies the strong influence of geometric structural issues on the folding process. The maximum predicting power is obtained when the topologic information contained in the cliquishness is used in combination with the contact order; this allows to reach a linear correlation as high as 0.74 with experimental folding rates recorder in 40 experimental measurements.

## IV Methods

### IV.1 Cross correlations

The linear correlation between two sets of data, $`\{x\}`$ and $`\{y\}`$ is obtained from the normalised scalar products of the covariations:

$$r=\frac{\underset{i}{}(x_i\overline{x})(y_i\overline{y})}{\sqrt{_i(x_i\overline{x})^2}\sqrt{_j(y_j\overline{y})^2}}$$
(4)
Without loss of generality, in the following we shall consider the sets of data to be with zero average and with unit norm, so that the expression of the correlation simplifies

$$r=\stackrel{}{x}\stackrel{}{y}$$
(5)
We now formulate the following problem. Two sets of data, $`\{x\}`$ and $`\{y\}`$ have linear correlation $`r_x`$ and $`r_y`$ respectively with a third (reference set), $`\{z\}`$. What is the maximum and minimum correlations we can expect between sets $`\{x\}`$ and $`\{y\}`$? We assume that $`r_x`$ and $`r_y`$ are positive since this condition can always be met by changing sign, if necessary, to the vector components.

The answer is easily found by decomposing $`\{x\}`$ and $`\{y\}`$ into their components parallel and orthogonal to $`\{z\}`$:

$$\stackrel{}{x}\stackrel{}{y}=b_{}c_{}+b_{}c_{}$$
(6)
Since $`{}_{}{}^{}c_{}^{}`$ is equal to $`r_xr_y`$, and hence is fixed, the maximum \[minimum\] correlation is found when $`b_{}`$ and $`c_{}`$ are \[anti\]parallel. Thus,

$$r_xr_y\sqrt{(1r_x^2)(1r_y^2)}rr_xr_y+\sqrt{(1r_x^2)(1r_y^2)}$$
(7)
Now we turn to a different, but related problem. How can we combine linearly $`\{x\}`$ and $`\{y\}`$, so to have the maximum correlation with $`\{z\}`$. The generic linear combination,

$$\stackrel{}{k}=\frac{\stackrel{}{x}+b\stackrel{}{y}}{\sqrt{1+b^2+2b\stackrel{}{x}\stackrel{}{y}}}$$
(8)
leads to the following correlations

$$\stackrel{}{k}\stackrel{}{z}=\frac{r_x+br_y}{\sqrt{1+b^2+2b\stackrel{}{x}\stackrel{}{y}}}$$
(9)
The maximum is achieved for

$$b=\frac{r_y\stackrel{}{x}\stackrel{}{y}r_x}{r_x\stackrel{}{x}\stackrel{}{y}r_y}$$
(10)
which yields

$$Max(\stackrel{}{k}\stackrel{}{z})=\sqrt{\frac{r_x^2+r_y^22\stackrel{}{x}\stackrel{}{y}r_xr_y}{1(\stackrel{}{x}\stackrel{}{y})^2}}$$
(11)
### IV.2 Generation of alternative compact structures

To generate the thirty randomly-collapsed structures used in the comparison of Fig. 8, we adopted a Monte Carlo technique. The length of the artificial proteins ranged uniformly in the interval 80-110. Starting from an open conformation, each structure was modified under the action of typical MC moves (single-bead, crankshaft, pivot) Sokal (1996). A newly generated modified configuration is accepted according to the ordinary Metropolis rule. The energy scoring function is composed of two terms: The first one contains a homopolymeric part that rewards the establishment of attractive interactions (cutoff of 6.0 Å) between any pair of non-consecutive residues. The second term is introduced to penalise structure realisations with radii of gyration larger than that found in naturally-occurring proteins with the same length. The Monte Carlo evaluation is embedded in a simulated annealing scheme Kirkpatrick et al. (1983) which allows to minimise efficiently the scoring function by slowly decreasing a temperature-like control parameter.

## V Acknowledgements

We are indebted with Amos Maritan for several illuminating discussions and with Fabio Cecconi and Alessandro Flammini for a careful reading of the manuscript. Support from INFM and MURST Cofin99 is acknowledged.

Figure captions

* Fig. 1. Correlation of cliquishness, relative and absolute contact order against folding rates of two-state folders. The values of the correlation coefficients are plotted as a function of the cutoff, $`d`$, used in the definition of the contact map.
* Fig. 2. Correlation of cliquishness, relative and absolute contact order against transition state placement of two-state folders. The values of the correlation coefficients are plotted as a function of the interaction cutoff, $`d`$.
* Fig. 3. Correlation of cliquishness, relative and absolute contact order against folding rates of two- and three-state folders. The values of the correlation coefficients are plotted as a function of the cutoff, $`d`$, used in the definition of the contact map.
* Fig. 4. Correlation of cliquishness, relative and absolute contact order against transition state placement of two- and three-state folders. The values of the correlation coefficients are plotted as a function of the interaction cutoff, $`d`$.
* Fig. 5. Scatter plot of cliquishness (left) and absolute contact order (right) versus folding rates of the 40 entries of Tables 1 and 2. The used values of $`d`$ are the optimal ones reported in the text. Filled circles, open squares and starred points denote proteins belonging to the $`\alpha `$, $`\alpha /\beta `$ and $`\beta `$ families, respectively.
* Fig. 6. Scatter plot of cliquishness (left) and absolute contact order (right) versus $`\theta _m`$ of the entries of Tables 1 and 2. The used values of $`d`$ are the optimal ones reported in the text. Filled circles, open squares and starred points denote proteins belonging to the $`\alpha `$, $`\alpha /\beta `$ and $`\beta `$ families, respectively.
* Fig. 7. Scatter plot of the logarithm of folding rates for the entries of Tables 1 and 2, against data from optimally combined cliquishness and contact order. The optimal linear superposition, see Methods, is obtained for $`b=0.7`$, ($`\{x\}`$ and $`\{y\}`$ being the cliquishness and contact order data respectively. Filled circles, open squares and starred points denote proteins belonging to the $`\alpha `$, $`\alpha /\beta `$ and $`\beta `$ families, respectively.
* Fig. 8. Scatter plot of average cliquishness versus absolute contact order, for randomly collapsed structures generated by stochastic numerical methods.