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.

Protein folding is one of the most studied biophysical problems1, and despite the fact that protein folding is a straightforward biophysical process2, up to now there is not a general agreement on how and why proteins fold3. Experimentally, protein folding kinetics is divided into two fundamental categories: Two-State (TS) folding and Multi-State (MS) folding. While Two-State kinetics can be considered as an “all-or-none” transition, Multi-State folding displays at least one or more intermediates. Measuring experimentally the type of protein kinetics is not an easy task3, and computational studies can help unraveling relevant mechanisms4. The classification of proteins in these two major groups and the related prediction of folding rates have been widely debated in recent years45678. Previous studies have focused on several different types of predictors91011, exploiting the main features of protein primary structures and protein contact map representations (for a review see ref. 3). The geometry of the native protein structure plays a relevant role to infer the value of the folding rate. For this task, different predictors have been proposed based on: structural topology measures such as contact order1213 and long range contact order14, clustering coefficient, characteristic path length and assortativity coefficient11, cliquishness15, chain length and amino acid composition910. These observables or combinations of them were usually evaluated by means of binary logistic regression (BLR) and support vector machine (SVM). In particular, SVM classifiers map the data into a higher dimensional feature space, that is usually not easily interpretable in terms of the original variables. Most predictors do not usually perform in the same way both for Two-State and Multi-State proteins, causing unbalanced value of sensitivity and sensibility according to the target of the analysis.

In this paper we focus on the problem of the discrimination of protein folding state, and not on the real-valued prediction of the folding rates. Our aim is to find physics-based, and easy interpretable observables, that can be related to the folding state classification. We propose novel observables based on the network properties of the native structure only, and we show that, together with a clear physical interpretation, they also predict with high performance if a protein behaves as Two-State or Multi-State by using simple discriminant methods. As done before by other authors161718, we represent the protein 3D structure as a contact map between amino-acid residues (Protein Contact Network PCN, see Fig. 1 for an example). The PCN is the adjacency matrix of a graph, whose links represent the contacts between residues. Our assumption is that the native PCN contains a clue of the protein folding kinetics. In this respect, we introduce three observables that should take into account that Multi-State proteins must be trapped into one or more intermediate states. First of all, we make the hypothesis that MS protein should have more configurational microstates to explore than TS proteins, and we implemented a measure of Network Entropy to quantify this aspect onto a combination of the full PCN with contact potentials as weights for the existing PCN links. Second, from a modified version of the PCN, that keeps only long-range contacts but preserving network connectivity, we evaluated the spectrum of the Laplacian matrix, since it has been shown that its vibrational properties can be used to model experimental data1920. Finally, in order to measure the folding cooperativity21, we evaluate the fraction of sequence separation (diagonals of the full PCN) that do not contain residue contact pairs. The rationale of this measure is that the more diffuse is the cooperation (most of the diagonal participate) the less probable is to be trapped in intermediate states. In order to keep these observables 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.

Figure 1

An example of different PCN representations for protein 1A6N, with MS folding kinetics and 151 C residues.

In Panel A the whole PCN is displayed. In Panel B, once calculated the number of diagonals b = 4 needed to break the protein network in more than one component, 3 backbone diagonals were removed. In Panel C the related Long-range Interaction Network (LIN) is shown1114 in which no backbone diagonal is present.

Results

As previously mentioned, we introduced three main groups of observables (see Methods), i.e.

Network Entropy-based ratio S, 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 λ, λ, λ, related to the highest-frequency vibrational modes of a modified PCN, in which the backbone is removed.

Inter-residue link density R0, 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 Machine13) 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 not significantly increased. The combinations of the largest Laplacian eigenvalues λ, λ, λ with the link density R0 produce the best performances of classification, with a top-score value given by the couple (λ, R0), 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 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 R0 (84.50% ± 1.30% correctly classified proteins, MCC = 0.67 ± 0.03). Classifying without cross-validation, i.e. using the entire set of (λ, R0) 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: (λ, R0) and (S, R0). 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.

Table 1

Classification performances of the newly defined observables and Matthews correlation coefficient MCC, based on quadratic discriminant analysis.

(a) Classification performances
%	λN	λN−1	λN−2	R0	SR
λN	76.6 ± 1.3	74.7 ± 1.4	70.7 ± 1.8	85.2 ± 1.4	78.4 ± 1.2
λN−1		76.7 ± 1.4	70.7 ± 1.4	88.3 ± 1.1	75.9 ± 2.3
λN−2			77.6 ± 1.1	87.3 ± 1.4	76.5 ± 1.9
R0				75.9 ± 1.2	84.5 ± 1.3
SR					80.4 ± 1.8
(b) Matthews correlation coefficient MCC
λN	0.57 ± 0.02	0.53 ± 0.02	0.46 ± 0.03	0.69 ± 0.03	0.60 ± 0.02
λN−1		0.58 ± 0.02	0.46 ± 0.03	0.76 ± 0.02	0.57 ± 0.04
λN−2			0.59 ± 0.02	0.74 ± 0.03	0.56 ± 0.04
R0				0.52 ± 0.02	0.67 ± 0.03
SR					0.59 ± 0.04

The tables show the performances of couples of observables, with the performance of the single observables along the diagonal; the best overall performance is bold-typed. The results presented are the average values of 10-fold cross-validation over 10000 instances and their standard deviation.

Figure 2

Scatterplots of two top-ranking classification couples of observables.

Panel A: λ and R0 (classification = 88.3% ± 1.1%, MCC = 0.76 ± 0.02). Panel B: S and R0 (classification = 84.5% ± 1.3%, MCC = 0.67 ± 0.03).

Since in previous studies9 it has been shown that the chain length N 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 as a variable for discrimination. In our dataset, the N 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 CO1213 and the long range contact order LRCO14. In the considered dataset the couple (N, 〈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 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: 81.4% ± 1.8, MCC = 0.61 ± 0.04).

Table 2

Classification performances of the observables used in literature and Matthews correlation coefficient MCC.

(a) Classification performances
%	NC	〈h〉	CO	LRCO
NC	78.2 ± 1.5	73.7 ± 2.1	78.4 ± 1.7	80.0 ± 1.1
〈h〉		57.3 ± 1.4	72.6 ± 2.8	62.7 ± 1.9
CO			68.4 ± 0.7	70.4 ± 2.1
LRCO				54.4 ± 2.4
(b) Matthews correlation coefficient MCC
NC	0.54 ± 0.03	0.44 ± 0.05	0.54 ± 0.04	0.57 ± 0.02
〈h〉		0.17 ± 0.03	0.44 ± 0.06	0.22 ± 0.04
CO			0.36 ± 0.01	0.44 ± 0.04
LRCO				0.14 ± 0.05

The tables show the performances of couples of observables, with the performance of the single observables along the diagonal; the best overall performance is bold-typed. The results presented are the average values of 10-fold cross-validation over 10000 instances and their standard deviation.

We also performed the same analysis on a different and more recent dataset7 (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 R0 with S had a 80.9% ± 1.3% ratio of correctly classified protein with MCC = 0.58 ± 0.03, while R0 for λ 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 R0 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, λ and R0 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 acoustics23). 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, 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, 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 Finkelstein24. 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 Cs 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 is the 3D Euclidean distance between the i and j residues. The protein contact network PCN is then obtained by choosing an upper threshold of 8 Å25262728:

In order to build our Entropy-based observables, we retrieved also data regarding amino acidic interactions, such as hydropathy indexes2930 and contact potentials, namely, 20 × 20 matrices describing the interactions between the 20 side-chains313233. 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 potentials3233 are very similar. For each protein, starting from the known residue sequence, we define the contact potential matrix P in which

where r is the amino acid residue corresponding to the i C, and the matrix P has the same size of the related PCN. Since P 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

The first observable introduced is associated with the-so called Entropy of a network ensemble3435. 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 ageing36. 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} 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 for two different ensembles, with a different number of constraints: in the first ensemble, we only fix the strength sequence {s} of the protein network (S), while in the second ensemble (S) we fix both strength sequence {s} and degree sequence {k}. The degree sequence {k} and the strength sequence {s} 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 π(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

where the average values of strength sequence and degree sequence over the network ensemble are enforced to match the real features of the selected protein, i.e. and . The network entropy observable S is defined as the ratio between the two entropies with a different number of constraints

with S ≥ S given the fewer number of constraints. The closer is the value of S to 1, the less relevant is the role of the degree sequence constraint {k}, and thus most of the information on possible PCN configurations is enclosed in the strength sequence {s} only. On the contrary, a large value of S implies that the given strength sequence is compatible with a larger number of degree sequences (corresponding to more possible PCNs). Thus, MS proteins could in principle have larger S values than TS proteins: having more stable (or metastable) configurations available could be reflected in a larger number of available configurations as measured by S.

Laplacian-based observables λ

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 visualization3738, and it has already been applied to characterize protein features39. Given an adjacency matrix A without self loops, we define the Laplacian operator as

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 membrane40. 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 state41, 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 PCN1114, 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 {λ}, the largest eigenvalues of L rescaled by the number of residues N. This rescaling was chosen because we observed a dependence of the largest eigenvalues on N, 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

contain all the links between residues with a sequence separation equal to d (ranging from 1 to N − 1) with respect to the protein backbone. The observable R0 is defined as the ratio between the number of d-diagonals without links and the number of residues N of the protein, i.e.

where δ is the Kronecker delta. For a given protein, a low value of R0 implies that the residues interact at many different levels of sequence separation (different values of d). On the contrary, a high value of R0 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 is the number of true positives, T the number of true negatives, F the number of false positives and F 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, 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 folding22. Each amino acid is associated with a hydropathy index h, 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 refers to the hydropathy index of residue i when the Kyte-Doolittle (KD) scale29 is considered. The average hydrophobicity 〈h〉 has been often coupled with N to classify the protein folding kinetics.

We also considered the classification power of structural topology measures such as contact order1213

where L is the total number of contacts for the given PCN, and long range contact order14

which were used before for TS/MS classification purposes.

Additional Information

How to cite this article: Menichetti, G. et al. Network measures for protein folding state discrimination. Sci. Rep. 6, 30367; doi: 10.1038/srep30367 (2016).