BACKGROUND
an important objective of molecular biochemistry is a detailed analysis of protein characteristics like functionality, stability, and dynamics.
this is a laborious and time consuming task due to the many aspects of protein function and the large spectrum of experimental methods required for their determination.
ideally, one would characterize experimentally the contribution of each individual amino acid residue, which is however not feasible for larger proteins.
this is why the biochemical assessment of proteins has to concentrate on a relatively small number of residues.
in enzymes, these are the residues directly involved in catalysis and substrate binding; resulting annotations can be found in dedicated databases like pdbsum  <cit> .
however, there are no equivalent databases available when one has to identify residues which are important for stability or other characteristics.
due to the enormous success of genome sequencing projects, the sequences of more than  <dig>  <dig> protein families  are known at date and thus, methods of computational biology are of utmost importance to support their characterization.
a large number of in silico approaches are at hand to identify important residues.
often, a family-specific multiple sequence alignment  is the main data source to elucidate the role of the residues; for latest reviews see refs.
<cit> .
most effective is the assessment of residue variation deduced from the corresponding msa columns.
the success of these analyses can be explained with the biochemical properties of the residues: for example, in most cases only one residue-type fulfills all critical requirements at catalytic sites, which prohibits a mutation.
accordingly, a strict residue conservation is a strong indicator signaling functionally important residues  <cit> .
in contrast, a prevalent but not exclusively found amino acid is often important for protein stability  <cit> , which similarly holds for ligand-binding sites  <cit> .
interestingly, these less conserved residue positions may bear a pattern indicative of dependencies in the occupancy of two or more positions.
the importance of these correlation signals and their consequences have long been realized  <cit> .
quite different approaches have been introduced to identify correlated residue pairs; see e.g.
refs.
<cit> .
unfortunately, these correlation signals, which are due to constraints imposed by the local environment of a residue, can be disturbed by neutral mutations.
if an msa contains sequences from many closely related species, neutral mutations in a predecessor may give rise to a strong correlation signal.
thus, the elimination of highly similar sequences improves the quality of correlation analysis  <cit> .
additionally, other approaches have been proposed to eliminate signals induced by a common evolutionary path of the proteins  <cit> .
all these methods for the analysis of correlation patterns are aimed at the identification of pairs of residues, which are functionally or structurally important.
more specific methods enable us to predict residue contacts.
for the latter application, transitive dependencies, which by definition interlink several pairs of residues, have to be eliminated as well  <cit> .
different approaches have proven applicable and these algorithms have been named global methods  <cit> .
among them are psicov  <cit> , dca  <cit> , and evfold  <cit> .
the common idea of global methods is to treat pairs of residues as mutually dependent entities and to minimize the effects of transitive covariation and phylogenetic noise.
in contrast, most algorithms like those described in refs.
<cit>  do not correct for transitive dependencies.
these approaches have been named local methods  <cit>  as they assume that pairs of residue positions are statistically independent of other pairs.
due to chaining effects, the identified residue positions constituting a pair, can be near to each other or far apart in the protein’s structure.
most of the local methods rely in one way or another on assessing the mutual information, which is commonly based on shannon’s entropy  <cit> .
thus, these local methods deduce a measure for mutual dependencies solely from the amino acid frequencies observed at the positions under study.
consequently, the biochemical properties of the residues are ignored, which may degrade the performance of the algorithm.
one of these local methods is the algorithm h2r  <cit> , which identifies in a first step mutual dependencies between pairs of residue positions and scores in a second step each residue position k by the conn-value, which is the number of significant pairwise correlations it is involved in.
mutagenesis studies with two enzymes demonstrated that positions with high conn-values have an increased probability of being important for enzyme function or stability  <cit> .
as we were interested to further improve performance of h2r in terms of specificity, we implemented h2rs, which additionally takes into account substitution frequencies for residue pairs.
moreover, h2rs determines a specific p-value for each analysis of a residue pair, which facilitates the selection of significant correlation signals.
to further standardize the analyses, h2rs normalizes the resulting conn-values to z-scores, which we named conz-values.
using a testbed consisting of  <dig> enzymes, we demonstrated in a comparison with the predecessor algorithm h2r and two alternative algorithms that a larger fraction of residues endowed by h2rs with high conz-values are located near ligand binding sites.
additionally, we studied in detail the predictions of h2r, h2rs, and the global method psicov for five well characterized enzymes.
it turned out that the outcome of local and global methods overlapped only marginally and that residues with high conz-values are functionally or structurally significant.
RESULTS
utilizing the von neumann entropy to improve the identification of correlated mutations
a classification or regression problem can be solved optimally by means of sophisticated classifiers like support vector machines, given that positive and negative examples are at hand during training.
however, there is no clear definition of a correlated mutation.
this is why we cannot model the positive cases and can only characterize as precisely as possible the standard situation.
thus, to create a null model, we can deduce mean substitution frequencies for residue pairs from a large number of samples by analyzing known proteins.
these substitution frequencies reflect the expected case and will allow us to identify more precisely deviations, which indicate mutual dependencies.
based on this argument, we anticipated an improvement in the identification of correlated mutations, if we additionally take into account the similarity of residue pairs together with their frequencies.
note that frequencies are the only source of information in the standard approach.
the algorithm h2r is based on shannon’s information theory  <cit>  and computes for each pair of residue positions k, l the term u  according to
uk,l=2hk+hl-hk,lhk+hl
here, h is the entropy of an individual column k
hk=-∑i=120paiklnpaik
and paik is the probability of amino acid ai at position k. the entropy h of two variables  k and l is
hk,l=-∑i,jpaik,ajllnpaik,ajl
and paik,ajl is the probability of the amino acid pair  at positions k and l. in this context, frequency values deduced from the columns of an msa served as estimates for probabilities.
due to normalization, u is a more reliable indicator of co-evolution than a raw mutual information value  <cit> .
as we were interested to improve specificity, we searched for an information theoretical concept allowing the integration of substitution frequencies determined for residue pairs.
the von neumann entropy  is a generalization of the classical shannon entropy and has been introduced in quantum statistical mechanics  <cit> .
in computational biology, the vne has been used successfully to characterize the conservation of individual residue positions  <cit> .
extending this concept to residue pairs, we aimed at a novel uvne term to replace u.
the core concept of the vne is the utilization of a so-called density matrix ρk,l, that is, a positive definite matrix whose trace  equals to  <dig>  ρk,l can be computed for each pair k, l according to:
ρk,l=pk,lapk,l
here, pk,l=diagp <dig> …,p <dig> and p1…p <dig> are the pairwise amino acid probabilities paik,ajl specified in formula .
these probabilities satisfy the normalization condition ∑i=1400pi= <dig>  a is a  <dig> ×  <dig> matrix that assesses the similarity of residue pairs and it is this matrix that allows us to model substitutions more precisely.
if a is equal to the identity matrix, then the vne is equal to the shannon entropy, that is, vne = h; see below.
based on ρk,l, the von neumann entropy vne can be calculated as
vnek,l=vneρk,l=-∑i=1400λilogλi
by means of the eigenvalues λi of ρk,l.
normalization analogous to formula , which reduces phylogenetic crosstalk, requires corresponding values vne and vne.
for their determination, we applied partial traces  <cit>  on ρk,l to deduce two density matrices ρkk,l and ρlk,l, which are specific for a pair of columns k, l. the elements of ρkk,l and ρlk,l were named si,j and ti,j, respectively, and were computed according to
si,j=∑u=120r20i-1+u,20j-1+u
and
ti,j=∑u=120r20u-1+i,20u-1+j
where ri,j denotes the appropriate entry in the density matrix ρk,l.
thus, this approach allows us to deduce all entropy terms from the density matrix ρk,l, which eliminates normalization problems.
we calculate the vneρmk,l for the residue positions m ∈ {k, l} analogously to equation  based on the eigenvalues λi of the  <dig> ×  <dig> matrix ρmk,l:
vneρmk,l=-∑i=120λilogλi
finally, we define the normalized uvne-value:
uvnek,l=vneρkk,l+vneρlk,l-vneρk,lvneρkk,l+vneρlk,l
computing these values is straightforward, if a matrix a is at hand.
computing a matrix a to assess the similarity of residue pairs
in the case of correlated mutations, the matrix a is a prerequisite to assess the similarity of residue pairs that occur in homologous proteins at corresponding positions.
to determine the  <dig> ×  <dig> values of a, we followed the concept introduced for the blosum approach to score the similarity of amino acid residues based on substitution frequencies  <cit> .
here, we extended this concept to pairs of residues, as similarly used in p2pmat  <cit> .
a precompiled and redundancy free set of  <dig>  <dig> protein 3d structures  <cit>  offered by the pisces server  <cit>  was used as a representative sample.
for each protein, the corresponding msa was taken from the hssp database  <cit>  to deduce pairwise substitution frequencies.
based on the 3d structure, those residue pairs k, l were identified which contacted each other in the protein.
the distances between the centers of any two heavy atoms belonging to one residue each were determined and alternatively the cut-offs  <dig>  Å and  <dig>  Å were chosen to select contacting pairs.
these values correspond to the interval of distances used during casp <dig> to identify contacts between residues and ligands  <cit> .
for these cut-offs, we deduced  <dig>  <dig>  <dig> and  <dig>  <dig>  <dig> contacts from  <dig>  <dig>  <dig> sequences, respectively.
then, substitution frequencies were determined by analyzing the corresponding columns of the msas; see figure  <dig> and methods.
the values of the two corresponding matrices a <dig>  and a <dig>  were normalized to affirm symmetry.
their comparison indicated highly similar values indicating that this distance is no critical parameter, which is in agreement with findings of casp <dig>  <cit> .
as we wanted to consider the larger number of contacts for the determination of the similarity values, we chose a = a <dig>  for all further computations.
this matrix is available as additional file  <dig>
a p-value for the strength of correlation signals deduced from a randomization test
our next goal was to introduce a universally applicable statistical measure for the strength of the pairwise correlations, and we opted for a randomization test.
here, the null hypothesis is that there is no dependency in the pairwise frequencies.
thus, we can assess the strength of each pairwise correlation by shuffling the content of the two columns k, l under study  <cit> .
as we shuffle the content column-wise, the entropy  of the two individual columns remains constant; however, we simultaneously degrade the putative correlation between the two residue positions.
then, we can compare the uvne value deduced from the unaltered combination of residue pairs with a distribution of uvne values resulting from many shuffling rounds.
thus, we can rate the correlation strength for this specific combination of residue pairs observed in columns k and l. consequently, if the uvne values are similarly large or surpass the uvne value, the correlation is statistically not significant.
on the other hand, if all uvne values are significantly lower, then this specific uvne value signals a pronounced dependency in the occupancy of the two residue positions, which indicates correlated mutations.
to compute this p-value efficiently, the number of randomized samples has to be minimized.
moreover, we need a statistical model which has to be valid, if the number of residue types is relatively small which may cause a skewed distribution.
the more conserved the residue positions are, the fewer pairwise frequencies occur and the more the distribution of pairwise frequencies deviates from a normal distribution; compare figure  <dig>  as we wanted to assess the extremeness of the uvne values, we selected a gumbel distribution  <cit>  for modeling.
this distribution is specified by only two parameters μ and β that can be determined in a straightforward manner; see methods and formulae 12– <dig>  to confirm that the gumbel distribution is a proper model, we determined histograms consisting of  <dig> uvne values each for all of  <dig>  <dig>  <dig> pairs of residue positions in our dataset.
prior to the computation of the next uvne value, columns were shuffled 100m times, where m is the number of sequences in the respective msa.
a kolmogorov smirnov test  <cit>  with α =  <dig>  confirmed that the distributions of these uvne values and the deduced gumbel distribution did not differ significantly for  <dig> % of all cases.
using the same dataset, we additionally made clear that the two parameters μ and β can be estimated with adequate precision after  <dig> instances of randomization.
thus, to compute a specific p-value for each residue pair, it is sufficient to compute  <dig> uvne values and to determine one value of the fitted cumulative gumbel distribution.
for a protein of length l, we apply this test n = l/ <dig> times, which suggests to introduce the bonferroni correction  <cit>  in order to reduce the number of false positive results caused by the frequent application of the test.
thus, a corrected cut-off c_o for the corresponding p-value p is
c_ok,l=μ-βloglog11-p/n.
c_o allows for a statistically meaningful and content specific selection of correlated residue positions.
μ and β are defined by formulae  and ; see methods.
for the identification of correlated mutations, a p-value p has to be selected beforehand.
then, all pairs of residue positions with uvne ≥ c_o are utilized to compute conn-values by counting the number of significantly correlated pairs k  is part of.
to further alleviate the comparison of different test sets, conn-values were transformed to z-scores conz; see formula .
an in silico testbed for the assessment of correlation methods
the ultimate validation of a correlated mutation is a biochemical experiment, which is frequently based on the replacement of residues by the standard amino acid alanine.
however, the detailed experimental analysis of a large number of mutations introduced in one protein like dihydrofolate reductase  <cit>  is still the exception.
this lack of reliable results impedes establishing a bona fide testbed for correlation methods and enforces the use of in silico surrogates.
it is known that many correlated mutations are in close proximity to functional sites  <cit> .
thus, a testbed has been created that consists of  <dig> enzymes whose structure and active site residues are well characterized  <cit> .
to assess the quality of correlation analysis, residue positions around functional sites have been counted as positives and all others as negatives  <cit> .
to broaden the statistical basis, we compiled a non-redundant dataset of  <dig> enzymes, whose functional sites, i.e.
catalytic and binding sites, are known and which are represented by a pdb structure and a corresponding msa in the hssp database; see materials.
to determine performance values,  <dig>  <dig> residues were classified and the distances between van der waals spheres were determined.
we regarded all  <dig> residues with a maximal distance of  <dig> Å to a functional site as positive cases and all other  <dig>  <dig> residues as negative cases.
the classification and the resulting performance depends on the chosen p-value and the cut-off for conz.
this is why we tested several combinations and summarized results in table  <dig>  for a p-value between 10- <dig> and 10- <dig> and a conz-threshold of  <dig> , the specificity was between  <dig>  and  <dig>  and precision was between  <dig>  and  <dig> .
for the p-value 10- <dig> and the conz-threshold of  <dig> , specificity was  <dig>  and precision  <dig> .
for p-values ≤ 10- <dig> and conz =  <dig>  the performance reached a plateau.
the comparison with the predecessor algorithm h2r made clear that the novel algorithm performed better: specificity and precision were up to 3% higher.
additionally, we analyzed the same dataset with the algorithms cmat  <cit>  and sca  <cit> , which predict pairs of correlated residue positions.
standalone versions as of february  <dig> were downloaded and applying the same criteria as above, performance was determined.
cmat was used with default parameters.
for sca, we selected three cut-off values  <dig> ,  <dig> , and  <dig> .
performance values were added to table  <dig>  cmat reached a specificity of  <dig>  and a precision of  <dig> .
for sca, the specificity increased from  <dig>  to  <dig> , and the precision from  <dig>  to  <dig> , for the cut-offs  <dig>  and  <dig> .
these results indicate that residue positions predicted by h2rs are more likely close to functional sites.
moreover, the number of false positives is lower, as indicated by the higher precision values determined for h2rs.
these numbers are a rough estimate of the algorithm’s performance due to the limitations of the in silico testbed.
however, all other alternative methods of performance evaluation  <cit>  are not applicable here: these are the analysis of simulated msas, the determination of the residues’ spatial distance or an assessment of free energy differences derived from double mutants.
in silico
testbed
for all programs, specificity and precision were deduced from the analysis of  <dig> enzymes with known catalytic and binding sites.
residues with a maximal distance of  <dig> Å to a functional site were regarded as positives.
all other residues were regarded as negatives.
h2r and cmat were used with default settings.
for h2rs, the cut-off was applied to the p-value.
for sca, three cut-off values were chosen.
an assessment of predicted coevolving residues in well-characterized enzymes
to evaluate performance of our algorithm in more detail, we analyzed the h2rs predictions for five well studied enzymes: three enzymes from tryptophan biosynthesis, named trpa, trpb, trpc, dihydrofolate reductase , and hexokinase .
trpa and trpb constitute the heteromeric tryptophan synthase complex, which catalyzes the final reaction of indole-3-glycerole phosphate and serine to tryptophan.
trpa cleaves indole-3-glycerol phosphate to glyceraldehyde-3-phosphate and indole, which is transported through a hydrophobic tunnel to the active center of trpb.
there, tryptophan is synthesized from serine and indole  <cit> .
for the localization of predicted residue positions, we utilized the 3d dataset with pdb id 1kfc, which is the trpa/trpb complex from salmonella typhimurium .
the enzyme indole-3-glycerol phosphate synthase  catalyzes the ring closure of an n-alkylated anthranilate to a 3-alkyl indole derivative, which is the fourth step in the tryptophan biosynthesis.
it adopts the widespread 8-barrel fold and has been studied in detail  <cit> .
here, we utilized the structure of trpc from sulfolobus solfataricus .
dhfr catalyzes the reduction of dihydrofolate to tetrahydrofolate via hydride transfer from nadph.
it has been found in most organisms and plays a critical role for cell proliferation and cell growth  <cit> .
we utilized the structure determined for dhfr from escherichia coli .
the hexokinase from schistosoma mansoni  is the first enzyme in the glycolytic pathway and catalyzes the transfer of a phosphoryl group to alpha-6-glucose .
the 3d crystal structure contains so <dig> anions in the catalytic cleft  <cit> .
smhk is the only enzyme of a larger set that has been analyzed previously by correlation analysis and for which the msa  was available online.
to generate smhk_cma, the authors have used a sophisticated protocol to merge several structure based msas  <cit> .
although local and global methods of correlation analysis have different objectives, we were interested to determine the overlap of their predictions.
this is why we also compared the outcome of h2rs and psicov  <cit> , which is a global method predicting residue contacts.
for psicov we analyzed the top l/ <dig> predictions, which is the recommended default for a protein sequence of length l. an msa was created for each enzyme by using delta-blast  <cit>  with the options max target threshold  <dig> and expect threshold 10- <dig>  the resulting sequences were realigned by means of mafft  <cit>  in linsi mode.
we were interested in an assessment of the most specific h2rs predictions.
this is why we chose the low cut-off 10- <dig> for the p-value and a conz-threshold of  <dig> .
to allow for a comparison, we also listed the conz-values for all residues predicted by h2r in table  <dig>  residues were regarded as functionally important, if they were close to a functional site specified in pdbsum  <cit> .
thus, all direct neighbors in the sequence were chosen and all residues with a 3d distance of maximally  <dig> Å .
for the enzymes sttrpa, sttrpb, sstrpc, ecdhfr, and smhk, h2r and h2r were used to identify important residue positions.
for these residues, annotation was deduced from literature.
the first column lists the name of the enzyme.
the second column gives the residue and its position.
the third column gives the conz-value deduced by h2rs from all uvne()-values based on a p-value of 10- <dig>  the column h2r lists mean conn-values resulting from  <dig> randomization tests.
the column psicov lists the number of contacting pairs the residue belonged to.
the last column lists the role of the residues, for details see results.
“nd” indicates that we did not find clues to the function of this residue.
sttrpa consists of  <dig> residues, and h2rs predicted two important residues, namely l <dig> and l <dig>  both residues are in close proximity to the substrate; see figure  <dig>  h2r predicted l <dig>  s <dig>  a <dig>  i <dig> and l <dig>  s <dig> stabilizes the inactive conformation of the active center  <cit> .
a <dig> and i <dig> are near the active site and l <dig> belongs to the trpa/trpb interface  <cit> .
l <dig> and l <dig> also belong to the  <dig> l/ <dig> predictions of psicov; see table  <dig>
for sttrpb, h2rs predicted  <dig> of the  <dig> residues as being important; see figure  <dig>  t <dig>  q <dig>  and v <dig> are in close proximity to the substrate binding residue k <dig>  <cit> .
c <dig> and f <dig> are located at the end of the hydrophobic tunnel  <cit>  and t <dig> and s <dig> are metal binding sites  <cit> .
g <dig> is important for the coordination of ion binding  <cit> , and s <dig> and p <dig> are in close proximity to the bound sodium ion.
m <dig> is in contact with f <dig> and s308; see above.
the role of the two residues c <dig> and q <dig> is unknown to us.
in contrast, h2r predicted five of these residues, namely c <dig>  q <dig>  p <dig>  f <dig>  s <dig>  and additionally p <dig>  g <dig>  l <dig>  and f <dig>  f <dig> is a metal binding site, g <dig> is near the binding site for the substrate and the function of p <dig> and l <dig> is unknown to us.
of the  <dig> h2rs predictions,  <dig> belong to the  <dig> l/ <dig> contacting residues predicted by psicov; see table  <dig>
for sstrpc, h2rs predicted  <dig> important positions; see figure  <dig>  v <dig> is near the active site.
i <dig> and l <dig> are catalytically important: after replacing each of these two residues by alanine, the activity of trpc dropped 30-fold  <cit> .
a <dig> is located next to the substrate binding site e <dig> and the catalytic residue s <dig>  <cit> ; s <dig> is known to be a phosphate binding site  <cit> .
the role of the two residues i <dig> and i <dig> is unknown to us.
h2r detected the phosphate binding site s <dig>  the catalytically important residue i <dig>  plus the residues a <dig>  y <dig>  m <dig>  n <dig>  and l <dig>  a <dig>  m <dig>  and n <dig> are near the active site.
the role of l <dig> is unknown however, the l187a mutation has no drastic effect on function and stability  <cit> .
the function of y <dig> is unknown to us.
all of the residue positions predicted by h2rs belonged to the  <dig> l/ <dig> contacting residue pairs predicted by psicov; see table  <dig>
for ecdhfr, h2rs predicted six important residue positions; see figure  <dig>  w <dig> and k <dig> are contacting the substrate, whereas a <dig> and m <dig> are in close proximity to the binding site a <dig> and the catalytic site i <dig>  respectively  <cit> .
h <dig> plays a significant role in the network of coupled motions required for a hydride transfer  <cit>  and a mutation of g <dig>  which lies in proximity of nadph, reduced the hydride transfer rate  <cit> .
the predecessor algorithm, h2r, identified g <dig>  h <dig>  plus d <dig>  whose function is unknown to us.
of the above sites, only a <dig> was an element of the  <dig> l/ <dig> predictions of psicov; see table  <dig>
smhk consists of a hk type- <dig>  and a hk type- <dig> domain ; see entry q <dig> of uniprot  <cit> .
h2rs identified  <dig> suspicious residues , which we number according to the pdbsum  <cit>  entry 1bdg.
a <dig>  c <dig>  and a <dig> are located at the very end of domain  <dig>  whereas c <dig> occurs at the very beginning of domain  <dig> and these four residues are flanking a ß-turn  <cit> .
k <dig> is a neighbor of q <dig> that binds glc, v <dig> is a neighbor of i <dig>  and of t <dig>   <cit> .
v <dig> is a neighbor of g <dig> and s <dig> that both bind so <dig> <cit> .
t <dig> is close to c <dig> .
only for two residues, namely t <dig> and v <dig>  their role is unknown to us; however both residues are located at the domain interface at a distance of not more than  <dig>  Å. h2r found c <dig> and additionally d <dig>  whose function is unknown to us.
<dig> of the h2rs predictions were in the  <dig> l/ <dig> predictions of psicov.
when utilizing the msa smhk_cma, h2rs predicted only three residues with a positive conz-value, which is given in brackets: k <dig> , t <dig> , and c <dig> .
t <dig> binds glc, and k <dig> is located next to the glc binding e <dig>  <cit> .
for c <dig>  see above.
interestingly, in the  <dig> sequences remaining in the msa after filtering, residue positions  <dig> and  <dig> were occupied in not more than 43% by cysteines, which form a disulfide bridge that stiffens the orientation of the two domains  <cit> .
alternatively, the following residue pairs were observed with more than 2% frequency: st , gv , sm , rt , hp , av  and ra .
these distinct pairwise combinations support nicely the idea of mutual dependencies and pairwise correlations.
although the number of cases is small, these well characterized proteins allow for a more realistic assessment of the prediction performance.
altogether, h2rs predicted  <dig> important residues and h2r  <dig>  respectively.
false positives were  <dig>  in the case of h2rs and  <dig>  in the case of h2r.
thus, the resulting precision is  <dig>  for h2rs and  <dig>  for h2r.
these results emphasize the relatively high specificity reached by computing conn-values and additionally suggest a considerable improvement for the novel algorithm.
discussion
h2rs is a major improvement over h2r
for all well-characterized enzymes studied in results, h2rs predicted a larger number and a higher fraction of residue positions for which we could rationalize an important role in function or stability.
here, we concentrated on the analysis of residues with a conz-value ≥  <dig> .
generally, this detailed analysis of five enzymes signals more precisely than the assessment of our in silico testbed the improved specificity of h2rs.
it was achieved i) by replacing shannon’s entropy by the von neumann entropy and ii) by integrating a more sensitive statistical approach that adapts to the composition of each pair of msa columns.
based on this dataset, we can expect a 10% increase in specificity to nearly 90%.
however, this improvement has to be paid with a much longer execution time: computing the von neumann entropy requires the determination of eigenvalues, which is time-consuming and the determination of p-values further increases the execution time by a factor of  <dig>  one way of accelerating the calculation of entropy values might be an application of the rényi entropy  <cit> , which is a generalization of the von neumann entropy.
for 0 < α ≠  <dig>  the α-rényi entropy is given by α-rek,l=11-αlog∑i=1400λiα and for α →  <dig>  we recover the neumann entropy vne.
interestingly, for α =  <dig>  the calculation of the α-rényi entropy does not require the eigenvalues of the matrix ρk,l but only the diagonal entries of the square of ρk,l, which drastically speeds up the computation.
however, it has not been tested yet whether the rényi entropy allows the adequate modeling of biological phenomena like residue substitutions.
global and local methods of correlation analysis complement each other
one goal in the design of h2r, which is a local method, was the identification of individual residue positions important for protein function or stability.
this is why we introduced the conn-value.
for two enzymes it has been shown that positions with high conn-values have an increased probability of being important for enzyme function or stability  <cit> .
the results presented here further confirm the high specificity to be gained with local methods, which is in agreement with data from the literature; see e.g.
refs.
<cit> .
the results obtained for smhk emphasize that not all correlated mutations are due to functional constraints:  <dig> of  <dig> residues with high conz-values were located at the domain interface and two of them  belong to a disulfide bond that interlinks the domains in some of the homologous proteins.
the other residue combinations observed at these two positions illustrate nicely that they were to a great extent occupied by unique residue pairs.
moreover, these findings emphasize a limitation of the in silico testbed.
structurally important residues often lay far apart from the catalytic center  <cit> .
as shown above, some bear a strong correlation signal and are identified by h2rs.
however, these hits are regarded as false positives and deteriorate the performance values deduced from the testbed.
whereas local methods consider transitive correlations as well, global methods aim at eliminating these dependencies.
the outcome of h2rs and the l/ <dig> predictions of the global method psicov overlapped only for  <dig> of  <dig> residue positions; see table  <dig>  this result can be explained by the scope of the methods: according to the desired function, global methods identify contacting residue pairs which are not necessarily enriched near functional sites.
using the msa smhk_cma, h2rs predicted only three residues known to be functionally important, albeit two with low conz-values.
using the same dataset, the algorithm comulator, which aims at identifying perturbations  <cit> , detected a network of six residue positions that surround the active site.
their occupancy almost perfectly separated the two main groups of glucokinases  <cit> .
in summary, these findings highlight the pros and cons of the different approaches and suggest that they supplement each other quite well.
msas have to be prepared carefully
a critical parameter of correlation analysis is the preparation of the input, i.e.
the msa.
for the prediction of intra-protein residue contacts, a strong correlation between the number of homologs and the prediction strength has been shown, which further increased, if orthologs and paralogs were included in the msa  <cit> .
for the sake of standardization, we used in all cases studied here the same methods of msa preparation without human intervention.
additionally we chose identical and very rigorous cut-offs for the identification of important residue positions.
this rigid protocol might be the reason for the considerably differing number of predictions: using the cut-off conz ≥  <dig>  and a p-values of 10- <dig>  h2rs predicted for sttrpa only  <dig>  but for sttrpb  <dig> important residue positions.
these differences suggest for the user an individual adjustment of the parameters for each protein family in order to optimize the benefit of correlation analysis.
CONCLUSIONS
the various global and local methods of correlation analysis have their own field of application and supplement each other.
we made plausible that residues in the vicinity of functional sites, which are a large portion of h2rs predictions, do not necessarily belong to residue pairs with the strongest global correlation signal.
the predictions of global methods can be assessed by the 3d distance of the involved residue pairs.
in contrast, the evaluation of local methods is more ambiguous.
due to the lack of a precise definition of a correlated mutation, it is difficult to specify positive cases.
this circumstance has drastic consequences and imposes restrictions to the design and the evaluation of algorithms.
with this in mind, we developed an algorithm that considers pairwise substitution frequencies and assesses the strength of the correlation signal statistically.
we made plausible that in silico testbeds only allow for a rough performance evaluation.
favorable is the detailed analysis of well characterized model systems, which is only feasible for a small number of cases.

 an important objective of molecular biochemistry is a detailed analysis of protein characteristics like functionality, stability, and dynamics. this is a laborious and time consuming task due to the many aspects of protein function and the large spectrum of experimental methods required for their determination. ideally, one would characterize experimentally the contribution of each individual amino acid residue, which is however not feasible for larger proteins. this is why the biochemical assessment of proteins has to concentrate on a relatively small number of residues. in enzymes, these are the residues directly involved in catalysis and substrate binding; resulting annotations can be found in dedicated databases like pdbsum. however, there are no equivalent databases available when one has to identify residues which are important for stability or other characteristics. due to the enormous success of genome sequencing projects, the sequences of more than protein families are known at date and thus, methods of computational biology are of utmost importance to support their characterization. a large number of in silico approaches are at hand to identify important residues. often, a family-specific multiple sequence alignment is the main data source to elucidate the role of the residues; for latest reviews see refs.. most effective is the assessment of residue variation deduced from the corresponding msa columns. the success of these analyses can be explained with the biochemical properties of the residues: for example, in most cases only one residue-type fulfills all critical requirements at catalytic sites, which prohibits a mutation. accordingly, a strict residue conservation is a strong indicator signaling functionally important residues. in contrast, a prevalent but not exclusively found amino acid is often important for protein stability, which similarly holds for ligand-binding sites. interestingly, these less conserved residue positions may bear a pattern indicative of dependencies in the occupancy of two or more positions. the importance of these correlation signals and their consequences have long been realized. quite different approaches have been introduced to identify correlated residue pairs; see e.g. refs.. unfortunately, these correlation signals, which are due to constraints imposed by the local environment of a residue, can be disturbed by neutral mutations. if an msa contains sequences from many closely related species, neutral mutations in a predecessor may give rise to a strong correlation signal. thus, the elimination of highly similar sequences improves the quality of correlation analysis. additionally, other approaches have been proposed to eliminate signals induced by a common evolutionary path of the proteins. all these methods for the analysis of correlation patterns are aimed at the identification of pairs of residues, which are functionally or structurally important. more specific methods enable us to predict residue contacts. for the latter application, transitive dependencies, which by definition interlink several pairs of residues, have to be eliminated as well. different approaches have proven applicable and these algorithms have been named global methods. among them are psicov, dca, and evfold. the common idea of global methods is to treat pairs of residues as mutually dependent entities and to minimize the effects of transitive covariation and phylogenetic noise. in contrast, most algorithms like those described in refs. do not correct for transitive dependencies. these approaches have been named local methods as they assume that pairs of residue positions are statistically independent of other pairs. due to chaining effects, the identified residue positions constituting a pair, can be near to each other or far apart in the protein’s structure. most of the local methods rely in one way or another on assessing the mutual information, which is commonly based on shannon’s entropy. thus, these local methods deduce a measure for mutual dependencies solely from the amino acid frequencies observed at the positions under study. consequently, the biochemical properties of the residues are ignored, which may degrade the performance of the algorithm. one of these local methods is the algorithm h2r, which identifies in a first step mutual dependencies between pairs of residue positions and scores in a second step each residue position k by the conn-value, which is the number of significant pairwise correlations it is involved in. mutagenesis studies with two enzymes demonstrated that positions with high conn-values have an increased probability of being important for enzyme function or stability. as we were interested to further improve performance of h2r in terms of specificity, we implemented h2rs, which additionally takes into account substitution frequencies for residue pairs. moreover, h2rs determines a specific p-value for each analysis of a residue pair, which facilitates the selection of significant correlation signals. to further standardize the analyses, h2rs normalizes the resulting conn-values to z-scores, which we named conz-values. using a testbed consisting of enzymes, we demonstrated in a comparison with the predecessor algorithm h2r and two alternative algorithms that a larger fraction of residues endowed by h2rs with high conz-values are located near ligand binding sites. additionally, we studied in detail the predictions of h2r, h2rs, and the global method psicov for five well characterized enzymes. it turned out that the outcome of local and global methods overlapped only marginally and that residues with high conz-values are functionally or structurally significant.  utilizing the von neumann entropy to improve the identification of correlated mutations a classification or regression problem can be solved optimally by means of sophisticated classifiers like support vector machines, given that positive and negative examples are at hand during training. however, there is no clear definition of a correlated mutation. this is why we cannot model the positive cases and can only characterize as precisely as possible the standard situation. thus, to create a null model, we can deduce mean substitution frequencies for residue pairs from a large number of samples by analyzing known proteins. these substitution frequencies reflect the expected case and will allow us to identify more precisely deviations, which indicate mutual dependencies. based on this argument, we anticipated an improvement in the identification of correlated mutations, if we additionally take into account the similarity of residue pairs together with their frequencies. note that frequencies are the only source of information in the standard approach. the algorithm h2r is based on shannon’s information theory and computes for each pair of residue positions k, l the term u according to uk,l=2hk+hl-hk,lhk+hl here, h is the entropy of an individual column k hk=-∑i=120paiklnpaik and paik is the probability of amino acid ai at position k. the entropy h of two variables k and l is hk,l=-∑i,jpaik,ajllnpaik,ajl and paik,ajl is the probability of the amino acid pair at positions k and l. in this context, frequency values deduced from the columns of an msa served as estimates for probabilities. due to normalization, u is a more reliable indicator of co-evolution than a raw mutual information value. as we were interested to improve specificity, we searched for an information theoretical concept allowing the integration of substitution frequencies determined for residue pairs. the von neumann entropy is a generalization of the classical shannon entropy and has been introduced in quantum statistical mechanics. in computational biology, the vne has been used successfully to characterize the conservation of individual residue positions. extending this concept to residue pairs, we aimed at a novel uvne term to replace u. the core concept of the vne is the utilization of a so-called density matrix ρk,l, that is, a positive definite matrix whose trace equals to ρk,l can be computed for each pair k, l according to: ρk,l=pk,lapk,l here, pk,l=diagp …,p and p1…p are the pairwise amino acid probabilities paik,ajl specified in formula. these probabilities satisfy the normalization condition ∑i=1400pi= a is a × matrix that assesses the similarity of residue pairs and it is this matrix that allows us to model substitutions more precisely. if a is equal to the identity matrix, then the vne is equal to the shannon entropy, that is, vne = h; see below. based on ρk,l, the von neumann entropy vne can be calculated as vnek,l=vneρk,l=-∑i=1400λilogλi by means of the eigenvalues λi of ρk,l. normalization analogous to formula, which reduces phylogenetic crosstalk, requires corresponding values vne and vne. for their determination, we applied partial traces on ρk,l to deduce two density matrices ρkk,l and ρlk,l, which are specific for a pair of columns k, l. the elements of ρkk,l and ρlk,l were named si,j and ti,j, respectively, and were computed according to si,j=∑u=120r20i-1+u,20j-1+u and ti,j=∑u=120r20u-1+i,20u-1+j where ri,j denotes the appropriate entry in the density matrix ρk,l. thus, this approach allows us to deduce all entropy terms from the density matrix ρk,l, which eliminates normalization problems. we calculate the vneρmk,l for the residue positions m ∈ {k, l} analogously to equation based on the eigenvalues λi of the × matrix ρmk,l: vneρmk,l=-∑i=120λilogλi finally, we define the normalized uvne-value: uvnek,l=vneρkk,l+vneρlk,l-vneρk,lvneρkk,l+vneρlk,l computing these values is straightforward, if a matrix a is at hand. computing a matrix a to assess the similarity of residue pairs in the case of correlated mutations, the matrix a is a prerequisite to assess the similarity of residue pairs that occur in homologous proteins at corresponding positions. to determine the × values of a, we followed the concept introduced for the blosum approach to score the similarity of amino acid residues based on substitution frequencies. here, we extended this concept to pairs of residues, as similarly used in p2pmat. a precompiled and redundancy free set of protein 3d structures offered by the pisces server was used as a representative sample. for each protein, the corresponding msa was taken from the hssp database to deduce pairwise substitution frequencies. based on the 3d structure, those residue pairs k, l were identified which contacted each other in the protein. the distances between the centers of any two heavy atoms belonging to one residue each were determined and alternatively the cut-offs Å and Å were chosen to select contacting pairs. these values correspond to the interval of distances used during casp to identify contacts between residues and ligands. for these cut-offs, we deduced and contacts from sequences, respectively. then, substitution frequencies were determined by analyzing the corresponding columns of the msas; see figure  and methods. the values of the two corresponding matrices a and a were normalized to affirm symmetry. their comparison indicated highly similar values indicating that this distance is no critical parameter, which is in agreement with findings of casp. as we wanted to consider the larger number of contacts for the determination of the similarity values, we chose a = a for all further computations. this matrix is available as additional file a p-value for the strength of correlation signals deduced from a randomization test our next goal was to introduce a universally applicable statistical measure for the strength of the pairwise correlations, and we opted for a randomization test. here, the null hypothesis is that there is no dependency in the pairwise frequencies. thus, we can assess the strength of each pairwise correlation by shuffling the content of the two columns k, l under study. as we shuffle the content column-wise, the entropy of the two individual columns remains constant; however, we simultaneously degrade the putative correlation between the two residue positions. then, we can compare the uvne value deduced from the unaltered combination of residue pairs with a distribution of uvne values resulting from many shuffling rounds. thus, we can rate the correlation strength for this specific combination of residue pairs observed in columns k and l. consequently, if the uvne values are similarly large or surpass the uvne value, the correlation is statistically not significant. on the other hand, if all uvne values are significantly lower, then this specific uvne value signals a pronounced dependency in the occupancy of the two residue positions, which indicates correlated mutations. to compute this p-value efficiently, the number of randomized samples has to be minimized. moreover, we need a statistical model which has to be valid, if the number of residue types is relatively small which may cause a skewed distribution. the more conserved the residue positions are, the fewer pairwise frequencies occur and the more the distribution of pairwise frequencies deviates from a normal distribution; compare figure  as we wanted to assess the extremeness of the uvne values, we selected a gumbel distribution for modeling. this distribution is specified by only two parameters μ and β that can be determined in a straightforward manner; see methods and formulae 12– to confirm that the gumbel distribution is a proper model, we determined histograms consisting of uvne values each for all of pairs of residue positions in our dataset. prior to the computation of the next uvne value, columns were shuffled 100m times, where m is the number of sequences in the respective msa. a kolmogorov smirnov test with α =  confirmed that the distributions of these uvne values and the deduced gumbel distribution did not differ significantly for % of all cases. using the same dataset, we additionally made clear that the two parameters μ and β can be estimated with adequate precision after instances of randomization. thus, to compute a specific p-value for each residue pair, it is sufficient to compute uvne values and to determine one value of the fitted cumulative gumbel distribution. for a protein of length l, we apply this test n = l/ times, which suggests to introduce the bonferroni correction in order to reduce the number of false positive results caused by the frequent application of the test. thus, a corrected cut-off c_o for the corresponding p-value p is c_ok,l=μ-βloglog11-p/n. c_o allows for a statistically meaningful and content specific selection of correlated residue positions. μ and β are defined by formulae and; see methods. for the identification of correlated mutations, a p-value p has to be selected beforehand. then, all pairs of residue positions with uvne ≥ c_o are utilized to compute conn-values by counting the number of significantly correlated pairs k is part of. to further alleviate the comparison of different test sets, conn-values were transformed to z-scores conz; see formula. an in silico testbed for the assessment of correlation methods the ultimate validation of a correlated mutation is a biochemical experiment, which is frequently based on the replacement of residues by the standard amino acid alanine. however, the detailed experimental analysis of a large number of mutations introduced in one protein like dihydrofolate reductase is still the exception. this lack of reliable results impedes establishing a bona fide testbed for correlation methods and enforces the use of in silico surrogates. it is known that many correlated mutations are in close proximity to functional sites. thus, a testbed has been created that consists of enzymes whose structure and active site residues are well characterized. to assess the quality of correlation analysis, residue positions around functional sites have been counted as positives and all others as negatives. to broaden the statistical basis, we compiled a non-redundant dataset of enzymes, whose functional sites, i.e. catalytic and binding sites, are known and which are represented by a pdb structure and a corresponding msa in the hssp database; see materials. to determine performance values, residues were classified and the distances between van der waals spheres were determined. we regarded all residues with a maximal distance of Å to a functional site as positive cases and all other residues as negative cases. the classification and the resulting performance depends on the chosen p-value and the cut-off for conz. this is why we tested several combinations and summarized results in table  for a p-value between 10- and 10- and a conz-threshold of, the specificity was between and and precision was between and. for the p-value 10- and the conz-threshold of, specificity was and precision. for p-values ≤ 10- and conz =  the performance reached a plateau. the comparison with the predecessor algorithm h2r made clear that the novel algorithm performed better: specificity and precision were up to 3% higher. additionally, we analyzed the same dataset with the algorithms cmat and sca, which predict pairs of correlated residue positions. standalone versions as of february were downloaded and applying the same criteria as above, performance was determined. cmat was used with default parameters. for sca, we selected three cut-off values,, and. performance values were added to table  cmat reached a specificity of and a precision of. for sca, the specificity increased from to, and the precision from to, for the cut-offs and. these results indicate that residue positions predicted by h2rs are more likely close to functional sites. moreover, the number of false positives is lower, as indicated by the higher precision values determined for h2rs. these numbers are a rough estimate of the algorithm’s performance due to the limitations of the in silico testbed. however, all other alternative methods of performance evaluation are not applicable here: these are the analysis of simulated msas, the determination of the residues’ spatial distance or an assessment of free energy differences derived from double mutants. in silico testbed for all programs, specificity and precision were deduced from the analysis of enzymes with known catalytic and binding sites. residues with a maximal distance of Å to a functional site were regarded as positives. all other residues were regarded as negatives. h2r and cmat were used with default settings. for h2rs, the cut-off was applied to the p-value. for sca, three cut-off values were chosen. an assessment of predicted coevolving residues in well-characterized enzymes to evaluate performance of our algorithm in more detail, we analyzed the h2rs predictions for five well studied enzymes: three enzymes from tryptophan biosynthesis, named trpa, trpb, trpc, dihydrofolate reductase, and hexokinase. trpa and trpb constitute the heteromeric tryptophan synthase complex, which catalyzes the final reaction of indole-3-glycerole phosphate and serine to tryptophan. trpa cleaves indole-3-glycerol phosphate to glyceraldehyde-3-phosphate and indole, which is transported through a hydrophobic tunnel to the active center of trpb. there, tryptophan is synthesized from serine and indole. for the localization of predicted residue positions, we utilized the 3d dataset with pdb id 1kfc, which is the trpa/trpb complex from salmonella typhimurium. the enzyme indole-3-glycerol phosphate synthase catalyzes the ring closure of an n-alkylated anthranilate to a 3-alkyl indole derivative, which is the fourth step in the tryptophan biosynthesis. it adopts the widespread 8-barrel fold and has been studied in detail. here, we utilized the structure of trpc from sulfolobus solfataricus. dhfr catalyzes the reduction of dihydrofolate to tetrahydrofolate via hydride transfer from nadph. it has been found in most organisms and plays a critical role for cell proliferation and cell growth. we utilized the structure determined for dhfr from escherichia coli. the hexokinase from schistosoma mansoni is the first enzyme in the glycolytic pathway and catalyzes the transfer of a phosphoryl group to alpha-6-glucose. the 3d crystal structure contains so anions in the catalytic cleft. smhk is the only enzyme of a larger set that has been analyzed previously by correlation analysis and for which the msa was available online. to generate smhk_cma, the authors have used a sophisticated protocol to merge several structure based msas. although local and global methods of correlation analysis have different objectives, we were interested to determine the overlap of their predictions. this is why we also compared the outcome of h2rs and psicov, which is a global method predicting residue contacts. for psicov we analyzed the top l/ predictions, which is the recommended default for a protein sequence of length l. an msa was created for each enzyme by using delta-blast with the options max target threshold and expect threshold 10- the resulting sequences were realigned by means of mafft in linsi mode. we were interested in an assessment of the most specific h2rs predictions. this is why we chose the low cut-off 10- for the p-value and a conz-threshold of. to allow for a comparison, we also listed the conz-values for all residues predicted by h2r in table  residues were regarded as functionally important, if they were close to a functional site specified in pdbsum. thus, all direct neighbors in the sequence were chosen and all residues with a 3d distance of maximally Å. for the enzymes sttrpa, sttrpb, sstrpc, ecdhfr, and smhk, h2r and h2r were used to identify important residue positions. for these residues, annotation was deduced from literature. the first column lists the name of the enzyme. the second column gives the residue and its position. the third column gives the conz-value deduced by h2rs from all uvne()-values based on a p-value of 10- the column h2r lists mean conn-values resulting from randomization tests. the column psicov lists the number of contacting pairs the residue belonged to. the last column lists the role of the residues, for details see results. “nd” indicates that we did not find clues to the function of this residue. sttrpa consists of residues, and h2rs predicted two important residues, namely l and l both residues are in close proximity to the substrate; see figure  h2r predicted l s a i and l s stabilizes the inactive conformation of the active center. a and i are near the active site and l belongs to the trpa/trpb interface. l and l also belong to the l/ predictions of psicov; see table  for sttrpb, h2rs predicted of the residues as being important; see figure  t q and v are in close proximity to the substrate binding residue k. c and f are located at the end of the hydrophobic tunnel and t and s are metal binding sites. g is important for the coordination of ion binding, and s and p are in close proximity to the bound sodium ion. m is in contact with f and s308; see above. the role of the two residues c and q is unknown to us. in contrast, h2r predicted five of these residues, namely c q p f s and additionally p g l and f f is a metal binding site, g is near the binding site for the substrate and the function of p and l is unknown to us. of the h2rs predictions, belong to the l/ contacting residues predicted by psicov; see table  for sstrpc, h2rs predicted important positions; see figure  v is near the active site. i and l are catalytically important: after replacing each of these two residues by alanine, the activity of trpc dropped 30-fold. a is located next to the substrate binding site e and the catalytic residue s; s is known to be a phosphate binding site. the role of the two residues i and i is unknown to us. h2r detected the phosphate binding site s the catalytically important residue i plus the residues a y m n and l a m and n are near the active site. the role of l is unknown however, the l187a mutation has no drastic effect on function and stability. the function of y is unknown to us. all of the residue positions predicted by h2rs belonged to the l/ contacting residue pairs predicted by psicov; see table  for ecdhfr, h2rs predicted six important residue positions; see figure  w and k are contacting the substrate, whereas a and m are in close proximity to the binding site a and the catalytic site i respectively. h plays a significant role in the network of coupled motions required for a hydride transfer and a mutation of g which lies in proximity of nadph, reduced the hydride transfer rate. the predecessor algorithm, h2r, identified g h plus d whose function is unknown to us. of the above sites, only a was an element of the l/ predictions of psicov; see table  smhk consists of a hk type- and a hk type- domain; see entry q of uniprot. h2rs identified suspicious residues, which we number according to the pdbsum entry 1bdg. a c and a are located at the very end of domain whereas c occurs at the very beginning of domain and these four residues are flanking a ß-turn. k is a neighbor of q that binds glc, v is a neighbor of i and of t. v is a neighbor of g and s that both bind so. t is close to c. only for two residues, namely t and v their role is unknown to us; however both residues are located at the domain interface at a distance of not more than Å. h2r found c and additionally d whose function is unknown to us. of the h2rs predictions were in the l/ predictions of psicov. when utilizing the msa smhk_cma, h2rs predicted only three residues with a positive conz-value, which is given in brackets: k, t, and c. t binds glc, and k is located next to the glc binding e. for c see above. interestingly, in the sequences remaining in the msa after filtering, residue positions and were occupied in not more than 43% by cysteines, which form a disulfide bridge that stiffens the orientation of the two domains. alternatively, the following residue pairs were observed with more than 2% frequency: st, gv, sm, rt, hp, av and ra. these distinct pairwise combinations support nicely the idea of mutual dependencies and pairwise correlations. although the number of cases is small, these well characterized proteins allow for a more realistic assessment of the prediction performance. altogether, h2rs predicted important residues and h2r respectively. false positives were in the case of h2rs and in the case of h2r. thus, the resulting precision is for h2rs and for h2r. these results emphasize the relatively high specificity reached by computing conn-values and additionally suggest a considerable improvement for the novel algorithm. discussion h2rs is a major improvement over h2r for all well-characterized enzymes studied in results, h2rs predicted a larger number and a higher fraction of residue positions for which we could rationalize an important role in function or stability. here, we concentrated on the analysis of residues with a conz-value ≥ . generally, this detailed analysis of five enzymes signals more precisely than the assessment of our in silico testbed the improved specificity of h2rs. it was achieved i) by replacing shannon’s entropy by the von neumann entropy and ii) by integrating a more sensitive statistical approach that adapts to the composition of each pair of msa columns. based on this dataset, we can expect a 10% increase in specificity to nearly 90%. however, this improvement has to be paid with a much longer execution time: computing the von neumann entropy requires the determination of eigenvalues, which is time-consuming and the determination of p-values further increases the execution time by a factor of one way of accelerating the calculation of entropy values might be an application of the rényi entropy, which is a generalization of the von neumann entropy. for 0 < α ≠  the α-rényi entropy is given by α-rek,l=11-αlog∑i=1400λiα and for α →  we recover the neumann entropy vne. interestingly, for α =  the calculation of the α-rényi entropy does not require the eigenvalues of the matrix ρk,l but only the diagonal entries of the square of ρk,l, which drastically speeds up the computation. however, it has not been tested yet whether the rényi entropy allows the adequate modeling of biological phenomena like residue substitutions. global and local methods of correlation analysis complement each other one goal in the design of h2r, which is a local method, was the identification of individual residue positions important for protein function or stability. this is why we introduced the conn-value. for two enzymes it has been shown that positions with high conn-values have an increased probability of being important for enzyme function or stability. the results presented here further confirm the high specificity to be gained with local methods, which is in agreement with data from the literature; see e.g. refs.. the results obtained for smhk emphasize that not all correlated mutations are due to functional constraints: of residues with high conz-values were located at the domain interface and two of them belong to a disulfide bond that interlinks the domains in some of the homologous proteins. the other residue combinations observed at these two positions illustrate nicely that they were to a great extent occupied by unique residue pairs. moreover, these findings emphasize a limitation of the in silico testbed. structurally important residues often lay far apart from the catalytic center. as shown above, some bear a strong correlation signal and are identified by h2rs. however, these hits are regarded as false positives and deteriorate the performance values deduced from the testbed. whereas local methods consider transitive correlations as well, global methods aim at eliminating these dependencies. the outcome of h2rs and the l/ predictions of the global method psicov overlapped only for of residue positions; see table  this result can be explained by the scope of the methods: according to the desired function, global methods identify contacting residue pairs which are not necessarily enriched near functional sites. using the msa smhk_cma, h2rs predicted only three residues known to be functionally important, albeit two with low conz-values. using the same dataset, the algorithm comulator, which aims at identifying perturbations, detected a network of six residue positions that surround the active site. their occupancy almost perfectly separated the two main groups of glucokinases. in summary, these findings highlight the pros and cons of the different approaches and suggest that they supplement each other quite well. msas have to be prepared carefully a critical parameter of correlation analysis is the preparation of the input, i.e. the msa. for the prediction of intra-protein residue contacts, a strong correlation between the number of homologs and the prediction strength has been shown, which further increased, if orthologs and paralogs were included in the msa. for the sake of standardization, we used in all cases studied here the same methods of msa preparation without human intervention. additionally we chose identical and very rigorous cut-offs for the identification of important residue positions. this rigid protocol might be the reason for the considerably differing number of predictions: using the cut-off conz ≥  and a p-values of 10- h2rs predicted for sttrpa only but for sttrpb important residue positions. these differences suggest for the user an individual adjustment of the parameters for each protein family in order to optimize the benefit of correlation analysis.  the various global and local methods of correlation analysis have their own field of application and supplement each other. we made plausible that residues in the vicinity of functional sites, which are a large portion of h2rs predictions, do not necessarily belong to residue pairs with the strongest global correlation signal. the predictions of global methods can be assessed by the 3d distance of the involved residue pairs. in contrast, the evaluation of local methods is more ambiguous. due to the lack of a precise definition of a correlated mutation, it is difficult to specify positive cases. this circumstance has drastic consequences and imposes restrictions to the design and the evaluation of algorithms. with this in mind, we developed an algorithm that considers pairwise substitution frequencies and assesses the strength of the correlation signal statistically. we made plausible that in silico testbeds only allow for a rough performance evaluation. favorable is the detailed analysis of well characterized model systems, which is only feasible for a small number of cases. 
comparing the outcome of the local method h2rs and of the global method psicov indicates that such methods supplement each other and have different scopes of application. we determined a by deducing substitution frequencies from contacting residue pairs observed in the homologs of proteins, whose structure is known. considering substitution frequencies for residue pairs by means of the von neumann entropy and a p-value improved the success rate in identifying important residue positions. to decide on significance, the algorithm calculates a p-value by performing a randomization test for each individual pair of residue positions. a common measure assessing these dependencies is the mutual information, which is based on shannon’s information theory that utilizes probabilities only. to compute the corresponding mutual information, a matrix a is required, which assesses the similarity of residue pairs. the identification of functionally important residue positions is an important task of computational biology. consequently, such approaches do not consider the similarity of residue pairs, which may degrade the algorithm’s performance. methods of correlation analysis allow for the identification of pairs of residue positions, whose occupancy is mutually dependent due to constraints imposed by protein structure or function. the outcome of h2rs and of a method that predicts contacting residue positions overlapped only marginally. the integration of proven statistical concepts and normalization allows for an easier comparison of results obtained with different proteins. to improve specificity of h2r, we developed a revised algorithm, named h2rs, which is based on the von neumann entropy. the analysis of a large in silico testbed demonstrated that specificity and precision were higher for h2rs than for h2r and two other methods of correlation analysis.
 the identification of functionally important residue positions is an important task of computational biology. methods of correlation analysis allow for the identification of pairs of residue positions, whose occupancy is mutually dependent due to constraints imposed by protein structure or function. a common measure assessing these dependencies is the mutual information, which is based on shannon’s information theory that utilizes probabilities only. consequently, such approaches do not consider the similarity of residue pairs, which may degrade the algorithm’s performance. one typical algorithm is h2r, which characterizes each individual residue position k by the conn-value, which is the number of significantly correlated pairs it belongs to.  to improve specificity of h2r, we developed a revised algorithm, named h2rs, which is based on the von neumann entropy. to compute the corresponding mutual information, a matrix a is required, which assesses the similarity of residue pairs. we determined a by deducing substitution frequencies from contacting residue pairs observed in the homologs of proteins, whose structure is known. in analogy to h2r, the enhanced algorithm computes a normalized conn-value. within the framework of h2rs, only statistically significant vne values were considered. to decide on significance, the algorithm calculates a p-value by performing a randomization test for each individual pair of residue positions. the analysis of a large in silico testbed demonstrated that specificity and precision were higher for h2rs than for h2r and two other methods of correlation analysis. the gain in prediction quality is further confirmed by a detailed assessment of five well-studied enzymes. the outcome of h2rs and of a method that predicts contacting residue positions overlapped only marginally. h2rs can be downloaded from http://www-bioinf.uni-regensburg.de.  considering substitution frequencies for residue pairs by means of the von neumann entropy and a p-value improved the success rate in identifying important residue positions. the integration of proven statistical concepts and normalization allows for an easier comparison of results obtained with different proteins. comparing the outcome of the local method h2rs and of the global method psicov indicates that such methods supplement each other and have different scopes of application. 