to date , most allosteric drugs 
 have been discovered through high - throughput 
 screening .
but growing databases of biomolecular structure and sequence 
 data , in conjunction with increases in computing power and improvements 
 in predictive algorithms , are enabling the rational de novo design 
 of allosteric drugs .
given the large number of published algorithms 
 for predicting allosteric mechanisms , it can be difficult to select 
 the most appropriate method for a given target .
this review serves 
 as an introduction for those who wish to use computational techniques 
 to develop allosteric drugs . after a broad overview of allosteric 
 drug discovery ,
first , we discuss bioinformatics and 
 molecular - dynamics methods to identify allosterically important sequence 
 positions .
second , we summarize the computational methods to predict 
 druggable pockets at these functionally relevant sites .
finally , we 
 describe how markov state models and topological analyses can tie 
 these single sequence sites to global protein function and dynamics .
allosteric drugs offer a number of advantages 
 that make them desirable 
 as drug candidates .
allosteric effectors , by definition , alter protein 
 activity by binding to a site distinct from the orthosteric pocket . 

one example of an allosteric system is fructose 1,6-bisphosphatase 
 ( shown in figure 2 ) .
because allosteric 
 sites are typically less evolutionarily conserved , allosteric drugs 
 can be highly selective , even among other members of the same protein 
 family . in some cases , allosteric sites are so unique among proteins that 
 an effector is said to have  absolute subtype specificity .
for 
 example , they can be active only in the presence of the endogenous 
 ligand , thus restricting their effect to certain tissues at certain 
 times , which may slow desensitization .
allosteric 
 effectors are generally saturable , meaning that they 
 have a maximal effect that does not necessarily correspond to complete 
 inhibition or activation .
for example , if the maximal effect is an 80% reduction in signaling , 
 overdosing will not fully eliminate an essential signal .
other advantages can include 
 noncompetitive inhibition ( i.e. , drug 
 activity can not be  overwhelmed  by high concentrations 
 of the endogenous ligand ) and pathway- or substrate - specific modulation , 
 which reduces unwanted activity by specifically targeting a single 
 protein function . for 
 example , if a protein is involved in multiple pathways , an allosteric 
 effector may impact the activity of each pathway differently depending 
 on the systems - biology context .
if a protein acts on multiple substrates , 
 the impact on activity may depend on the biological context .
despite many potential advantages of allosteric therapeutics , it 
 has been challenging to identify predictive approaches to discovering 
 allosteric drugs . in recent decades ,
the pharmaceutical industry has 
 favored more traditional targets for three primary reasons : the relative 
 ease of assay development around orthosteric sites ; access to high - throughput , 
 high - resolution x - ray crystallography ; and advances in ligand- and 
 receptor - based computational methods to optimize ligand - binding affinity 
 at a substrate - competitive site .
this structure - based approach is 
 thought to significantly reduce the time and cost of hit - to - lead and 
 lead - to - drug development by reducing the number of compounds that 
 need be synthesized .
comparing 
 computer - aided drug discovery ( cadd ) and high - throughput screening 
 ( hts ) reported that the two methods had hit rates of 35 and 0.021% , 
 respectively .
in contrast , allosteric 
 drugs are uniquely challenging from a rational 
 drug - design perspective . because experimental assays typically measure 
 orthosteric function rather than ligand binding at the allosteric 
 site , efficient development of allosteric drugs requires that the 
 complex structure  activity relationships ( sars ) governing both 
 binding affinity and allosteric activity be considered simultaneously .
further , allosteric sites are 
 less likely to be evolutionarily conserved . while this enables increased 
 subtype specificity , it also increases the chances of evolved resistance and can complicate testing in 
 evolutionarily distant animal models .
additionally , allosteric effectors are particularly susceptible 
 to  mode switching ,  where relatively minor chemical 
 changes can drastically affect ligand efficacy .
structurally similar drug metabolites , therefore , may have varying 
 and unpredictable distributions and allosteric effects .
optimizing allosteric modes of action requires methods that are 
 very different from those used in orthosteric drug discovery .
while drug designers may desire to target a single protein function , 
 an allosteric effector may also alter other functions , hindering a 
 full mechanistic understanding of the pharmacology .
also , the benefits of spatiotemporal specificity are lost 
 if the distribution of the endogenous ligand changes with progression 
 of the disease state .
finally , assessment 
 of the limited number of known allosteric pockets indicates that they 
 are generally shallow and present flat 
 sars .
these structural features similarly 
 challenge existing rational drug - discovery paradigms and the general 
 practice of developing selective compounds by optimizing affinity . despite these challenges
, allosteric drug discovery has gained 
 momentum recently due to a number of developments .
first , several allosteric drugs across a broad range of 
 pharmacological target classes have been rationally designed , encouraging pursuit of others , as evidenced by the number of allosteric 
 drugs currently in clinical trials .
the 
 recent elucidation of new membrane protein crystal structures for 
 gpcrs and ion channels have assisted in the structure - based design approach to these successes . 
 finally , advances in our understanding of allosteric mechanisms have 
 supported development of additional rational design strategies ( see 
 below ) .
our understanding of allosteric mechanisms has advanced 
 considerably 
 since the initial conception of monod , wyman , and changeux .
this revised view supports newly established and emerging computational 
 advances that comprehensively map conformational landscapes and predict 
 communication between allosteric and orthosteric sites .
for example , 
 the physical mechanisms of allostery generally alter the entropic 
 and enthalpic factors that define the conformational landscape and , 
 therefore , govern protein function .
the observed correlation between allosteric modulation 
 and protein structural dynamics is varied : major conformational rearrangements 
 occur in some cases , as compared with subtle shifts in conformational 
 populations in others .
an excellent metaphor for these phenomena 
 is kornev and taylor s classification of  domino  
 versus  violin  models of allosteric signal transduction ( figure 1 ) .
further , allosteric signals are transmitted through a range 
 of structural motifs , from rigid core regions to flexible linkers .
allostery may occur through essential residues 
 along a single allosteric path or through 
 many weak pathways connecting one site to another , acting in concert .
as such , it is not surprising that many of the 
 allostery - prediction methods discussed in this review ( some of which 
 do not use structure / geometry information at all ) in practice may 
 identify noncontiguous groups of residues as being allosterically 
 linked .
such predictions should not immediately be assumed to be wrong 
 but rather may indicate a non-domino  model of allostery . when a 
 small molecule binds to the allosteric site of a protein , 
 information is transferred through the protein molecule to its active 
 site .
the first 
 mechanism , here defined as the  domino model  , is a 
 sequential set of events propagating linearly from the allosteric 
 site to the active site .
binding of the effector triggers local structural 
 changes that sequentially propagate via a single pathway to the active 
 site .
it was suggested that this mechanism is applicable for the pdz 
 domain family , g protein - coupled receptors , 
 the chymotrypsin class of serine proteases , and hemoglobin .
the second mechanism , defined here conceptually 
 as a  violin model  , is based on vibration pattern changes 
 inside the protein . in a violin
its pitch can be changed by a slight 
 movement of the violin player s finger on the fingerboard . 

information about the finger movement is , thus , transferred throughout 
 the whole body of the violin with no specific pathway for the signal 
 transduction . by analogy ,
protein allosteric site is a fingerboard 
 of the protein and a small signaling molecule is the player s 
 finger .
if a protein is in a particular vibration mode , it is possible 
 to suggest that binding a small effector molecule to a specific site 
 can change this mode .
the signal , thus , will be spread throughout 
 the whole protein including its active site
. the  domino model  
 is a reliable way to transfer information in a macro world , but on 
 a molecular level , with significant thermal motions of the protein , 
 this mechanism will be prone to random triggering of the domino chain 
 reaction , creating noise in the signaling system . thermal motions 
 in the case of the  violin model
in fact , the permanent motion of the molecule is a prerequisite 
 for this mechanism . reproduced with permission from ref ( 45 ) .
recent work has also revealed that protein allostery is not 
 merely 
 a transition between two discrete protein conformations , as initially 
 thought , but rather a shift in the equilibrium populations of many 
 conformations , induced by effector binding .
it is becoming increasingly clear that the kinetics of these transitions 
 define the mechanisms of allostery .
empowered 
 with these new understandings and advances in molecular simulation 
 ( in terms of speed of calculation and improved methodologies ) , the 
 era of allosteric drug discovery is now on the cusp of radical advancement .
orthosteric and allosteric pockets ( yellow and red , respectively ) 
 are bound to an endogenous ligand and an allosteric effector , respectively . 

note that the allosteric site is distant from the orthosteric site 
 such that there is no overlap between the bound poses of the allosteric 
 and orthosteric ligands .
despite the distance between them , the allosteric 
 effector measurably modifies the enzymatic activity at the orthosteric 
 site .
so - called tried - and - true  design principles  
 are still 
 being developed .
however , a few general principles have begun to emerge . 
 for example , many argue that it is insufficient to design a ligand 
 that merely binds to an allosteric site ; rather , the effector must 
 make contact with certain key binding - pocket atoms to have the desired 
 effect .
these key atoms can often be identified 
 through mutagenesis experiments and crystallographic studies of other 
 allosteric ligands .
a promising set of design principles is 
 encapsulated in the  allo - network 
strategy , a rational 
 approach that adopts two simultaneous but orthogonal approaches to 
 ligand design . on the protein - structure 
 level ,
the primary focus is to target a single protein function or 
 an interaction with a single partner . on the signaling - pathway level , 
 the  allo - network 
strategy suggests targeting less - connected 
 upstream proteins instead of the more direct , though potentially highly 
 connected , signaling proteins themselves .
when applied to early stage 
 design , the allo - network method is predicted to increase the likelihood 
 that a given allosteric effector will proceed through the drug - approval 
 process .
several 
 published examples for recently approved allosteric drugs 
 serve to illustrate the current state of the art for emerging allosteric 
 drug design principles .
they also represent the significant advances 
 that have been made to utilize structure - based methods for challenging 
 druggable sites such as protein  protein interfaces and for 
 membrane proteins such as ion channels .
the available details of the 
 discovery and optimization of these compounds do not include the methods 
 discussed in this review ; however , they highlight where these predictive 
 techniques could contribute to the allosteric design process . in 2011 , gilmartin et al . of glaxosmithkline reported the discovery 
 of a pharmacokinetically optimized allosteric mek inhibitor , gsk1120212 .
the first generation inhibitor was discovered 
 by high throughput screening , and the 
 subsequent ternary crystal structure showed the allosteric pocket 
 to be adjacent to the orthosteric atp - bound site . by 2012
there were 14 allosteric mek1/2 inhibitors in clinical 
 trials , because it was recognized that 
 an inhibitor developed for this allosteric pocket afforded two very 
 unique opportunities to avoid adverse clinical effects : the high doses 
 required to compete against 1 mm cellular atp concentrations and inhibition 
 of closely related atp - binding sites in other kinases .
the unique 
 efficacy properties of gsk1120212 highlight both the opportunity and 
 challenge of allosteric drug design .
report that , 
 although some other mek allosteric drugs demonstrate inhibition of 
 the erk1/2 pathway in vitro , this has not translated into efficacy 
 in patients .
gsk1120212 has since been 
 approved in the united states under the name trametinib for treatment 
 of metastatic melanoma caused by the v600e mutation . in 2012 , 

these inhibitors produce an allosteric effect by binding 
 at an interdomain interface and stabilizing a preexisting autoinhibited 
 state of the protein .
the original discovery of the allosteric site 
 was accomplished using a fragment - based hts technique , followed by 
 optimization using x - ray crystallography and structure - based sar . 

the authors discuss their experience with a few confounding factors 
 in the allosteric design process , namely the need to use the full - length 
 protein in their screening construct to observe the exerted allosteric 
 effect , and the subsequent directed evolution study of resistance 
 mutations that could occur at the allosteric site .
hackos et 
 al . of genentech published on the discovery of positive 
 allosteric modulators ( pams ) for glun2a - containing nmda receptors 
 in 2016 .
pams are allosteric ligands 
 which increase the effect of the endogenous signaling molecule and 
 do not cause a change in its absence .
the authors note that the validation of this allosteric 
 site was reinforced by its similarity to an analogous allosteric site 
 in ampa receptors , but that the nmda receptor site has elements of 
 asymmetry that the ampa receptor site did not . in comparing two similar 
 compounds , gne-6901 and gne-8324 ,
the authors make comments that indicate 
 evidence of mode switching or a shallow sar landscape , and they further 
 characterize the details of the allosteric mechanism using mutagenesis 
 experiments . in summary
these examples demonstrate allosteric 
 drug discovery 
 can be successful at protein sites often considered to be undruggable . 

it is apparent that these successes can be further built upon through 
 computational methods that allow for rational rather than serendipitous 
 hts discovery of new allosteric binding sites and a deeper understanding 
 of allosteric mechanisms that overcome design challenges such as mode 
 switching .
protein - sequence analysis 
 is a useful tool to detect and characterize allosteric pathways and 
 pockets . here , we classify sequence - based methods into two groups : 
 ( 1 )  single site 
methods , which produce a list of individual 
 functional sequence positions ; ( 2 )  coupled site 
methods , 
 which produce a list of groups comprised of two or more sequence positions 
 that appear to be functionally linked based on their coevolution .
these 
 challenges include how to select and aggregate clean , relevant sequences 
 as input ; interpret the output ; and integrate sequence - analysis results 
 with other forms of data . determining the biological meaning of a 
 strong signal is also problematic .
while many analysis methods identify 
 evolutionarily important residues , the specific biological role of 
 these residues can not be inferred without additional knowledge . for 
 example
, it is difficult to determine , based on sequence alone , whether 
 an evolutionarily significant residue plays an allosteric role , or 
 whether its role is related to another essential process ( e.g. , substrate 
 binding , maintaining protein structure , etc . ) .
most techniques require many sequences to establish statistical significance . 
 to obtain the required number of sequences , researchers often lower 
 the stringency of their search parameters , resulting in alignments 
 that contain sequences with lower similarity or incomplete coverage 
 of the original query .
while some analysis methods manage to detect 
 meaningful coevolution over a wide range of multiple sequence alignment 
 ( msa ) conservation and noise levels , others are more susceptible to 
 messy data . for a more complete discussion 
 of these topics and
how they affect coevolution analysis methods , 
 readers are directed to an excellent recent review by de juan et al . by our definition , single - site evolutionary 
 analysis methods return 
 a list of predicted functional sequence positions but do not suggest 
 specific linkages between sites .
once a researcher has constructed 
 an msa , the conservation or phylogenetic relevance of each column 
 can be used to infer the evolutionary importance of each sequence 
 position .
this importance is sometimes a hallmark of thermodynamically 
 critical residues that participate in allostery .
though single - site 
 methods only return a list of single high - scoring sequence positions , 
 the inner workings of some single - site methods are based on the aggregate 
 or correlated behaviors of multiple sequence positions ( e.g. , to determine 
 baseline residue probabilities within a multiple - sequence alignment 
 or construct a phylogenetic tree ) .
single - site methods for detecting 
 allostery are advantageous because they lack much of the noise often 
 associated with correlation analysis .
these analyses are also appealing 
 because of their simplicity : there are usually fewer parameters to 
 set , and the results can be visualized directly by highlighting key 
 residues on a three - dimensional ( 3d ) protein structure .
shannon 
 entropy , one of the simplest nontrivial sequence - analysis metrics , was used widely in early works to identify conserved 
 sequence positions for drug - design or mutagenesis experiments .
similar in form to thermodynamic entropy from 
 statistical mechanics , shannon entropy measures the population diversity 
 of residues at a single msa position .
it is also central to mutual 
 information ( mi ) , a popular coupled - pair metric .
the mi of two sequence 
 positions is defined as the sum of the individual position entropies , 
 minus the entropy of the positions considered jointly .
while we do 
 not cover the mathematical details of these methods here , interested 
 readers are directed to previous articles on these topics .
shannon entropy does not consider amino acid similarity ( e.g. , 
 in the shannon entropy framework , a leucine - to - isoleucine mutation 
 is considered mathematically equivalent to a leucine - to - arginine mutation ) . 

other entropy measures , such as the relative shannon entropy ( also 
 called the kullback  leibler divergence ( kld ) ) and the von neumann entropy , attempt to 
 overcome this limitation and , as a result , may be more useful in the 
 search for allosteric sites
. relative shannon entropy / kld accounts 
 for some measure of the protein s chemical environment by considering 
 each mutation with respect to the background amino acid frequencies 
 calculated from the msa
. this analysis may be particularly useful 
 when searching for allosteric sites in proteins that reside in membranes 
 or other noncytosolic compartments , where background residue probabilities 
 or mutational preferences may be biased due to different biochemical 
 contexts .
in contrast , von neumann entropy , a concept borrowed from 
 quantum statistical mechanics , is calculated using amino acid similarity 
 matrices .
identifying an optimal amino acid similarity metric is nontrivial 
 and may well depend on the nature of the system ( e.g. , in a well - packed 
 protein , residue size may be a sensitive metric , whereas surface - site 
 comparison may require the user to prioritize charge ) . in a recent 
 publication describing these types of entropy ,
johansson and toh explored 
 how the two metrics can be mixed to detect enzyme active sites with 
 maximum sensitivity .
zhang et al . 
 constructed a variety of new analysis methods in 2008 
 by combining shannon or von neumann entropy , phylogenetic tree structure , 
 and a novel gap - treatment approach . in 

benchmarking their method , they compared their results to evolutionary 
 trace and consurf ( discussed in greater detail below ) .
two of their 
 hybrid approaches outperformed all other techniques in detecting significant 
 residues across a variety of proteins : the improved zoom method , which 
 incorporates a tree breakdown of subalignments , and the physiochemical 
 similarity zoom method , which extends the improved zoom method with 
 von neumann entropy and tree - branch - size normalization .
lichtarge , 
 bourne , and cohen pioneered the evolutionary trace ( et ) method .
the 
 approach has become quite popular , largely because the algorithm is 
 intuitive and its results are readily visualizable .
et aligns a number of sequences and constructs a phylogenetic 
 tree , and then monitors the conservation of sequence positions at 
 major tree branching points . by slicing the tree at different similarity 
 cutoffs ,
the evolutionary significance of these sequence positions is implied 
 by their conservation in the sequences beyond the next branch . in 
 their first paper , the authors demonstrated 
 that et can detect functionally important sites in sh2 , sh3 , and dna - binding 
 domains .
work has since been published on et validation , parameter 
 optimization , and best - use practices . in a method often referred to as  difference - et ,  

the user runs et on two related proteins and considers differences 
 in the high - ranking residues and their scores .
the sequence positions 
 with strongly varying scores may suggest specificity determinants 
 or differences in allosteric and/or orthosteric mechanisms .
notably , 
 difference - et has been used in the study of gpcr specificity . to better account for varying rates of evolution in different 
 subtrees 
 and correlated mutations , in 2004 mihalek et al .
this work also introduces 
 the zoom et method ( not related to improved zoom , above ) , which adds 
 higher weight to sequences that are more similar to the protein of 
 interest . in the introductory work
, they used real - valued and zoom 
 et to detect the functional residues in a kinase domain , and then 
 compared the performance of both methods to regular ( integer - valued ) 
 et and entropy .
given unpruned sequence data sets , the real - valued 
 et and zoom et methods outperformed the others by a significant margin . 

in contrast , integer - valued et prevailed in most respects when pruned 
 data was available .
an automated web server is available to perform 
 real - valued et calculations , generate reports , and visualize results 
 ( http://mammoth.bcm.tmc.edu/etserver.html ) . in 2008 , merkl et al
. 
 introduced a method called h2r that serves as a segue between single - site 
 and coupled - site approaches .
h2r generates 
 a mutual - information matrix for an msa , and then discards all but 
 the strongest detected coupled pairs . for each sequence position
k , the method returns conn(k ) , the number 
 of top - ranked pairs that include k. initial work 
 proved that h2r can successfully detect functionally significant residues 
 across a range of proteins .
more recently , h2rs , an improved version 
 of h2r , has been released .
this method 
 modifies the original by using von neumann instead of shannon entropy 
 and performing more detailed checks for statistical significance . 

h2rs is available as a web server and a stand - alone application at http://www - bioinf.uni - regensburg.de/. second - order sequence analysis detects residue 
 pairs that mutate 
 in concert more frequently than would be expected given random genetic 
 events .
coevolving residue pairs are assumed to be functionally linked , 
 often because they serve essential roles in allostery or structural 
 integrity .
the immediate output of second - order allostery analysis 
 is a list of residue pairs with associated correlation strengths . 

combining these individual pairwise correlations into a single picture 
 of the entire protein is a separate task . on the most basic level , 
 the strongest correlations that include a residue or site of interest 
 can suggest thermodynamic coupling to other sites , possibly related 
 to allostery .
more complex analyses use hierarchical clustering or 
 principal component analysis to analyze these linkages and uncover 
 strongly linked networks of coevolving residues .
several 
 simple yet reliable residue - coupling analyses have maintained a presence 
 in the literature over the past decades .
these basic approaches are 
 advantageous because they are easier to understand and have been shown 
 to score consistently well in a wide range of tests .
however , they 
 may fail to detect correlations in more complex cases . though more complex methods exist , many of these basic methods still 
 appear as analysis options in coevolution - detection software packages 
 and web servers . in this review ,
mutual information 
 ( mi ) is one of the most straightforward and long - lived coupling metrics . 

the mi between two sequence positions is defined as the sum of the 
 shannon entropies of both positions , minus their joint entropy .
due 
 to its simplicity and favorable mathematical properties , mi analysis 
 is the basis for a number of more complex coevolution methods .
for example , uncorrelated pairs 
 of high - entropy sequence positions are likely to have a higher mi 
 than uncorrelated pairs of low - entropy positions . to compensate for this and other shortcomings ,
further , methods to estimate baseline values for correlation ( e.g. , 
 resampling or sequence shuffling ) can improve mi analysis .
another relatively direct 
 coevolution metric , the mclachlan - based substitution correlation ( mcbasc ) , looks for similar patterns of variation in the 
 columns of an msa , weighting for residue similarity using the mclachlan 
 scoring matrix .
analogous methods can 
 be constructed using different substitution matrices , but mcbasc continues 
 to be a popular choice in the literature . in 2002 ,
kass and horovitz analyzed 
 the groel complex using a chi - squared test to detect significant residue 
 coevolution in an msa .
the analysis suggested intra- and interchain 
 contact pairs and has continued to appear in the literature under 
 the name 
method developed by lockless 
 and ranganathan is perhaps the most widely used sequence - based method 
 for allostery prediction .
sca draws an 
 analogy to statistical physics by calculating a  coupling energy  
 between each sequence - position pair .
the original sca method computes 
 a conservation value for each sequence position i in an msa , applies one of several types of perturbation to another 
 position j ( depending on the sca version ) , and finally recalculates the conservation 
 at position i for the sequences that satisfy the 
 perturbation . by calculating the change in individual and joint conservation 
 over a variety of perturbations
, sca establishes a  coupling 
 energy  that indicates the evolutionary coupling of positions i and j. the output of the sca method 
 is an n  n matrix of coupling 
 energies , where n is the number of sequence positions 
 in the alignment . in early work , the researchers manually identified 
 strongly coupled residue pairs that included one functional member 
 ( per experiment ) .
more recent versions of sca have grouped this matrix 
 into meaningful clusters of coevolving residues using hierarchical 
 or spectral clustering .
refinements 
 of sca have achieved improved statistical properties 
 by resampling the original distribution . in a 2011 paper , sca
was effectively used to engineer a light - sensitive 
 lov2 domain onto the surface of dhfr at a location that sca had identified 
 as energetically linked to the enzyme active site .
some variants of 
 the resulting protein chimera were found to have acquired light - dependent 
 activity .
further work has used 
 sca to design artificial ww domains . in 2011 , 
 an sca analysis of antigen 85c from mycobacterium 
 tuberculosis suggested new sites that potentially could be 
 exploited in drug design . a number of 
 projects have also demonstrated how sca can be used to target mutations 
 that affect protein function . inspired by earlier work on the sequence correlation entropy 
 ( sce ) 
 method , dima and thirumalai published 
 an sca variant in 2006 .
this variant 
 controls for specific protein composition by calculating the background 
 probability that a given amino acid will be present at a random sequence 
 position .
this probability is determined by considering only the sequences 
 being analyzed , as opposed to all sequences in the swiss - prot database .
further , they borrowed a coupled two - way clustering procedure from 
 gene - sequence analysis to define the sectors . in validating this method
, the authors analyzed the pdz , 
 gpcr , and lectin families of proteins and were able to quantitatively 
 predict functional residues , which were in agreement with experimental 
 findings .
 
 method in which correlation is established by excluding certain sequences 
 from an msa and monitoring how entropies change .
another popular perturbation - based 
 method was published in 2004 by dekker et al .
this method , explicit likelihood of subset covariation 
 ( elsc ) , relies on similar principles but returns correlation scores 
 in the form of probabilities rather than statistical energies .
elsc 
 was shown to be superior to sca in contact prediction when tested 
 on a range of protein families .
it has since been implemented on web 
 servers and has been a popular benchmark method 
 in the literature . in 2009 , 

weigt et al . proposed a mutual - information - based method called direct 
 coupling analysis ( dca ) that disentangles directly interacting residues 
 from large networks of indirectly coupled sequence positions .
while this method is typically used in structure 
 prediction to identify spatially adjacent sequence positions , it may 
 find application in the study of short - range allosteric interactions . 

a more efficient implementation of the dca method , known as  mean 
 field 
both introductory papers show that dca is a robust predictor of 
 both intra- and interprotein contacts and that it can hint at the 
 existence of unobserved protein conformations .
related work has shown 
 that dca can be used in conjunction with structural models to generate 
 predictive models of protein complexes , determine the sequence positions that contribute to protein - interaction 
 specificity , and describe the conformational 
 ensembles of proteins in crystallographic or near - crystallographic 
 states .
a web server and software package 
 are available to perform dca analysis at http://dca.rice.edu/portal/dca/home .
psicov is another popular 
 contact - detection method that may find productive use in the study 
 of allostery .
mathematically , psicov 
 relies on an estimated inverse of the msa covariance matrix , which 
 acts as a matrix of correlations between all sequence - position pairs 
 that inherently controls for the variations in all other positions . 

the code has been published online at http://bioinfadmin.cs.ucl.ac.uk / downloads / psicov/. recurrence quantification analysis ( rqa ) , another second - order sequence - analysis 
 technique , is best used when much is already known about the mechanism 
 under investigation ( e.g. , physiochemical amino acid properties such 
 as charge or hydrophobicity are known to drive the allostery ) .
rqa 
 itself is a general method in nonlinear dynamics : in the context of protein sequences , it considers a scalar - value 
 vector that represents some property of a given sequence . in introductory 
 work by zbilut et al .
, the method was 
 used to properly classify 56 tem-1/-lactamase mutants with 
 impaired function based on their hydrophobicity profiles .
further 
 rqa work used hydrophobicity scores to classify proteins as allosteric 
 or nonallosteric , study p53 mutants , and reveal interaction partners in viral - envelope 
 proteins . in 2005 ,
colafranceschi 
 et al . investigated the effect of choosing different physicochemical 
 amino acid descriptors and changing the numerical parameters of the 
 rqa algorithm .
more recently , a comparison 
 method based on rqa measurements , known as cross - rqa , effectively 
 detected protein allostery .
interested 
 readers are directed to a review by zbilut - webber , which provides 
 examples of rqa applied to a range of computational biology problems .
some work 
 has been done to competitively benchmark the performance of these 
 methods . in 2004 ,
fodor and aldrich compared omes , mi , sca , and mcbasc 
 in a variety of tests . in short , the study found that performance 
 is largely dependent on the way that different methods determine background 
 residue probabilities and handle positional conservation .
a follow - up study investigated how effectively 
 coevolution analysis finds thermodynamically linked residue pairs . in general , spatially contiguous linked pairs 
 were detected , but long - range couplings did not agree with experiments . in 2010 ,
brown and brown introduced a new pair - scoring method , 
 called z - scored - product normalized mutual information ( znmi ) , and 
 compared it to the accuracy and reproducibility of mi , two versions 
 of sca , omes , and elsc .
the authors presented 
 a thorough meta analysis of method performance and the impact of input - parameter 
 selection .
though none of the tools tested was particularly powerful , 
 znmi was the most robust prediction tool .
brown and brown also found 
 that the use of multiple subalignments produced more accurate and 
 reproducible results .
a comparative analysis of sca and dca 
 revealed that the top 35 
  sectons  found via spectral clustering of the dca matrix 
 corresponded to pairs , triplets , and quadruplets of spatially contiguous 
 residues .
in contrast , a similar analysis 
 of the sca matrix produced spatially adjacent clusters of many residues 
 each .
these different results validate the stated goals of each method : 
 dca aims to find contacting pairs , whereas sca aims to find potentially 
 distant groups that are thermodynamically linked in a certain function . in 2014 , pel et al .
investigated seven coevolution analysis 
 methods to find the hallmark covarying pairs in gpcr alignments .
they considered three variants of mi , mclachlan 
 based substitution correlation , sca , elsc , and omes .
omes and elsc 
 were the most robust methods for finding the residues responsible 
 for subfamily divergence .
their study tests 
 omes , two variants of mi , 
 sca , psicov , and di , and finds that psicov and di are best at identifying 
 contacting residues .
omes and mip excel at removing false positives 
 from the lists of predicted contacts .
while the authors focused on 
 detecting inter- and intramolecular contacts , their analysis also 
 provided useful insights to guide the productive use of each method . 

for example , all methods benefit from repeatedly shuffling the msa 
 and rerunning the analyses in order to provide a baseline and remove 
 false positives .
finally , the authors found that the consensus of 
 di and psicov provides a more robust prediction of contacting residues 
 than any single prediction method alone .
the software used to perform 
 this analysis is available through the prody evol program ( http://prody.csb.pitt.edu/evol/ ) . in the course of introducing 
 new types of mi analysis ( dbzpx2 , 
 dgbzpx2 , and nbzpx2 ) and evaluating the effectiveness of msa simulation 
 ( a topic beyond the scope of this paper ) , ackerman et al . in 2012 

new  methods ( the zpx2 family , dca , 
 and log(r ) ( not discussed here ) ) were significantly 
 superior to the  old  methods ( omes , mcbasc , elsc , and 
 sca ) .
contact 
 map webviewer ( cmweb ) automatically 
 constructs an msa and visualizes a variety of coevolution analyses : 
 mutual information , sca , elsc , omes , and an early method presented 
 by gbel et al .
the same server 
 can also compare the results of these methods to user - uploaded data 
 ( e.g. , results the user obtained using some other type of analysis ) . 

the cmweb server can be accessed at http://cmweb.enzim.hu/. the coevolution analysis of protein 
 residues server hosted by the gerstein lab ( http://coevolution.gersteinlab.org/coevolution/ ) can perform a large number of the coupled - site analyses presented 
 in this review , including sca , elsc , mi , and mcbasc - type methods employing 
 different scoring matrices . the server can also validate the results 
 of these methods in predicting residue distances in a crystal structure .
mistic ( mutual information server to infer coevolution ) is an automated 
 web server that accepts user - submitted msas or collects them from 
 pfam .
mistic uses a corrected form of 
 mi to infer coevolving pairs and offers several analysis methods that 
 combine structure and coevolution .
it 
 can be accessed at http://mistic.leloir.org.ar/. caps ( coevolution analysis using protein sequences ) is a 
 unique 
 algorithm that combines phylogenetic , 3d , and msa data to predict 
 coupled sequence positions .
versions 1 and 2 are 
 hosted on web servers at http://bioinf.gen.tcd.ie/caps/ and http://caps.tcd.ie/ , respectively .
the interprotein - correlated mutations server ( i - coms , http://i-coms.leloir.org.ar/index.php ) focuses on detecting contacts at protein  protein interfaces , 
 though it can also return intrachain correlations .
the server automatically builds alignments ; performs mi , 
 dca , or psicov analysis ; generates visualizations of the results ; 
 and allows users to download data taken from various points in the 
 data - collection and analysis workflow . in 2012 ,
they employed an automated 
 workflow to control for various types of noise in sequence alignments , 
 using the msa sequence profile to establish prior knowledge about 
 the protein .
while the authors only studied a few mi variants , they 
 stressed that their profile - based method could be extended to more 
 complex analysis techniques .
their approach , correlated mutation analysis 
 tool ( cmat ) , is available on a web server at http://binfolab12.kaist.ac.kr / cmat/. access to consurf , a single - site detection method similar 
 to et , 
 is available at http://consurf.tau.ac.il/. the prody python package , 
 referred to above , can compute a variety of coevolution metrics . in 2014 ,
released the latest version of the powerful 
 bio3d r package for protein structure and sequence analysis .
while it focuses on structural analysis , the 
 package can compute shannon entropy and offers useful functions to 
 create and manipulate sequence alignments . also in 2014 , li et al .
published the cormut software package for r , 
 which computes mi , a metric called the 
jaccard index ,  and the conditional selection pressure metric ka / ks.table 1 summarizes the selected coevolution web servers / software 
 packages and their capabilities discussed here .
allosteric 
 effectors communicate with orthosteric pockets by altering protein 
 dynamics , either through large - scale structural changes or through 
 more subtle changes in correlated residue motions . a proper understanding 
 of the allosteric mechanism is impossible unless one fully appreciates 
 the underlying dynamic conformational flexibility and energetic landscape .
computational methods to characterize receptor flexibility include 
 molecular dynamics ( md ) and monte carlo simulations .
briefly , md simulations 
 represent molecular systems as beads ( e.g. , atoms ) connected by springs 
 ( e.g. , bonds ) .
newton s equations of motion are numerically 
 integrated to propagate the dynamics of this classically formulated 
 system .
simulation time is typically discretized into steps of 1 or 2 fs , 
 as required to accurately simulate the fastest atomic motions ( i.e. , 
 hydrogen bond stretching ) . at each step , the potential energy of the 
 system is calculated by considering the energies associated with bonded 
 atoms ( e.g. , due to bond stiffness , angle bending , torsion rotation , 
 etc . ) and nonbonded atom pairs ( e.g. , due to electrostatic and van 
 der waals forces ) .
the forces on each atom , calculated from the negative 
 gradient of the potential energy , are then used to update the atomic 
 positions at each step . in a related method , called 
( amd ) , the underlying potential energy landscape is modified by raising 
 the energy wells .
these modifications facilitate transitions between 
 states that do not normally occur on the time scales of traditional 
 md simulations .
interested readers are directed to recent reviews 
 that describe md simulations in the context of drug discovery .
in contrast , monte carlo based methods take a stochastic approach 
 to conformational sampling . at each monte carlo 
 step ,
a markov - chain procedure is used to generate a slightly modified 
 system conformation .
this new conformation is then accepted or rejected 
 at random , biased by the relative potential energy of the new conformation 
 as compared to that of the previous step ( the so - called 
as these simulations 
 are stochastic , they do not typically explore conformational space 
 via the same time - dependent path traced in an md simulation .
indeed , 
 a key limitation of the monte carlo method is that the time information 
 connecting the states is lost .
nevertheless , monte carlo simulations 
 are widely used because they tend to more efficiently sample system 
 conformations that are statistically representative of the equilibrium 
 state .
as mentioned 
 above , allosteric signals can be transmitted through 
 large - scale conformational changes or subtle shifts in the correlated 
 motions / interactions of individual residues .
accessing the time scales required to observe many large - scale conformational 
 changes is presently computationally intractable using unbiased simulation - based 
 methods , but these simulations are well - suited 
 to the study of allosteric signals that are transmitted via fast , 
 local fluctuations ( i.e. , nanosecond - scale changes in the coordinated 
 motions of adjacent residues ) .
distilling the dynamics sampled 
 by these simulation - based methods into a simple network representation 
 of protein motions often facilitates allosteric analysis . rather than tracking the position of every atom
, each protein constituent 
 ( e.g. , amino acid ) is represented as a single node ( located , for example , 
 at the residue center of mass ) .
each pair of nodes is connected by 
 an edge whose length is inversely proportional to the degree of interdependence 
 between their motions , such that nodes connected by short edges are 
 highly interdependent .
different metrics of  interdependence  
 have been used , including metrics based on correlated motions ( e.g. , 
 mutual information , fluctuations in atomic positions , etc . ) or the number of specific noncovalent 
 interactions .
once a network representation 
 of protein dynamics has been constructed , 
 various network - analysis techniques ( described below ) can identify 
 important pathways of communication between distant residues that 
 may contribute to allostery .
interested readers are directed to a 
 recent review by atilgan et al . for more detailed information .
a number 
 of simulation - analysis techniques consider the interdependence of 
 individual residues along an allosteric path or paths . in 2008 , bradley 
 et al . analyzed molecular dynamics simulations of escherichia 
 coli nikr , a transcriptional repressor of the nikabcde nickel 
 permease , using a novel network - space clustering approach .
residue interdependence 
 were considered , based on fluctuations in residue motions and the 
 number of polar contacts between residue - pair members .
first identified sets of residues whose motions and/or contact 
 counts were correlated with residues in both allosteric and orthosteric 
 pockets .
( upgma ) , a hierarchical 
 agglomerative approach that groups residue sets according to the degree 
 of similarity in their correlation patterns .
the authors ultimately 
 identified residue sets likely to mediate communication between the 
 nikr allosteric ni - binding pocket and the orthosteric 
 dna - binding pocket .
others have pursued orthogonal methods to 
 represent and analyze protein dynamics in silico . in 2009 , mcclendon 
 et al .
introduced mutinf , a novel md - analysis technique for identifying 
 statistically significant correlated motions . unlike other techniques
, mutinf uses an internal coordinate 
 scheme . rather than considering the positional coordinates of residue 
 atoms directly ,
it quantifies correlated motions using a mutual - information 
 metric derived from residue torsion angles .
this approach makes mutinf 
 particularly well - suited to study allosteric signals transmitted through 
 the gearlike correlated motions of side - chain rotamers . using multiple 
 short simulations ,
the authors were able to identify statistically 
 significant correlations . applying the technique to interleukin-2 
 provided a viable theory of allosteric transmission involving a population - shift - mediated 
 signal transmitted between the allosteric and orthosteric pockets 
 via a hydrophobic core .
leibler 
 divergence metric to improve the sensitivity of mutinf analysis . a paper by van wart et al .
, published 
 in 2012 , studied imidazole 
 glycerol phosphate synthase ( igps ) , an important enzyme in the histidine 
 biosynthetic pathway of plants , fungi , and microbes . following a molecular dynamics simulation of the enzyme
, 
 they built a representative network of residue interactions as inferred 
 by correlated motions .
the edges connecting physically distant protein 
 residues were ignored to emphasize correlations resulting from immediate 
 physical interactions .
the single shortest path ( in network space ) 
 between the allosteric and orthosteric sites was then identified . 

two of the residues along this path had been shown previously by experiment 
 to be involved in allostery .
van wart et al . later expanded 
 on their work by creating the freely 
 available weighted implementation of suboptimal paths ( wisp ) program .
wisp can identify not only the single optimal 
 path connecting allosteric and orthosteric sites , but also multiple 
 other near - optimal paths . while allosteric signaling may occur through 
 a single optimal path in some cases , for many proteins it is likely 
 the summed ( or perhaps synergistic ) effect of signaling through multiple 
 near - optimal pathways . while several efficient algorithms exist 
 to calculate the single 
 optimal path between two nodes in a network space ( e.g. , the floyd 
warshall 
 algorithm , dijkstra s algorithm , etc . ) , calculating near - optimal 
 pathways is far more computation intensive .
first , 
 the edges between physically distant nodes are removed regardless 
 of the interdependence of their motions , as in van wart s original 
 paper .
second , all nodes that can not possibly participate in optimal 
 or near - optimal pathways of allosteric communication are deleted . 
 for each node
, dijkstra s algorithm is used to calculate the 
 shortest path connecting the source node , the node being considered 
 ( ni ) , and the sink node . 

if this shortest path is still too long ( in network space ) to be an 
 optimal or near - optimal pathway , no other path passing through ni can be optimal or near - optimal . 
 therefore ,
ni is deleted 
 from the network . after thus greatly simplifying the network , it becomes 
 computationally tractable to calculate notable pathways of allosteric 
 communication between source and sink residues using a recursive , 
 bidirectional approach . in 2014 , levine et al .
introduced a 
 novel analysis framework , called 
 n - body information theory ( nbit ) , that uses multivariate mutual information 
 to measure residue interdependence .
this technique is unique in that 
 it considers not only the interdependence ( mutual information ) between 
 two residues , but also the mutual information between groups of up to n residues .
specifically , a metric called 
  coordination information  is used to measure the degree 
 to which the correlated motions of a set of residues are shared with 
 another residue not in the set ( i.e. , the  contribution of 
 a site to all possible n - body correlations with another 
 site  ) .
the authors use this technique to analyze md simulations 
 of the bacterial transporter leut and identify many of the same residues 
 previously found through experiment .
residue 
 interdependence to consider how highly interconnected ( and often rigid ) communities of residues interact with one another . in this 
 view
, allostery is achieved when largely isolated communities of highly 
 interconnected residue nodes (  small - world networks  ) 
 engage in long - range interactions .
studied the allosteric mechanism of prfar binding to the igps 
 heterodimer . after simulating the system 
 in both the absence and presence of the allosteric modulator , they 
 calculated the mutual information between the positional vectors of 
 all residue  residue pairs . ignoring the correlations between 
 nodes that were not frequently physically adjacent ( 5.0  
 for at least 75% of the simulation ) , the authors next used the girvan  newman 
 algorithm to identify communities of highly interconnected nodes . 

this analysis provides a powerful illustration of how an allosteric 
 modulator might manipulate entire protein  modules  
 to transmit a signal .
supported by nmr experiments , the authors ultimately 
 concluded that prfar binding decreases the correlated motions at sidel 
 of the heterodimer and induces conformational changes in sider that 
 ultimately impact the hydrogen - bond network near the active site .
sethi et al . used a similar technique in 2008 
 to study allostery in trna : protein complexes .
the use of accelerated md ( amd ) simulations in this 
 same mutual - information / community - analysis workflow has also proved 
 fruitful .
amd simulations are particularly helpful when the allosteric 
 mechanism being studied occurs on time scales that are longer than 
 those accessible through traditional md . for example , in 2012 , gasper 
 et al . used a mutual - information / community - analysis technique to study 
 -thrombin .
one occurs on slower time scales that might not have been accessible 
 were it not for the amd .
others have used amd to similarly study allosteric 
 effects in the m2 muscarinic receptor . in 2014 ,
blacklock and verkhivker used a similar technique 
 to analyze 
 md trajectories of the hsp90 chaperone in various functional states . to verify the existence of communities in this 
 system , they first calculated force constants for each residue , corresponding 
 to the energy cost of residue displacement over the course of the 
 equilibrium simulation . as residues with high force constants are 
 typically rigid , community boundaries ( e.g. , the boundary between 
 a rigid module and a flexible interdomain hinge site ) often correspond 
 to locations where these constants change abruptly .
in fact , the  hot 
 spot  regions did correspond to residues known experimentally 
 to be important hinge sites or sites of dimerization .
to further 
 characterize this community - based organization , the 
 authors next used protein - structure network analysis to build a network 
 representation of hsp90 dynamics .
they calculated cross - correlation 
 matrices from the md simulations and identified residues with high 
 degrees of  centrality  ( the number of interacting residues ) 
 and  betweenness  ( the frequency with which a given 
 node lies along the shortest path in network space between two other 
 nodes ) . based on this analysis , the authors concluded that hsp allosteric 
 communication is mediated by highly stable central nodes , surrounded 
 by residues with high degrees of betweenness that may  shield  
 the more critical residues from the random brownian / thermal motions 
 common to all proteins .
one 
 straightforward approach is to consider any identified pocket other 
 than the orthosteric to be allosteric .
however , coupling 
 these algorithms with techniques that account for pocket dynamics 
 is often useful .
some allosteric pocket conformations are extremely 
 rare in the absence of the allosteric effector itself .
cryptic 
 pockets  are not readily evident in static apo- or ligand - bound 
 crystal structures , but they may be predictively sampled with molecular 
 dynamics simulations .
computational methods for pocket detection can be classified as 
 geometry , knowledge , and energy based .
we describe these three types 
 of methods in sections 3.3.2 , 3.3.3 , and 3.3.4 .
geometry - based methods identify pockets by considering only the positions 
 of the receptor atoms themselves .
these methods are well - suited to 
 situations where speed of computation is critical ( e.g. , when pocket 
 detection is applied to conformations extracted from entire md trajectories ) 
 or where pockets are particularly well - defined .
on the other hand , 
 accuracy in pocket detection suffers when pockets are shallow , partially 
 collapsed , or highly flexible . and
pocket ranking by simple geometric 
 metrics ( e.g. , volume ) does not always correlate strongly with ligand 
 binding affinities or druggability because geometry - based methods 
 do not consider the chemical properties of a given pocket .
they can be classified into three subcategories : sphere 
 based , -shape based , and grid 
 based .
sphere - based methods , such as pass , phecom , pocasa , and surfnet , first cover 
 protein surfaces with spheres .
each sphere is evaluated and eliminated 
 if judged unlikely to occupy a surface cavity .
in contrast , 
 methods such as apropos , cast , castp , and fpocket use -shape schemes 
 to identify protein pockets . in brief ,
an -shape is like a 
 convex hull , except that it potentially follows the surface of the 
 source points more closely than a convex hull .
the degree to which 
 it deviates from the convex hull is controlled by the  value 
 ( for more details , see edelsbrunner et al . ) .
apropos compares multiple -shape representations of a protein 
 at different values of . pockets 
 are identified by subtracting  higher - resolution 
in contrast , cast generates 
 a voronoi diagram based on protein - atom locations , removes any voronoi 
 edges and vertices that fall completely outside the receptor , and 
 identifies pockets as collections of empty triangles ( for details , 
 see ref ( 211 ) ) .
fpocket 
 uses a voronoi tessellation to calculate receptor -spheres 
 ( i.e. , spheres that include four protein atoms on their surface but 
 no protein atoms in their interior ) . when clustered , these spheres tend to congregate at potential pocket 
 sites .
grid - based geometric methods for pocket detection include 
 ghecom , ligsite , ligsite , pocket , pocket - picker , pocketfinder , q - sitefinder , and povme / povme pocket i d .
these techniques consider 
 points spaced along a grid that encompasses the protein receptor . 

each point is evaluated for the likelihood of pocket occupancy , and 
 high - scoring points are clustered to identify binding pockets , whether 
 orthosteric or allosteric .
they include the presence of steric 
 clashes with the protein , detected protein  solvent  protein 
 events , buriedness indices , predicted methyl 
probe interaction energies , etc . as a representative example of the geometry - based class , 

consider 
 povme / povme pocket i d , a grid - based pocket - finding algorithm .
povme pocket i d identifies surface and internal cavities by first 
 creating a low - resolution 3d grid of equidistant points that encompasses 
 the entire protein .
points are removed if they are physically near 
 any protein heavy atom or fall outside the convex hull defined by 
 the protein -carbons .
the volumes corresponding to the remaining 
 low - resolution points , which tend to encompass protein pockets , are 
 then replaced with smaller , higher - resolution 3d grids , and the same 
 point - eliminating process is repeated .
finally , any remaining high - resolution 
 points that have fewer than a user - defined number of neighbors are 
 eliminated iteratively until no such points remain , and stretches 
 of contiguous points are grouped together into distinct pockets .
the locations of these identified pockets can be fed into the povme 
 program , which calculates pocket metrics such as volume and shape . 

the latest version of povme can efficiently analyze entire md trajectories , 
 allowing the user to monitor pocket changes over time and identify 
 sampled conformations that are geometrically distinct .
this approach 
 is useful for identifying cryptic pockets that only occasionally manifest 
 themselves over the course of an md simulation .
knowledge - based algorithms 
 draw on large structural and genomic databases 
 to determine binding - pocket locations based on sequence and overall 
 structure rather than pocket shape .
algorithms such as 3dligandsite , findsite , and 
 pocketome infer these locations by querying 
 existing databases for ligand - bound proteins that are structurally 
 similar to the target protein .
these algorithms work particularly 
 well when the target protein has a highly conserved binding site . 

other methods , such as frpred , use sequence conservation , in conjunction 
 with supporting methods like surface - residue - prediction algorithms , 
 to identify surface residues that are likely biologically important . finally , some hybrid methods , such as the evolutionary - trace 
 method , conseq , and concavity , combine sequence 
 and structural homology to identify likely binding pockets .
the section single - site evolutionary analysis methods describes 
 methods like these in further detail .
energy - based 
 methods identify binding sites by evaluating whether or not a given 
 protein region interacts favorably with docked small organic probes . 

this additional information 
 does come at a cost ; energy - based methods are more computationally 
 demanding and so do not scale well , limiting their applicability when 
 analyzing large - scale structural data sets or md trajectories .
methods such as sitemap and sitehound use single probes ( typically methane or water 
 molecules ) to assess the likelihood of small - molecule / protein binding . 

in contrast , methods such as ftmap , grid , and mcss consider multiple chemically diverse probes 
 and so identify not only pockets , but also druggable hot spots that 
 are more likely to function as orthosteric or allosteric binding sites .
ftmap is the computational equivalent of the multiple 
 solvent crystal structures ( mscs ) technique . rather than superimposing multiple crystal structures obtained by 
 soaking a protein in six to eight different organic solvents , ftmap 
 computationally docks 16 distinct organic probes into a user - provided 
 protein model .
the docked positions of these probes are then clustered , 
 and the clusters are ranked by their boltzmann - averaged energies . 
 in an alternate implementation called ftsite ,
adjacent , mutually accessible clusters are combined into 
  sites  that are then ranked by the number of nonbonded 
 contacts between the protein and the corresponding probes .
both these 
 techniques can be applied to multiple structures , whether derived 
 from x - ray crystallography or simulation , to potentially identify 
 cryptic sites that are not always evident when single structures are 
 considered .
for example , ivetac et al . used it to detect five potentially druggable 
 sites on 1ar and 2ar .
it has also been used to evaluate the druggability of p53 cryptic 
 pockets , which were subsequently used to prospectively discover novel 
 reactivation compounds . over the past 10 years , 
 markov state models ( msms )
have been used increasingly to analyze 
 molecular dynamics simulations , suggesting that they are already promising 
 tools to study allostery in a drug - discovery context . in brief , a markov state model is a stochastic kinetic model that 
 describes the probability of transitioning between discrete states 
 at a fixed time interval .
these models are required to have the markovian property ( i.e. , 
 that the probability of transitioning between discrete states is independent 
 of previous transitions ) . by clustering protein structures extracted 
 from an md trajectory , it is possible to identify discrete conformational 
 states for use in msms .
transitions between conformational clusters observed over the course 
 of an md trajectory are tallied , and the msm is then built from the 
 transition probabilities between these distinct clusters .
as they 
 are built solely on these transition probabilities , msms can draw 
 upon multiple trajectories , enabling efficient sampling of the entire 
 conformational space through many independent simulations carried 
 out in parallel .
msms can be built with the software package msmbuilder , 
 available at http://msmbuilder.org / latest/. emma , available at http://www.emma-project.org/latest/ , is another useful building tool .
because independent trajectories 
 are aggregated as a postprocessing step , the simulations can be carried 
 out on a variety of computational resources , including independent 
 desktop machines , local clusters , or high - performance computing resources . 

used tens of thousands 
 of trajectories carried out on google s exascale cloud computing 
 resources . by subsequently  knitting together  the simulation 
 data within a msm framework , the researchers elucidated the activation 
 pathway of gpcr 2ar .
msm / md analysis provides access to the thermodynamic , kinetic , 
 and structural characteristics of the protein conformational ensemble 
 ( i.e. , a robust description of the free - energy landscape of the protein ) .
the thermodynamics of the various conformational 
 states can be calculated from the equilibrium distribution .
it is 
 also possible to resolve the transition kinetics between individual 
 states , the concerted or principal protein motions , metastable states , 
 and the transition pathways between discrete states .
lastly , the md trajectories 
 also provide representative receptor snapshots for use in structure - based 
 drug design .
md - based 
 msms were originally developed to study protein folding , but they can also be used to study molecular phenomena that are 
 directly applicable to allostery and drug discovery , including cryptic 
 site identification , protein  ligand interactions , and protein 
 function .
msms have been used to 
 study binding - site conformational heterogeneity , revealing conformations 
 that are not readily evident in any crystal structure . in 2012 ,
bowman 
 and co - workers used msms to identify and validate cryptic allosteric 
 sites in tem-1 -lactamase , interleukin-2 , and rnase h. similarly , in a recent study of gpcr conformational transitions , 
 kohlhoff et al . identified intermediate conformations that they used 
 successfully in a subsequent virtual - screening campaign .
shukla and co - workers also identified a key 
 conformational intermediate in the src kinase activation pathway that 
 could be used as a target for drug discovery . while traditional molecular dynamics simulations have been 
 used to identify cryptic allosteric sites and elucidate binding - site 
 conformational diversity , markov state models can characterize the 
 probability and dynamics of these conformations , permitting rational 
 design against dynamic pockets .
protein 
 binding that are critical to both allosteric and orthosteric structure - based 
 drug discovery .
used msms derived from multiple independent 
 md simulations to accurately determine the kon and binding free energy of benzamidine to trypsin .
used msms 
 in conjunction with metadynamics , a variant of molecular dynamics , 
 to determine the koff for the same system .
while this approach is currently too computationally 
 demanding for virtual screening , it has identified transitory sites 
 important for binding and shed new light on the interactions that 
 govern ligand - binding kinetics .
msms also provide useful insights 
 into how a protein conformational 
 ensemble begets function and transmits allosteric signals .
most studies 
 have focused on the transition between two functional states . in this 
 context , msms have been used to study src kinase , the 2 adrenergic receptor ( 2ar ) , and bacterial response 
 regulator ntrc .
recently used 
 this same technique to study allostery in the cyclic - nucleotide binding 
 domain ( cbd ) of protein kinase a. by 
 building and simulating models of the cbd , with and without bound 
 camp , the study provided insight into how ligand binding modulates 
 the protein conformational ensemble . as it is challenging to 
 determine the impact of an effector molecule 
 on the conformational ensemble , most msm studies to date have not 
 examined allostery directly .
this is best shown in our study of the pka cbd , where 
 the effector molecule camp significantly modulates the kinetics of 
 the transition from the active to the inactive states of the protein , 
 while simultaneously leaving the reverse transition between the inactive 
 and active states unmodified .
this modulation 
 seems to be regulated by transient interactions between the effector 
 molecule and the protein .
second , recent msm studies suggest that 
 there is not one but many pathways between functional states . studies 
 of 2ar , ntrc , and the pka cbd each revealed multiple 
 transition pathways that arose from transient interactions within 
 the system , suggesting a probabilistic model of allostery in addition 
 to one that describes allostery in terms of concerted pathways
. studies like these can provide useful insights 
 applicable to structure - based drug discovery targeting transitional 
 conformational states .
ligand 
 interactions and allostery , but building these models is computationally 
 demanding , limiting their use as screening tools .
this weakness aside , 
 msms do provide a strong foundation for rational drug discovery and 
 enhanced understanding of allosteric mechanisms .
like msms , energy - landscape 
 theories were originally developed to study protein folding .
but they 
 too can be applied to the study of allosteric mechanisms and protein 
 dynamical modes . per statistical mechanics ,
the underlying free - energy 
 landscape of a system determines the probability of observing a particular 
 configuration .
any topological analysis aimed at understanding allostery 
 must be based on a thermodynamic understanding of the impact that 
 effector binding has on this energy landscape .
it is helpful 
 to view allosteric binding as an event that perturbs the conformational 
 exchange , or folding pathways , of a protein .
the principles of statistical 
 mechanics applied to polymers in solution have provided much insight 
 into the protein - folding process .
as proteins fold , the number of 
 native contacts increases , and the entropy of the polymer decreases . 

folding is impossible when the required conformational transitions 
 are occluded by large free - energy barriers or wells .
for example , 
 a large local minimum along the folding pathway could force a protein 
 to undergo a glass transition , leading to a kinetically trapped state . 

( i.e. , the folding pathway must be funnel shaped 
 and relatively smooth ) . in terms of allostery and energy - landscape theory ,
a positive 
 allosteric modulator can be seen as a compound that increases the 
 likelihood of the binding state by deepening the binding - state well 
 or decreasing the depth of competing wells . on the other hand
, a negative 
 allosteric modulator may shift the location of the binding - site well , 
 shallow out the well , or increase the depth of competing wells .
using 
 these guiding principles , heuristics can be designed to identify each 
 of the above scenarios .
have developed a heuristic 
 to pick out locally 
 frustrated residues . while the overall folding process 
 obeys the principle of minimum frustration
, this does not preclude 
 the possibility of having small local barriers that must be overcome 
 to fold . to detect these local frustrations ,
the authors systematically 
 perturbed the amino acid sequence and compared the perturbed total 
 energy with that of the native state .
this comparison provides the 
 effective stabilization energy associated with each amino acid pair . 

if a given native - state interaction is highly stabilizing ( i.e. , the 
 interaction energy is favorable compared to mutants ) , it is said to 
 be  minimally frustrated .  on the other hand , if a pair 
 of residues is sufficiently destabilizing ( i.e. , the interaction energy 
 is weak compared with that of the mutants ) , it is  frustrated .  

frustrated residues can be thought of as weak points in the global 
 structure , contributing to barriers that may introduce slow degrees 
 of freedom . as such
jenik et al . have implemented this heuristic as a web server http://www.frustratometer.tk/. weinkam et al .
used a dual - topology 
 go model to survey the free - energy 
 landscape of the 
 calmodulin - gfp ca sensor protein , the maltose binding 
 protein , and the csl transcription factor .
they introduced a truncated 
 gaussian distance term to the soft - sphere atom overlap term implemented 
 in modeller , resulting in either one basin corresponding to the bound / unbound 
 distances or two corresponding to both , depending on the pairwise 
 distance and cutoff .
subsequent md simulations of the resulting model 
 were analyzed using their pseudocorrelation map algorithm .
they found that their model reproduces allosteric 
 motions and sheds insight into the macroscopic mechanism .
protein dynamics 
 and subsequent perturbations from binding events can also be studied 
 using normal - mode analysis ( nma ) .
nma assumes that the 
 potential energy landscape in the vicinity of a minimized atomic structure 
 is approximately harmonic .
this simplifying assumption allows for 
 diagonalization of the hessian matrix . solving the eigenvalue problem 
 produces the eigenvectors ( movement direction ) and eigenvalues ( vibrational 
 frequencies ) .
nma is known to reproduce the slow degrees of freedom 
 of protein motion well . applying additional 
 constraints to simulate
the impact of an effector on binding - site 
 residues can yield new fundamental modes that may provide insight 
 into the allosteric mechanism .
ming and wall applied normal - mode 
 analysis to investigate the effect of tri - n - acetyl - d - glucosamine upon lysozyme .
they 
 assumed that greater perturbation to the conformational distribution 
 corresponds to increased likelihood that a binding event will occur 
 at a regulatory or active site .
the degree of perturbation is monitored 
 by comparing the distributions predicted by nma of the protein with 
 and without the ligand bound , via the kullback 
leibler divergence . 
 in this work , ming and wall generated many arbitrary poses and investigated 
 the allosteric potential , defined by the kld .
the use of normal - mode analysis removes the need to integrate the 
 equations of motion .
however , for the harmonic approximations to be 
 accurate , the starting structure needs to be in a local minima .
an 
 initial molecular mechanics minimization can ensure that proteins 
 of interest are near what might be seen in their native environments 
 ( e.g. , a gpcr embedded in a lipid bilayer with heterogeneous composition ) . 

however , alternate options include simplified interaction hamiltonians 
 such as those offered by elastic network models ( enms ) . the complexity 
 and computational cost of simulating detailed atomic potentials inspired 
 development of simpler enms in which uniform harmonic potentials are 
 used to model all interactions .
a popular enm is the 
 gaussian network model ( gnm ) , wherein neighboring residues are connected 
 by virtual springs to create a network / graph of interactions . here 

the interatomic potential energy , u , of the system 
 can be expressed as1where  is the uniform 
 force constant of the harmonic springs , r is 
 the n - dimensional column vector of the instantaneous 
 fluctuations of specific atoms , and  is the n  n kirchoff ( connectivity ) matrix of residue 
 interactions given by2where rc is a residue cutoff distance .
ming and wall have developed 
 a backbone - enhanced elastic network model ( benm ) in which the interactions 
 between connected residues are scaled by an additional factor . using 
 an objective function defined by the kld of the marginal distribution 
 of all - atom and c atomic coordinates
, the authors 
 showed that benm can reproduce the atomistic mean - squared displacement 
 for bovine pancreatic trypsin inhibitor ( bpti ) . to consider the effects of allostery
, the authors again 
 used the kld of local conformational distributions before and after 
 ligand binding . to locate allosteric 
 binding sites , su et al .
all neighboring residues ( represented by their 
 corresponding c atoms ) were connected by virtual 
 harmonic springs of equivalent strength .
a spherical probe then sampled 
 the protein surface to locate the potential effector interaction points . 

additional springs were attached between the effector point and the 
 local residues in the protein  effector model . the free - energy 
 difference ( g ) between the protein 
regions with large g values were assumed 
 to be important for the allosteric mechanism . using this method , su 
 et al . studied the hsp70 and the glua2 ampa receptors .
used a similar coarse - grained network - interaction 
 model to predict the location of allosteric sites . to identify potential pockets for effector / substrate binding , 
 the authors generated a 
residue interaction graph  
 ( rig ) , with residues as nodes and interatomic contacts as edges .
local closeness ,  defined as the 
 number of total neighboring residues weighted by 1/r , where r is the inter - residue distance . regions with high local closeness were likely 
 to be in pockets and , thus , were putative binding sites .
the extent 
 to which each binding site was coupled to the intrinsic protein motions 
 was then explored by estimating the strain ( so - called  binding 
 leverage  ) on ligand 
a ligand that bound to a site with high binding 
 leverage could restrict protein dynamics .
the binding leverage describes the extent to which changes 
 in binding 
 sites and protein conformational degrees of freedom are coupled ; it 
 does not specifically predict allosterically linked sites . to predict 
 allostery , mitternacht and berezovsky introduced  leverage 
 coupling ,  which selects for sites that are strongly coupled 
 to the same degree of freedom .
the mode 
 associated with the most binding leverage is considered most relevant , 
 based on the assumption that there is one dominant motion promoting 
 allostery .
over the past decade , advances in computing power and predictive 
 algorithms coupled with the increased availability of structural and 
 biochemical data have revealed new opportunities for rational design 
 of allosteric drugs .
the emergence of novel computational approaches 
 to describe and predict allosteric phenomena across a range of scales , 
 from the coordinated atomic movements in a single receptor molecule 
 to complex allosteric signaling networks , is ushering in a new era 
 wherein computational methods can be used to prospectively predict , 
 discover , and characterize allosteric sites and effector molecules . 
 within the context of a drug - discovery program , such approaches hold 
 the potential for developing drugs with increased specificity and 
 selectivity , as well as the ability to gain new and more comprehensive 
 understanding of old targets .
the convergence of advances in ( i ) theoretical 
 msm - based frameworks and msm building software , ( ii ) community md 
 codes that can achieve > 100 ns / day sampling for realistic sized 
 systems 
 on single gaming / commodity gpu processors , and ( iii ) pocket and druggable 
 site - detection algorithms now make it possible for researchers even 
 in industrial settings , with fast - paced timelines and stringent quality 
 standards , to apply these approaches to drug targets already in their 
 arsenals .
the application of these methods to kinases and gpcrs seems 
 particularly worthwhile , given the existence of assays and structural 
 data , and the challenges faced by existing drug candidates in the 
 clinic .
 allosteric drug development holds 
 promise for delivering medicines 
 that are more selective and less toxic than those that target orthosteric 
 sites . to date , the discovery of allosteric binding sites and lead 
 compounds has been mostly serendipitous , achieved through high - throughput 
 screening . over the past decade , 
 structural data has become more readily 
 available for larger protein systems and more membrane protein classes 
 ( e.g. , gpcrs and ion channels ) , which are common allosteric drug targets . 
 in parallel , 
 improved simulation methods now provide better atomistic 
 understanding of the protein dynamics and cooperative motions that 
 are critical to allosteric mechanisms . as a result of these advances , 
 
 the field of predictive allosteric drug development is now on the 
 cusp of a new era of rational structure - based computational methods . 
 here 
 , we review algorithms that predict allosteric sites based on 
 sequence data and molecular dynamics simulations , describe tools that 
 assess the druggability of these pockets , and discuss how markov state 
 models and topology analyses provide insight into the relationship 
 between protein dynamics and allosteric drug binding . in each section , 
 we first provide an overview of the various method classes before 
 describing relevant algorithms and software packages . 
Review
Motivation and Organization
Introduction
Computational Methods for
Studying Allostery
Energy Landscape and Topological Analyses
Conclusions