BACKGROUND
researchers carrying out structure-based drug design  are constantly looking to improve the modelling of protein conformational change and its relationship to ligand binding.
it is well known that protein-target conformational flexibility can lead to problems in, for instance, small molecule binding-mode prediction and structure-activity relationship interpretation
<cit> .
it has been said that a lack of appreciation of the dynamics of macromolecular complexation is holding back progress in virtual screening
<cit> .
changes in protein conformation, when experimentally observed, can lead to the discovery of highly prized cryptic binding sites
<cit>  and allosteric pockets
<cit> .
for the discovery of protein-protein interaction inhibitors, where there is no small endogenous small molecule to mimic, the treatment proteins as flexible entities is especially important
<cit> .
nmr is perhaps the experimental technique most able to report on protein structure and dynamics in solution.
however, multiple crystal structures of the same protein can provide much information on the conformational alternatives adopted by protein structures when interacting with other molecules
<cit> .
x-ray crystal structures are the basis of the majority of sbdd projects, and drug companies often amass hundreds of crystal structures of the same protein with different ligands bound over the course of a drug discovery project.
it has been estimated that the pharmaceutical industry as a whole currently solves  <dig>  macromolecular crystal structures each year
<cit> .
normally only one set of model coordinates is provided per structure solution but many similar models can be found that fit the experimental data equally well
<cit> .
single sets of structure coordinates can also be used to generate ensembles using conformational sampling techniques
<cit>  or simulations
<cit> .
taken together, experimental and computer generated ensembles should provide a fuller picture of a protein molecule’s true nature
<cit> .
many authors have encouraged the use of protein-structure ensembles in drug discovery
.
here we would like to distinguish this approach from more traditional  sbdd by the coining the expression “ensemble-based drug discovery” .
it is still not straightforward with existing molecular modelling and bioinformatics tools to make full use of the large and ever increasing number of structures in some drug target families in the pdb
<cit>  and in proprietary collections.
software for the analysis of structural ensembles is readily available of course.
root mean squared deviation  of equivalent atoms, after optimal superposition, is the most widely used measure of pairwise protein structure similarity.
similarly root mean squared fluctuation  is used to measure residue positional variability.
principal component analysis , or essential dynamics, is a very effective way of distilling the most important motions from molecular dynamics  trajectories
<cit> , and structures derived from x-ray crystallography
<cit>  and nmr
<cit> .
many programs provide functionality for their calculation, for example gromacs analysis modules
<cit> , bio3d
<cit>  and dynamite
<cit> .
wordom
<cit>  is a package for the analysis of md simulations which also provides pca, along with many other analysis techniques.
many software packages designed for the analysis of md snapshots require a consistent set of atoms and residues, making their use on crystal structures awkward.
comparisons of the essential dynamics produced by md with principal components derived from crystal and nmr structures have been used to validate the results of the former approach
<cit> .
prody
<cit>  facilitates the comparison of crystal structure and md trajectory pc’s with normal modes calculated from a single structure.
other programs allow pairwise comparison of protein structures in order to identify hinge regions and interdomain motions.
these include molmovdb
<cit> , dyndom
<cit> , flexprot
<cit>  and fatcat
<cit> .
dihedral angle pca
<cit>  is a less common approach but can be used to identity hinge regions in proteins
<cit> .
another technique employed is distance difference matrices, for example using struster
<cit> .
as well as the coordinates of a protein, one can also study the pockets formed by its surface and how they change as the protein conformation changes.
these pockets are of particular interest to those involved in drug discovery.
mdpocket is designed for the analysis of pockets in md simulations using fpocket
<cit> .
provar
<cit>  uses calculated surface properties, generated by a range of programs, assigned to individual residues allowing comparison between the results and within ensembles of structures without superposition.
the techniques developed here are designed to complement the existing approaches described above, especially those employing rmsd and coordinate pca that rely upon molecular superposition.
molecular superposition becomes problematic when large conformational changes occur in a protein, involving multiple domains and rmsd can violate the triangle inequality rule
<cit> .
in such cases it is useful to compare the curvature and torsion of a spline fitted to cα atoms, an approach that was first used for the comparison of protein backbone conformations in 1978
<cit> .
although the uses of this and related approaches  remain relatively uncommon, curvature and torsion of cα splines have been utilised to find the conserved core of homologous proteins
<cit> , to analyse secondary structure motifs
<cit>  and to assess structural alignments of distantly related proteins
<cit> .
this sort of analysis is very efficient and scales well with the size and number of structures in an ensemble.
in the work described here, cα-spline curvature and torsion, together with side-chain conformation, intermolecular interaction fingerprints and pocket properties, are represented as per-residue descriptors and grouped by alignment position in a way that is analogous to correlated mutation analysis sequence analysis techniques
<cit> .
summary statistics and inter-residue relationships calculated in this way are mapped onto representative 3d structures to aid visual interpretation of the results.
all the methodologies are programmed in a purpose built python package called polyphony, in which the “plug-in” architecture allows other descriptors of protein structure generated by 3rd party programs, for example fpocket
<cit> , ncont
<cit> , credo
<cit>  and piccolo
<cit> , to be added with very little effort.
matplotlib
<cit> , pymol
<cit> , jalview
<cit> , and ete
<cit>  are used for visualisation of the results.
below various novel metrics and algorithms are described for the statistical analyses of local geometry of proteins, compared across ensembles of structures.
these include a tailor-made variance measure that is used to distinguish random thermal-like fluctuations from significant conformational changes.
in addition, existing techniques such as pca are applied in new ways to differential geometry descriptors of protein conformation.
a novel approach to using the output from fpocket
<cit>  to identify distinct and cryptic pockets from an ensemble of structures is described.
due to the scalability and automation of the general approach taken it can be applied to large sets of experimental structures or calculated conformations.
these advantages are illustrated below by a comparison of conformational commonalities within evolutionary related proteins and by comparing snapshots from an md simulation with x-ray crystal structures of the same protein.
the ultimate aim of this work is to extend the repertoire of tools available to those wanting to do ensemble-based drug discovery.
however, it is hoped that the methods will also be used to discover fundamental mechanisms in the makeup of protein machines.
all the known mechanistic changes in cdk <dig> and related structures were discovered afresh but in addition, more subtle, evolutionarily conserved allosteric changes were revealed.
RESULTS
the initial emphasis of polyphony is on the analysis of structures from x-ray crystallography but it can be used on nmr structures and snapshots along the trajectories of protein simulations.
cyclin-dependent kinase  <dig>  is often used to validate computational approaches which treat protein-ligand complexes as flexible entities  due to the large number of crystal structures deposited in the pdb
<cit>  and to the conformational changes that can be observed in these structures.
it is a protein kinase involved in cell cycle progression and inhibitors of its function have been designed with the purpose of exploring their utility as anticancer agents
<cit> .
kinase activity is switched on by binding of cyclin and phosphorylation by cdk activating kinase
<cit> .
cdk <dig> is inhibited by the binding of cip and ink <dig> family proteins
<cit> .
the binding of cyclin a or cyclin e results in conformational change in which the pstaire region  rotates and swings in so that the glu  <dig> takes its place as part of the catalytic triad
<cit> .
the n and c terminal lobes also move closer together, pivoting about the hinge region, and the t  loop moves to expose the catalytic cleft.
phosphorylation of thr  <dig> further organises the substrate binding site by causing additional changes to the t-loop
<cit> .
the glycine rich loop is very flexible and shows structural heterogeneity in cdk <dig> structures as well as in protein kinases in general.
the cmgc insert  is a flexible region in the c-terminal lobe.
see figure
<dig> for the locations of the substructures mentioned above.figure  <dig>
the structure of the cdk <dig> - cyclin complex .
this structure is phosphorylated at thr  <dig>  and bound to cyclin a <dig>  and adp .
the pstaire region is shown in orange, the t-loop in yellow, the glycine rich loop in magenta, the cgmc insert in purple, the αef-αf loop in cyan and the hinge region in red.
analysis of cdk <dig> crystal structures
a polyphony script was used to download the structures of the 95% identity sequence cluster containing pdb code 1hck, chain a from the rcsb website
<cit> .
at the time of analysis this cluster contained  <dig> structures and  <dig> chains.
the vast majority were solved with the p <dig> space group and over half had unit cell dimensions of around  <dig>   <dig>  72 Å. however, there were examples of  <dig> other space groups.
starting models with pdb codes were cited for  <dig> of the structures.
of these,  <dig> started directly or indirectly from 1hck or 1hcl
<cit> ,  <dig> from 1fin
<cit> ,  <dig> from 1qmz
<cit>  and  <dig> from 1b39
<cit> .
there were no other significant groupings of structures reported to start from the same model.
within grouping  average rmsds to their respective starting model range from  <dig>  to  <dig>  Å .
overall  average rmsds to these starting structures range from  <dig>  to  <dig>  Å .
the ramachandran plot redrawn
a ramachandran plot for all  <dig> chains combined  was plotted and separate secondary structure regions colour coded.
figure
2b shows the equivalent plot of curvature against torsion with colours derived from the ramachandran plot.
it can be seen that in this latter plot the residues with the same secondary structure are roughly in the same area of the graph, although these areas are less distinct.figure  <dig>
the ramachandran plot compared to the plot of residue curvature against torsion.
ramachandran plot for all residues in all chains.
colours delineate the regions specified in
<cit>  plot of curvature against torsion for each residue in all chains.
colours denote the region the ramachandran plot shown in .
contour levels are at fibonacci numbers due to the large number of overlapping points.
local conformational variability
figure
3a shows how conformational variability for each cdk <dig> residue calculated from spline κ/τ compares to rmsf.
the differences between the two measures increase as variability increases.
rmsf is calculated after fitting a sliding window of  <dig> contiguous residues to equivalents in a reference structure.
this fit becomes arbitrary at high degrees of conformational divergence.
figure
3b shows the conformational variability for each residue calculated from φ/ψ torsion angles and spline κ/τ.
variation is calculated in different ways for the two descriptors  but it is perhaps surprising that they show no correlation.
this difference is probably due to the fact that φ/ψ torsion is a more local measure of conformation than κ/τ, which is calculated for a spline whose shape is dependent on the conformation of neighbouring residues.
changes in φ/ψ torsion angles could be neutralised by compensatory changes in the conformation of neighbouring residues.figure  <dig>
cα-trace κ/τ variation compared to two other metrics of backbone variation  root mean squared fluctuation  in cα position after least squares fitting a sliding window of  <dig> residues to a reference structure and  variation in φ/ψ torsion angle .
labels show alignment position numbers indicating that the largest differences lie largely in the t-loop .
subgroup comparisons
in order to divide proteins into monomeric and  cyclin bound forms, the chains were clustered using the tanimoto similarity of the per residue protein-protein interaction fingerprint
<cit>  of contacts extracted from the piccolo database
<cit> .
figures
4a and b below show curvature and torsion values respectively for residues in the region around the conserved dfg motif which lies at the start of the t-loop.
the plots are colour coded according to the presence  or absence  of a biologically relevant protein-protein interaction.
they clearly show the difference in conformation in this region that occurs upon protein binding as shown in 3d in figure
<dig> figure  <dig>
plot of curvature and torsion of cdk <dig> chains as a function of sequence for residues n <dig> to a <dig> in the region of the dfg loop.
monomeric chains are coloured blue and protein-bound chains, red.
curvature and  torsion as a function of sequence.
pymol  cα trace for residues n <dig> to a <dig> in the region of the dfg loop.
monomeric chains are coloured blue and protein-bound chains, red.
note that κ/τ values are not calculated for residues within  <dig> residues of the termini and the frayed ends adjacent to disordered residues .
conformational variability and average temperature factor
of course loading all  <dig> chains into a molecular graphics program, aligning them, colouring them etc.
is slow, memory intensive, and the results can be messy.
instead, summary views can be generated with polyphony, where properties of the whole ensemble are projected onto one or more representative structures in pymol.
figure
<dig> shows two such summary views.
in figure
6a the per residue relative variability in backbone and side-chain conformation are shown.
this can be compared to the relative average normalised cα temperature factor in figure
6b.
temperature factor normalisation is done using the per protein mean and standard deviation, calculated after the removal of outlier residues, following the procedure of smith et al.
<cit> .
relative backbone conformational variability and average temperature factor are grossly similar in distribution along the sequence.
however the former is a more local property showing sharper variations.
residues that stand out as having particularly variable backbone conformations are gly  <dig>  the last glycine in the glycine-rich loop, phe  <dig> of the dfg motif and val  <dig> near the end of the t-loop.
the side-chains of tyrosines  <dig> and  <dig> also stand out as conformationally variable.
hinge residues, in which a small conformational change leads to a large shift in a subdomain, are not easily spotted in this analysis.figure  <dig>
conformational variability vs. average temperature factor.
pymol putty cartoon of 3pxz
<cit>  chain a showing summaries over  <dig> cdk <dig> chains of  relative variability in conformation of backbone κ/τ  and sidechains : a grey colour indicates disordered residues in greater than 50% of structures, and  average cα temperature factor .
a red colour and thick tube indicates the highest variability and highest average temperature factor.
the labels highlight particularly conformationally variable residues.
the influence of crystal contacts
viewed in a 2d plot , it is clear that, in the c-terminal half of the protein, after the t-loop, the conformational variability and average temperature factors are higher in the protein bound structures of cdk <dig>  as expected, parts of the t-loop and the loop immediately preceding the pstaire helix are more ordered in the bound form.
enticingly, this presents a picture of a loss of entropy at the protein-protein interface in the t-loop upon binding, which is compensated for by a gain in entropy elsewhere.
however, unnatural crystal contacts are also protein-protein interactions and must be taken into consideration.
protein-protein interactions can result in lower temperature factors at the interface
<cit> .
for instance, residues 20– <dig> have lower average temperature factors in the bound form but probably as a result of conserved crystal contact with another cyclin molecule in the asymmetric unit.
barrett and nobel
<cit> , using molecular dynamics simulations in solution, found that cyclin bound cdk <dig> becomes more rigid on phosphorylation, except for the cmgc insert region which becomes more mobile.figure  <dig>
comparison of conformational variability, average temperature factor and presence of crystal contacts in cdk <dig>  alignment position is labelled on the x axis.
monomeric chains are coloured blue and protein bound chains, red.
conformational variability using curvature and torsion  average normalised temperature factor  crystal contacts as calculated using ccp4
<cit>  ncont.
the value plotted is the average proportion of structures with at least one contact within 5 Å. the location of the pstaire helix is highlighted in orange and the t-loop is highlighted in yellow.
hinges and loop flips
once groups or clusters of structures are defined, they can be compared in terms of conformation and also interactions.
shown in figure
<dig> are the variances  in curvature and torsion between bound and unbound cdk <dig> chains.
interestingly, variances in the curvature and/or torsion at individual residues are sharply defined and not smoothed out due to spline fitting.
it is clear from figure
<dig>  and from looking at a 3d representation of the data using the polyphony pymol api, that hinges are highlighted by this analysis.
in particular there is a big change between monomeric and protein bound structure in τ at residue gln  <dig> of the hinge region, the κ and τ at phe  <dig> of the dfg motif, and a definite signal at asn  <dig> at the end of the pstaire helix.
in addition many residues light up in the t-loop and the cyclin binding region.
this ability to highlight what might be called micro-hinges, which include those between secondary structure elements, as well as between domains, is an advantage of this approach.
once noise, or random thermal motion, is reduced in this way it becomes apparent that some changes in backbone conformation are compensated for by opposite changes at a neighbouring residue.
these paired peaks are the signals produced by loop regions that occupy different but conserved conformations in the groups of proteins .
examples that can be seen in figure
<dig> are the loops centred at glu  <dig> and thr  <dig>  the cα trace emerges from the conserved core and bends one way but must then bend back again to return to the conserved core.
this sort of conformational change is distinct from a hinge because it only affects local structure.
it is worth noting here that, because curvature is always positive , a curve to the left , which is then compensated by a curve to the right , is not apparent from these plots.
curvature and torsion must be taken together to observe such differences.figure  <dig>
plot of alignment position  against the group variance
s
i
for monomeric chains  and protein bound  cdk <dig>  curvature
s
i
and  torsion
s
i
positions of the highlighted residues in the structure of cdk <dig>
principal components analysis
in order to facilitate a principal components analysis , a full matrix of data was generated .
the pca used the non-linear iterative partial least squares  algorithm implemented in pychem
<cit> .
the resulting scores plots for backbone and side-chain conformation are shown in figures
9a and b respectively.
it can be seen that the first principal component  divides monomeric and protein-bound structures.
in the backbone conformation pca  indicates a high contribution from the 72– <dig> ten loop tip but also cys  <dig>  gln  <dig> of the hinge region and asn  <dig> at the end of the pstaire helix.
in the side-chain analysis , the conformation arg  <dig> is the major discriminating factor.
it forms a salt bridge with gln  <dig> on cyclin binding, seemingly helping to lock the pstaire helix in place .
the charged residues in this salt bridge are not conserved amongst cdk isoforms.figure  <dig>
pca scores plot for cdk <dig>  monomeric chains are coloured blue and protein bound chains red.
backbone conformation using curvature and torsion.
the variance explained by pc <dig> is 45% and 51% by pc <dig> and pc <dig> combined.
side-chain conformation using the position of a terminal atom in cartesian coordinates after transformation of the residue to a reference origin .
pc <dig> explains 24% and pc <dig> plus pc <dig> explains 35% of the variance.
labels show the pdb codes of outliers.
pca discriminatory side-chain conformations.
salt bridge formation between arg  <dig> and glu  <dig> in cyclin bound structures .
tyr  <dig> in monomeric , protein-bound but unphosphorylated , protein-bound and phosphorylated at thr  <dig> .
the rest of the t-loop is cut away for clarity.
a structure that stands out as a red dot amongst the blues in figure
9a, i.e.
protein bound structure in amongst the monomers is 1buh
<cit> .
this structure is a complex of cdk <dig>  not with cyclin, but with cell cycle-regulatory protein ckshs <dig>  which binds at a different site on the c-terminal lobe of cdk <dig>  the chain b of 1fq1
<cit>  also stands out in this way in the side-chain conformation pca plot .
this is the structure of kinase-associated phosphatase  in complex with phospho-cdk <dig> .
kap binds to a different region on the surface of cdk <dig> than cyclin, but one that overlaps with that of ckshs <dig>  another structure that stands out is 1w98
<cit> , which is a structure of p-cdk <dig> in complex with a truncated cyclin e <dig>  structures 3pxz, 3pxq, 3py <dig>  and 3pxf
<cit>  lie in the centre of both the backbone and side-chain pca plots.
these are monomeric structures that are located closer to the region occupied by the protein-bound structures.
interestingly, they come from a series of structures in complex with allosteric inhibitors that cause conformational changes to the pstaire helix.
pc <dig> in figure
9a and b both subdivide the protein-bound cdk <dig> structures into two groups.
the 179– <dig> yys motif, which forms part of the αef-αf loop , contributes most to backbone conformation pc <dig> in the figure
9a.
this loop changes shape on phosphorylation of thr  <dig> and is known to be coupled to the activation loops in a variety of kinases
<cit> .
this loop lies underneath the t-loop, which is excluded from the pca because so many structures  are disordered in this region.
the side-chain conformation of tyr  <dig>  in the αef-αf loop, more or less single-handedly separates structures on the pc <dig> axis in figure
9b.
figure
10b shows its different conformations in monomeric, cyclin bound and cyclin bound plus phosphorylated cdk <dig>
the first two principal components explain a relatively low amount of the variance, especially for side-chain conformations .
unlike the usual cartesian coordinate pca, we are measuring local changes that act in concert rather than big blocks of secondary structure that move together.
many local conformational states are conserved or change in an uncorrelated manner and therefore do not strongly influence the first principal components.
correlated conformational changes
in addition to pca, there is the capability within polyphony of identifying correlated residue conformational differences between individual pairs of residues.
this can be done for backbone or side-chain conformations.
the κ/τ and x, y, z correlation coefficients, respectively, are averaged for each residue pair.
lines are drawn between the most highly correlated residues in an analogous way to young et al.
<cit>  and the covariance web produced by dynamite
<cit> .
these published methods use cartesian coordinates after structure superposition, whilst here the correlations are calculated for local conformations.
only correlations between residues with a separation of  <dig> residues in the protein sequence are shown in figure
<dig> but this cut-off is under user control in polyphony.figure  <dig>
correlated conformational change.
the  <dig> most highly correlated pairs of residues separated in sequence by at least  <dig> amino acids.
backbone.
the correlation coefficients range from  <dig>  between asn  <dig> at the end of the pstaire helix and phe  <dig> of the dfg motif, to  <dig>   side-chains.
leu  <dig> and leu  <dig> in the core of cdk <dig> switch rotamers on protein binding.
most of the cdk <dig> main chain correlations found in this way are consistent with the results of the pca and known mechanisms of action but also highlight previously unpublished interrelationships between residues.
for instance, phe  <dig> of the dfg motif is correlated with asn  <dig>  identified above as a hinge residue at the end of the pstaire region .
cys  <dig> in the αef-αf loop and gln  <dig>  which is identified above as the pivotal residue in the hinge region, are also linked to these two residues via a network of correlations.
surprisingly leu  <dig> is also paired with these residues, albeit with a lower correlation coefficient  but this probably due to crystal contacts with this residue which are much more common in monomeric compared to protein-bound structures .
similarly the side-chain conformations of his  <dig> in the c-lobe are correlated with asp  <dig> in the n-lobe .
this is almost certainly an experimental artefact, due to crystal contacts .
more interestingly, there seems to be a repacking of part of the hydrophobic core on cyclin binding.
the neighbouring residues leu  <dig> and leu  <dig> change conformation in a correlated way  .
pocket analysis
the pocket detection program fpocket
<cit>  was chosen for integration with polyphony because it is an open source program and it provides a druggability score for each pocket.
however, it is possible to integrate other programs because of polyphony’s plug-in architecture.
as employed by provar
<cit> , the assignment of pocket attributes to residues facilitates the comparison of structures of the same or homologous proteins, given a sequence alignment.
like provar, the fraction of structures in which a particular residue  is found to be part of a pocket are reported, here expressed as a percentage.
again like provar, these percentages are mapped onto the coordinates of a representative structure as colour-coded surfaces.
in polyphony there is a novel methodology which shows the user distinct  predicted druggable pockets in the structure in which they occur .
in this way cryptic pockets can easily be identified and visualised.
figure
<dig> shows the results of this procedure applied to the cdk <dig> structures.
the highest ranked pocket  was an atp binding pocket which is connected to a water filled cavity behind the pstaire helix in 2c5y
<cit> .
this cavity would not exist when the pstaire helix moves in on binding to cyclin.
thus this cavity is a potential binding site for allosteric inhibitors of cyclin binding.
in fact, fragments that bind in this cavity were recently discovered in later published structures
<cit> .
the pocket selected second , which lies under the cmgc insert also shows potential as a binding site for allosteric modulators.
the region  is known to be flexible from simulations and analysis of crystallographic temperature factors
<cit> .
it forms part of the binding site for cell cycle–regulatory protein ckshs1
<cit>  and kinase associated phosphatase
<cit> .
there are no small molecule ligands bound to this pocket in the cdk <dig> crystal structures but there are in p38α
<cit>  and jnk-1
<cit> .
the third pocket identified  is close by in the c-lobe.figure  <dig>
location of fpocket-predicted druggable binding sites selected by the polyphony distinct pocket selection procedure .
the blue pocket was ranked highest and is the atp binding site in 2c5y
<cit>  chain a. the second pocket is coloured red and is from 1oi9
<cit>  chain c. the third in green is from 2r3l
<cit>  chain a.  the blue pocket from  viewed from above, overlaid with a recently discovered allosteric ligand from 3pxz
<cit> .
the ligand from 2c5y is shown in blue sticks.
intermolecular interactions
atomic scale interactions are extracted from three different sources 1/ the piccolo database
<cit>  of predicted biologically relevant protein-protein interactions 2/ the credo database
<cit>  of protein - small molecule interactions and 3/ crystal contacts calculated using the ccp4
<cit>  program ncont.
this program is not open source but is available without cost to academic and non-profit institutions.
the two data sources are mysql databases created in-house.
piccolo is currently being updated monthly and can be downloaded and installed locally for free.
the original credo was frozen in april  <dig> and replacement by a new postgresql database
<cit> .
residue-based structure interaction fingerprints
<cit>  where extracted from piccolo and credo and analysed in various ways including clustering  using tanimoto similarity.
two papers have explored the differences in inhibitor binding in cdk <dig> atp binding site between active and inactive forms, by comparing two
<cit>  and twelve
<cit>  structures.
a comparison of credo derived sifts for the  <dig> protein-bound and monomeric cdk <dig> chains with ligands bound, revealed that the set of residues interacting with small molecule ligands were almost identical, except for glu  <dig>  this residue is part of the pstaire helix and the conserved catalytic triad and swings into the atp binding site on activation
<cit> .
none of the  <dig> ligands bound to monomeric cdk <dig> contact this residue , whereas  <dig> of the  <dig>  bound to protein complexes do.
sequence alignment view
sequence alignment viewers are, of course, a very useful way of visualising similarities and differences in the amino acid sequences of related proteins.
in many cases the simple single letter amino acid codes are augmented with colours to indicate a particular property of each residue such as hydrophobicity or known features such as the occurrence of a post translation modifications.
this idea has been extended to include descriptions of the residue conformation or environment with the folded protein
<cit> .
here we are concerned, in the main, with multiple structures of the same protein.
thus, a plain sequence alignment is not very informative.
however once annotations, such as colour coded b-factors and protein-protein contact counts, are added, the sequence alignment view becomes a convenient way to visualise the similarities and in conformation and interactions across a large number of structures.
jalview
<cit>  was chosen for this purpose because users can import colour-coded residue features from a text file.
it also allows the import of newick format trees, which can then be used to sort sequences.
polyphony can create jalview feature and newick tree files for any of the residue-based properties.
figure
<dig> shows what this looks like.figure  <dig>
a sequence alignment sorted and annotated by structural and conformational features.
dendrogram and  sequence alignment view of cdk <dig> structures shown in jalview
<cit> .
the alignment is ordered by the tree, generated by hierarchical clustering using backbone curvature-and-torsion pearson similarity.
the red vertical line in the dendrogram cuts the chains into  <dig> large clusters, which correspond to monomeric , cyclin bound  and cyclin bound and phosphorylated .
amino acids in the alignment are colour-coded by b-factor  and protein-protein contacts from piccolo
<cit>  .
white gaps indicate disordered residues.
application to the analysis of homologous proteins
even when only a handful of structures of a given protein are available they can be usefully compared with those of homologous proteins.
differences between the structures of homologous proteins include architectural as well as conformational changes.
to separate out these two sources of variation, in polyphony, conformational variance is calculated separately by protein.
then the locations of intra-protein variance are compared across homologous proteins.
below this type of analysis is illustrated on members of the cmgc family of protein kinases
<cit> .
a structural alignment was taken from the homstrad database
<cit> .
figure
<dig> shows the aligned variance plots for cdk <dig>  erk <dig>  p38α and jnk <dig>  the regions of known conformational change are highlighted.
the glycine rich loop has high values in all proteins.
a peak is visible in the hinge region, particularly in cdk <dig> and p38α.
the dfg motif has a peak in all proteins except erk <dig>  the activation loop has a high peak in the cdks and erk <dig>  there also peaks in the consensus plot in lesser-known regions, particularly in the αef-αf loop region.figure  <dig>
plot of the group variance
s
for cgmc protein kinases.
groups were clusters of very closely related structures .
cdk <dig>  n =  <dig> in  <dig> groups  cdk <dig>   <dig> in  <dig>  erk <dig>   <dig> in  <dig>  p38α,  <dig> in  <dig>  jnk <dig>   <dig> in  <dig>  consensus .
the colours indicate the glycine rich loop , the hinge region , the dfg motif , the t-loop  and αef-αf loop .
tyr  <dig>  in cdk <dig> interacts with the phosphorylated thr  <dig> and has distinct conformations in monomeric, cyclin-bound and cyclin-bound plus phosphorylated forms .
in p38α the αef-αf peak occurs at met  <dig>  in contrast to cdk <dig> tyr  <dig>  this residue is not in close proximity to the phosphorylated residues thr  <dig> and tyr  <dig>  however, its backbone conformation is strongly correlated with residues 180– <dig> .
it takes two main conformations, one where the side-chain is surface exposed and one where it is buried under the cmgc insert.
the insert itself shifts as result and is  <dig> or 3 Å closer in to the rest of the protein in the latter conformation.
as mentioned in the pocket analysis section above, this region has recently become the focus of inhibitor discovery and design in jnk- <dig> and p38α
<cit> .
in erk <dig> it was recently discovered that the cgmc insert, which is a helix-turn-helix motif, is involved with nuclear localisation
<cit>  and even dna binding
<cit> .
our analysis reveals what appears to be an evolutionary conserved link between phosphorylation and the conformation and dynamics of this functionally important part of these proteins.
use on computed conformations
one of the aims of this project was to create new ways to compare experimental structural ensembles with those generated by computational methods.
as shown in figure
<dig>  variability of backbone curvature and torsion provides a more local measure of conformational change than average b-factor.
figure
<dig> shows variability over  <dig> p38α x-ray chains and a md simulation downloaded from the model database
<cit> .
the simulation was of the structure 1a9u
<cit>  for 12 ns.
for this analysis  <dig> snapshots separated by 50 ps, covering the last 9 ns of the simulation, were extracted as a multimodel pdb format file.
many of the conformational changes seen amongst the crystal structures are likely to occur on a much longer timescale than that covered by the simulation.
however there is some qualitative similarity in the variability observed in these structures and those generated by the simulation.
for instance, the change in conformation at the key hinge residue met  <dig> is reproduced, indicating that this is a due to a high frequency motion.
there is also clearly some dynamics in the αef-αf loop centred on met  <dig>  which is highlighted in the previous section as being a region of significant conformation heterogeneity.
this residue remains buried under the cgmc insert during this short simulation but the αef-αf loop appears to stretch out as the cgmc moves away from the c-terminal lobe as a whole.
there are also discrepancies between the crystal structure conformational difference and the simulated dynamics.
for instance, in the region 240– <dig> there is partial unravelling of the first helix of the cgmc insert in the simulation, whereas it is well conserved between crystal structures.
this could be due to the influence of crystal contacts or other reasons.
work is under way to further assess the use of polyphony for the analysis of md trajectories.figure  <dig>
backbone conformational variability for x-ray structures  and an md simulation of p <dig> .
CONCLUSIONS
the suite of new methods for the analysis of ensembles of protein structures is selected on the premise that there are advantages to comparing structures without superimposing them.
the approach taken is analogous to sequence analysis techniques and starts from a sequence alignment.
the alignment itself is trivial because the sequences involved are almost identical.
the amino acid types in the columns in the alignment do not differ  but the conformations of the residues do.
backbone conformational changes can be those that have very little effect on the shape of the protein, such as those near the termini or at the end of flexible loops.
others are hinging motions that change the relationship between domains.
one difference between these two extremes of variation is that the former tend to be randomly distributed and the latter tend to be consistent across multiple structures.
another difference is that the latter tend to coincide with environmental changes, such as a binding event.
when a binding event is a protein-protein interaction with a large interface, correlated conformational changes tend to occur at the interface and allosterically.
all of these effects can be detected by statistical analysis of all the local conformational differences within an ensemble without the need to superimpose the structures and calculate euclidean distances in cartesian space.
once discovered in this way, the consequences of these significant conformational changes can be observed by visual comparison of the structures and the use of existing analytical techniques.
the obvious drawback with this approach is the requirement for a large sample size.
one must also be wary of systematic experimental artefacts such a conserved non-biological crystal contacts.
these issues can be addressed by complementing x-ray structures with nmr derived and computationally generated ensembles.
however, there is great benefit not only in solving the crystal structures of new proteins but also in solving the structure of same protein multiple times, especially when co-crystallised with new binding partners.
tools and methods for superposition-independent statistical analysis of protein structure ensembles were developed.
the general approach, and the individual methodologies within it, were validated by the rediscovery of the published findings of the many authors who have compared crystal structures of cdk <dig> since jeffrey et al.
in 1995
<cit> .
the major conformational changes that occur on cyclin binding, such as the movement of the pstaire helix and the opening and closing of the gap between the n and c-terminal lobes, were detected via symptomatic changes in hinge residues.
in addition, more subtle changes were also detected and found to be conserved in other cmgc kinases.
these include correlated changes in the αef-αf loop linking phosphorylation sites on the activation loop to the cmgc insert.
this information provides a further clue to the role of this region whose importance is only recently beginning to be revealed and targeted with small molecule ligands in jnk <dig> and p38α.
thus far no structures of cdk <dig> with small molecules bound in this region have been published.
the pocket analysis above reveals a pocket can exist there and it is predicted to be druggable, illustrating the utility of ensemble-based drug discovery.

 researchers carrying out structure-based drug design are constantly looking to improve the modelling of protein conformational change and its relationship to ligand binding. it is well known that protein-target conformational flexibility can lead to problems in, for instance, small molecule binding-mode prediction and structure-activity relationship interpretation. it has been said that a lack of appreciation of the dynamics of macromolecular complexation is holding back progress in virtual screening. changes in protein conformation, when experimentally observed, can lead to the discovery of highly prized cryptic binding sites and allosteric pockets. for the discovery of protein-protein interaction inhibitors, where there is no small endogenous small molecule to mimic, the treatment proteins as flexible entities is especially important. nmr is perhaps the experimental technique most able to report on protein structure and dynamics in solution. however, multiple crystal structures of the same protein can provide much information on the conformational alternatives adopted by protein structures when interacting with other molecules. x-ray crystal structures are the basis of the majority of sbdd projects, and drug companies often amass hundreds of crystal structures of the same protein with different ligands bound over the course of a drug discovery project. it has been estimated that the pharmaceutical industry as a whole currently solves macromolecular crystal structures each year. normally only one set of model coordinates is provided per structure solution but many similar models can be found that fit the experimental data equally well. single sets of structure coordinates can also be used to generate ensembles using conformational sampling techniques or simulations. taken together, experimental and computer generated ensembles should provide a fuller picture of a protein molecule’s true nature. many authors have encouraged the use of protein-structure ensembles in drug discovery. here we would like to distinguish this approach from more traditional sbdd by the coining the expression “ensemble-based drug discovery”. it is still not straightforward with existing molecular modelling and bioinformatics tools to make full use of the large and ever increasing number of structures in some drug target families in the pdb and in proprietary collections. software for the analysis of structural ensembles is readily available of course. root mean squared deviation of equivalent atoms, after optimal superposition, is the most widely used measure of pairwise protein structure similarity. similarly root mean squared fluctuation is used to measure residue positional variability. principal component analysis, or essential dynamics, is a very effective way of distilling the most important motions from molecular dynamics trajectories, and structures derived from x-ray crystallography and nmr. many programs provide functionality for their calculation, for example gromacs analysis modules, bio3d and dynamite. wordom is a package for the analysis of md simulations which also provides pca, along with many other analysis techniques. many software packages designed for the analysis of md snapshots require a consistent set of atoms and residues, making their use on crystal structures awkward. comparisons of the essential dynamics produced by md with principal components derived from crystal and nmr structures have been used to validate the results of the former approach. prody facilitates the comparison of crystal structure and md trajectory pc’s with normal modes calculated from a single structure. other programs allow pairwise comparison of protein structures in order to identify hinge regions and interdomain motions. these include molmovdb, dyndom, flexprot and fatcat. dihedral angle pca is a less common approach but can be used to identity hinge regions in proteins. another technique employed is distance difference matrices, for example using struster. as well as the coordinates of a protein, one can also study the pockets formed by its surface and how they change as the protein conformation changes. these pockets are of particular interest to those involved in drug discovery. mdpocket is designed for the analysis of pockets in md simulations using fpocket. provar uses calculated surface properties, generated by a range of programs, assigned to individual residues allowing comparison between the results and within ensembles of structures without superposition. the techniques developed here are designed to complement the existing approaches described above, especially those employing rmsd and coordinate pca that rely upon molecular superposition. molecular superposition becomes problematic when large conformational changes occur in a protein, involving multiple domains and rmsd can violate the triangle inequality rule. in such cases it is useful to compare the curvature and torsion of a spline fitted to cα atoms, an approach that was first used for the comparison of protein backbone conformations in 1978. although the uses of this and related approaches remain relatively uncommon, curvature and torsion of cα splines have been utilised to find the conserved core of homologous proteins, to analyse secondary structure motifs and to assess structural alignments of distantly related proteins. this sort of analysis is very efficient and scales well with the size and number of structures in an ensemble. in the work described here, cα-spline curvature and torsion, together with side-chain conformation, intermolecular interaction fingerprints and pocket properties, are represented as per-residue descriptors and grouped by alignment position in a way that is analogous to correlated mutation analysis sequence analysis techniques. summary statistics and inter-residue relationships calculated in this way are mapped onto representative 3d structures to aid visual interpretation of the results. all the methodologies are programmed in a purpose built python package called polyphony, in which the “plug-in” architecture allows other descriptors of protein structure generated by 3rd party programs, for example fpocket, ncont, credo and piccolo, to be added with very little effort. matplotlib, pymol, jalview, and ete are used for visualisation of the results. below various novel metrics and algorithms are described for the statistical analyses of local geometry of proteins, compared across ensembles of structures. these include a tailor-made variance measure that is used to distinguish random thermal-like fluctuations from significant conformational changes. in addition, existing techniques such as pca are applied in new ways to differential geometry descriptors of protein conformation. a novel approach to using the output from fpocket to identify distinct and cryptic pockets from an ensemble of structures is described. due to the scalability and automation of the general approach taken it can be applied to large sets of experimental structures or calculated conformations. these advantages are illustrated below by a comparison of conformational commonalities within evolutionary related proteins and by comparing snapshots from an md simulation with x-ray crystal structures of the same protein. the ultimate aim of this work is to extend the repertoire of tools available to those wanting to do ensemble-based drug discovery. however, it is hoped that the methods will also be used to discover fundamental mechanisms in the makeup of protein machines. all the known mechanistic changes in cdk and related structures were discovered afresh but in addition, more subtle, evolutionarily conserved allosteric changes were revealed.  the initial emphasis of polyphony is on the analysis of structures from x-ray crystallography but it can be used on nmr structures and snapshots along the trajectories of protein simulations. cyclin-dependent kinase is often used to validate computational approaches which treat protein-ligand complexes as flexible entities due to the large number of crystal structures deposited in the pdb and to the conformational changes that can be observed in these structures. it is a protein kinase involved in cell cycle progression and inhibitors of its function have been designed with the purpose of exploring their utility as anticancer agents. kinase activity is switched on by binding of cyclin and phosphorylation by cdk activating kinase. cdk is inhibited by the binding of cip and ink family proteins. the binding of cyclin a or cyclin e results in conformational change in which the pstaire region rotates and swings in so that the glu takes its place as part of the catalytic triad. the n and c terminal lobes also move closer together, pivoting about the hinge region, and the t loop moves to expose the catalytic cleft. phosphorylation of thr further organises the substrate binding site by causing additional changes to the t-loop. the glycine rich loop is very flexible and shows structural heterogeneity in cdk structures as well as in protein kinases in general. the cmgc insert is a flexible region in the c-terminal lobe. see figure  for the locations of the substructures mentioned above.figure the structure of the cdk - cyclin complex. this structure is phosphorylated at thr and bound to cyclin a and adp. the pstaire region is shown in orange, the t-loop in yellow, the glycine rich loop in magenta, the cgmc insert in purple, the αef-αf loop in cyan and the hinge region in red. analysis of cdk crystal structures a polyphony script was used to download the structures of the 95% identity sequence cluster containing pdb code 1hck, chain a from the rcsb website. at the time of analysis this cluster contained structures and chains. the vast majority were solved with the p space group and over half had unit cell dimensions of around 72 Å. however, there were examples of other space groups. starting models with pdb codes were cited for of the structures. of these, started directly or indirectly from 1hck or 1hcl, from 1fin, from 1qmz and from 1b39. there were no other significant groupings of structures reported to start from the same model. within grouping average rmsds to their respective starting model range from to  Å. overall average rmsds to these starting structures range from to  Å. the ramachandran plot redrawn a ramachandran plot for all chains combined was plotted and separate secondary structure regions colour coded. figure  2b shows the equivalent plot of curvature against torsion with colours derived from the ramachandran plot. it can be seen that in this latter plot the residues with the same secondary structure are roughly in the same area of the graph, although these areas are less distinct.figure the ramachandran plot compared to the plot of residue curvature against torsion. ramachandran plot for all residues in all chains. colours delineate the regions specified in plot of curvature against torsion for each residue in all chains. colours denote the region the ramachandran plot shown in. contour levels are at fibonacci numbers due to the large number of overlapping points. local conformational variability figure  3a shows how conformational variability for each cdk residue calculated from spline κ/τ compares to rmsf. the differences between the two measures increase as variability increases. rmsf is calculated after fitting a sliding window of contiguous residues to equivalents in a reference structure. this fit becomes arbitrary at high degrees of conformational divergence. figure  3b shows the conformational variability for each residue calculated from φ/ψ torsion angles and spline κ/τ. variation is calculated in different ways for the two descriptors but it is perhaps surprising that they show no correlation. this difference is probably due to the fact that φ/ψ torsion is a more local measure of conformation than κ/τ, which is calculated for a spline whose shape is dependent on the conformation of neighbouring residues. changes in φ/ψ torsion angles could be neutralised by compensatory changes in the conformation of neighbouring residues.figure cα-trace κ/τ variation compared to two other metrics of backbone variation root mean squared fluctuation in cα position after least squares fitting a sliding window of residues to a reference structure and variation in φ/ψ torsion angle. labels show alignment position numbers indicating that the largest differences lie largely in the t-loop. subgroup comparisons in order to divide proteins into monomeric and cyclin bound forms, the chains were clustered using the tanimoto similarity of the per residue protein-protein interaction fingerprint of contacts extracted from the piccolo database. figures  4a and b below show curvature and torsion values respectively for residues in the region around the conserved dfg motif which lies at the start of the t-loop. the plots are colour coded according to the presence or absence of a biologically relevant protein-protein interaction. they clearly show the difference in conformation in this region that occurs upon protein binding as shown in 3d in figure  figure plot of curvature and torsion of cdk chains as a function of sequence for residues n to a in the region of the dfg loop. monomeric chains are coloured blue and protein-bound chains, red. curvature and torsion as a function of sequence. pymol cα trace for residues n to a in the region of the dfg loop. monomeric chains are coloured blue and protein-bound chains, red. note that κ/τ values are not calculated for residues within residues of the termini and the frayed ends adjacent to disordered residues. conformational variability and average temperature factor of course loading all chains into a molecular graphics program, aligning them, colouring them etc. is slow, memory intensive, and the results can be messy. instead, summary views can be generated with polyphony, where properties of the whole ensemble are projected onto one or more representative structures in pymol. figure  shows two such summary views. in figure  6a the per residue relative variability in backbone and side-chain conformation are shown. this can be compared to the relative average normalised cα temperature factor in figure  6b. temperature factor normalisation is done using the per protein mean and standard deviation, calculated after the removal of outlier residues, following the procedure of smith et al.. relative backbone conformational variability and average temperature factor are grossly similar in distribution along the sequence. however the former is a more local property showing sharper variations. residues that stand out as having particularly variable backbone conformations are gly the last glycine in the glycine-rich loop, phe of the dfg motif and val near the end of the t-loop. the side-chains of tyrosines and also stand out as conformationally variable. hinge residues, in which a small conformational change leads to a large shift in a subdomain, are not easily spotted in this analysis.figure conformational variability vs. average temperature factor. pymol putty cartoon of 3pxz chain a showing summaries over cdk chains of relative variability in conformation of backbone κ/τ and sidechains: a grey colour indicates disordered residues in greater than 50% of structures, and average cα temperature factor. a red colour and thick tube indicates the highest variability and highest average temperature factor. the labels highlight particularly conformationally variable residues. the influence of crystal contacts viewed in a 2d plot, it is clear that, in the c-terminal half of the protein, after the t-loop, the conformational variability and average temperature factors are higher in the protein bound structures of cdk as expected, parts of the t-loop and the loop immediately preceding the pstaire helix are more ordered in the bound form. enticingly, this presents a picture of a loss of entropy at the protein-protein interface in the t-loop upon binding, which is compensated for by a gain in entropy elsewhere. however, unnatural crystal contacts are also protein-protein interactions and must be taken into consideration. protein-protein interactions can result in lower temperature factors at the interface. for instance, residues 20– have lower average temperature factors in the bound form but probably as a result of conserved crystal contact with another cyclin molecule in the asymmetric unit. barrett and nobel, using molecular dynamics simulations in solution, found that cyclin bound cdk becomes more rigid on phosphorylation, except for the cmgc insert region which becomes more mobile.figure comparison of conformational variability, average temperature factor and presence of crystal contacts in cdk alignment position is labelled on the x axis. monomeric chains are coloured blue and protein bound chains, red. conformational variability using curvature and torsion average normalised temperature factor crystal contacts as calculated using ccp4 ncont. the value plotted is the average proportion of structures with at least one contact within 5 Å. the location of the pstaire helix is highlighted in orange and the t-loop is highlighted in yellow. hinges and loop flips once groups or clusters of structures are defined, they can be compared in terms of conformation and also interactions. shown in figure  are the variances in curvature and torsion between bound and unbound cdk chains. interestingly, variances in the curvature and/or torsion at individual residues are sharply defined and not smoothed out due to spline fitting. it is clear from figure  and from looking at a 3d representation of the data using the polyphony pymol api, that hinges are highlighted by this analysis. in particular there is a big change between monomeric and protein bound structure in τ at residue gln of the hinge region, the κ and τ at phe of the dfg motif, and a definite signal at asn at the end of the pstaire helix. in addition many residues light up in the t-loop and the cyclin binding region. this ability to highlight what might be called micro-hinges, which include those between secondary structure elements, as well as between domains, is an advantage of this approach. once noise, or random thermal motion, is reduced in this way it becomes apparent that some changes in backbone conformation are compensated for by opposite changes at a neighbouring residue. these paired peaks are the signals produced by loop regions that occupy different but conserved conformations in the groups of proteins. examples that can be seen in figure  are the loops centred at glu and thr the cα trace emerges from the conserved core and bends one way but must then bend back again to return to the conserved core. this sort of conformational change is distinct from a hinge because it only affects local structure. it is worth noting here that, because curvature is always positive, a curve to the left, which is then compensated by a curve to the right, is not apparent from these plots. curvature and torsion must be taken together to observe such differences.figure plot of alignment position against the group variance s i for monomeric chains and protein bound cdk curvature s i and torsion s i positions of the highlighted residues in the structure of cdk principal components analysis in order to facilitate a principal components analysis, a full matrix of data was generated. the pca used the non-linear iterative partial least squares algorithm implemented in pychem. the resulting scores plots for backbone and side-chain conformation are shown in figures  9a and b respectively. it can be seen that the first principal component divides monomeric and protein-bound structures. in the backbone conformation pca indicates a high contribution from the 72– ten loop tip but also cys gln of the hinge region and asn at the end of the pstaire helix. in the side-chain analysis, the conformation arg is the major discriminating factor. it forms a salt bridge with gln on cyclin binding, seemingly helping to lock the pstaire helix in place. the charged residues in this salt bridge are not conserved amongst cdk isoforms.figure pca scores plot for cdk monomeric chains are coloured blue and protein bound chains red. backbone conformation using curvature and torsion. the variance explained by pc is 45% and 51% by pc and pc combined. side-chain conformation using the position of a terminal atom in cartesian coordinates after transformation of the residue to a reference origin. pc explains 24% and pc plus pc explains 35% of the variance. labels show the pdb codes of outliers. pca discriminatory side-chain conformations. salt bridge formation between arg and glu in cyclin bound structures. tyr in monomeric, protein-bound but unphosphorylated, protein-bound and phosphorylated at thr. the rest of the t-loop is cut away for clarity. a structure that stands out as a red dot amongst the blues in figure  9a, i.e. protein bound structure in amongst the monomers is 1buh. this structure is a complex of cdk not with cyclin, but with cell cycle-regulatory protein ckshs which binds at a different site on the c-terminal lobe of cdk the chain b of 1fq1 also stands out in this way in the side-chain conformation pca plot. this is the structure of kinase-associated phosphatase in complex with phospho-cdk. kap binds to a different region on the surface of cdk than cyclin, but one that overlaps with that of ckshs another structure that stands out is 1w98, which is a structure of p-cdk in complex with a truncated cyclin e structures 3pxz, 3pxq, 3py and 3pxf lie in the centre of both the backbone and side-chain pca plots. these are monomeric structures that are located closer to the region occupied by the protein-bound structures. interestingly, they come from a series of structures in complex with allosteric inhibitors that cause conformational changes to the pstaire helix. pc in figure  9a and b both subdivide the protein-bound cdk structures into two groups. the 179– yys motif, which forms part of the αef-αf loop, contributes most to backbone conformation pc in the figure  9a. this loop changes shape on phosphorylation of thr and is known to be coupled to the activation loops in a variety of kinases. this loop lies underneath the t-loop, which is excluded from the pca because so many structures are disordered in this region. the side-chain conformation of tyr in the αef-αf loop, more or less single-handedly separates structures on the pc axis in figure  9b. figure  10b shows its different conformations in monomeric, cyclin bound and cyclin bound plus phosphorylated cdk the first two principal components explain a relatively low amount of the variance, especially for side-chain conformations. unlike the usual cartesian coordinate pca, we are measuring local changes that act in concert rather than big blocks of secondary structure that move together. many local conformational states are conserved or change in an uncorrelated manner and therefore do not strongly influence the first principal components. correlated conformational changes in addition to pca, there is the capability within polyphony of identifying correlated residue conformational differences between individual pairs of residues. this can be done for backbone or side-chain conformations. the κ/τ and x, y, z correlation coefficients, respectively, are averaged for each residue pair. lines are drawn between the most highly correlated residues in an analogous way to young et al. and the covariance web produced by dynamite. these published methods use cartesian coordinates after structure superposition, whilst here the correlations are calculated for local conformations. only correlations between residues with a separation of residues in the protein sequence are shown in figure  but this cut-off is under user control in polyphony.figure correlated conformational change. the most highly correlated pairs of residues separated in sequence by at least amino acids. backbone. the correlation coefficients range from between asn at the end of the pstaire helix and phe of the dfg motif, to side-chains. leu and leu in the core of cdk switch rotamers on protein binding. most of the cdk main chain correlations found in this way are consistent with the results of the pca and known mechanisms of action but also highlight previously unpublished interrelationships between residues. for instance, phe of the dfg motif is correlated with asn identified above as a hinge residue at the end of the pstaire region. cys in the αef-αf loop and gln which is identified above as the pivotal residue in the hinge region, are also linked to these two residues via a network of correlations. surprisingly leu is also paired with these residues, albeit with a lower correlation coefficient but this probably due to crystal contacts with this residue which are much more common in monomeric compared to protein-bound structures. similarly the side-chain conformations of his in the c-lobe are correlated with asp in the n-lobe. this is almost certainly an experimental artefact, due to crystal contacts. more interestingly, there seems to be a repacking of part of the hydrophobic core on cyclin binding. the neighbouring residues leu and leu change conformation in a correlated way. pocket analysis the pocket detection program fpocket was chosen for integration with polyphony because it is an open source program and it provides a druggability score for each pocket. however, it is possible to integrate other programs because of polyphony’s plug-in architecture. as employed by provar, the assignment of pocket attributes to residues facilitates the comparison of structures of the same or homologous proteins, given a sequence alignment. like provar, the fraction of structures in which a particular residue is found to be part of a pocket are reported, here expressed as a percentage. again like provar, these percentages are mapped onto the coordinates of a representative structure as colour-coded surfaces. in polyphony there is a novel methodology which shows the user distinct predicted druggable pockets in the structure in which they occur. in this way cryptic pockets can easily be identified and visualised. figure  shows the results of this procedure applied to the cdk structures. the highest ranked pocket was an atp binding pocket which is connected to a water filled cavity behind the pstaire helix in 2c5y. this cavity would not exist when the pstaire helix moves in on binding to cyclin. thus this cavity is a potential binding site for allosteric inhibitors of cyclin binding. in fact, fragments that bind in this cavity were recently discovered in later published structures. the pocket selected second, which lies under the cmgc insert also shows potential as a binding site for allosteric modulators. the region is known to be flexible from simulations and analysis of crystallographic temperature factors. it forms part of the binding site for cell cycle–regulatory protein ckshs1 and kinase associated phosphatase. there are no small molecule ligands bound to this pocket in the cdk crystal structures but there are in p38α and jnk-1. the third pocket identified is close by in the c-lobe.figure location of fpocket-predicted druggable binding sites selected by the polyphony distinct pocket selection procedure. the blue pocket was ranked highest and is the atp binding site in 2c5y chain a. the second pocket is coloured red and is from 1oi9 chain c. the third in green is from 2r3l chain a. the blue pocket from viewed from above, overlaid with a recently discovered allosteric ligand from 3pxz. the ligand from 2c5y is shown in blue sticks. intermolecular interactions atomic scale interactions are extracted from three different sources 1/ the piccolo database of predicted biologically relevant protein-protein interactions 2/ the credo database of protein - small molecule interactions and 3/ crystal contacts calculated using the ccp4 program ncont. this program is not open source but is available without cost to academic and non-profit institutions. the two data sources are mysql databases created in-house. piccolo is currently being updated monthly and can be downloaded and installed locally for free. the original credo was frozen in april and replacement by a new postgresql database. residue-based structure interaction fingerprints where extracted from piccolo and credo and analysed in various ways including clustering using tanimoto similarity. two papers have explored the differences in inhibitor binding in cdk atp binding site between active and inactive forms, by comparing two and twelve structures. a comparison of credo derived sifts for the protein-bound and monomeric cdk chains with ligands bound, revealed that the set of residues interacting with small molecule ligands were almost identical, except for glu this residue is part of the pstaire helix and the conserved catalytic triad and swings into the atp binding site on activation. none of the ligands bound to monomeric cdk contact this residue, whereas of the bound to protein complexes do. sequence alignment view sequence alignment viewers are, of course, a very useful way of visualising similarities and differences in the amino acid sequences of related proteins. in many cases the simple single letter amino acid codes are augmented with colours to indicate a particular property of each residue such as hydrophobicity or known features such as the occurrence of a post translation modifications. this idea has been extended to include descriptions of the residue conformation or environment with the folded protein. here we are concerned, in the main, with multiple structures of the same protein. thus, a plain sequence alignment is not very informative. however once annotations, such as colour coded b-factors and protein-protein contact counts, are added, the sequence alignment view becomes a convenient way to visualise the similarities and in conformation and interactions across a large number of structures. jalview was chosen for this purpose because users can import colour-coded residue features from a text file. it also allows the import of newick format trees, which can then be used to sort sequences. polyphony can create jalview feature and newick tree files for any of the residue-based properties. figure  shows what this looks like.figure a sequence alignment sorted and annotated by structural and conformational features. dendrogram and sequence alignment view of cdk structures shown in jalview. the alignment is ordered by the tree, generated by hierarchical clustering using backbone curvature-and-torsion pearson similarity. the red vertical line in the dendrogram cuts the chains into large clusters, which correspond to monomeric, cyclin bound and cyclin bound and phosphorylated. amino acids in the alignment are colour-coded by b-factor and protein-protein contacts from piccolo. white gaps indicate disordered residues. application to the analysis of homologous proteins even when only a handful of structures of a given protein are available they can be usefully compared with those of homologous proteins. differences between the structures of homologous proteins include architectural as well as conformational changes. to separate out these two sources of variation, in polyphony, conformational variance is calculated separately by protein. then the locations of intra-protein variance are compared across homologous proteins. below this type of analysis is illustrated on members of the cmgc family of protein kinases. a structural alignment was taken from the homstrad database. figure  shows the aligned variance plots for cdk erk p38α and jnk the regions of known conformational change are highlighted. the glycine rich loop has high values in all proteins. a peak is visible in the hinge region, particularly in cdk and p38α. the dfg motif has a peak in all proteins except erk the activation loop has a high peak in the cdks and erk there also peaks in the consensus plot in lesser-known regions, particularly in the αef-αf loop region.figure plot of the group variance s for cgmc protein kinases. groups were clusters of very closely related structures. cdk n =  in groups cdk in erk in p38α, in jnk in consensus. the colours indicate the glycine rich loop, the hinge region, the dfg motif, the t-loop and αef-αf loop. tyr in cdk interacts with the phosphorylated thr and has distinct conformations in monomeric, cyclin-bound and cyclin-bound plus phosphorylated forms. in p38α the αef-αf peak occurs at met in contrast to cdk tyr this residue is not in close proximity to the phosphorylated residues thr and tyr however, its backbone conformation is strongly correlated with residues 180–. it takes two main conformations, one where the side-chain is surface exposed and one where it is buried under the cmgc insert. the insert itself shifts as result and is or 3 Å closer in to the rest of the protein in the latter conformation. as mentioned in the pocket analysis section above, this region has recently become the focus of inhibitor discovery and design in jnk- and p38α. in erk it was recently discovered that the cgmc insert, which is a helix-turn-helix motif, is involved with nuclear localisation and even dna binding. our analysis reveals what appears to be an evolutionary conserved link between phosphorylation and the conformation and dynamics of this functionally important part of these proteins. use on computed conformations one of the aims of this project was to create new ways to compare experimental structural ensembles with those generated by computational methods. as shown in figure  variability of backbone curvature and torsion provides a more local measure of conformational change than average b-factor. figure  shows variability over p38α x-ray chains and a md simulation downloaded from the model database. the simulation was of the structure 1a9u for 12 ns. for this analysis snapshots separated by 50 ps, covering the last 9 ns of the simulation, were extracted as a multimodel pdb format file. many of the conformational changes seen amongst the crystal structures are likely to occur on a much longer timescale than that covered by the simulation. however there is some qualitative similarity in the variability observed in these structures and those generated by the simulation. for instance, the change in conformation at the key hinge residue met is reproduced, indicating that this is a due to a high frequency motion. there is also clearly some dynamics in the αef-αf loop centred on met which is highlighted in the previous section as being a region of significant conformation heterogeneity. this residue remains buried under the cgmc insert during this short simulation but the αef-αf loop appears to stretch out as the cgmc moves away from the c-terminal lobe as a whole. there are also discrepancies between the crystal structure conformational difference and the simulated dynamics. for instance, in the region 240– there is partial unravelling of the first helix of the cgmc insert in the simulation, whereas it is well conserved between crystal structures. this could be due to the influence of crystal contacts or other reasons. work is under way to further assess the use of polyphony for the analysis of md trajectories.figure backbone conformational variability for x-ray structures and an md simulation of p.  the suite of new methods for the analysis of ensembles of protein structures is selected on the premise that there are advantages to comparing structures without superimposing them. the approach taken is analogous to sequence analysis techniques and starts from a sequence alignment. the alignment itself is trivial because the sequences involved are almost identical. the amino acid types in the columns in the alignment do not differ but the conformations of the residues do. backbone conformational changes can be those that have very little effect on the shape of the protein, such as those near the termini or at the end of flexible loops. others are hinging motions that change the relationship between domains. one difference between these two extremes of variation is that the former tend to be randomly distributed and the latter tend to be consistent across multiple structures. another difference is that the latter tend to coincide with environmental changes, such as a binding event. when a binding event is a protein-protein interaction with a large interface, correlated conformational changes tend to occur at the interface and allosterically. all of these effects can be detected by statistical analysis of all the local conformational differences within an ensemble without the need to superimpose the structures and calculate euclidean distances in cartesian space. once discovered in this way, the consequences of these significant conformational changes can be observed by visual comparison of the structures and the use of existing analytical techniques. the obvious drawback with this approach is the requirement for a large sample size. one must also be wary of systematic experimental artefacts such a conserved non-biological crystal contacts. these issues can be addressed by complementing x-ray structures with nmr derived and computationally generated ensembles. however, there is great benefit not only in solving the crystal structures of new proteins but also in solving the structure of same protein multiple times, especially when co-crystallised with new binding partners. tools and methods for superposition-independent statistical analysis of protein structure ensembles were developed. the general approach, and the individual methodologies within it, were validated by the rediscovery of the published findings of the many authors who have compared crystal structures of cdk since jeffrey et al. in 1995. the major conformational changes that occur on cyclin binding, such as the movement of the pstaire helix and the opening and closing of the gap between the n and c-terminal lobes, were detected via symptomatic changes in hinge residues. in addition, more subtle changes were also detected and found to be conserved in other cmgc kinases. these include correlated changes in the αef-αf loop linking phosphorylation sites on the activation loop to the cmgc insert. this information provides a further clue to the role of this region whose importance is only recently beginning to be revealed and targeted with small molecule ligands in jnk and p38α. thus far no structures of cdk with small molecules bound in this region have been published. the pocket analysis above reveals a pocket can exist there and it is predicted to be druggable, illustrating the utility of ensemble-based drug discovery. 
furthermore, an evolutionarily-conserved conformational link from the location of this pocket, via the αef-αf loop, to phosphorylation sites on the activation loop is discovered. a transient but druggable allosteric pocket in cdk is predicted to occur under the cmgc insert. in addition, molecular dynamics simulations are increasingly being used to further explore the conformational landscape of these complexes. the methodologies are encoded in a python package called polyphony, which is released as open source to accompany this paper. novel methodologies are described for the analysis of multiple structures of a protein. structure-based drug design is an iterative process, following cycles of structural biology, computer-aided design, synthetic chemistry and bioassay. in favorable circumstances, this process can lead to the structures of hundreds of protein-ligand crystal structures. currently, methods capable of the analysis of ensembles of crystal structures and md trajectories are limited and usually rely upon least squares superposition of coordinates. the approaches are tested on crystal structures of cdk and other cmgc protein kinases and a simulation of p38α.
 structure-based drug design is an iterative process, following cycles of structural biology, computer-aided design, synthetic chemistry and bioassay. in favorable circumstances, this process can lead to the structures of hundreds of protein-ligand crystal structures. in addition, molecular dynamics simulations are increasingly being used to further explore the conformational landscape of these complexes. currently, methods capable of the analysis of ensembles of crystal structures and md trajectories are limited and usually rely upon least squares superposition of coordinates.  novel methodologies are described for the analysis of multiple structures of a protein. statistical approaches that rely upon residue equivalence, but not superposition, are developed. tasks that can be performed include the identification of hinge regions, allosteric conformational changes and transient binding sites. the approaches are tested on crystal structures of cdk and other cmgc protein kinases and a simulation of p38α. known interaction - conformational change relationships are highlighted but also new ones are revealed. a transient but druggable allosteric pocket in cdk is predicted to occur under the cmgc insert. furthermore, an evolutionarily-conserved conformational link from the location of this pocket, via the αef-αf loop, to phosphorylation sites on the activation loop is discovered.  new methodologies are described and validated for the superimposition independent conformational analysis of large collections of structures or simulation snapshots of the same protein. the methodologies are encoded in a python package called polyphony, which is released as open source to accompany this paper. issue-copyright-statement© the author 2014 