allosteric regulation is a key mechanism 
 whereby proteins respond 
 to environmental stimuli that modulate their activity .
classic models of allostery ( e.g. , the mwc and knf models ) suggest that a binding 
 event at an allosteric site induces substantial conformational changes 
 in the primary catalytic site .
however , allostery has since been observed 
 in the absence of large - scale conformational changes , suggesting that subtle alterations in protein dynamics can induce 
 a population shift in the conformational ensemble without substantially 
 changing the mean conformation of the protein .
recent advances in both correlated - residue clustering and dynamical 
 network analysis have helped computationally quantify allosteric states .
dynamical network models of allostery often focus on the single 
 most direct path of residues leading from the allosteric to the primary 
 active site .
however , few researchers have considered the state changes 
 of slightly longer ( suboptimal ) allosteric pathways
. the statistical 
 distribution of these additional pathways may be useful for locating 
 accessible residues that , if disrupted via pharmacological or mutational 
 means , could modulate the allosteric regulation of important drug 
 targets . in this paper
, we introduce weighted implementation 
 of suboptimal 
 paths ( wisp ) , a tool that compliments current dynamical network models 
 of allostery by rapidly calculating the primary communicating path 
 between two residues as well as the slightly longer suboptimal paths . 

we illustrate the utility of the wisp method using the biological 
 system hish - hisf , a well - characterized glutamine amidotransferase 
 enzyme . to facilitate the broader adoption 
 of this method , we have also created a wisp plugin for the popular 
 visual molecular dynamics ( vmd ) package .
wisp has been specifically tested on several operating systems , 
 using several versions of python , numpy , scipy , and networkx ( table 1 ) .
wisp has been 
 tested on a number 
 of operating systems , using various versions of numpy , scipy , and 
 networkx .
we note that installation of necessary packages under windows 
 was difficult ; however , the command - line version of the program was 
 successfully executed after installing the appropriate dependencies 
 using the activepython software package .
as input , wisp 
 accepts an aligned molecular dynamics trajectory in the common multiframe 
 pdb format .
trajectory postprocessing 
 is necessary prior to wisp analysis , as most trajectories are not 
 initially aligned or pdb formatted .
the freely available visual molecular 
 dynamics ( vmd ) software package can be used to perform the necessary 
 alignment and conversion .
wisp , similar to 
 other dynamical network analysis tools , is based on the dynamic interdependence among protein constituents 
 ( e.g. , amino acids ) .
a protein system is first simplified by representing 
 each constituent as a single node .
for example , depending on user - specified 
 wisp parameters , an amino acid can be represented by a node positioned 
 at the residue center of mass , the side - chain center of mass , the 
 backbone center of mass , or the  carbon . as a default
the interdependence among nodes 
 is represented as a connecting edge with an associated numeric value 
 that reflects its strength .
there are numerous methods for describing 
 the interdependence among nodes in a protein network .
typically , this 
 interdependence is represented by a matrix c with 
 values corresponding to the weights of each edge . by default
, wisp 
 generates an n matrix c by calculating the correlated motion among node  node pairs 
 as shown in eqs 1 and 2:12where n is the number 
 of 
 nodes , i and j are indices corresponding 
 to individual nodes , ri(t ) is the location of node i at 
 time t , and cij is the matrix element at position ( i , j ) . the absolute value of cij is larger when the motions of two nodes are highly 
 correlated 
 or anticorrelated . in order to compute signaling pathways , it is useful 
 to construct a matrix where the opposite is true , i.e. , where small 
 values indicate highly correlated or anticorrelated motions .
consequently , 
 the correlation matrix is functionalized according to eq 3 , as outlined in previous works.3 as a point of clarification , each wij can be thought of as a  distance  
 in functionalized correlation space . throughout the remainder of this 

paper , concepts like length and distance will refer to spans in this 
 space , unless specifically described as
 we further note that , while wisp s 
 default functionalized correlation matrix is generally useful , any 
 user - specified matrix that defines signaling strength as inversely 
 proportional to edge length can be used . in order to improve 
 the speed of subsequent path - finding steps , 
 the complexity of the functionalized correlation matrix w must be reduced . to this end
first , a 
 contact - map matrix mcontact is used to 
 separate entries in w that represent pairs of physically 
 distant residues from those that represent adjacent residues . by default
, mcontact is constructed using pcutoff , a user - specified cartesian cutoff distance that 
 represents physical proximity .
the average location of each 
 atom over the course of the aligned molecular dynamics trajectory 
 is first calculated , followed by a pairwise cartesian distance comparison . 

two nodes are considered to be in physical contact if the average 
 locations of any of their associated residue atoms come within pcutoff of one another .
mcontact entries are set to zero for all node  node pairs 
 that are not in physical contact .
a simplified , functionalized correlation 
 matrix wsimp is then constructed by multiplying w and mcontact element - wise . 

the entries of wsimp that equal zero represent 
 node  node interactions that are subsequently ignored .
second , to further reduce the complexity of the functionalized 
 correlation matrix w , a pruning algorithm identifies 
 nodes that only participate in pathways having lengths in network 
 space that are greater than another cutoff ( dcutoff ) . as the ultimate goal is to identify suboptimal paths 
 with lengths less than dcutoff
, these 
 nodes can be effectively discarded as well . to identify these nodes , 

an fnp 
 is the optimal pathway between two user specified nodes na and nb that is forced to pass through a given third node ni .
for any two fixed nodes na and nb , each third node ni is associated with a single fnp .
the set 
 of all fnps can therefore be generated by iterating over all the nodes , ni , of the system .
to 
 calculate an fnp , dijkstra s algorithm , included in networkx , is first used to identify the optimal paths 
 between na  ni and nb  ni , respectively .
the fnp has a length equal to the 
 sum of these two constituent paths .
any path between na and nb that passes through ni must have a length equal to or greater than that 
 of the associated fnp .
consequently , if the length of the fnp is greater 
 than dcutoff , all entries in wsimp associated with ni are set to zero , so that ni is effectively ignored .
having generated wsimp , we are now ready to search for both the 
 single optimal and multiple suboptimal paths between na and nb .
the optimal path is fairly easy to identify using 
 dijkstra s algorithm , mentioned above .
in contrast , identifying 
 all suboptimal paths is difficult because the number of possible pathways 
 between na and nb grows rapidly as the total 
 number of nodes increases . to identify suboptimal paths , a recursive
simultaneous searches start from na and nb ( figure 1 , in blue 
 and red , respectively ) and recursively traverse the nodes of the dynamical 
 network .
the recursive algorithm ignores the connections / correlations 
 between nodes that are physically distant ( figure 1 , gray lines ) .
additionally , nodes eliminated using the fnp 
 technique described above are likewise ignored ( figure 1 , gray circles ) , resulting in substantial speedups . as soon 
 as any of the lengthening paths grows longer than dcutoff ,
that branch of the recursion is killed ( figure 1 , red  x  ) .
simultaneous searches start 
 from na and nb ( blue and red , respectively ) and recursively 
 traverse the nodes of the dynamical network .
connections / correlations 
 between nodes that are physically distant are ignored ( gray lines ) . 

nodes eliminated using the fnp technique are also ignored ( gray circles ) . 
 as soon as any of the lengthening paths grows too long , that branch 
 of the recursion is killed ( red  x  ) . at each recursive 
 step , all branches originating from na and nb are compared for common nodes ( asterisk ) . if a common node exists , 
 the two paths are joined . if the length of this composite path is 
 sufficiently short , a suboptimal path has been identified . at each recursive step , all branches originating 
 from na and nb are compared for common nodes ( figure 1 , the node marked with an asterisk ) . if a common 
 node exists , 
 the two paths are joined at this node . if the length of this composite 
 path is less than dcutoff , a suboptimal 
 path has been identified .
as wisp has been developed to take advantage 
 of multiple processors , running the program on a multicore system 
 can lead to further speedups beyond the software optimizations described 
 above .
the program output is a directory containing 
 multiple files , including the specific w and mcontact matrices used .
the primary output file 
 is a tcl script that , when loaded into vmd , draws three - dimensional 
 splines representative of the optimal and suboptimal paths .
user defined 
 parameters control the relationship between spline thickness , color , 
 opacity , and path length .
useful information is also given as comments 
 in the tcl file , including path lengths and participating protein 
 residues .
in addition to the command - line 
 program , we have also developed a visual molecular dynamics ( vmd ) plugin and tcl - based gui for easy preparation 
 and visualization of wisp results .
the main window of the wisp 
 gui ( figure 2 ) allows the user to specify the 
 molecular trajectory as well as the allosteric - signal source and sink 
 residues .
several additional window interfaces allow the user to modify 
 more advanced program options if needed .
all options available through 
 the wisp command - line interface are available to users of the gui .
the 
 gui is used to visualize the allosteric pathways between leu50:hisf 
 and glu180:hish . in the main window
( top left ) , the user selects the 
 relevant molecule and which residues to use as the source and sink . 

the user may also select to load the visualization into vmd upon job 
 completion .
the setting option windows ( left and bottom right ) allow 
 the user to specify additional wisp arguments .
once satisfied with the run specifications , the user may 
 click 
 the  run wisp 
the plugin loads the visualization of the 
 allosteric pathways into the main vmd window , where the appearance 
 can be further modified according to the user s preferences .
the molecular dynamics simulations 
 of hish - hisf used in the current study have been described previously . in brief , a model of the hish 
hisf apo 
 dimer was prepared from the 1gpw crystal 
 structure ( thermotoga maritima ) . to 
 generate the corresponding holo structure , the 1ox5 crystal structure ( saccharomyces cerevisiae ) , which contains a cocrystallized prfar allosteric effector molecule , 

was aligned to the apo model , effectively positioning prfar within 
 the 1gpw : hisf allosteric site .
the aligned 1ox5 prfar was then merged 
 with the 1gpw - based apo model to yield the corresponding holo structure . 

following solvation with tip3p water molecules and 1 ns of harmonic 
 constrained equilibration , 20 ns of production dynamics with a 2 fs 
 time step were run for both the apo and holo systems using namd , the charmm27 force field , and the same prfar parametrization used previously .
consequently , one natural strategy for rational drug design is to 
 impede or agonize protein function via allosteric modulation .
classic 
 views of allostery suggest that the binding of an effector molecule 
 at an allosteric site induces large conformational shifts that alter 
 the activity of the primary site .
however , as allostery is not necessarily 
 limited to large shifts , this reasoning does not explain some examples 
 of regulation at a distance .
recently showed that significant backbone deformations 
 are not required for an allosteric effect ; rather , in the absence 
 of large conformational changes , subtle shifts in local dynamics driven 
 by entropic effects govern certain types 
 of allostery .
quasi - harmonic analysis ( e.g. , like that used 
 by software packages 
 such as carma to calculate entropy ) is commonly 
 used to build dynamical network models that quantify signaling pathways 
 among protein constituents .
optimal and suboptimal pathways are calculated 
 that connect protein constituents believed to be important for allostery 
 ( i.e. ,  sources  and  sinks  ) . an optimal 
 pathway is the shortest distance traversed between source and sink 
 along weighted edges ( e.g. , as determined by correlated motions ) , 
 and suboptimal pathways are those closest in length to , but not including , 
 the optimal path .
existing tools can compute optimal and suboptimal 
 pathways between residues ; however , these 
 programs lack the speed required to compute more than 50 suboptimal 
 pathways within a reasonable amount of time ( several hours or days ) . 
 as statistics related to suboptimal pathways may provide important 
 insights that can not be gleaned from the single optimal pathway , faster 
 algorithmic advances must be made .
wisp is designed to facilitate 
 the calculation of hundreds of suboptimal 
 pathways in minutes , thereby permitting fast and robust statistical 
 analysis of biological systems modeled as dynamical networks . for 
 example
, using a modern workstation with 24 cores , we recently used 
 a 20  000-frame trajectory to identify 750 pathways .
wisp loaded 
 and analyzed the trajectory , generated the functionalized correlation 
 matrix , and identified the 750 pathways in 21 min and 52 s. when the 
 calculation was repeated using a copy of the functionalized correlation 
 matrix saved from the first run , the 750 pathways were identified 
 in only 5 min and 44 s. to demonstrate the utility of the wisp 
 algorithm , we used it to 
 study hish - hisf , a multidomain globular protein known to exhibit allostery . 

the activity of hish - hisf , which regulates the fifth step of the histidine 
 biosynthetic pathway in plants , fungi , and microbes , is substantially 
 altered by the allosteric effector n1-[(5-phosphoribulosyl)-formimino]-5-aminoimidazole-4-carboxamide 
 ribonucleotide ( prfar ) .
guided by previous 
 work , we investigated the suboptimal 
 pathways between residues leu50:hisf and glu180:hish using 20 ns molecular 
 dynamics simulations of both apo and holo hish - hisf . a total 
 of 700 pathways ( figure 3 ) between 
 leu50:hisf and glu180:hish
were calculated using wisp s default 
 correlation ( eqs 13 ) and contact - map matrices , described in the materials 
 and methods .
had only the two optimal pathways ( apo vs holo ) 
 been considered , we would have concluded that communication between 
 the allosteric and primary site is fundamentally different in the 
 presence and absence of the prfar effector molecule ( figures 3 and 4 ) .
the optimal pathway 
 between leu50:hisf and glu180:hish in the apo state was leu50:hisf 
  phe49:hisf  phe77:hisf 
the 700 shortest paths 
 between 
 leu50:hisf and glu180:hish , shown as red splines , derived from ( a ) 
 the apo trajectory and ( b ) the holo trajectory .
wisp allows the user 
 to choose between a number of graphical settings to better visualize 
 signaling among nodes .
a histogram of 
 the 700 path lengths associated with the apo and holo trajectories 
 is shown .
the path distribution is largely shifted to the left for the holo 
 ( allosteric ) state .
this shift likely results from a more coherent 
 signal in the holo simulation , indicating a possible decrease in the 
 entropy along the pathways due to prfar binding .
however , when we considered multiple suboptimal paths , it 
 became 
 apparent that allosteric signaling may be far more intricate .
the 
 optimal path in the apo simulation is the shortest suboptimal path 
 in the holo simulation ( top 0.3% ) , and the optimal path in the holo 
 simulation is the 13th shortest suboptimal path in the apo simulation 
 ( top 2.0% ) . in light of this multipathway analysis , the idea that 
 prfar binding fundamentally alters a solitary line of communication 
 between the allosteric and primary site becomes less tenable .
rather , 
 the binding of the effector molecule likely has small effects on multiple 
 pathways , both optimal and suboptimal , that when taken together yield 
 a substantial allosteric effect .
the lengths of the two optimal pathways of the two systems 
 did not differ substantially ( apo , 2.97 ; holo , 2.84 ) .
consequently , 
 had only these two pathways been considered , some might have mistakenly 
 concluded that the allosteric consequences of prfar binding are minor . 

in contrast , when hundreds of suboptimal paths were also considered , 
 a large prfar - dependent shift in communication between the allosteric 
 and primary site became apparent . to demonstrate this shift , we generated 
 a histogram of all path lengths for both the holo and apo simulations 
 ( figure 4 ) .
the distribution derived from the 
 holo trajectory is substantially skewed toward shorter path lengths , 
 suggesting that the motions of the residues connecting the allosteric 
 and primary sites are more tightly correlated when prfar is bound . 

an overall  dynamical tightening  and loss of entropy 
 along the pathways may therefore explain the allosteric signal .
to identify protein residues critical for allosteric transmission , 
 we counted the number of times each residue appeared in any of the 
 700 paths associated with the apo and holo trajectories , respectively 
 ( i.e. , the degeneracy of each node , figure 5 ) .
notably , a number 
 of residues had large effector - molecule - dependent shifts in degeneracy , 
 i.e. , hisf : leu47 ( shifts down ) , val69 ( shifts up ) , ala70 ( shifts 
 up ) , ile73 ( shifts up ) , asp74 ( shifts up ) , pro76 ( shifts down ) , and 
 ala97 ( shifts down ) and hish : lys181 ( slight shift down ) ( table 2 ) .
importantly , these residues , which may be crucial 
 for the regulation of protein activity , did not all appear in the 
 optimal apo and holo paths and so would not have been identified had 
 the suboptimal paths been ignored .
previous studies in evolutionary 
 conservation have shown that hisf : leu47 , val69 , ala70 , and ile73 
 are partially or strongly conserved and hisf : pro76 and ala97 and 
 hish : lys181 are strictly conserved across the entire glutamine amidotransferase 
 family of enzymes .
hisf : asp74 is not 
 conserved , but this amino acid is still predicted to play a role in 
 allostery .
compounds that target ( i.e. , 
 selectively bind ) these critical residues may serve as useful precursors 
 to future allosteric - modulating small molecules .
the total number of times 
 a given residue participates in any of the 700 paths ( i.e. , node degeneracy ) 
 is shown for ( a ) hisf and ( b ) hish .
green indicates the holo state , 
 blue indicates the apo state , and cyan indicates an overlap .
note 
 that leu50:hisf and glu180:hish are present in all 700 paths . we note that our decision to specifically 
 analyze the 700 shortest 
 paths between leu50:hisf and glu180:hish was arbitrary . in order to 
 better
assess the minimum number of paths required to reliably predict 
 node degeneracy , we analyzed the holo trajectory by varying the number 
 of paths considered and calculating the degeneracy of selected residues / nodes 
 implicated in the allosteric mechanism ( figure 6 ) .
we note that the degeneracy of these nodes had largely converged 
 by 350 paths .
a similar result was obtained when the apo simulation 
 was analyzed ( data not shown ) .
given that the relative importance 
 of suboptimal paths in determining the competency of an allosteric 
 signal is likely highly system dependent , we do not necessarily recommend 
 this exact number of paths for all analyses .
however , we are hopeful 
 that this general benchmark will help guide future researchers in 
 their efforts .
normalized node degeneracy as a function of the number 
 of suboptimal 
 paths calculated ( holo simulation ) .
the degeneracy of selected residues / nodes 
 as a function of paths searched was calculated and normalized by dividing 
 by the total number of paths .
similar results were obtained when 
 the apo simulation was analyzed ( data not shown ) .
we present wisp , a program that rapidly calculates 
 both optimal 
 and suboptimal communication pathways between distinct protein residues . 

the program is available as a vmd plugin or a standalone command - line 
 script .
wisp outputs path members and lengths that can be subsequently 
 used in the analysis of path distributions , node degeneracy , and other 
 metrics of interest to scientists studying the molecular mechanisms 
 of allostery the utility of our program was presented by performing 
 a dynamical 
 analysis of the hish - hisf protein . in our test case ,
allosteric modulation 
 was likely the result of subtle changes in multiple suboptimal pathways 
 rather than large changes in a single optimal path .
additionally , 
 we showed that prfar binding causes a large shift toward shorter path 
 lengths ( i.e. , more correlated motions ) in 700 communication pathways 
 between residues hisf : leu50 and hish : glu180 .
this shift explains the 
 strong allosteric effects of the prfar modulator ( figure 4 ) .
remarkably , the significant shift in collective 
 correlated dynamics occurred even at relatively short ( tens of nanoseconds ) 
 time scales , suggesting that the allosteric signal is rapidly transmitted . 

the multiple suboptimal pathways are dominated by a few select residues , 
 as indicated by the shift in node degeneracy between the apo and holo 
 states ( figure 5 and table 2 ) .
wisp has been successfully tested on a number of 
 platforms ( table 1 ) .
we are hopeful that the 
 program will be a useful 
 tool for the computational - biology community . a numerical representation of 
 the same data from figure 5 .
the comparison 
 between the apo and holo states suggests that certain residues are 
 more sensitive to the allosteric effector prfar than others ( shaded 
 columns ) .
 allostery 
 can occur by way of subtle cooperation among protein 
 residues ( e.g. , amino acids ) even in the absence of large conformational 
 shifts . 
 dynamical network analysis has been used to model this cooperation , 
 helping to computationally explain how binding to an allosteric site 
 can impact the behavior of a primary site many ngstroms away . 
 
 traditionally , computational efforts have focused on the most optimal 
 path of correlated motions leading from the allosteric to the primary 
 active site . 
 we present a program called weighted implementation of 
 suboptimal paths ( wisp ) capable of rapidly identifying additional 
 suboptimal pathways that may also play important roles in the transmission 
 of allosteric signals . aside from providing signal redundancy , suboptimal 
 paths traverse residues that , if disrupted through pharmacological 
 or mutational means , could modulate the allosteric regulation of important 
 drug targets . to demonstrate the utility of our program , we present 
 a case study describing the allostery of hish - hisf , an amidotransferase 
 from t. maritima thermotiga . 
 wisp and 
 its vmd - based graphical user interface ( gui ) can be downloaded from http://nbcr.ucsd.edu/wisp . 
Introduction
Materials and Methods
Results/Discussion
Conclusion