text:
# The 𝐵𝑒𝑝𝑝𝑜SAX High Energy Large Area Survey (HELLAS) - VI. The radio properties

## 1 Introduction

Radio follow-up of X-ray sources have played an important role in the optical identification program of the X-ray sources by providing position accurate to $``$ 1′′. This was particular true for the X-ray sources detected during the $`Einstein`$, ASCA and $`Beppo`$SAX X-ray missions for which the typical positional error is a circle of about 1 arcmin radius. A radio detection within the X-ray error box gives fundamental information on the position of the X-ray sources. In fact, since the majority of the bright extragalactic X-ray sources are associated with AGN, a radio source within the X-ray error box is physically associated to the X-ray source with a very high probability. For this reason, many of the optical identification programs of the X-ray sources have made an intensive use of the radio data (see for example Stocke et al. 1991 for the identification of the X-ray sources of the $`Einstein`$ Extended Medium Sensitivity Survey (EMSS) and Akiyama et al. 2000 for the identification of the X-ray sources in the ASCA Large Sky Survey (LSS) ).

Moreover, the radio data coupled with the optical and X-ray photometry allow us to compute the broad-band two point spectral indices $`\alpha _{ro}`$ and $`\alpha _{ox}`$ providing valuable information on the nature of the X-ray source population even in the absence of optical spectroscopy (Stocke et al. 1991).

A complete radio follow-up of X-ray selected samples is also an important tool to study the differences between AGNs that are strong radio-sources (radio-loud, RL) and those that are radio-quiet (RQ). Although the two classes have similar spectral index distributions (SEDs) outside the radio band (Elvis et al. 1994), their luminosity functions show differences in all the bands in which they have been studied (La Franca et al. 1994). In the optical band, using the PG sample of optically selected AGN, Padovani (1993) has shown that the shapes of the luminosity functions for RL and RQ are different. Similar results have been obtained by Della Ceca et al. 1994 and Ciliegi et al. 1995 studying the X-ray luminosity function (XLF) of RL and RQ separately. In particular Della Ceca et al. (1994) found a flattening of the XLF of the RL sample for $`L_x10^{44.5}`$ erg s-1. As a result the expected fraction of RL AGNs is a function of the X-ray flux limit in X-ray surveys. They predict that this fraction is $``$ 13 per cent for $`f_x(0.33.5keV)2\times 10^{13}`$ $`\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ and decreases to $``$2.5 per cent for $`f_x(0.33.5keV)2\times 10^{15}`$ $`\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ .

Radio follow-ups of X-ray selected samples seems to confirm this prediction. In fact, while shallow X-ray sample like the $`ASCA`$ Large Sky Survey (Akiyama et al. 2000) with a flux limit of $`f_x(210keV)1\times 10^{13}`$ $`\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ shows a fraction of RL around 10 per cent, deep ($`f_x(0.52.0keV)510\times 10^{15}`$ $`\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ ) ROSAT samples show a fraction of RL between 2 and 4 per cent (Ciliegi et al. 1995, de Ruiter et al. 1997, Zamorani et al. 1999).

In this paper we report the results of the radio follow-up of all the 147 X-ray sources detected by the $`Beppo`$SAX-MECS instrument in the framework of the High Energy LLarge Area Survey (HELLAS). This survey has observed about 85 deg2 of the sky in the 5-10 keV band down to a flux of 4-5 $`\times `$ 10-14 $`\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ . The whole survey and catalogue is described by Fiore et al. (2001), while the synthesis models for the X-ray background and the correlation with the soft X-rays have been investigated by Comastri et al. (2001) and Vignali et al. (2001). Finally, the spectroscopic identification of the HELLAS sources and the study of their evolution have been presented by La Franca et al. (2002, hereafter LF02).

## 2 Radio Observations

The 20 HELLAS sources with a declination further south than $``$40 deg have been observed with the Australia Telescope Compact Array (ATCA) while the 127 sources with DEC$`>`$40 deg have been observed with the Very Large Array (VLA). For these latter sources a complete covering at 20 cm down to the 5 $`\sigma `$ flux limit of 2.5 mJy is already available with the NRAO/VLA Sky Survey (NVSS, Condon et al. 1998) while the FIRST survey (Faint Images of the Radio Sky at Twenty centimeters, White et al. 1997) is available only for 27 HELLAS sources (5 $`\sigma `$ limit of $``$ 1 mJy).

In order to obtain information also on the radio spectral properties of the HELLAS sources we adopted the following strategy. All the 147 HELLAS sources have been observed at 6 cm down to a 1 $`\sigma `$ flux limit of $``$ 0.10-0.25 mJy. For the 20 HELLAS sources observed with the ATCA, we take advantage of the fact that the 6 and 3 cm receivers of the ATCA share a common feed-horn and we observed simultaneously also at 3 cm, obtaining a 3cm flux limit of $``$ 0.22 mJy (1 $`\sigma `$ level).

The wavelength of 6 cm and the flux limit reached (0.10-0.25 mJy) can be considered a good compromise between a deep radio survey and the necessity of avoiding strong contamination from spurious radio sources within the X-ray error box. The expected number of 6 cm sources is in fact N(S)=(0.42$`\pm `$0.05)(S/30)-1.18±0.19 where N(S) is the number of sources per arcmin2 with a flux density $`>`$S $`\mu `$Jy (Fomalont et al. 1991). Considering for the X-ray error box a circle of 1 arcmin radius (see below) and a 3 $`\sigma `$ limit of 0.3 mJy, we expect that only 0.1 radio source lies just by chance within the HELLAS error box. Going deeper in the radio flux (reaching, for example, a 3 $`\sigma `$ limit of 0.05 mJy) will increase the number of chance coincidence within the HELLAS error box to 0.7. On the other hand, at 20 cm the situation is more critical, since at 0.3 mJy the number of chance coincidence expected within an HELLAS error box is $``$ 0.3.

### 2.1 The ATCA observations

The ATCA observations of the 20 HELLAS sources were performed on June 1999. They were made with the ATCA simultaneously at two different frequencies : 4.848 and 8.640 GHz (referred to as 6 and 3 cm in the rest of the paper). The synthesized beam (full width at half power) is 2 arcsec at 6 cm and 1 arcsec at 3 cm. The primary flux density calibrator was PKS B1934$``$638, whose flux densities at different frequencies are incorporated directly in the calibration software. The data were calibrated and reduced using the ATCA reduction package MIRIAD (Multi-channel Image Reconstruction Image Analysis and Display). For each field a 512$`\times `$512 pixel image was constructed, with a pixel size of 1 arcsec at 6 cm and of 0.3 arcsec at 3 cm. The minimum root mean square (rms) noise obtained in each field is $``$ 0.25 mJy at 6 cm and $``$ 0.22 mJy at 3 cm (1 $`\sigma `$ level).

### 2.2 The VLA observations

The VLA observations were performed on 11 April 2000 at 4.885 GHz (6 cm) in C configuration. With this configuration and frequency, the synthesized beam size is $``$ 4 arcsec. All the data were analyzed with the NRAO AIPS reduction package. The data were calibrated using 3C286 as primary flux density calibrator. As for the ATCA data, for each field a 512$`\times `$512 pixel image ($``$ 8.5$`\times `$8.5 arcmin2) was constructed with a pixel size of 1 arcsec. For the majority of the fields the 1 $`\sigma `$ noise obtained in the central area is comparable to the expected one ($``$ 0.1 mJy). There are however some fields in which the noise is slightly higher due to problems during data acquisitions or due to the presence of a nearby strong radio source.

## 3 X-ray/radio associations

Using our 6 cm maps, we searched for radio sources within the X-ray error box of all the 147 HELLAS sources published in Fiore et al. 2001. The HELLAS error box has been assumed to be of 90 arcsec, in order to be absolutely conservative in the cross-correlation process, even though it must be noted that on average the $`Beppo`$SAX position of the HELLAS sources are better defined than 90 arcsec (see Appendix I of Fiore et al. 2001 for a detailed discussion of the position accuracy of the HELLAS sources). For the low number of high-Galactic latitude fields with neither a target nor a known X-ray source in the same field of view the error box has been assumed to be of 2 arcmin since no correction to the astrometry was possible (Vignali 2001).

For all the radio sources detected within the HELLAS error box we have computed the probability of random association with the X-ray source. Assuming that the radio sources belong to a Poissonian distributed population of sources,

$`P_{XR}=1e^{N(S)\pi d^2}`$
gives the probability to have a random association within a distance $`d`$ (distance between the X-ray position and the radio position of the possible counterpart) with a radio source having a flux greater than S. N(S) is the expected number of source with flux greater than that of the possible counterpart (see Section 2).

As first step, P$`{}_{XR}{}^{}<`$0.01 was chosen as a convenient threshold to accept a radio/X-ray association. We have 32 radio/X-ray associations with P$`{}_{XR}{}^{}<`$0.01. However, we note that this threshold represents only a starting point in the radio X-ray association process since further and more tightening constraints on the most likely associations come from the position of the optical counterparts of the HELLAS sources reported in LF02. Using a subsample of 118 HELLAS sources (enclosed in a region with DEC$`<`$79 deg and outside 5$`{}_{}{}^{h}<`$RA$`<6.5^h`$ and 17$`{}_{}{}^{h}<`$RA$`<20^h`$) they report the optical identification of 61 HELLAS : 37 Type 1 AGN, 9 Type 2 AGN, 5 narrow emission-line galaxies (ELG), 6 Cluster, 2 BL Lac, 1 Radio Galaxy and 1 Star. Of these 61 sources, 24 have a radio counterpart within 5 arcsec and with P$`{}_{RO}{}^{}<`$ 0.0002 (PRO has been calculated as PXR using the distance between the radio and optical position as value for $`d`$). Assuming that all the 24 radio/optical associations are real associations (as suggested by their low PRO values), we analyzed their PXR values. For 19 radio/optical associations we found P$`{}_{XR}{}^{}<`$0.01, 4 have 0.01$`<`$P$`{}_{XR}{}^{}<`$0.1 and one has P$`{}_{XR}{}^{}=`$0.217. The large value of PXR for the 5 sources with P$`{}_{XR}{}^{}>`$0.01 is due to the fact that for these sources the distance between the radio/optical position and the X-ray position is greater than 60 arcsec (see LF02 for a discussion on the procedure of the optical identification of the HELLAS sources). In the light of these results, we increase our threshold for the radio X-ray associations to P$`{}_{XR}{}^{}=`$0.1 ($`i.e.`$ we expect that about 10 per cent of the proposed identifications may be spurious positional coincidences). With this new threshold we found 53 X-ray/radio associations. We do not have X-ray sources with two (or more) radio counterparts both with P$`{}_{XR}{}^{}<`$0.1.

### 3.1 Optical counterparts of the X-ray/radio associations

Starting from the radio position of the 53 X-ray/radio associations, we searched for optical counterparts within 5 arcsec from the radio position using the optical positions of the 61 HELLAS sources identified by LF02, the USNO-A2.0 111http://archive.eso.org/skycat/servers/usnoa optical catalogue, the APM 222http://www.ast.cam.ac.uk/$``$apmcat/ optical catalogue and the Nasa Extragalactic Database (NED). As said above, 24 X-ray/radio associations have been identified with sources in LF02 (10 Type 1 AGN, 4 Type 2 AGN, 2 BL LAC, 3 Clusters, 4 ELGs and 1 Radio galaxy), 1 has been identified with a z=0.708 Radio galaxy in the Lockman Hole using NED (see table 2 source 116 in Lehmann et al. 2000 for a description of this source), 13 have an optical (R band) identification in the USNO and/or APM catalogue while 15 X-ray/radio associations do not have an optical identification brighter than R=20. The list of all the 147 HELLAS sources with their most probable optical and radio counterparts is shown in Table 1. The first three columns show the X-ray position, the offset in arcmin in the $`Beppo`$SAX field and the 5-10 keV X-ray flux in units of 10-13 $`\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ (from Fiore et al. 2001). Columns IV, V and VI give the R band magnitude (from LF02 if not otherwise specified), spectroscopic identification (1 for Type I AGN, 2 for Type II AGN, 3 for BL LAC, 4 for Cluster, 5 for ELG and 6 for Radio galaxies) and redshift from LF02. Columns VII, VIII, and IX give the radio flux at 3cm, 6 cm and 20 cm, respectively. All the 20 cm data are from the NVSS catalogue version 2.17 if not otherwise specified. For the X-ray sources without radio counterparts, we report at 3 and 6 cm the 3 $`\sigma `$ upper limit of our maps while at 20 cm we report the limit of the radio catalogue used (NVSS or FIRST). Column X gives the radio spectral index calculated using the available frequencies (mainly 6 and 20 cm). Column XI gives the position of the radio sources obtained from our 6 cm map (except when the source is detected only at 20 cm). The typical positional error of 6 cm sources is below 1 arcsec. Columns XII and XIII give the distance between the radio and the X-ray position and the probability $`P_{XR}`$ of a random association. Columns XIV and XV give the distance between the radio and the optical position and the probability $`P_{RO}`$ of a random association. Finally columns XVI and XVII give the two point spectral indices $`\alpha _{ro}`$ and $`\alpha _{ox}`$ defined as $`\alpha _{ro}=log(F_{5GHz}/F_{2500})/5.38`$ and $`\alpha _{ox}=log(F_{2500}/F_{2keV})/2.605`$ where $`F_{5GHz}`$ is the radio flux at 6 cm obtained from the present observation or by extrapolating the 20 cm flux density with a power-law slope with $`\alpha _r`$=0.7, $`F_{2500}`$ is the optical flux at 2500Å derived from the R magnitude with an optical slope of $`\alpha _o`$=0.5 while $`F_{2keV}`$ is the monochromatic X-ray flux at 2 keV derived from the measured $`Beppo`$SAX 5-10 keV flux and assuming an energy spectral index $`\alpha _x`$=0.6 (Fiore et al. 2001).

#### 3.1.1 Comparison with the optical identification reported in LF02

As said above, 24 of the 61 sources with an optical identification reported in LF02 have a radio counterpart within 5 arcsec. In addition to this information, it is interesting to note that no other radio sources have been found within their X-ray error box, as well as no radio sources have been detected in the X-ray error box of the remaining 37 sources. Moreover, among the 13 X-ray sources classified as empty fields down to R=21 in LF02, only two (0134$``$3006 and 1304$``$0533) have a radio counterpart (see Table 1). These two radio sources have not an optical identification brighter than R=21 within 5 arcsec from their radio position, in agreement with their classification as empty field.

These results, under the hypothesis that a radio source within the X-ray error box is physically associated to the X-ray source with a very high probability, strengthen the goodness of the optical identification reported in LF02.

### 3.2 Radio-Loud and Radio-Quiet Classification

The two point spectral index $`\alpha _{ro}`$ has been commonly used to discriminate between radio-loud and radio-quiet objects. We were able to calculate the $`\alpha _{ro}`$ value for a total of 84 sources: 55 HELLAS sources identified in LF02 (we did not calculated the two point spectral indicies for the 6 clusters due to the indetermination of their optical magnitude), for the z=0.708 Radio galaxy in the Lockman Hole (see above), for the 13 X-ray/radio associations that have been identified with an USNO and/or an APM source and for the 15 X-ray/radio associations without an optical counterpart (in the latter case using an optical upper limit of R=20). The $`\alpha _{ro}`$ distribution of these 84 sources is shown in Figure 1 (sources detected both in radio and optical band: shaded histogram; 3$`\sigma `$ upper limits (X-ray sources with an optical counterpart but without a radio counterpart): open solid line histogram; 3$`\sigma `$ lower limits (X-ray sources with a radio counterpart but without an optical counterpart): open dashed line histogram).

The values of $`\alpha _{ro}`$ used in literature to discriminate between radio-loud and radio-quiet objects are 0.20 (see Elvis et al. 1994, Giommi, Menna & Padovani 1999) or 0.35 (Zamorani et al. 1981, Della Ceca et al. 1994). To choose between these two values, we applied to our $`\alpha _{ro}`$ distribution (Figure 1) the Lee statistic (see Fitchett & Merrit 1988 and Lee 1979 for details). The Lee statistic is a test for a bimodality of a distribution. The location of the maximum corresponds to the value of the variable which best separates the two subsamples. The results of the Lee statistic are shown in Figure 2. We find that the Lee parameter reaches it maximum for $`\alpha _{ro}`$0.35 - 0.36. This results has been obtained using only the distribution of the sources detected both in radio and optical band (shaded histogram in Figure 1). However similar results are obtained treating the upper and lower limits as detections. Consequently the value of $`\alpha _{ro}`$=0.35 has been used to separate the RL from the RQ objects.

We have 18 RL sources with $`\alpha _{ro}0.35`$ (11 have been detected both in the radio and optical bands and 7 have a lower limit on $`\alpha _{ro}`$ greater than 0.35), 58 RQ sources (24 detected both in the radio and optical bands, 34 with an $`\alpha _{ro}`$ upper limit lower than 0.35) while for 8 sources it was impossible to give a reliable classification due to the fact that they have an $`\alpha _{ro}`$ lower limit lower than 0.35. However, given the fact that their $`\alpha _{ro}`$ lower limit are close to the threshold of 0.35 (the lowest value is $`\alpha _{ro}>`$0.23, see Table 1 and the open dashed line histogram in Figure 1) we assumed these 8 sources to be RL objects. Over 84 sources analyzed we have therefore 26 RL and 58 RQ objects. For the remaining 63 HELLAS sources it was impossible to calculate the two point spectral index $`\alpha _{ro}`$ (or at least a lower or an upper limit) due to the lack of a radio and optical identification, $`i.e.`$ for these sources only an X-ray detection and a radio upper limit is available. However also the upper limit on their radio-to-X-ray ratio can be used to obtain useful informations on the radio nature of these source. In Figure 3 (top panel) we report the distribution of the radio-to-X-ray ratio for the RQ and RL objects in the atlas of quasars published by Elvis et al. (1994) while in the the other two panels we report the distribution of the 84 HELLAS sources for which it was possible to obtain a classification as RQ and RL object (middle panel) and the distribution of the radio-to-X-ray ratio upper limits for the 63 HELLAS sources without a classification in RL and RQ (bottom panel). Although the atlas of 47 quasars from Elvis et al. (1994) is an heterogeneous sample of objects selected mainly from the PG (Schmidt and Green 1983), 3C and Parkes catalogues, it gives a useful reference for the radio-to-X-ray ratio values for RQ and RL objects. As shown in figure 3 the RQ and RL objects show a significantly different distribution of their radio-to-X-ray ratio. A Kolmogorov-Smirnov (KS) test of the radio-to-X-ray ratio distributions for the HELLAS RQ and RL (reported in the middle panel of figure 3) shows that the two distributions are drawn from the same population with a probability P$`{}_{KS}{}^{}<5\times 10^4.`$ The same result has been obtained performing a KS test between the distribution of the RL and that of the 63 sources for which only and X-ray detection and a radio upper limit is available. Viceversa, a KS test between this latter and the distribution of the HELLAS RQ objects indicate that the hypothesis that the two distribution are drawn from the same population can not be excluded (PKS=0.69). It is therefore reasonable to assume that all the 63 sources for which only and X-ray detection and a radio upper limit is available are radio-quiet objects. Under this hypothesis, the fraction of radio-loud objects in the HELLAS sample is $``$18 per cent (26 RL over a total of 147 sources).

## 4 Discussion

### 4.1 Radio, Optical and X-ray properties

To examine the general properties of the radio/X-ray associations, in Figure 4 we show the 6 cm radio flux versus the R-band magnitude for the 53 radio/X-ray associations of the HELLAS sample. Superimposed are the lines corresponding to constant values for the two point spectral index $`\alpha _{ro}`$. The dot dashed line shows the value of $`\alpha _{ro}`$=0.35 used to separate RL and RQ objects. In Figure 5 we plot the radio versus the hard (5-10 keV) X-ray luminosities. The latter has been corrected for absorption using the hardeness ratio reported in Fiore et al. (2001).

It is clear from Figures 4 and 5 that the radio-loud objects are associated mainly with Type 1 AGNs with L$`{}_{510keV}{}^{}\genfrac{}{}{0pt}{}{_>}{^{}}10^{44}`$ erg/s, while all the Type 2 AGNs and Emission Line Galaxies are radio-quiet sources with $`\alpha _{ro}`$ $`\genfrac{}{}{0pt}{}{_<}{^{}}`$ 0.15 and L$`{}_{510keV}{}^{}\genfrac{}{}{0pt}{}{_<}{^{}}10^{44}`$ erg/s.

The lack of low luminosity radio-loud sources is well explained by the dependence of the radio-loud fraction to the X-ray and optical luminosity. The fraction of the radio-loud AGN changes, in fact, from $``$18% at L$`{}_{x}{}^{}4\times 10^{44}`$ erg/s to $``$2% at L$`{}_{x}{}^{}2\times 10^{43}`$ in X-ray selected samples (Della Ceca et al. 1994) and from $``$20%-50% at M$`{}_{B}{}^{}`$24.5 to $``$7%-8% at fainter absolute magnitude in optically selected samples (Padovani 1993, La Franca et al. 1994, Goldschmidt et al. 1999).

The lack of high X-ray luminosity Type 2 AGN is also a well know effect. Only few examples have been reported prior the $`Chandra`$ mission (see, for example, Almaini et al. 1995, Boyle et al. 1998, Franceschini et al. 2000). However, even in the deeper $`Chandra`$ and XMM surveys only a couples of sources considered as the prototype of Type 2 QSO (L$`{}_{210}{}^{}>10^{44}`$, $`N_H>10^{23}`$ cm-2, optical emission lines with FWHM $`\genfrac{}{}{0pt}{}{_<}{^{}}`$ 2000 km/s) have been detected (Norman et al. 2002, Stern et al. 2002). Their space density will be better constrained by the forthcoming deep X-ray ($`Chandra`$ and XMM) and far-infrared surveys (SIRTF).

Finally, in Figure 6 we show the $`\alpha _{ox}`$-$`\alpha _{ro}`$ plot for the 84 HELLAS sources for which we were able to calculate $`\alpha _{ox}`$ and $`\alpha _{ro}`$. From this figure it is evident that the RQ Type 1 AGNs and the Type 2 AGNs occupy the same region in the plane $`\alpha _{ox}`$-$`\alpha _{ro}`$ and then the two classes cannot be distinguished only on the basis of their multiband photometric data. On the other hand, the region with $`\alpha _{ox}`$ $`\genfrac{}{}{0pt}{}{_<}{^{}}`$ 1.5 and $`\alpha _{ro}\genfrac{}{}{0pt}{}{_>}{^{}}`$ 0.35 is occupied only by Type 1 AGN. We would therefore conclude that all the unidentified sources with $`\alpha _{ro}`$ $`\genfrac{}{}{0pt}{}{_>}{^{}}`$ 0.30 - 0.35 and $`\alpha _{ox}\genfrac{}{}{0pt}{}{_<}{^{}}`$ 1.5 are likely to be Type 1 AGN.

### 4.2 Comparison with other X-ray samples

#### 4.2.1  The $`EinsteinObservatory`$, ROSAT and ASCA surveys

Previous radio follow-up of X-ray selected samples have shown different percentages of radio/X-ray associations. Using 20 cm observations of ROSAT X-ray selected sample, Ciliegi et al. (1995), de Ruiter et al. (1997) and Zamorani et al. (1999) found only $``$ 10% of radio/X-ray coincidences, while Stocke et al. (1991) found a significantly higher fraction (about 27%) for the X-ray sources in the $`Einstein`$ Extended Medium Sensitivity Survey (EMSS). For the HELLAS sources we have 53 radio/X-ray associations over a total of 147 sources. From the statistical analysis of the radio/X-ray associations (see Section 3), we estimate that about 10% ($`i.e.`$ $``$ 5 sources) of the proposed associations may be spurious and therefore the percentage with the “correct” radio/X-ray associations is in the range $``$33-36%. This value is comparable to that obtained by Akiyama et al. (2000) in the course of the $`ASCA`$ Large Sky Survey ($``$ 35%), one of the few large area samples selected in the hard X-ray band and with a complete follow-up in the radio and optical bands.

In Table 2 we report a summary of different X-ray surveys for which a complete radio follow-up is available. All these surveys are “AGN dominated” in the sense that the large majority of their identification ($``$75%) are associated with AGNs. For each survey we report the X-ray and radio limits, the radio-to-X-ray limit ratio fr/f2keV assuming an X-ray spectral index $`\alpha _x`$=1.0, the percentage of radio/X associations (number of radio/X-ray association over the total number of X-ray sources in the survey) and the fraction of Radio-Loud objects in the radio/X-ray associations. As clearly shown in the table, the percentage of radio/X-ray associations is a function, as expected, of fr/f2keV: the lower this ratio is ($`i.e`$ deeper radio data in comparison to the X-ray flux limit) the higher is the fraction of X-ray objects with a radio counterpart. The top panel of Figure 3 suggests (as proposed also by Akiyama et al. 2000) that this change in the percentage of radio/X-ray associations is due to the fact that when fr/f2keV becomes very low we start seeing also radio-quiet AGNs. To check if this is the case, in Figure 7 we plot the 6 cm radio flux as function of the two point spectral index $`\alpha _{ro}`$ for the HELLAS sources and for the 12 radio/X-ray associations found in the ASCA LSS (using a radio spectral index $`\alpha _r`$=0.7 to convert the radio fluxes from 20 cm to 6 cm). The solid line shows the limit beyond which a source can not be detected in the optical band due to the magnitude limit of the optical data. The line in Figure 7 has been drawn using R=21. As clearly shown in figure, going deeper in the radio data we start to see the radio-quiet population ($`\alpha _{ro}<`$0.35). This population is almost completely absent at radio flux level greater than 10-20 mJy but becomes the dominant population at radio fluxes lower than $``$1 mJy. From Figure 7 it is also evident that deep optical data are needed in order to detect in the optical band radio-loud sources with radio fluxes lower than a few mJy.

Moreover it is interesting to note that the 3 X-ray surveys with the highest percentage of radio/X-ray associations (EMSS, LSS and HELLAS) are all selected in the medium or hard X-ray band (see Table 2), while the 3 X-ray surveys with the lowest percentage of radio/X-ray associations (Lockman, Marano and CRSS) are all selected in the soft X-ray band. In addition to the lowest radio-to-X-ray flux ratio available in the EMSS, LSS and HELLAS surveys, the higher percentage of radio/X-ray associations could also be due to the selection in the hard X-ray band. In fact, since the radio emission is transparent against obscuring material, obscured AGNs are detected in an unbiased way in the radio wavelength as well as in the hard X-ray band, but they could be missed in samples selected in the soft X-ray band where obscuration effects are much more relevant. To test this hypothesis we calculate the fraction of radio/X-ray associations in the HELLAS and LSS surveys using a radio flux limits of 3.7 mJy and 3.2 mJy respectively, in order to have an fr/f2keV similar to that of the soft X-ray selected samples (see Table 2). We found a percentage of radio/X-ray associations around 20%, i.e. still higher than that found in soft X-ray selected samples. In conclusion the higher fraction of radio/X-ray associations found in the EMSS, LSS and HELLAS surveys in comparison to the soft X-ray selected samples is probably due to a combination of two effects: deeper radio data in comparison to the X-ray flux limit and the selection in harder X-ray bands.

#### 4.2.2 The deep $`Chandra`$ surveys

Recently new deep X-ray samples have been obtained thanks to the revolutionary capabilities of the $`Chandra`$ X-ray satellite. The deep surveys (Brandt et al. 2001, Rosati et al. 2002) have reached flux limits at least two orders of magnitude lower than the fainter samples listed in Table 2. Extensive programs of optical identifications showed that most of the sources in these samples are not AGN. The bulk of the optical identifications in these deep X-ray surveys is associated with normal galaxies, the majority of which show strong emission lines indicative of star formation activity (Barger et al. 2001, Brandt et al. 2001,Bauer et al. 2002). However, the most surprising finding is the discovery of a sizeable number of relatively bright X-ray sources spectroscopically identified with early-type “normal” galaxies without any signature of nuclear activity in the optical spectra (Fiore et al. 2000, Barger et al. 2001). These optically normal X-ray luminous galaxies were thought to be relatively rare, unusual objects, but during the recent deep X-ray surveys obtained with $`Chandra`$ and XMM we saw a remarkable increase of their number (Comastri et al. 2002).

The radio properties of these new deep X-ray samples have been studied using very deep VLA radio survey (down to few micro Jy at 1.4 GHz) obtained in the SSA13 field (Barger et al. 2001) and in the HDF-N (Bauer et al. 2002). These studies showed a different behavior of these new samples in comparison to the “AGN-dominated” samples. In fact, in the new deep X-ray surveys a very high percentage of radio/X-ray associations (up to 80% in the SSA13 field, see Barger et al. 2001) has been found, with the highest matching fraction found among the emission line galaxies, which are comprised of apparently normal and starburst galaxies at redshift of z$``$0.1-1.3 thought to be undergoing recent star formation. The nature of the radio/X-ray associations is therefore completely different between that found in the “AGN dominated” samples reported in Table 2 and that found in the new deep X-ray surveys. While the former are associated with an AGN activity, the latter are mainly associated with a star forming activity.

### 4.3 The radio spectral index

Using all the radio data available at 3, 6 and 20 cm we were able to calculate the radio spectral index (or an upper/lower limit when the source is detected in only one radio band) for a total of 50 HELLAS sources. Since the majority of the spectral indices have been calculated using the 6 cm fluxes obtained with our observations and the 20 cm fluxes reported in the public catalogues (NVSS and FIRST), their values should be treated as approximate since the images used to measure the fluxes were not matched in resolution. We also caution that variability among the radio AGN is a legitimate concern and could lead to inaccurate spectral index estimates since the two radio bands were observed several years apart.

All the radio spectral index $`\alpha _r`$ are reported in column X of Table 1, while in Figure 8 we plot $`\alpha _r`$ as function of the two point spectral index $`\alpha _{ro}`$. The horizontal dashed line denotes the typical spectral index, $`\alpha _r`$=0.5, separating steep and flat spectrum radio sources. Steep spectral indices ($`\alpha _r>`$0.5) often indicate radio emission from star formation or lobe-dominated AGN, while flat ones ($`\alpha _r<`$0.5) often indicate emission from core-dominated AGN (e.g. Kellermann and Owen 1988).

A statistical analysis using both the measured values and the limits to $`\alpha _r`$ has been performed with the software package ASURV, which implements the methods described by Feigelson & Nelson (1985) and Isobe, Feigelson & Nelson (1986). We divided our sources according to their optical classification studying the mean radio spectral index of sources with broad emission lines (AGN1, open circles in Figure 8) and with narrow emission lines (AGN2 plus Emission Line Galaxies, filled symbols in Figure 8). We found that while the Type 1 AGNs have a mean radio spectral index of $`\alpha _r=0.25\pm `$0.10, the Type 2 AGNs plus the Emission Line Galaxies have a steeper radio spectral index $`\alpha _r=0.69\pm `$0.11. This result is in agreement with the idea that the radio emission from Type 1 AGNs is core dominated with a flat radio spectral index due to the self-absorption process. On the other hand, in Type 2 AGN and Emission Line Galaxies the nonthermal radio emission (either due to a compact nucleus or to the integrated emission arising from the supernovae remnants in a starburst region) takes place on larger physical scale without self-absorption. For the case of the star formation these scales correspond approximately to 0.1-1.0 Kpc as observed in the local ELGs population (Condon 1989, Condon 1992) while for the core dominated sources these scales correspond to few parsec as well-know from many different radio surveys (see Nagar et al. 2000,2002 and Giovannini et al. 2001 for the most recent works). Under this hypothesis, we expect a more extended radio emission from the HELLAS Type 2 AGNs and ELGs in comparison to the radio emission from the HELLAS Type 1 AGNs. Sub-arcsecond radio observations are needed to test this hypothesis.

### 4.4 Radio Quiet and Radio Loud AGNs in the HELLAS spectroscopically identified subsample

Starting from the 61 spectroscopically identified sample of 118 HELLAS sources published by LF02 we selected a subsample of 46 AGNs (37 Type 1 + 9 Type 2) brighter than 5 $`\times `$ 10-14 $`\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ in the 5-10 keV band and brighter than R=21 and R=19 (for Type 1 and Type 2 AGNs) in the optical band. Radio data are available for all the 46 HELLAS AGNs. Using the two point spectral index $`\alpha _{ro}`$=0.35 we divided the sources in RL and RQ. We have 6 ($``$ 13%) RL sources, all classified as Type 1 AGN.

This fraction of RL AGNs is well in agreement with the predictions of Della Ceca et al. (1994). Using the de-evolved X-ray Luminosity function for RL and RQ AGNs they predicted a fraction of RL AGNs of $``$10% for an X-ray flux limit of $`5\times 10^{14}`$ $`\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ (see their Figure 11).

## 5 Summary and Conclusion

We present the results of a complete radio follow-up down to a 6 cm flux limit of about 0.3 mJy (3 $`\sigma `$) obtained with the VLA and ATCA radio telescopes of all the 147 X-ray sources detected in the $`Beppo`$SAX HELLAS survey. Our major conclusions are as follows.

$``$ From a statistical analysis based on the X-ray, optical and radio position, we found 53 X-ray/radio associations over a total of 147 sources. We estimated that about 10% ($`i.e.`$ $``$ 5 sources) of the proposed associations may be spurious and therefore the percentage of the “real” X-ray/radio associations is in the range $``$33-36%. This percentage is in agreement to that obtained in other X-ray surveys selected in the hard X-ray band and is significantly greater than the percentage found in surveys selected in the soft X-ray band. This is probably due to a combination of two effects: the availability of deeper radio data in comparison to the X-ray flux limit and a selection in harder X-ray band. Obscured AGNs are in fact missed in soft X-ray selected samples, while they are detected in the radio and hard X-ray bands since both these wavelengths are transparent against obscuring material. This change in the percentage and nature of the radio/X-ray associations that we start to see in the HELLAS survey is much more evident in the new deep X-ray $`Chandra`$ samples. Radio observations at the micro Jy level of the $`Chandra`$ samples showed, in fact, that the radio/X-ray associations in these samples are dominated by normal galaxies and starburst at z$``$ 0.1-1.3.

$``$ The value of $`\alpha _{ro}=0.35`$ used to separate the radio-loud and radio-quiet sources has been estimated using the Lee statistic. Of the 53 HELLAS radio/X-ray associations, 26 have been classified as radio-loud sources. From the analysis of the ratio between the radio and the X-ray fluxes, we have assumed that all the 63 sources for which it was impossible to calculate the two point spectral index $`\alpha _{ro}`$ (and than to obtain a classification in RQ or RL sources) are likely to be radio-quiet objects. Under this hypothesis, the fraction of radio-loud objects in the HELLAS sample is $``$18 per cent.

$``$ The analysis of the multiband photometric data based on the two point spectral indices $`\alpha _{ro}`$ and $`\alpha _{ox}`$ has shown that the identified radio-loud objects are associated mainly with Type 1 AGN with L$`{}_{510keV}{}^{}\genfrac{}{}{0pt}{}{_>}{^{}}10^{44}`$ erg/s, while the identified Type 2 AGNs and Emission Line Galaxies are radio-quiet sources with L$`{}_{510keV}{}^{}\genfrac{}{}{0pt}{}{_<}{^{}}10^{44}`$ erg/s. RQ Type 1 AGNs and Type 2 AGNs cannot be distinguished only on the basis of their photometric data. On the basis of these results we would conclude that all the unidentified sources with $`\alpha _{ro}`$ $`\genfrac{}{}{0pt}{}{_>}{^{}}`$ 0.30 - 0.35 and $`\alpha _{ox}\genfrac{}{}{0pt}{}{_<}{^{}}`$ 1.5 are likely to be Type 1 AGN.

$``$ The analysis of radio spectral index has shown that Type 1 AGNs have a mean spectral index flatter than Type 2 AGNs and Emission Line Galaxies. This result is in agreement with the idea that the core-dominated radio emission from Type 1 AGNs is self-absorbed, while in AGN2 and Emission Line Galaxies the radio emission take place on larger physical scale, without self-absorption.

$``$ Finally, using a subsample of 46 AGNs (37 Type 1 and 9 Type 2, see section 4.4) we studied the fraction of the RL and RQ AGNs. We have 6 ($``$ 13%) RL sources , all classified as Type 1 AGN. This fraction of RL AGNs is well in agreement with the results of Della Ceca et al. 1994 which predicted a fraction of RL AGNs of $``$10% for an X-ray flux limit of $`5\times 10^{14}`$ $`\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ .

Acknowledgements

This paper is based on observations collected at the Very Large Array (VLA) Radio Telescope and at the Australia Telescope Compact Array (ATCA). The VLA is a facility of National Radio Astronomy Observatory (NRAO) which is operated by Associated Universities, Inc., under cooperative agreement with the National Science Foundation. The ATCA is part of the Australia Telescope which is funded by the Commonwealth of Australia for operation as a National Facility managed by CSIRO. This work was supported by the Italian Ministry for University and Research (MURST) under grant COFIN-00-02-36 and by Italian Space Agency ASI contracts I/R/113/01 and I/R/073/01. CV also acknowledges the NASA LTSA grant NAG5-8107 for financial support.


format: text

text:
# On Design of Optimal Nonlinear Kernel Potential Function for Protein Folding and Protein Design

## 1 Introduction

Potential function plays critical roles in computational studies of protein folding, protein structure prediction, and protein sequence design . A variety of empirical potential functions have been developed for coarse-grained models of proteins, where amino acid residues are not represented at atomic details. One prominent class of potentials are knowledge-based potentials derived from statistical analysis of database of protein structures . In this class of potentials, the interaction between a pair of residues are estimated from its relative frequency in database when compared with a reference state or a null model. This approach has been successfully applied in fold recognition, in threading, and in many other studies . The drawback of this class of potential is that there are several conceptual difficulties. These include the neglect of chain connectivity in the reference state, and the problematic implicit assumption of Boltzmann distribution . An alternative approach is to find a set of parameters such that the potential functions are optimized by some criterion, e.g., maximized energy difference between native conformation and a set of alternative (or decoy) conformations . This approach has been shown to be very effective in recognizing native structures from alternative conformations  and in folding membrane proteins . However, if a large number of native protein structures are to be simultaneously stabilized against a large number of decoy conformations, no such potential functions can be found .

There are three key steps in developing effective empirical potential function using optimization: (1) the functional form of the potential, (2) the generation of a large set of decoys for discrimination, and (3) the optimization techniques. The initial step of choosing an appropriate functional form so far has been straightforward. Empirical pairwise potentials are usually all in the form of weighted linear sum of interacting residue pairs (see reference  for an exception). In this functional form, the weight coefficients are the parameters of the potential, which are optimized for discrimination. The same functional form is also used in statistical potential, where the weight coefficients are derived from database statistics. For the task of decoy generation, an efficient method to obtain millions of decoys is gapless threading, i.e., a decoy conformation can be obtained by mounting the sequence of the native protein to different parts of the conformation of a larger protein of an unrelated structure . Alternatively, a large number of challenging decoys can be generated by using either a chain growth method  or a protein structure refinement method , both are technically more complex. The optimization techniques that have been used include perceptron learning and linear programming , and analytical solution has also been obtained when assumptions about the distribution of data are made . The objectives of optimization are often maximization of energy gap between native conformation and the average of decoy conformation, or energy gap between native and decoys with lowest energy, or the $`z`$-score of the native structure .

In this study, we develop an alternative formulation of protein potential function in the form of mixture of nonlinear Gaussian kernel functions. We also use a different optimization technique based on quadratic programming. Instead of maximizing the energy gap, here a function related to bounds of expected classification errors is optimized . In addition, we explore the relationship between protein folding and protein design. We study a simplified version of the protein folding problem. Our goal is to identify the native protein structure from an ensemble of alternative or decoy structures for a given amino acid sequence. We also study a simplified version of the protein design problem. Our goal of protein sequence design is to identify a protein sequence that is most compatible with a given three-dimensional coarse-grained structure. In this study, we do not address the problem of how to generate candidate conformation or candidate sequence by searching either the conformation space or the sequence space.

Experimentation with the nonlinear function developed shows that it can discriminate simultaneous 440 native proteins against 14 million structure decoys generated by gapless threading. Results of similar experiments with perceptron learning was negative, as reported in literature . We also test our potential function for protein design. Using the same training set of 440 proteins and a corresponding set of 14 millions sequence decoys, we succeed in developing a nonlinear function that correctly classifies all 440 native sequences. In contrast, we cannot obtain a weighted linear sum potential using the state-of-the-art interior point solver of linear programming method as reported in , such that it is capable of classifying perfectly all the native sequences. We also perform blind tests for both native structure and native sequence recognition. Taking 194 proteins unrelated to the 440 training set proteins, the nonlinear potential achieves a success rate of 92.8% and 98.4% in sequence design and in structure recognition, respectively. Both results compare favorably with optimal linear potential and statistical potential.

The rest of the paper is organized as follows. We first describe theory and model of linear and nonlinear function, including the kernel model and the optimization technique. We then explain details of computation. We further describe experimental results of learning and results of blind test. We conclude with discussion.

## 2 Theory and Models

### Modeling Protein Folding Potential.

To model protein computationally, we first need a method to describe its geometric shape and its sequence of amino acid residues. Frequently, a protein is represented by a $`d`$-dimensional vector $`𝒄R^d`$. For example, we can represent a protein as a vector $`𝒄R^d`$, $`d=20`$, by measuring the solvent accessible surface areas of each of the 20 residue types. A method that is widely used is to count nonbonded contacts of various types of amino acid residue pairs in a protein structure. In this case, the count vector $`𝒄R^d,d=210`$, is used as the protein descriptor. Once the structural conformation of a protein $`𝒔`$ and its amino acid sequence $`𝒂`$ is given, the protein description $`f:(𝒔,𝒂)R^d`$ will fully determine the $`d`$-dimensional vector $`𝒄`$. In the case of contact vector, $`f`$ corresponds to the mapping provided by specific contact definition, e.g., two residues are in contact if their distance is below a specific cut-off threshold distance.

Based on the classical experiments of Anfinsen , a fundamental requirement for protein folding is that the native structure $`𝒔_N`$ with native amino acid sequence $`𝒂_N`$ must have an energy $`H(f(𝒔_N,𝒂_N))`$ that is the lowest among a set of alternative structure called decoys $`𝒟=\{(𝒔_D,𝒂_N)\}`$, where the native amino acid sequence $`𝒂_N`$ takes a decoy conformation $`𝒔_D`$ that is different from the native conformation $`𝒔_N`$ :

$$H(f(𝒔_N,𝒂_N))<H(f(𝒔_D,𝒂_N))\text{for all }(𝒔_D,𝒂_N)𝒟$$

Sometimes we can further require that the energy difference must be greater than a constant $`b>0`$:

$$H(f(𝒔_N,𝒂_N))+b<H(f(𝒔_D,𝒂_N))\text{for all }(𝒔_D,𝒂_N)𝒟$$

A widely used functional form for protein potential function $`H`$ is the weighted linear sum of pairwise contacts . The linear sum energy $`H`$ is then:

$$H(f(𝒔,𝒂))=H(𝒄)=𝒘𝒄.$$
(1)
As soon as the weight vector $`𝒘`$ is specified, the potential function is fully defined. For such linear potentials, the basic requirement for protein potential is then:

$$𝒘(𝒄_N𝒄_D)<0$$

or

$$𝒘(𝒄_N𝒄_D)+b<0$$

if we require that the energy difference between a native structure and a decoy must be greater than a real value $`b`$.

### Two Geometric Views of Linear Protein Folding Potentials.

There is a natural geometric view of the inequality requirement for weighted linear sum potentials. A useful observation is that each of the inequalities divides the space of $`R^d`$ into two halfs separated by a hyperplane. The hyperplane is defined by the normal vector $`(𝒄_N𝒄_D)`$ and its distance $`b/𝒄_N𝒄_D`$ from the origin. The weight vector $`𝒘`$ must be located in the half-space opposite to the direction of the normal vector $`(𝒄_N𝒄_D)`$ (Fig 1a). This half-space can be written as $`𝒘(𝒄_N𝒄_D)+b<0`$.

When there are many inequalities to be satisfied simultaneously, the intersection of the half-spaces forms a convex polyhedron . If the weight vector is located in the polyhedron, all the inequalities are satisfied. Potentials with such weight vector $`𝒘`$ can discriminate the native protein from the set of all decoys.

For each native protein $`i`$, there is one convex polyhedron $`𝒫_i`$ formed by the set of inequalities associated with its decoys. If our potential can discriminate simultaneously $`n`$ structures from a union of sets of decoys, the weight vector $`𝒘`$ must be located in a smaller convex polyhedron $`𝒫`$ that is the intersection of the $`n`$ convex polyhedra:

$$𝒘𝒫=\underset{i=1}{\overset{n}{}}𝒫_i.$$

There is yet another geometric view of the inequality requirements. The relationship $`𝒘(𝒄_N𝒄_D)+b<0`$ for all decoys and native protein structures can be regarded as a requirement that all points $`\{𝒄_N𝒄_D\}`$ are located on one side of a hyperplane, which is defined by its normal vector $`𝒘`$ and its distance $`b/𝒘`$ to the origin (Fig 1b). We can show that such a hyperplane exists if the origin is not contained within the convex hull of the set of points $`\{𝒄_N𝒄_D\}`$ (see Appendix).

The second geometric view is dual and mathematically equivalent to the first geometric view. In the first view, a point $`𝒄_N𝒄_D`$ determined by the structure-decoy pair $`c_N=(𝒔_N,𝒂_N)`$ and $`c_D=(𝒔_D,𝒂_N)`$ corresponds to a hyperplane representing an inequality, a solution weight vector $`𝒘`$ corresponds to a point located in the final convex polyhedron. In the second view, each structure-decoy pair is represented as a point $`𝒄_N𝒄_D`$ in $`R^d`$, and the solution weight vector $`𝒘`$ is represented by a hyperplane separating all the points $`𝒞=\{𝒄_N𝒄_D\}`$ from the origin.

### Optimal Linear Potentials.

Several optimization methods have been applied to find the weight vector $`𝒘`$. The Rosenblantt perceptron method works by iteratively updating an initial weight vector $`𝒘_0`$ . Starting with a random vector, e.g., $`𝒘_0=\mathrm{𝟎}`$, we test each native protein and its decoy structure. Whenever the relationship $`𝒘(𝒄_N𝒄_D)+b<0`$ is violated, we update $`𝒘`$ by adding to it a scaled vector $`\eta (𝒄_N𝒄_D)`$. The final weight vector is therefore a linear combination of protein and decoy count vectors:

$$𝒘=\eta (𝒄_N𝒄_D)=\underset{N𝒩}{}\alpha _N𝒄_N\underset{D𝒟}{}\alpha _D𝒄_D.$$
(2)
Here $`𝒩`$ is the set of native proteins, and $`𝒟`$ is the set of decoys. The set of coefficients $`\{\alpha _N\}\{\alpha _D\}`$ gives a dual form representation of the weight vector $`𝒘`$ as an expansion of the training examples, including both native and decoy structures.

If the final convex polyhedron $`𝒫`$ is non-empty, there can be infinite number of choices of $`𝒘`$, all with perfect discrimination. But how do we find a weight vector $`𝒘`$ that is optimal? This depends on the criterion for optimality. For example, one can choose the weight vector $`𝒘`$ that minimizes the variance of energy gaps between decoys and natives: $`\mathrm{arg}_𝒘\mathrm{min}\frac{1}{|𝒟|}\left(𝒘(c_Nc_D)\right)^2\left[\frac{1}{|𝒟|}_D\left(𝒘(𝒄_N𝒄_D)\right)\right]^2`$ as used in reference , or minimizing the $`Z`$-score of a large set of native proteins, or minimizing the $`Z`$-score of the native protein and an ensemble of decoys , or maximizing the ratio $`R`$ between the width of the distribution of the energy and the average energy difference between the native state and the unfolded ones . A series of important works using perceptron learning and other optimization techniques  showed that effective linear sum potentials can be obtained.

Here we describe yet another optimality criterion according to the second geometric view. We can choose the hyperplane $`(𝒘,b)`$ that separates the points $`\{𝒄_N𝒄_D\}`$ with the largest distance to the origin. Intuitively, we want to characterize proteins with a region defined by the training set points $`\{𝒄_N𝒄_D\}`$. It is desirable to define this region such that a new unseen point drawn from the same protein distribution as $`\{𝒄_N𝒄_D\}`$ will have a high probability to fall within the defined region, and non-protein points following a different distribution, which is assumed to be centered around the origin when no a priori information is available, will have a high probability to fall outside the defined region. In this case, we are more interested in modeling the region or support of the distribution of protein data, rather than estimating its density distribution function. For linear potential, regions are half-spaces defined by hyperplanes, and the optimal hyperplane $`(𝒘,b)`$ is the one with maximal distance to the origin. This is related to the novelty detection problem and single-class support vector machine studied in statistical learning theory . In our case, any non-protein points will need to be detected as outliers from the protein distribution characterized by $`\{𝒄_N𝒄_D\}`$. Among all linear functions derived from the same set of native proteins and decoys, the optimal weight vector $`𝒘`$ is likely to have the least amount of mislabellings. The optimal weight vector $`𝒘`$ can therefore be found by solving the following primal quadratic programming problem:

Minimize $`{\displaystyle \frac{1}{2}}𝒘^2`$
subject to $`𝒘(𝒄_N𝒄_D)+b<0\text{ for all }N𝒩\text{ and }D𝒟.`$
The solution maximizes the distance $`b/𝒘`$ of the plane $`(𝒘,b)`$ to the origin. The dual form of the same quadratic programming problem can be written as :

Minimize $`{\displaystyle \underset{i,j}{}}\alpha _i\alpha _j𝒙_i𝒙_j`$
subject to $`\alpha _i0\text{ and }{\displaystyle \underset{i}{}}\alpha _i=1`$
where $`𝒙_i,𝒙_j\{(𝒄_N𝒄_D)\}`$.

### Nonlinear Potential.

However, it is possible that no such weight vector $`𝒘`$ exists, i.e., the final convex polyhedron $`𝒫=_{i=1}^n𝒫_i`$ may be an empty set. First, for a specific native protein $`i`$, there may be severe restriction from some inequality constraints, which makes $`𝒫_i`$ an empty set. Some decoys are very difficult to discriminate due to perhaps deficiency in protein representation. In these cases, it is impossible to adjust the weight vector so the native structure has a lower energy than the decoy. Second, even if a weight vector $`𝒘`$ can be found for each native protein, i.e., $`𝒘`$ is contained in a nonempty polyhedron, it is still possible that the intersection of $`n`$ polyhedra is an empty set, i.e., no weight vector can be found that can stabilize all native proteins against the decoys simultaneously. Computationally, the question whether a solution weight vector $`𝒘`$ exists can be answered unambiguously in polynomial time , and recent studies using millions of decoys strongly suggest that when the number of native protein structure reaches 300–400, no such weight vector can be found . When perfect discrimination is impossible, a technique that minimizes the percentage of unsatisfied inequalities or the error rate was developed in reference .

A fundamental reason for this failure is that the functional form of linear sum of pairwise interaction is too simplistic. It has been suggested that higher order interactions such as three-body or four-body contacts should be incorporated . Functions with polynomial terms using upto 6 degree of Chebyshev expansion has also been used to represent pairwise interactions .

Here we propose an alternative approach. At this time we still limit ourselves to pairwise contact interactions. We introduce a nonlinear potential function analogous to the dual form of the linear function in Equation (2), which takes the following form:

$$H(f(𝒔,𝒂))=H(𝒄)=\underset{D𝒟}{}\alpha _DK(𝒄,𝒄_D)\underset{N𝒩}{}\alpha _NK(𝒄,𝒄_N)$$
(3)
where $`\alpha _D0`$ and $`\alpha _N0`$ are parameters of the potential function to be determined, and $`𝒄_D`$ is a contact vector of decoy $`D`$ in the set of decoys $`𝒟=\{(𝒔_𝒟,𝒂_𝒩)\}`$, $`𝒄_N`$ a contact vector of native structure $`N`$ in the set of native training proteins $`𝒩=\{(𝒔_N,𝒂_N)\}`$. The difference of this functional form from linear function in Equation (2) is that a kernel function $`K(𝒙,𝒚)`$ replaces the linear term. A convenient kernel function $`K`$ is:

$$K(𝒙,𝒚)=e^{𝒙𝒚^2/2\sigma ^2}.$$

Intuitively, the potential surface has smooth Gaussian hills of height $`\alpha _D`$ centered on the location $`𝒄_D`$ of decoy structure $`D`$, and has smooth Gaussian cones of depth $`\alpha _N`$ centered on the location $`𝒄_N`$ of native structures $`N`$. Ideally, the value of the potential function will be $`1`$ for contact vectors $`𝒄_N`$ of native proteins, and will be $`+1`$ for contact vectors $`𝒄_D`$ of decoys.

### Optimal Nonlinear Potential.

To obtain the nonlinear potential, our goal is to find a set of parameters $`\{\alpha _D,\alpha _N\}`$ such that $`H(f(𝒔_N,𝒂_N))`$ has energy value close to $`1`$ for native proteins, and the decoys have energy values close to $`+1`$. There are many different choices of $`\{\alpha _D,\alpha _N\}`$. We use an optimality criterion originally developed in statistical learning theory . First, we note that we have implicitly mapped each structure and decoy from $`R^{210}`$ through the kernel function of $`K(𝒙,𝒚)=e^{𝒙𝒚^2/2\sigma ^2}`$ to another space with dimension as high as tens of millions. Second, we then find the hyperplane of the largest margin distance separating proteins and decoys in the space transformed by the nonlinear kernel. That is, we search for a hyperplane with equal and maximal distance to the closest native proteins and the closest decoys. Such a hyperplane can be found by obtaining the parameters $`\{\alpha _D\}`$ and $`\{\alpha _N\}`$ from solving the following Lagrange dual form of quadratic programming problem:

Maximize $`{\displaystyle \underset{i𝒩𝒟}{}}\alpha _i{\displaystyle \frac{1}{2}}{\displaystyle \underset{i,j𝒩𝒟}{}}y_iy_j\alpha _i\alpha _je^{𝒄_i𝒄_j^2/2\sigma ^2}`$
subject to $`0\alpha _iC`$
where $`C`$ is a regularizing constant that limits the influence of each misclassified conformation , and $`y_i=+1`$ if $`i`$ is a native protein, and $`y_i=1`$ if $`i`$ is a decoy. These parameters lead to optimal classification of unseen test set proteins against decoys .

### Modeling Sequence Design Potential.

For protein sequence design, we assume that native sequence is more stable on the native conformation than on a different structure of another protein with low sequence identity. We use the method of gapless sequence threading to generate sequence decoys by mounting a sequence fragment from a different protein of a larger size of an unrelated structure to the native structure . Because all native contacts are retained, we find that such sequence decoys are quite challenging.

For protein design, we seek potential functions that allows the search and identification of sequences most compatible with a specific given coarse-grain three-dimensional structure. We use a model analogous to the Anfisen experiments. We require that the native amino acid sequence $`𝒂_N`$ mounted on the native structure $`𝒔_N`$ has the lowest energy compared to a set of unrelated alternative sequences $`𝒟=\{𝒔_N,𝒂_D\}`$ mounted on the same native protein structure $`𝒔_N`$:

$$H(f(𝒔_N,𝒂_N))<H(f(𝒔_N,𝒂_D))\text{for all }𝒂_D𝒟$$

Equivalently, the native sequence will have the highest probability to fit into the specified native structure. This is the same principle described in . Much work has been done using linear design function of sum of contact pairs in the form of $`H(f(𝒔,𝒂))=H(𝒄)=𝒘𝒄`$ . The discussion of the two geometric views of the weight vector of linear potential and the optimality criterion also applies to the sequence design problem.

We now explore nonlinear potential function for sequence design using the following functional form:

$$H(f(𝒔,𝒂))=H(𝒄)=\underset{D𝒟}{}\alpha _DK(𝒄,𝒄_D)\underset{N𝒩}{}\alpha _NK(𝒄,𝒄_N)$$
(4)
where $`\alpha _D0`$ and $`\alpha _N0`$ are coefficients to be determined, and $`𝒄_D=f(𝒔_N,𝒂_D)`$ is the contact vector of a decoy sequence $`𝒂_D`$ mounted on its native protein structure $`𝒔_N`$, and $`𝒄_N=f(𝒔_N,𝒂_N)`$ is the contact vector of a native sequence $`𝒂_N`$ from the set of native training proteins $`𝒩`$ mounted on the native structure $`𝒔_N`$. The only difference from nonlinear folding potential of Equation (3) is that here $`𝒟`$ is a set of sequence decoys mounted on native protein structures, rather than a set of structure decoys. Again, we use kernel function $`K(𝒙,𝒚)=e^{𝒙𝒚^2/2\sigma ^2}`$.

The optimal parameters $`\{\alpha _N\}`$ and $`\{\alpha _D\}`$ are obtained similarly by solving the following Lagrange dual convex quadratic programming problem:

Maximize $`{\displaystyle \underset{i𝒩𝒟}{}}\alpha _i{\displaystyle \frac{1}{2}}{\displaystyle \underset{i,j𝒩𝒟}{}}y_iy_j\alpha _i\alpha _je^{𝒄_i𝒄_j^2/2\sigma ^2}`$
subject to $`0\alpha _iC`$
where $`C`$ is a regularizing constant , and $`y_i=+1`$ if $`i`$ is a native protein, and $`y_i=1`$ if $`i`$ is a decoy.

## 3 Computational Methods

### Alpha Contact Maps.

Because protein molecules are formed by thousands of atoms, their shapes are complex. In this study we use the count vector of pairwise contact interactions derived from the edge simplices of the alpha shape of a protein structure. Edge simplices in the alpha shape represent nearest neighbor interactions that are in physical contacts. They encode precisely the same contact information as a subset of the edges in the Voronoi diagram of the protein molecule. These Voronoi edges are shared by two interacting atoms from different residues, but intersect with the body of the molecule modeled as the union of atom balls. We refer to references  for further theoretical and computational details.

### Generating Structure Decoys and Sequence Decoys by Threading.

Maiorov and Crippen used the gapless threading method for generating a large number of structure decoys . In this method, the sequence of a smaller protein $`𝒂_N`$ is threaded through the structure of an unrelated larger protein and takes the conformation $`𝒔_D`$ of a fragment of the larger protein . Along the way, the sequence of the smaller protein can take the conformations of many fragments of the larger protein, each provides a structure decoy. With this approach, we generate for each native protein $`\{(𝒔_N,𝒂_N)\}`$ a set of structure decoys $`\{(𝒔_D,𝒂_N)\}`$ (Fig 2a).

We can generate sequence decoys in an analogous way, as already suggested in . In this case, we thread the sequence of a larger protein through the structure of a smaller protein, and obtain sequence decoys by mounting a fragment of the native sequence from the large protein to the full structure of the small protein. We therefore have for each native protein $`(𝒔_N,𝒂_N)`$ a set of sequence decoys $`(𝒔_N,𝒂_D)`$ (Fig 2b).

### Protein Data.

Following reference , we use protein structures contained in the Whatif database  in this study. Whatif database contains a representative set of sequence-unique protein structures generated from X-ray crystallography. Structures selected for this study all have pairwise sequence identity $`<30\%`$, R-factor $`<0.21`$, and resolution $`<2.1`$. Whatif database contains less structures than Pdbselect because the R-factor and resolution criteria are more stringent . Nevertheless, it provides a good representative set of currently all known protein structures.

We use a list of 456 proteins compiled from the 1998 release (Whatif98) of the Whatif database , which was kindly provided by Dr. Vendruscolo. There are 192 proteins with multiple chains in this dataset. Some of them have extensive interchain contacts. For these proteins, it is possible that their conformations may be different if there are no interchain contacts present. We use the criterion of Contact Ratio to remove proteins that have extensive interchain contacts. Contact Ratio is defined here as the number of interchain contacts divided by the total number of contacts a chain makes. For example, protein 1ept has four chains A,B,C, and D. The intra chain contact number of chain B is 397. Contacts between chain A and chain B is 178, between B and C is 220, between B and other heteroatoms is 11. The Contact Ratio of chain B is therefore $`(178+220+11)/(397+178+220+11)=51\%`$. Thirteen protein chains are removed because they all have Contact Ratio $`>30\%`$. We further remove three proteins because each has $`>10\%`$ of residues missing with no coordinates in the Protein Data Bank file. The remaining set of 440 proteins are then used as training set for developing both folding and design potential functions. Using the sequence and structure threading method described earlier, we generated a set of 14 080 766 sequence decoys and a set of structure decoys of the same size.

### Learning Linear Potential.

For comparison, we develop optimal linear potential following the method and computational procedure described in reference. We apply the interior point method as implemented in BPMD by Mészáros  to search for a weight vector $`𝒘`$. We use two different optimization criteria as described in reference . The first is:

Identify $`𝒘`$
subject to $`𝒘(𝒄_N𝒄_D)<ϵ\text{ and }|𝒘_i|10`$
where $`𝒘_i`$ denotes the $`i`$-th component of weight vector $`𝒘`$, and $`ϵ=1\times 10^6`$. Let $`𝒞=\{c_Nc_D\}`$, and $`|𝒞|`$ the number of decoys. The second optimization criterion is:

Minimize $`\mathrm{min}{\displaystyle \frac{1}{|𝒞|}}{\displaystyle \left(𝒘(c_Nc_D)\right)^2}\left[{\displaystyle \frac{1}{|𝒞|}}{\displaystyle \left(𝒘(𝒄_N𝒄_D)\right)}\right]^2`$
subject to $`𝒘(𝒄_N𝒄_D)<ϵ`$
### Learning Nonlinear Kernel Potential.

We use SVMlight (http://svmlight.joachims.org/)  with Gaussian kernels and a training set of 440 native proteins plus 14 080 766 decoys to obtain the optimized parameter $`\{\alpha _N,\alpha _D\}`$. The regularization constant $`C`$ takes default value, which is estimated from the training set $`𝒩𝒟`$ as implemented in SVMlight:

$$C=|𝒩𝒟|^2/\left[\underset{𝒙𝒩𝒟}{}\sqrt{K(𝒙,𝒙)2K(𝒙,\mathrm{𝟎})+K(\mathrm{𝟎},\mathrm{𝟎})}\right]^2.$$
(5)
Since we cannot load all 14 millions decoys into computer memory simultaneously, we use three heuristic strategies for training. In the first method, the total 14 080 766 decoys are divided randomly into 34 subsets. The training process starts with all native proteins and the first decoy set. After the training process converges, we select the decoy structures that have $`\alpha _i0`$ (these decoys are called support vectors ). From the second iteration on, we combined the native proteins, the new decoy set, and all decoys that appeared as support vectors in any of the previous training steps. These form the new training set for the next iteration of learning. This process is continued until all decoys have been exhausted in training.

In the second method, each of the 34 decoy subsets was combined with the native set in turn to form a training set. The decoy support vectors are then selected from each of the 34 training processes. All these decoy support vectors are combined, along with native proteins for a final training process. The results are taken as the optimized protein potential.

The third method is similar to the procedure reported in . We first randomly selected a subset of decoys that fits into the computer memory. For example, we pick every 51st decoy from the list of 14 million decoys. This leads to an initial training set of 276 095 decoys and 440 native proteins. An initial protein potential is then obtained after learning. Next the energies for all 14 million decoys and all 440 native proteins are evaluated. Three decoy sets were collected based on the evaluation results: the first set of decoys contains the violating decoys which have lower energy than the native structures; the second set contains decoys with the lowest absolute energy, and the third set contains decoy support vectors identified in previous training process. The union of these three subsets of decoys are then combined with the 440 native protein as the training set for the next iteration of learning. This process is repeated until the energy difference to native protein for all decoys are greater than 0.0. Using this strategy, the number of iterations typically is between 2 and 10. During the training process, we set the cost factor $`j`$ in SVMlight to 120, which is the factor by which training errors on native proteins outweighs errors on decoys.

The value of $`\sigma ^2`$ for the Gaussian kernel $`K(𝒙,𝒚)=e^{𝒙𝒚^2/2\sigma ^2}`$ is chosen by experimentation. We find that if the value of $`\sigma ^2`$ is too large, no $`\{\alpha _N,\alpha _D\}`$ can be found that can perfectly classifies the 440 training proteins and their decoys, i.e., the problem is unlearnable. If the value of $`\sigma ^2`$ is too small, the performance in blind-test will deteriorate because of overfitting. For the 440 native proteins and the 14 080 766 sequence decoys, the value of $`\sigma ^2`$ is chosen between 250.0 (for training method 1) and 138.9 (for training method 3). The final folding potential is obtained with $`\sigma ^2=227.3`$ and the final design potential is obtained with $`\sigma ^2=138.9`$, both derived using the third training method.

## 4 Results

### Linear Folding Potentials from Structure Decoys.

Using perceptron learning and a large number of decoys generated by gapless threading, Vendruscolo et al showed that if the number of native proteins exceeds 270, it is impossible to find parameters $`𝒘`$ for potential function $`H(𝒔,𝒂)=𝒘𝒄`$, such that all native structures have lower energies than decoys . That is, no $`𝒘`$ can be found such that $`𝒘𝒄_N<𝒘𝒄_D`$ holds for all decoys.

Tobi et al showed that pairwise contact potential $`H(𝒔,𝒂)=𝒘𝒄`$ can be found to distinguish a different set of 572 proteins from a set of 28 213 009 structure decoys generated by gapless threading . In this study, two residues are defined to be in contact if the geometric centers of their side chains are at a distance between 2.0 Å and 6.4 Å. To search for the optimal weight vector $`𝒘`$, the authors used linear programming solver based on interior point method as implemented in BPMD by Mészáros .

We succeeded in reproducing the results of Tobi et al  and found a $`𝒘`$ vector that enables perfect discrimination of the same set of 572 proteins used in  against 28 261 307 structure decoys we generated by gapless threading. We used the same contact definition based on Euclidean distance between geometric centers of side chains, as well as the same optimization criterion described earlier in the Methods section. The values of the 210 pairwise contact potentials (data not shown) are similar although not identical to the values listed in .

### Linear Design Potentials from Sequence Decoys.

Structure decoys generated by gapless threading are obtained by taking a fragment of the structure of a large protein such that it contains exactly the same number of amino acid residues as the native protein. However, contacts in the resulting fragment structure often contain only a fraction of the total contacts in native protein.

Design decoys are generated by taking a fragment sequence of a larger protein and thread it onto the native conformation of a smaller protein. As a result, although the contact count vector $`𝒄`$ may be very different, the number of contacts in the design decoy is exactly the same as that in native protein. For decoys generated by gapless threading, design decoys therefore are far more challenge to discriminate then structure decoys (Fig 3).

Our experimentation confirms this observation. After generating 14 080 766 sequence design decoys for the 440 proteins in the training set, we apply the same interior point method to search for an optimal $`𝒘`$ that can discriminate native sequences from decoy sequences. That is, we search for parameters $`𝒘`$ for $`H(𝒔,𝒂)=𝒘𝒄`$, such that $`𝒘𝒄_N<𝒘𝒄_D`$ for all sequences. However, we fail to find a feasible solution for the weight vector $`𝒘`$. This indicates that no $`𝒘`$ exists capable of discriminating perfectly 440 native sequences from the 14 million decoy sequences. We repeated the same experiment using the set of 572 native proteins from reference  and 28 261 307 sequence decoys. The result is also negative.

### Learning Nonlinear Kernel potential.

To overcome the problems associated with linear potentials, we use the same set of 440 native proteins and 14 million decoys to obtain nonlinear kernel folding and design potentials. In both cases, we succeeded in finding a function in the form of Equation (3) that can discriminate all 440 native proteins from 14 million decoys.

Unlike statistical potentials where each native protein structure in the database contribute to the empirical potential, only a subset of native structures contribute and have $`\alpha _N0`$. In addition, a small fraction of decoys also contribute to the potential function. Table 1 list the details of the potential, including the numbers of native proteins and decoys that participate in Equation (3). These number represent about $`60\%`$ of native proteins and $`<0.1\%`$ of decoys from the original training data.

### Discrimination Tests for Folding Potential.

Blind test in discriminating native proteins from decoys using an independent test set is essential to assess the effectiveness of folding potentials. To construct a test set, we first take the entries in Whatif99 database that are not present in Whatif98. After eliminating proteins with chain length less than 46 residues, we obtain a set of 201 proteins. These proteins all have $`<`$ 30% sequence identities with any other sequence in either the training set or the test set proteins. Since 139 of the 201 test proteins have multiple chains, we use the same criteria applied in training set selection to exclude 7 proteins with $`>30\%`$ Contact Ratio or with $`>10\%`$ residues missing coordinates in the PDB files. This leaves a smaller set of test proteins of 194 proteins. Using gapless threading, we generate a sets of 3 096 019 structure decoys from the set of 201 proteins. This is a superset of the decoy set generated using 194 proteins.

For comparison, we also test the discrimination results of the optimal linear potential taken as reported in reference , as well as the statistical potential developed by Miyazawa and Jernigan. Here we use the contact definition reported in , that is, two residues are declared to be in contact if the geometric centers of their side chains are within a distance of 2.0 – 6.4 Å.

To test nonlinear folding potential functions for discriminating native proteins from structure decoys in both the 194 and the 201 test sets, we take the structure $`𝒔`$ from the conformation-sequence pair $`(𝒔,𝒂_N)`$ with the lowest energy as the predicted structure of the native sequence. If it is not the native structure $`𝒔_N`$, the discrimination failed and the folding potential does not work for this protein. The results of discriminating native structures using nonlinear folding potential are summarized in Table 2. There are 3 and 8 misclassified native structures for the 194 set and 201 set, respectively. These correspond to a failure rate of 1.5% and 4.0%, respectively. The optimal nonlinear kernel folding potential performs better than the optimal linear potential based on calculation using potential values taken as reported in reference  (failure rates 3.6% and 6.5% for the 194 set and 201 set, respectively). Consistent with previous reports , statistical potential has about $`43.8\%`$ (81 out of 194) and $`43.2\%`$ (87 out of 201) failure rates for the 194 set and the 201 set, respectively.

### Discrimination Tests for Design Potential.

We use the same 194 set and 201 set of natives proteins and generate a set of 3 096 019 sequence decoys for testing the design potential. We take the sequence $`𝒂`$ from the conformation-sequence pair $`(𝒔_N,𝒂)`$ for a protein with the lowest energy as the predicted sequence. If it is not the native sequence $`𝒂_N`$, the discrimination failed and the design potential does not work for this protein.

Sequence decoys obtained by gapless threading are quite challenging, since all native contacts of the protein structures are maintained. No linear design potential function can be found using linear programming method that would succeed in the challenging task of identifying all 440 native sequences in the training set.

We succeeded in obtaining nonlinear design potential capable of discriminating all of the 440 native sequences (Table 3). It also works well in the test set. It succeeded in correctly identifying 92.8% (180 out of 194) of native sequences in the independent test set. This compares favorably with results obtained using optimal linear folding potential taken as reported in , which succeeded in identifying 80.9% (157 out of 194) of the test set, and the Miyazawa-Jernigan statistical potential (success rate 58.2%, 113 out of 194).

### Discrimination Test Using Challenging Decoys.

For the protein potentials derived with simple decoys generated by gapless threading, a more challenging test is to discriminate native proteins from an ensemble of explicitly generated three dimensional decoy structures with a significant number of near-native conformations . Here we evaluate the performance of nonlinear folding and design potential using three decoy sets from the database “Decoys ’R’ Us: the 4state\_reduced set, the Lattice\_ssfit set, and the lmsd set. We compare our results in performance with results reported in literature using optimal linear potential  and statistical potential  (Table 4). For the 4state\_reduced set of decoys, nonlinear design potential has the best performance in terms of identifying the native structure. The only misclassified protein 1sn3 has three disulfide bonds, which are not modeled explicitly in the protein description of a vector in $`R^{210}`$. The correlation of root mean square distance of conformations to the native structure and energy value in the 4state set are shown in Fig 4. For proteins 1r69, 2cro, 3icb, and 1ctf, the correlation coefficient between the rmsd values and the energy values are good. We note that the nonlinear potential is obtained by learning from training decoys that are obtained by gapless threading. It is likely that nonlinear kernel potential can be further improved if more realistic structure decoys are included in training.

### Nonlinear Potential Function for Folding and Design.

Because nonlinear potential is expressed as a kernel expansion of native protein and decoys, structure decoys and sequence decoys in general lead to different protein functions. For example, the contact count vectors $`𝒄`$ can be very different for a sequence decoy of a protein and a structure decoy of the same protein. The potential surface defined by the folding potential and the design potential therefore may be different. There are 268 out of 440 native proteins participating in folding potential function with $`\alpha `$ value ranging from 0.01 - 28.04, of which 185 ($`65\%`$) are between 0.01 and 2.00. There are 347 out of 440 native proteins participating in design potential function, with $`\alpha `$ value ranging from 0.02 to 110.64, of which 269 ($`77\%`$) are between 0.02 and 2.00. Therefore, the majority of the native proteins have similar $`\alpha `$ values for both folding and design potentials. Fig 5 shows the difference $`\mathrm{\Delta }\alpha _i`$ of the coefficient $`\alpha _i`$ for protein $`i`$ appearing in both folding potential and design potential. For the majority of the native proteins, $`\mathrm{\Delta }\alpha _i`$ values are small. That is, most native proteins contribute similarly in design potential and in folding potential. This is expected, because the main differences between the two potentials are due to differences in decoy sets. The native proteins that contribute most differently to folding potential and design potential are small proteins. For example, there are 34 native proteins whose $`|\mathrm{\Delta }\alpha |>5.0`$, and they all come from the top 100 smallest protein chains. Table 5 listed 20 proteins with highest $`\alpha `$ values for both kernel folding and design potentials. They are all small proteins and there are 15 out of the 20 proteins that appear both in folding and design potentials. It is possible that the energy values by kernel folding potential and by kernel design potential may be similar for many structure-sequence pairs $`(𝒔,𝒂)`$. Figure 6a shows that the test 194 proteins have similar energy values by the kernel folding and kernel design potentials.

We also compare the energy value of the potential functions for each of the 210 unit vector $`𝒄_1=\{1,0,\mathrm{},0\}^T,\mathrm{}`$, and $`𝒄_{210}=\{0,\mathrm{},1\}^T`$. We normalize these values so $`\mathrm{max}H(𝒄_i)=1`$ for both potentials (Fig 7a). There is strong correlation ($`R=0.91`$) for folding and design potentials.

However, other methods reveal that kernel folding and design potentials are different. One method is to compare the energy values of a subset of decoy structures that are challenging. That is, we compare energies of decoys with $`\alpha _i0`$. Fig 7b shows that for decoys appearing in the design potentials, there is little correlation in energy values calculated by design potential and by folding potential. Similarly, there is no correlation between energy values calculated by folding potential and by design potentials for the set of structure decoys entering the design potential function (Fig 7c). It seems that although the values of $`\alpha _N`$s are similar for the majority of the native proteins, design potential and folding potential can give very different energy values for some conformations. This suggests that the overall fitness for design and folding potential may be somewhat different.

### Comparison with Other Potentials.

Potential functions derived from kernel models have the form of $`H(𝒄)=_{d𝒟}\alpha _dK(𝒄,𝒄_d)_{N𝒩}\alpha _nK(𝒄,𝒄_n)`$, which cannot be directly compared with other potential functions of the form $`H(𝒄)=𝒘𝒄`$. Here we follow the same approach as above and compare the energy values of the unit vectors, the native proteins, and the decoys using different potential functions. The energy values by nonlinear folding and design potentials have very little correlation with energy values evaluated using either optimal linear folding potential  or statistical potential  (Fig 6c-f). It seems nonlinear potentials developed in this work contains information absent in other potential functions.

## 5 Discussion

A basic requirement for computational studies of protein folding and protein design is an effective potential function, which allows the search and the identification of native structures and native sequences. Current empirical potential functions are based on weighted linear sum of pairwise interactions. A major flaw of such potentials is that they cannot recognize the native structures of a large number of proteins from alternative decoy conformations .

There are several routes towards improving empirical potential functions. One approach is to introduce higher order interactions, where three-body or four-body interactions are explicitly incorporated in the potential function . The effectiveness of this approach is likely to be assessed in the future with large scale validation studies similar to those carried out for pairwise potentials . A different approach to improve empirical potential function is to introduce nonlinear terms. Recently, Fain et al uses sums of Chebyshev polynomials upto order 6 for hydrophobic burial and each type of pairwise interactions .

In this work, we propose a different framework for developing empirical nonlinear protein potential functions. Instead of using nonlinear function for each term of pairwise interactions, we use a set of simple Gaussian kernel functions located at both native proteins and decoys as the basis set. We regard the decoys as equivalent to the reference state or null model used in statistical potential. The expansion coefficients $`\{\alpha _N\},N𝒩`$ and $`\{\alpha _D\},D𝒟`$ of these Gaussian kernels determines the protein function. As long as the native proteins and decoys are represented as unique vectors $`𝒄R^d`$, the Gram matrix of the kernel function is full-rank. Therefore, the kernel function effectively maps the protein space into a high dimensional space in which effective discrimination with a hyperplane is easier to obtain. The optimization criterion here is not $`Z`$-score, rather we search for the hyperplane in the transformed high dimensional space with maximal separation distance between the native protein vectors and the decoy vectors. This choice of optimality criterion is firmed rooted in a large body of studies in statistical learning theory, where expected number of errors in classification of unseen future test data is minimized probabilistically by balancing the minimization of the training error (or empirical risk) and the control of the capacity of specific types of functions of potential function .

This approach is general and flexible, and can accommodate other protein representation, as long as the final descriptor of protein and decoy is a $`d`$-dimensional vector $`𝒄`$. In addition, different forms of nonlinear functions can be designed using different kernel functions, such as polynomial kernel and sigmoidal kernels. It is also possible to adopt different optimality criterion, for example, by minimizing the margin distance expressed in 1-norm instead of the standard 2-norm Euclidean distance.

A useful observation obtained from this study is that sequence decoys obtained from gapless threading are quite challenging. In fact, we found that no linear potential function exists that can discriminate a training set of 440 native sequence from sequence decoys generated by gapless threading. The success of nonlinear potential in perfect discrimination of this training set native sequences and its good performance in identifying the native sequences in an unrelated test set of 194 proteins indicate that nonlinear kernel potential is a general strategy for developing effective potential function for protein sequence design.

It is informative to examine the three misclassified proteins by the kernel folding potential (1bx7, 1hta, and 3erd). Hirustasin 1bx7 contain five disulfide bonds, which are not modeled explicitly by the protein description. 1hta (histone Hmfa) exists as a tetramer in complex with DNA under the physiological condition. Its nature structure may not be the same as that of a lone chain. The two terminals of this protein are rather flexible, and their conformations are not easy to determine. 3erd.a (estrogen receptor $`\alpha `$ ligand-binding domain) has extensive contacts with ligand. These unmodeled interactions may alter protein conformation. Among the 14 native sequences misclassified by the kernel design potential (1a73.a, 1bd8, 1bea, 1bm8, 1bn8.a, 1bvy.f, 1cku.a, 1dpt.a, 1hta, 1mro.c, 1ops, 1qav.a, 1ubp.b, 3ezm.a), several have extensive interchain interactions, although the contact ratio is below the rather arbitrary threshold of 30%: 1a73.a has Contact Ratio of $`23\%`$, 1mor.c has $`24\%`$, 1upb.b has $`19\%`$, and 1qav.a has $`13\%`$. It is possible that the substantial contacts with other chains would alter the confirmation of the protein. Amylase inhibitor 1bea contains 5 disulfide bond not explicitly modeled. 1cku.a (electron transfer protein) contains an iron/sulfur cluster, which covalently bind to four Cys residues and prevent them from forming 2 disulfide bonds. These covalnt bonds are not reflected in the protein description. 1bvf (oxidoreductase) is complexed with a heme and an FMN group. The conformations of 1cku.a and 1bvf.f may be different upon removing of these functionally important hetero groups. Altogether, there are some rationalization for 10 of the 16 misclassified proteins.

Our goal in this study is to explore an alternative formulation of potential function and assess the effectiveness of this new approach with experimental data. The nonlinear potential functions obtained in this study should be further improved. For example, unlike the study of optimal linear potential , where explicitly generated three-dimensional decoys structures are used in training, we used only structure decoys generated by threading. The test results using the 4state\_reduced set and the lattice\_ssfit are comparable or better with other residue-based potential (see Fig 4 and Table 4). It is likely that further incorporation of explicit three-dimensional decoy structures in the training set would improve the protein potential.

In summary, we show in this study an alternative formulation of protein function using a mixture of Gaussian kernels. We demonstrate that this formulation can lead to effective folding potential and design potential that perform well in independent tests. For protein sequence design where challenging decoys are available from gapless threading, nonlinear kernel potential can have perfect classification in the training set of 440 proteins, while linear potential and statistical potential failed. It also performs better in an independent test set of 194 proteins and reduce the misclassification to 40% of that of optimal linear potential. Our results suggest that more sophisticated functional form other than the simple weighted sum of contact pairs can be useful for studying protein folding and protein design. This approach can be generalized for any other protein representation, e.g., with descriptors for explicit hydrogen bond and higher order interactions.

## 6 Acknowledgment

We thank Dr. Hao Li for very helpful comments on protein design, Dr. Yang Dai for discussion on linear programming, and Dr. Mészáros for help in using the BPMD package. Dr. Michele Vendruscolo kindly provided the list of the set of 456 proteins. This work is supported by funding from National Science Foundation CAREER DBI0133856, DBI0078270, and American Chemical Society/Petroleum Research Fund 35616-G7.

## 7 Appendix

###### Lemma 1

For a potential function in the form of weighted linear sum of interactions, a decoy always has energy values higher than the native structure by at least an amount of $`b>0`$, i.e.,

$$𝒘(𝒄_D𝒄_N)>b\text{ for all }\{(c_Dc_N)|D𝒟\text{ and }N𝒩\}$$
(6)
if and only if the origin $`\mathrm{𝟎}`$ is not contained within the convex hull of the set of points $`\{(𝐜_D𝐜_N)|D𝒟\text{ and }N𝒩\}`$.

Proof:  Suppose that the origin $`\mathrm{𝟎}`$ is contained within the convex hull $`𝒜=\text{conv}(\{𝒄_D𝒄_N\})`$ of $`\{𝒄_D𝒄_N\}`$ and Equation (6) holds. By the definition of convexity, any point inside or on the convex hull $`𝒜`$ can be expressed as convex combination of points on the convex hull. Specifically, we have:

$$\mathrm{𝟎}=\underset{(c_Dc_N)𝒜}{}\lambda _{c_Dc_N}(𝒄_D𝒄_N),\text{ and }\lambda _{c_Dc_N}=1,\lambda _{c_Dc_N}>0.$$

That is, we have the following contradiction:

$$0=𝒘\mathrm{𝟎}=𝒘\underset{c_Dc_N}{}\lambda _{c_Dc_N}(𝒄_D𝒄_N)=\underset{c_Dc_N}{}\lambda _{(c_D,c_N)}𝒘(𝒄_D𝒄_w)>\underset{c_Dc_N}{}\lambda _{c_Dc_N}b=b.$$

Because the convex hull can be defined as the intersection of half hyperplanes derived from the inequalities, if a half hyperplane has a distance $`b>0`$ to the origin, all points contained within the convex hull will be on the other side of the hyperplane . Therefore, $`𝒘(𝒄_D𝒄_N)>b`$ will hold for all $`\{(𝒄_D𝒄_N)\}`$.