Network measures for protein folding state discrimination

Proteins fold using a two-state or multi-state kinetic mechanisms, but up to now there is not a first-principle model to explain this different behavior. We exploit the network properties of protein structures by introducing novel observables to address the problem of classifying the different types of folding kinetics. These observables display a plain physical meaning, in terms of vibrational modes, possible configurations compatible with the native protein structure, and folding cooperativity. The relevance of these observables is supported by a classification performance up to 90%, even with simple classifiers such as discriminant analysis.

as independent as possible from the protein size, they were accordingly rescaled by a function of residue chain length. In this paper we show that these observables perform very well even with a simple discriminant classifier, that allows to give a intuitive biophysical interpretation to our results.

Results
As previously mentioned, we introduced three main groups of observables (see Methods), i.e.
• Network Entropy-based ratio S R , a measure related to the possible configurations consistent with the native protein structure PCN combined with the contact potential matrix P (see Methods). • Laplacian-based observables λ N , λ N−1 , λ N−2 , related to the highest-frequency vibrational modes of a modified PCN, in which the backbone is removed. • Inter-residue link density R 0 , measuring the fraction of sequence separation (diagonals of the full PCN) that do not contain residue contact pairs.
All the defined observables were considered as features for classification of TS and MS proteins by Fisher discriminant analysis (see Methods).
Classification with combinations of the novel observables is very high and balanced between the two classes, with performances up to 88% and Matthews Correlation Coefficient MCC = 0.76 (see Methods). Previous results, with a more complex classifier (nonlinear Support Vector Machine 13 ) that does not allow a simple interpretation of the results as in our case, were around 80%. The details of all performances of the analyzed signatures, both with single observables and their combinations in couples, are shown in Table 1. We also considered higher-dimensional signatures (with combinations of 3 and 4 observables) but the performance was In Panel C the related Long-range Interaction Network (LIN) is shown 11,14 in which no backbone diagonal is present. not significantly increased. The combinations of the largest Laplacian eigenvalues λ N , λ N−1 , λ N−2 with the link density R 0 produce the best performances of classification, with a top-score value given by the couple (λ N−1 , R 0 ), with 88.33% ± 1.10% correctly classified proteins, and a highly homogeneous performance on both classes (87.20% ± 2.21% MS and 89.07% ± 1.01% TS correctly classified proteins, MCC = 0.76 ± 0.02). Considering the full PCN matrix, the classification performance of the Laplacian eigenvalue was reduced (see Methods). The entropy ratio S R is the best single classifier (80.36% ± 1.81% correctly classified proteins, MCC = 0.59 ± 0.04), and it also has a very high performance in combination with R 0 (84.50% ± 1.30% correctly classified proteins, MCC = 0.67 ± 0.03). Classifying without cross-validation, i.e. using the entire set of (λ N−1 , R 0 ) features, we obtain a performance of 90.48%, with 92.00% for MS proteins and 89.47% for TS proteins. In Fig. 2 we show the scatter-plot for two top-scoring couples: (λ N−1 , R 0 ) and (S R , R 0 ). It appears that MS and TS proteins are almost linearly separated in this parameter space, and this may allow a simple interpretation in terms of the observables, as stated in the Discussion section. Since in previous studies 9 it has been shown that the chain length N C is a good classifier of folding classes, we rescaled our observables in order to keep them as much independent as possible from protein length. Moreover, as a comparison for classification performance, we used N C as a variable for discrimination. In our dataset, the N C parameter correctly classifies 78.23% ± 1.52 proteins, with a large unbalance between correctly classified proteins from the two classes: 57.61% ± 3.04 for MS proteins and 91.80% ± 1.29 for TS proteins (MCC = 0.54 ± 0.03). We also evaluated the classification power of some measures typically used in literature, i.e. the average hydrophobicity value 〈 h〉 22 , the contact order CO 12,13 and the long range contact order LRCO 14 . In the considered dataset the couple (N C , 〈 h〉 ) correctly classifies 73.71% ± 2.15 proteins, with MCC = 0.44 ± 0.05. Both the structural topology measures perform poorly: CO correctly classifies 68.42% ± 0.65 proteins with MCC = 0.36 ± 0.01 while LRCO guesses right 54.38% ± 2.41 proteins with MCC = 0.14 ± 0.05. The best performing couple of observables (N C and LRCO) reaches a performance of 80% ± 1, with a high heterogeneity of performance for TS and MS proteins (MCC = 0.57 ± 0.02). A complete summary of the performances of classical measures (both singularly and in couples) can be found in Table 2. We also performed Discriminant Analysis with all possible combinations of new and classical observables, but the results did not outperform the performances obtained with the new observables only (CO and S R : 81.4% ± 1.8, MCC = 0.61 ± 0.04).
We also performed the same analysis on a different and more recent dataset 7 (http://kineticdb.protres.ru/ db/index.pl), containing 85 proteins for which the PDB structure is available. In this dataset we found the same top-performing variables: the combination of R 0 with S R had a 80.9% ± 1.3% ratio of correctly classified protein with MCC = 0.58 ± 0.03, while R 0 for λ N−2 we found 79.7% ± 1.1% and MCC = 0.59 ± 0.02.

Discussion
Trying to deduce properties of the proteins from their structure is still an open challenge: in this paper we propose novel observables, based on the network properties of PCN, that allow the discrimination between proteins with a different folding dynamics (TS proteins that present only two configurations (folded/unfolded) and MS proteins with a richer landscape of stable and metastable states) by looking only at information on their native state structure. Since previous work used the number of protein residues as a discriminating variable, but this may have some bias as we will discuss later, we only considered size-independent observables, by rescaling their value with an appropriate function of protein size. Thanks to this processing our analysis resulted more performing and robust, in particular in the "grey region" of short MS and long TS proteins. The performance of these observables is achieved by simple Fisher Discriminant Analysis, that allows a plain physical interpretation, in terms of vibrational modes, possible configurations compatible with the native protein structure, and folding cooperativity. The observable R 0 simply counts the density of inter-residue distances in which there are no contacts, but nonetheless it results to be very powerful for this classification purposes though being independent on protein size. This means that the information contained in the PCN bands (the diagonals of the related adjacency matrix containing links between d-neighboring nodes) is very relevant, and possibly more complex measures could be developed based on this information. Another class of observables that we have introduced is based on the eigenvalues of the Laplacian operator applied to PCNs. In analogy with the physical Laplacian operator (acting on Euclidean space) the eigenvalues and eigenvectors can be put in relation with the main vibrational modes of the network and their respective frequencies. We remark that for Laplacian observables it was important to emphasize the role of long-range residue contacts, that effectively characterize the protein 3D structure, by removing the protein backbone with a procedure that preserves the network connectivity as a unique component: based on PCN properties, this filtering of non-relevant links is protein-specific, differently from the more general definition of long-range interaction commonly used, with a unique threshold for all proteins. Our Discriminant Analysis showed that in general TS and MS proteins are better classified by larger Laplacian eigenvalues, corresponding to high-frequency vibrating modes, at difference with small eigenvalues such as the Fiedler number. From the best performing couple of observables, λ N−1 and R 0 see Fig. 2, we deduce that TS proteins have larger values of fast-vibrating frequencies, and a larger number of inter-residue contacts (as can be seen in Fig. 2A). An interesting remark is that the vibrating modes associated with large eigenvalues tend to be more localized in specific residue chain regions (such as focusing modes in optics and whispering modes in acoustics 23 ). It seems thus that the vibrating dynamics associated with specific regions of the residue may have a relevant role in these folding processes. The other observable we introduced, based on the concept of network ensemble, depends on an estimation of the size of network ensembles (from a canonical Statistical Mechanics point of view) that share common constraints (in our case the degree and the strength sequence of the PCN). As expected by the physical meaning of entropy, that counts the number of "microstates" corresponding to a "macrostate" characterized by some fixed constraints, TS proteins show a smaller value of S R , that can be interpreted as a smaller number of topological configurations available to the related networks. A high number of available PCN states, given a fixed degree and strength sequence, is thus very likely associated with MS proteins, with more intermediate states during the folding process.
We remark that some known parameters used for TS-MS classification, such as the residue number N C , are "extensive" variables that depend on protein size, thus their performances might be biased from the fact that many MS proteins are "long" and many TS proteins are "short" (e.g. considering protein size of the used database), and might not reflect real physical properties of the studied proteins. Since our analysis has a better performance of at least 10% with respect to these observables, this means that protein size is very relevant but not crucial to characterize MS and TS classes.
In conclusion, the high classification performance achieved, together with a direct physical interpretation, indicate that the newly introduced network-based observables can be relevant for a better comprehension of protein folding processes.

Methods
Experimental data. The proteins to be classified as TS or MS are obtained from the manually-annotated dataset curated by Ivankov and Finkelstein 24 . The dataset consists of 63 proteins, 25 of which are classified as Multi-State and 38 as Two-State. The protein structures are taken from the Protein Data Bank (www.rcsb.org). We model the protein structure with its alpha carbon (C α ) trace. We collapse the entire protein structures into related contact matrices between the C α s of the residues. Contact matrices represent a common way of modeling proteins, that guarantees a good representation of the complex relationship between structure and function of proteins, while cutting out the redundant information embedded in the whole 3D structure. Contact matrices are essentially networks in which the role of nodes is played by residues and edges or "contacts" depend upon a notion of "distance" between each couple of residues. The position of an entire amino acid is usually collapsed into the corresponding C α and the ordering of nodes is physically justified by the primary structure of the protein,  i.e. the protein backbone. The backbone is composed by residues that are in sequence and whose distance ranges 3-4 Å, the so-called "peptide bond". Once obtained the C α spatial distribution, the contact matrix D is considered, where each element d ij is the 3D Euclidean distance between the i th and j th residues. The protein contact network PCN is then obtained by choosing an upper threshold of 8 Å 25-28 : In order to build our Entropy-based observables, we retrieved also data regarding amino acidic interactions, such as hydropathy indexes 29,30 and contact potentials, namely, 20 × 20 matrices describing the interactions between the 20 side-chains [31][32][33] . Each element of the contact potential represents the interaction strength between a pair of amino acids at contact. In this paper we provide results only for the contact potential matrix M, as described in ref. 31, since the results obtained with other potentials 32,33 are very similar. For each protein, starting from the known residue sequence, we define the contact potential matrix P in which where r i is the amino acid residue corresponding to the i th C α , and the matrix P has the same size of the related PCN. Since P ij can assume negative values, in order to have only positive weights (necessary for the calculation of network entropy) we shifted their values in each PCN to have the smallest weight equal to one. For the Entropy calculations only, we use the Hadamard (element-wise) product of PCN and P.
Entropy-based measure S R . The first observable introduced is associated with the-so called Entropy of a network ensemble 34,35 . Network entropy is related to the logarithm of the number of typical networks (in our case the possible PCNs) that satisfy some given constraints based on node and link features of a real network instance (the studied protein). We hypothesize that the network structure of the protein native state retains information related to the protein folding process (such as the possible intermediate states that could be represented as non-native PCNs). It has been recently applied in a biological context, as a measure of the "parameter space" available to the cell (in terms of gene expression profile or clonal diversity) and it allowed to successfully characterize different cell states related to different cancer stages or to physiological ageing 36 . In our approach, each protein is considered as an undirected weighted network, in which we integrate the information on the topological structure given by protein contacts with the information on residue interactions given by {P ij } as weights for the existing PCN links (related to the contact potential matrix M described in ref. 31 and explained in details in the "Experimental Data" Section, Eq. 2). For each protein, we calculated the Network Entropy S BS for two different ensembles, with a different number of constraints: in the first ensemble, we only fix the strength sequence {s i } of the protein network (S s ), while in the second ensemble (S ks ) we fix both strength sequence {s i } and degree sequence {k i }. The degree sequence {k i } and the strength sequence {s i } are respectively defined as the number of contacts and the weighted sum of contacts of the nodes in the network.
Network Entropy can be generally defined as where weights w are positive and π ij (w) is the probability to observe weight w between residue i and j. The constraints previously defined for the calculation of maximum network entropy are written as Laplacian-based observables λ i . The Laplacian operator L on networks is a positive semi-definite operator that plays a major role in the study of diffusion processes on networks, in node clustering and network visualization 37,38 , and it has already been applied to characterize protein features 39 . Given an adjacency matrix A without self loops, we define the Laplacian operator as Scientific RepoRts | 6:30367 | DOI: 10.1038/srep30367 Remarkably, in case of a N-lattice network, the eigenvalue problem for the Laplacian operator can be put in analogy with the discretization of an N-dimensional elastic membrane 40 . With this analogy in mind, the eigenvalues of the Laplacian matrix can be associated with the oscillating frequencies (harmonics) of the vibrating modes on the membrane, with the largest eigenvalues corresponding to the highest frequencies.
Since we suppose that the relevant information on folding kinetics can be contained in the long-range contacts of the native folded state 41 , we decided to partially remove the backbone contacts from the original PCN. In more detail, for each protein we evaluated b, the number of d-diagonals (the set of links between nodes at a distance of d residues along the backbone, see Eq. 8) needed to break the protein network in more than one component. Then, b − 1 diagonals were removed from the PCN. We verified that with the full PCN our classification performance was significantly reduced by 2-3% on average, justifying our choice to emphasize the role of long-range connections with respect to the protein backbone. We remark that this procedure is specific for each protein, i.e. the parameter b depends on the PCN of each protein, and moreover, once removed b − 1 diagonals, the PCN is still connected, thus generating a unique eigenvalue spectrum for L. Also other authors considered a reduced PCN 11,14 , but they used a unique threshold to define long-range interactions, taking only inter-residue distances d > 12 independently from protein size and structure (see Fig. 1). Once the laplacian spectrum of this modified version of PCN was computed for each protein, we considered as observables {λ i }, the largest eigenvalues of L rescaled by the number of residues N C . This rescaling was chosen because we observed a dependence of the largest eigenvalues on N C , so our observables could in principle be linearly dependent on the number of residues. According to the vibrational interpretation of the Laplacian, these eigenvalues represent the highest vibrational frequencies associated with the long-range structure of the protein.
Inter-residue link density R 0 . Each PCN is an adjacency matrix, in which the d-diagonals ij contain all the links between residues with a sequence separation equal to d (ranging from 1 to N C − 1) with respect to the protein backbone. The observable R 0 is defined as the ratio between the number of d-diagonals without links and the number of residues N C of the protein, i.e.
where δ is the Kronecker delta. For a given protein, a low value of R 0 implies that the residues interact at many different levels of sequence separation (different values of d). On the contrary, a high value of R 0 indicates that, in such protein structure, the residue interactions are more localized and show less cooperativity.
Statistical analysis for classification. Fisher Discriminant Analysis, a robust classifier that allows plain interpretations of the classifying parameters due to the simple boundaries separating the classes, was applied to single observables and to their combinations (i.e. couples, triplets and quadruplets of observables). A 10-fold cross-validation with 10000 re-samplings was used to assess the performance of our classifiers, that will be described by the average value over the re-samplings and by the standard deviation as confidence interval. Given the presence of homologous proteins, in each partition of the 10-fold cross-validation all the homologous proteins were kept together, to reduce the risk of overestimating the classifier performance. In order to characterize the homogeneity of the classification performance over both TS and MS classes, we consider the Matthews Correlation Coefficient MCC, defined as where T P is the number of true positives, T N the number of true negatives, F P the number of false positives and F N the number of false negatives. A coefficient of + 1 defines a perfect prediction, 0 is nothing more than a random prediction, while − 1 reflects total disagreement between prediction and observation.
Classical observables. We investigated the classification power of some classical observables, with the main purpose to have a comparison with the introduced measures. The statistical analysis and the estimation of the classification performances were the same as for the novel observables (see previous Subsection).
We analyzed the performance of the chain length N C , the average hydrophobicity value 〈 h〉 , the contact order CO and the long range contact order LRCO.
The parameter h is chosen since hydrophobic force has always been indicated as one of the major drivers for protein folding 22 . Each amino acid is associated with a hydropathy index h i , a number representing the hydrophobic or hydrophilic properties of its side-chain, thus each protein can be associated with an average hydrophobicity value 〈 h〉 : where h i KD refers to the hydropathy index of residue i when the Kyte-Doolittle (KD) scale 29 is considered. The average hydrophobicity 〈 h〉 has been often coupled with N C to classify the protein folding kinetics.
We also considered the classification power of structural topology measures such as contact order 12