BACKGROUND
conformational changes in proteins can take place in a wide variety of ways, not all of which have been formally classified.
one important class of motions is shear, in which stacked side chains of the protein can slide without losing contact.
in this work we focus on the largest class, domain hinge bending, in which one structural domain of the protein moves relative to another domain about a hinge which connects the two  <cit> .
such motions typically involve the slowest degrees of freedom of that protein and so are difficult to predict by existing methods.
the prediction of ligand binding motions of the protein receptor has considerable potential applications in protein-protein and protein-ligand docking.
many methods can predict the side chain rearrangements required for docking  <cit>  but these assume that the large scale conformation is already nearly correct.
thus there is a need for a method that will put the receptor in the correct large scale conformation which can be a productive starting point  <cit> .
much work has been done in this area.
molecular dynamics   <cit>  explicitly computes the dynamical trajectory of molecules modeled as point masses connected by linear and nonlinear springs and can be used to predict conformational change, but usually only small- or moderate-scale domain motions can be reproduced  <cit>  with many biologically relevant motions remaining out of reach  <cit> .
accordingly several methods used md to account for the fast fluctuations of proteins in drug docking by first computing the protein trajectory using md  <cit> .
one limitation of such techniques is that they may not escape the vicinity of an initial conformation, even in a time span experimentally known to be sufficient for conformational change  <cit> .
althaus et al created a combinatorial tree of side-chain rotamers which they explored using a branch-and-cut algorithm,  <cit>  without varying the backbone conformation.
sandak et al.
created a flexible-receptor docking code which articulates the protein at a hinge point, but leaves the two resulting domains rigid  <cit> .
this method suffered from the opposite problem: it could generate large scale protein motions, but had no way of dealing with even small side chain rearrangements, a weakness leading to failure  <cit> .
the described methods are good at either treating the side-chain flexibility, or the large scale conformational changes, but not both simultaneously.
conformation explorer uses sandak et al.
's idea of moving domains about a hinge point to generate large scale conformational change, but also includes equilibration steps which permit relaxation and adjustment of all atoms.
normal modes have also been used by many authors to predict the conformational changes of proteins  <cit> .
comparison of the atomic coordinates of homologous pairs of proteins shows that the lowest order modes are most involved in conformational change,  <cit>  but also that multiple modes are needed to accurately represent the motion  <cit> .
it is possible to determine the correct combination of normal modes that will reproduce a desired motion, but this requires knowledge of at least a few interatomic distance constraints for the final structure  <cit> .
in a different approach, a docked protein-ligand complex was displaced along the lowest-frequency normal mode directions to minimize non-bonded energy terms in an md force field  <cit> .
however a normal mode expansion assumes a quadratic potential and so is accurate only for small fluctuations about an equilibrium structure; therefore the method cannot be used to predict larger scale conformational changes such as we treat in this work.
the method of lindahl et al.
gains improvements of  <dig>  to  <dig> Å for several proteins;  <cit>  our method recapitulates much larger conformational changes as we will show.
maiorov and abagyan  <cit>  rigidified all protein bonds except those in the interdomain linker and interface using internal coordinate modeling, and then used the biased probability monte carlo protocol to generate potential alternate conformations of the protein.
the method succeeded in generating a large number of alternate conformations, and some of these were somewhat similar to alternate conformations known crystallographically.
however without referring to the known alternate conformations, it was impossible to determine which of the many predicted structures was thermodynamically plausible.
further, many energy evaluations and minimizations were expended in evaluating generated conformers which were later discarded.
lastly, it was not easy to know how long a thorough exploration of conformation space would take, and no clear way to restrict the search to a given region of interest.
our method is similar in several ways to maiorov et al.
's, but also addresses these limitations.
in more recent work, de groot et al.
<cit>  showed they could find the holo conformations of several ligand-binding proteins.
the method relies on tconcoord,  <cit>  which determines flexible regions by analyzing hydrogen bonding networks.
once these are known, an ensemble of plausible structures is generated.
an interative process involving docking, md refinement, and filtering by radius of gyration then generates holo structures.
however the radius of gyration must be provided by experiments, because their fitness function does not take into account the receptor conformation.
our work resembles de groot's in the use of docking and md refinement, but differs in that the hinge location is determined by the user , and in that it requires no experimental information about the holo complex.
our method is also computationally cheaper.
key to our method is the conformational sampling and the use of a fitness function f which includes terms to discriminate the holo conformation from the generated ensemble.
the mentioned f deserves motivation.
there are numerous potentials which pick out the correct ligand-binding pose in a given receptor  <cit> .
however these typically have no terms to discriminate favorable from unfavorable receptor conformations.
other force fields have been trained and tested against ensembles of crystallographically obtained protein structures, but it is not clear how they would perform when the holo structure is unavailable  <cit> .
to overcome these limitations f includes terms  which help discriminate predicted holo complexes from a generated ensemble.
in this work we show that the conformation explorer can predictably, controllably, and economically generate multiple alternate conformations of proteins without experimental constraints  <cit>  on the holo complex.
the conformations can be equilibrated to ensure they are free of steric strain   <cit> .
we demonstrate that as an application of the method, we can predict the ligand bound, or holo, conformation of a protein given only the ligand and the coordinates of the protein in the apo state.
the exploration of phase space is performed in such a way as to minimize f, which includes terms measuring the radius of gyration rg, stability ef, and docked energy gd.
we will formally define all of these terms in the methods section.
methods
as mentioned, the method involves the identification of the hinge and rotations about that hinge, followed by equilibrations and docking runs.
the direction of the rotations and the choice of intermediate structures to which the rotations are applied, are subject to optimization.
we describe how to perform these steps.
determining and recording the hinge location
the hinge can be detected in several ways:
theoretically
hinge predictors such as flexoracle,  <cit>  tlsmd,  <cit>  hingemaster,  <cit>  and others  <cit>  can be used to predict the hinge location if at least one structure is known.
the hingemaster web server makes the use of the first three easy  <cit>  and was used in this work.
crystallographically
if the protein has been crystallized in two different conformations, these can be inspected visually to determine the hinge location.
this is usually the best and most accurate method  <cit> .
by other experimental means
proteolysis,  <cit>  hydrogen/deuterium exchange,  <cit>  optical trapping,  <cit>  fluorescent resonance energy transfer ,  <cit>  nuclear magnetic resonance,  <cit>  and many other techniques can be used to determine the hinge location.
assignment of domains and calculation of centers of mass
our method relies on the identification of a "stationary" domain s , a mobile domain m , and a linker or hinge region l .
the linker can be single, double, or triple stranded, and the two domains can be continuous or discontinuous  <cit> .
each residue is assigned to s, m, or l. the centers of mass   <cit>  of s, m, and l are labelled xs, xm, and xl, respectively.
these domain and com definitions are used in the subsequent preparation and manipulation of the structure.
preparation of protein and ligand
the conformation explorer at this time can handle only single protein chains.
therefore all additional peptides, ligands, metals, water, and dissolved ions are removed.
the ligand of interest is docked to the structure using autodock.
the entire receptor is in the grid map; we do not prejudice the docking with binding cleft information.
the simulation is not sensitive to this initial docking, since there are subsequent re-docking steps as we will explain.
it is important for the ligand to be docked before the equilibration , so that the latter does not lead to side chains filling in the binding cleft and and otherwise blocking the ligand from re-docking  <cit> .
the protein is then put into standard orientation.
this is the convention that xl coincide with the origin, xm lie along the z-axis, and xs lie in the -y part of the yz plane .
the coordinates of the ligand are required in pdb format, as is the gromacs .itp  and the autodock .pdbq  files.
the latter two can be generated from the first using the prodrg program  <cit> .
definition of rotations
the three euler rotations mentioned earlier are the space-centered rotation of the m domain about the z-axis, rz, about the x-axis, rx, and about the y-axis, ry, applied in that order.
.
the θz, θx, and θy are the angles of the corresponding rotations.
the ligand is not rotated, i.e.
remains stationary with respect to s.
equilibration
the preceding rotation step invariably results in unphysical bond lengths and angles in the boundary between l and m, and often in steric clashes between m and the rest of the protein and/or ligand.
a molecular dynamics equilibration is performed using the gromacs mdrun program for  <dig> time steps .
we solvated using a water box with a neutralizing atmosphere of cl- and na+ ions.
we found that  <dig> time steps  were sufficient to relieve the most significant steric strains, and  <dig> time steps were sufficient to allow an exponential decrease in enthalpy to level off , indicating a better equilibrated structure.
it was not sufficient time to allow substantial domain motions.
equilibration significantly increased the accuracy of ligand-binding prediction, compared to what we got with no equilibration.
angle calculation
we calculate the angular position of m of a generated structure in θyθxθz space by determining the ryrxrz, rotation that would have to be applied to m of the starting structure in standard orientation, to obtain a structure similar to the generated one.
this is done by first structurally aligning the generated structure with the starting structure by minimizing the root-mean-square deviation  between the s domains of the two structures.
the θy, θx, and θz angles are then computed.
this calculation is performed immediately following each equilibration step above.
note that in the course of the equilibration the angular position of m will change.
we refer to this as drift.
re-docking
at the end of each equilibration the ligand is removed from the structure file and then docked again to the protein using autodock  <dig>  <cit> .
the docking code reports multiple poses  with corresponding free energies of binding; we record the lowest of the latter as the gd, as we discuss below.
the corresponding ligand coordinates are then appended to the protein structure file.
terms of f: gd, rg, ef
the fitness function f we devised to discriminate the holo structure from the many conformations generated includes four quantities which separately have strengths and weaknesses, but which together form a predictor.
in this section we describe the terms, which are combined in a weighted sum to create f.
gd, as mentioned, was computed using autodock at the end of each equilibration.
this quantity was found early on to have significant ability to discriminate between near-holo structures and decoys.
however in some cases an unnatural cleft was generated by the mentioned rotations, to which the ligand docked with high affinity.
such unnatural conformations often have large radii of gyration and can be filtered out on that basis.
we also found that this measure did not vary smoothly with srmsd   <cit> , thus no gradient was available along which to minimize.
for these reasons we added terms to f as follows.
rg typically decreases during the cleft closure often associated with ligand binding, as has been observed by small angle x-ray scattering ,  <cit>  and as we will further demonstrate here.
this quantity is the distance rg from the center of mass at which all the mass of the protein's cα atoms can be concentrated to result in the same moment of inertia:  <cit>
rg=∑imri2m
where i counts over all atoms, ri is the distance  from the atom to the center of mass, and m is the total mass.
it also varies smoothly with srmsd, and so helps deal with the noise issues mentioned above and leads to improved convergence.
however rg alone is not a good predictive measure because it is trivially possible to minimize it with an unstable protein structure consisting of distorted interpenetrating domains.
this problem led to the introduction of the next quantity into f.
energy of folding, or ef, can be used to discriminate energetically feasible conformers from those with excessive interpenetration or unnatural orientation of domains.
this quantity is estimated using foldx.
the latter is a force field empirically fitted by guerois et al.
using a database of mutationally induced changes in protein stability  <cit> .
however since a wide variety of conformers, including both the holo and the apo, are stable, this measure alone cannot find the holo structure in an ensemble.
the point of computing this quantity is again to exclude unphysical structures.
figure of merit: srmsd
once the alternate conformers have been generated by rotation, equilibration, and re-docking, we are faced with the problem of determining how distant these are from the known target conformer, in our case the protein from the co-crystallized complex.
for this purpose maiorov et al.
used a "static" root-mean-square deviation , defined as the rmsd of the cα atoms in m, given that the cα atoms of s are optimally superimposed.
we use the same measure in the current work.
five-fold leave-one-out cross-validation of f
we computed gd, ef, and rg <dig> for every structure in each of the five generated ensembles .
we then divided the proteins into a single test protein and a training set consisting of the remaining four proteins.
we did this five times, choosing a different test protein each time.
each time, f was fitted to the training set as follows.
we first define f as:
f=xλ
where x is a matrix with one row for each structure in the training set and one column for each of the quantities gd, ef, and rg <dig>  the element in a given row and column is the physical quantity corresponding to that column, computed on the structure corresponding to that row.
λ is a three-element column vector of weights to be fitted as follows.
we assert that:
xλ≈y
where y is a column vector with each row containing the quantity srmsd, computed for the structure corresponding to that row.
rows in y correspond to rows in x.
λ is then given by:  <cit>
λ=-1xty
five different λ's were computed in five-fold leave-one-out bootstrapping,  <cit>  with one λ used for each protein, except where noted below.
holo, predicted, apo, and infeasible structures
for each protein studied, we discuss three protein-ligand configurations .
the holo structure is the crystallographically determined structure of the enzyme bound with ligand, which we use as a gold standard, and which would presumably be unknown in a practical situation.
the apo structure is the protein structure which would be known in practice and which is used as a starting point in our analysis.
the generated ensemble is the set of all structures generated by rotating m, equilibrating and docking.
the predicted structure is the conformer chosen from the generated ensemble as the best approximation of the holo structure.
in the ligand docking example explained here, it is selected based on the described f. an infeasible structure is one that could not be generated for one of two reasons:  the rotation was effected but the steric strains generated were too large and so the equilibration step failed to converge, or  the equilibration converged but the drift left the structure not in the targeted cell  but in a neighboring one, even after two successive attempts at rotation and equilibration.
in both cases the rotated but unequilibrated structure in the target cell was marked as infeasible.
this marked the corresponding cell as being on the boundary of the accessible conformation space and prevented further attempts to generate conformers in it.
minimization of f by line search algorithm
we minimized f using the line search algorithm, which successively varies each of the angles  over its full range.
after each angle is varied, the conformer that minimized f over the full range of that angle is used as a starting point for the next minimization.
the sampling of θy, θx, θz characteristic of this method can be appreciated graphically in figure  <dig> .
the line search minimization can be considered converged when no ry, rx, or rz rotation from the lowest-scoring conformer can further reduce the value of f. the algorithm is summarized in additional file  <dig>  in the following sections we describe the proteins we used for training, testing, and applying f.
glutamine binding protein
the motion of glnbp  <cit>  as it binds glutamine involves large-displacement domain hinge bending, experimentally determined to take ~5ns  <cit>   <dig> ns molecular dynamics simulations  of the apo structure failed to result in domain closure, possibly in part because the dynamics of the apo structure were computed with no ligand present  <cit> .
the apo or starting structure used,  <cit>  and the holo structure  <cit>  both originate in e. coli.
for glnbp flexoracle  <cit>  found a hinge at residues 86- <dig> and 182- <dig>  we selected residues  <dig> , <dig>  and  <dig> as the hinge location for the purpose of generating the rotations, but since the flexure occurs at the boundary between l and m, this adjustment makes little difference.
biotin carboxylase
acetyl-coa carboxylase , found in all animals, plants, and bacteria, catalyzes the carboxylation of acetyl-coa to malonyl-coa, the first committed step of fatty acid synthesis.
the first half-reaction is the formation of carboxybiotin which is catalyzed by the biotin carboxylase  subunit  <cit> .
pyruvate carboxylase  is found in many eukaryotes and some prokaryotes.
it plays a role in gluconeogenesis, mediating the carboxylation of pyruvate to oxaloacetate.
it has three functional domains, of which biotin carboxylase  is one.
the half-reaction catalyzed by bc is common to acc and pc, although the second half-reaction catalyzed by a different subunit differs from enzyme to enzyme depending on the substrate  <cit> .
upon binding atp, m rotates approximately 45° with respect to s. the large scale of this change and the existence of structural domains makes bc a good test case for our method.
we used the bc subunit of aquifex aeolicus pc  <cit>  as the starting or apo structure.
the ligand-free bc subunit of e. coli acc has been crystallized , but many residues are disordered and so the structure of the m domain is not entirely clear  <cit> .
no atp-bound structure of a. aeolicus bc is available, and so we used atp-bound bc from e. coli acc  <cit>  as our holo structure.
the sequence differences between the starting and holo structures had an impact on the results as we will discuss.
flexoracle  <cit>  predicts a hinge at residues 86- <dig> and 182-185; this can be used to choose the location of l.
ribose binding protein
rbp is a member of a large family of bacterial periplasmic binding proteins,  <cit>  and displays clear domain hinge bending motion  <cit> .
it has a triple stranded hinge, with a small fragment of the n-terminus crossing over onto the m-domain.
it binds ribose, a small, oxygen-rich ligand which we will discuss in more depth later.
adenylate kinase
adenylate kinase is a popular example of a domain hinge bending protein .
even though the holo and starting structures come from the same organism, they are from different compartments and have significant structural differences.
the holo structure was extracted from the mitochondrial intermembrane space, while the starting was taken from the mitochondrial matrix .
it naturally binds mg-atp and amp  <cit> , but it has been crys-tallized with so <dig> .
apo
holo
the five proteins  were each submitted to ce.
for each protein, the parameters of f were fitted using the results of the other four proteins.
in a preliminary round we computed closure of rbp binding its natural ligand, ribose, and of adk binding so <dig>  but srmsd of the predicted structures was high .
in the final round of calculations these two proteins were submitted with atp as a ligand , resulting in much lower srmsd.
glnbp and bc were submitted with their natural ligands, and mura with an antibiotic.
the motion of adk involves not two domains but three.
here we predicted the motion of the lid domain  with respect to the core domain .
since the ampbd domain  moves separately, for computing srmsd we aligned based only on the core domain.
mura
mura, which has been cyrstallized bound to the antibiotic t <dig>  <cit> .
the peculiarity of this ligand is that it binds to the open conformation of mura, rather than the closed, and thus interferes with the mechanism of motion  <cit> .
in this case if the predicted bound conformation is correct, it should not differ significantly from the open structure by the measure of srmsd.
the starting structure used was taken from the complex of mura with an analogue of its natural ligand  <cit> .
abc transporter periplasmic ligand-binding protein
abc-pbp is a very sparsely characterized protein.
its natural ligand is unknown, and few other details are available about its function  <cit> .
it has been crystallized in the apo form .
these characteristics make the protein a typical case in which a motion predictor like ce would be useful.
in this work we make a blind prediction of its holo conformation.
RESULTS
natural vs. non-natural ligands
in a first fitting, we docked each receptor with the ligand it was crystallized with in holo .
we observed that for two of the proteins , the ligand penetrated deeply into one of the two domains rather than binding at the cleft between the two.
observing that in both cases the ligand was small and oxygen-rich , we hypothesized:
<dig>  that the docking is not effective for any small and oxygen-rich ligand.
<dig>  that the precise choice of ligand is unimportant for the objective of predicting closed structures.
we ran ce on glnbp with non-natural  ligands .
we found that the three highest srmsd's were obtained for so <dig>  ribose, and oxalic acid, the three small oxygen-rich ligands.
in the three cases the ligand penetrated one of the two domains rather than binding in the cleft.
.
finding that atp yielded the lowest srmsd's, we ran ce on all five proteins using atp.
we found that in the cases of rbp and adk, the srmsd was lower than with their natural  ligands.
these two results support hypothesis .
we also found that for glnbp and mura, the results were slightly worse using atp than the natural ligand, though only by an average of  <dig> Å. thus the specific ligand matters, but only slightly.
in our final result we use atp in place of so <dig> and ribose, but make no such substitution when the natural ligand is not small and oxygen-rich .
the three small, oxygen-rich ligands  led to the worst results, because the ligands bound deep inside one of the two domains rather than at the cleft.
other ligands led to results comparable to those of glutamine.
glutamine binding protein
we ran ce on this protein to demonstrate the online results browser, and to rationalize the three terms of the fitness function.
the structure with lowest gd is shown in inset  <dig> e, and is also qualitatively similar to the holo.
the structure with lowest gyration radius is shown in inset  <dig> f. the holo tends to have a smaller rg than the apo structure.
rg cannot be used alone in f since it is trivially possible to minimize this quantity with a compact structure which is unstable, has significantly distorted domains, and does not resemble the holo.
the gyration radius appears squared for the physical reasons mentioned earlier.
the structure with most favorable ef, inset  <dig> g, is quite close to the starting structure.
this should not be surprising, since such structures have gone through a minimum of structural manipulation and so have low steric strain.
biotin carboxylase
the significant sequence and structural differences between the starting and holo structures, which come from different organisms, made comparison on the basis of srmsd somewhat difficult.
the structure that minimizes f  bears clear similarities to the holo in relative domain orientation.
ribose binding protein
in figure  <dig>  the structure that minimizes f with atp as a ligand shows good agreement with the holo.
the natural ligand led to poorer performance as mentioned.
adenylate kinase
the starting structure had 17Å, while the predicted structure  had  <dig> Å srmsd  - almost  <dig> Å lower.
the qualitative agreement was excellent, as can be appreciated in figure  <dig>  on the other hand docking so <dig>  a ligand it has been crystallized with, led to poor results.
mura
we have now shown that the conformation explorer can predict the bound conformation for four cases in which a large scale closing motion is required.
but what if the ligand binds to the open conformation?
one such case is that of mura.
the predicted structure is more closed than the holo, indicating that f is biased towards closed structures .
abc-pbp
the above five proteins are well characterized, with a known holo complex, unlike cases of scientific interest for which ce is useful.
in order to demonstrate a practical application, we submitted abc transporter perplasmic ligand binding protein  to ce.
since the natural ligand is unknown, we used atp.
the predicted structure is shown in figure  <dig>  we await future experimental validation of the predicted coordinates, which can be downloaded from molmovdb.
fitness function
as mentioned in generating the predicted structure for each of glnbp, bc, adk, rbp, and mura, we trained on the remaining four proteins.
the final fitness function, trained on all five proteins, is:
f= <dig> ⋅gd+ <dig> ⋅ef+ <dig> ⋅rg <dig>
computer time
computational cost ranged from about  <dig> hours  to  <dig> hours  on a single processor.
for comparison, the method of seeliger et al.
<cit>  requires about four days on a 50-node cluster.
web server
the described tool can be used through our online server at http://molmovdb.org.
the user must first provide the structural coordinates of the receptor and ligand, as described in figure  <dig>  he/she must then select the hinge points.
for this we make use of the hinge annotation tool on the database of macromolecular motions http://molmovdb.org morph page, as described in separate work  <cit> .
briefly, it consists of arrow buttons which control a moving window of two residues, highlighted in the jmol viewer .
the selected residue numbers can also be queried using the "?"
button.
once the desired hinge location has been highlighted, the "submit" button causes it to be recorded.
the user's job is queued upon submission.
each time an angle is varied over its full range by the line search minimizer, the user is emailed a progress notice.
this process can be tracked and inspected using a viewer similar to that described in figure  <dig>  a few iterations later the job will converge and the user will receive a link to a viewer that animates the ligand binding motion.
links to structural coordinate files of the trajectory and the predicted holo structure are also provided.
discussion
prediction of protein conformational change
for rbp, adk, glnbp, and bc, large scale motions of the m domain are required for binding.
for three of those, we successfully predicted the closed state within 6Å srmsd.
for bc, the predicted srmsd was higher, but still about 4Å lower than for the starting structure, and with qualitative similarities to the holo .
ligands binding to the open conformation
mura is a special case in which the ligand binds the open form.
our predicted structure was closed, with srmsd higher than for its starting structure, by about 4Å.
this is mostly because f, as a result of being trained on rbp, adk, glnbp, and bc for this protein, is explicitly biased towards closed structures through its gyration_radius term.
if the docked_energy term were used alone in place of f, srmsd would be much better, though still  <dig> Å higher than for the starting structure .
small, oxygen-rich ligands
the ligands so <dig> and ribose docked deep inside one of the two domains of adk and rbp , leading to poor srmsd.
so <dig>  ribose, and oxalic acid behaved similarly with glnbp, leading to the highest srmsd compared with other ligands .
this may be due to an insufficient desolvation penalty in the docking force field, or to the lack of explicit ions in the docking, at least for this class of ligands.
the above issue goes to the root of the utility of docking in predicting closed structures.
we observed that the docked ligands in the predicted structures had relatively few receptor contacts in common with those in experimentally observed holo structures.
and yet each of our five-fold leave-one-out studies verified the utility of the docked_energy term in f, for the purpose of detecting closed structures.
it may be that the docking force field is biased to favor closed structures, in which the ligand can maximize the number of contacts in a non-specific manner.
the cleft may also have non-specific ligand-binding properties  <cit> .
this would explain the success of using various non-natural ligands to predict the closed form of glnbp .
it would also explain why ce performed well on rbp when atp was used instead of its natural ligand .
ce's performance on glnbp, bc, and mura was somewhat worse  when atp was substituted for the natural ligand, indicating at least some predicted interactions are specific.
CONCLUSIONS
we conclude that for many hinge bending proteins it is possible to generate conformers similar to the closed structure given the apo.
we note that due to shortcomings in the docking method, any small, oxygen-rich ligands should be replaced, and atp is a good alternative.
the hinge location was predicted using an existing hinge prediction server, and iterative rotations, equilibrations, and docking runs resulted in prediction of closed structures similar to those known crystallographically, in four of the five cases.
the computational cost was moderate, permitting practical implementation on a single processor.
authors' contributions
sf conceived the idea, did the work, and wrote the manuscript.
mg refined the idea, directed the work and provided critical reading and guidance on the manuscript.
both authors have read and approved the work.
supplementary material
additional file 1
implementation of the line search algorithm.
starting from the apo structure, we generate conformations in the pitch, yaw, and roll rotational directions.
after exploring in each direction, the conformer that minimizes f is the starting point for exploration in the next direction.
the algorithm is converged when no rotation is possible in any direction that further minimizes f. * a particular direction is exhaustively explored when one conformer has been generated or attempted every 15°  in that direction, holding the other two direction angles constant .
click here for file
acknowledgements
we greatly appreciate help from robert bjornson in running an early version of conformation explorer on yale's bulldog cluster.
william grenawitzke wrote a script to dock ligands to generated protein conformers using autodock.
qingping xu and johan Åqvist provided useful advice.
trilok shahi provided editorial help.
sf gratefully acknowledges support from the swedish research council through the essence project grant .

 conformational changes in proteins can take place in a wide variety of ways, not all of which have been formally classified. one important class of motions is shear, in which stacked side chains of the protein can slide without losing contact. in this work we focus on the largest class, domain hinge bending, in which one structural domain of the protein moves relative to another domain about a hinge which connects the two. such motions typically involve the slowest degrees of freedom of that protein and so are difficult to predict by existing methods. the prediction of ligand binding motions of the protein receptor has considerable potential applications in protein-protein and protein-ligand docking. many methods can predict the side chain rearrangements required for docking but these assume that the large scale conformation is already nearly correct. thus there is a need for a method that will put the receptor in the correct large scale conformation which can be a productive starting point. much work has been done in this area. molecular dynamics explicitly computes the dynamical trajectory of molecules modeled as point masses connected by linear and nonlinear springs and can be used to predict conformational change, but usually only small- or moderate-scale domain motions can be reproduced with many biologically relevant motions remaining out of reach. accordingly several methods used md to account for the fast fluctuations of proteins in drug docking by first computing the protein trajectory using md. one limitation of such techniques is that they may not escape the vicinity of an initial conformation, even in a time span experimentally known to be sufficient for conformational change. althaus et al created a combinatorial tree of side-chain rotamers which they explored using a branch-and-cut algorithm, without varying the backbone conformation. sandak et al. created a flexible-receptor docking code which articulates the protein at a hinge point, but leaves the two resulting domains rigid. this method suffered from the opposite problem: it could generate large scale protein motions, but had no way of dealing with even small side chain rearrangements, a weakness leading to failure. the described methods are good at either treating the side-chain flexibility, or the large scale conformational changes, but not both simultaneously. conformation explorer uses sandak et al.'s idea of moving domains about a hinge point to generate large scale conformational change, but also includes equilibration steps which permit relaxation and adjustment of all atoms. normal modes have also been used by many authors to predict the conformational changes of proteins. comparison of the atomic coordinates of homologous pairs of proteins shows that the lowest order modes are most involved in conformational change, but also that multiple modes are needed to accurately represent the motion. it is possible to determine the correct combination of normal modes that will reproduce a desired motion, but this requires knowledge of at least a few interatomic distance constraints for the final structure. in a different approach, a docked protein-ligand complex was displaced along the lowest-frequency normal mode directions to minimize non-bonded energy terms in an md force field. however a normal mode expansion assumes a quadratic potential and so is accurate only for small fluctuations about an equilibrium structure; therefore the method cannot be used to predict larger scale conformational changes such as we treat in this work. the method of lindahl et al. gains improvements of to Å for several proteins; our method recapitulates much larger conformational changes as we will show. maiorov and abagyan rigidified all protein bonds except those in the interdomain linker and interface using internal coordinate modeling, and then used the biased probability monte carlo protocol to generate potential alternate conformations of the protein. the method succeeded in generating a large number of alternate conformations, and some of these were somewhat similar to alternate conformations known crystallographically. however without referring to the known alternate conformations, it was impossible to determine which of the many predicted structures was thermodynamically plausible. further, many energy evaluations and minimizations were expended in evaluating generated conformers which were later discarded. lastly, it was not easy to know how long a thorough exploration of conformation space would take, and no clear way to restrict the search to a given region of interest. our method is similar in several ways to maiorov et al.'s, but also addresses these limitations. in more recent work, de groot et al. showed they could find the holo conformations of several ligand-binding proteins. the method relies on tconcoord, which determines flexible regions by analyzing hydrogen bonding networks. once these are known, an ensemble of plausible structures is generated. an interative process involving docking, md refinement, and filtering by radius of gyration then generates holo structures. however the radius of gyration must be provided by experiments, because their fitness function does not take into account the receptor conformation. our work resembles de groot's in the use of docking and md refinement, but differs in that the hinge location is determined by the user, and in that it requires no experimental information about the holo complex. our method is also computationally cheaper. key to our method is the conformational sampling and the use of a fitness function f which includes terms to discriminate the holo conformation from the generated ensemble. the mentioned f deserves motivation. there are numerous potentials which pick out the correct ligand-binding pose in a given receptor. however these typically have no terms to discriminate favorable from unfavorable receptor conformations. other force fields have been trained and tested against ensembles of crystallographically obtained protein structures, but it is not clear how they would perform when the holo structure is unavailable. to overcome these limitations f includes terms which help discriminate predicted holo complexes from a generated ensemble. in this work we show that the conformation explorer can predictably, controllably, and economically generate multiple alternate conformations of proteins without experimental constraints on the holo complex. the conformations can be equilibrated to ensure they are free of steric strain. we demonstrate that as an application of the method, we can predict the ligand bound, or holo, conformation of a protein given only the ligand and the coordinates of the protein in the apo state. the exploration of phase space is performed in such a way as to minimize f, which includes terms measuring the radius of gyration rg, stability ef, and docked energy gd. we will formally define all of these terms in the methods section. methods as mentioned, the method involves the identification of the hinge and rotations about that hinge, followed by equilibrations and docking runs. the direction of the rotations and the choice of intermediate structures to which the rotations are applied, are subject to optimization. we describe how to perform these steps. determining and recording the hinge location the hinge can be detected in several ways: theoretically hinge predictors such as flexoracle, tlsmd, hingemaster, and others can be used to predict the hinge location if at least one structure is known. the hingemaster web server makes the use of the first three easy and was used in this work. crystallographically if the protein has been crystallized in two different conformations, these can be inspected visually to determine the hinge location. this is usually the best and most accurate method. by other experimental means proteolysis, hydrogen/deuterium exchange, optical trapping, fluorescent resonance energy transfer, nuclear magnetic resonance, and many other techniques can be used to determine the hinge location. assignment of domains and calculation of centers of mass our method relies on the identification of a "stationary" domain s, a mobile domain m, and a linker or hinge region l. the linker can be single, double, or triple stranded, and the two domains can be continuous or discontinuous. each residue is assigned to s, m, or l. the centers of mass of s, m, and l are labelled xs, xm, and xl, respectively. these domain and com definitions are used in the subsequent preparation and manipulation of the structure. preparation of protein and ligand the conformation explorer at this time can handle only single protein chains. therefore all additional peptides, ligands, metals, water, and dissolved ions are removed. the ligand of interest is docked to the structure using autodock. the entire receptor is in the grid map; we do not prejudice the docking with binding cleft information. the simulation is not sensitive to this initial docking, since there are subsequent re-docking steps as we will explain. it is important for the ligand to be docked before the equilibration, so that the latter does not lead to side chains filling in the binding cleft and and otherwise blocking the ligand from re-docking. the protein is then put into standard orientation. this is the convention that xl coincide with the origin, xm lie along the z-axis, and xs lie in the -y part of the yz plane. the coordinates of the ligand are required in pdb format, as is the gromacs.itp and the autodock.pdbq files. the latter two can be generated from the first using the prodrg program. definition of rotations the three euler rotations mentioned earlier are the space-centered rotation of the m domain about the z-axis, rz, about the x-axis, rx, and about the y-axis, ry, applied in that order.. the θz, θx, and θy are the angles of the corresponding rotations. the ligand is not rotated, i.e. remains stationary with respect to s. equilibration the preceding rotation step invariably results in unphysical bond lengths and angles in the boundary between l and m, and often in steric clashes between m and the rest of the protein and/or ligand. a molecular dynamics equilibration is performed using the gromacs mdrun program for time steps. we solvated using a water box with a neutralizing atmosphere of cl- and na+ ions. we found that time steps were sufficient to relieve the most significant steric strains, and time steps were sufficient to allow an exponential decrease in enthalpy to level off, indicating a better equilibrated structure. it was not sufficient time to allow substantial domain motions. equilibration significantly increased the accuracy of ligand-binding prediction, compared to what we got with no equilibration. angle calculation we calculate the angular position of m of a generated structure in θyθxθz space by determining the ryrxrz, rotation that would have to be applied to m of the starting structure in standard orientation, to obtain a structure similar to the generated one. this is done by first structurally aligning the generated structure with the starting structure by minimizing the root-mean-square deviation between the s domains of the two structures. the θy, θx, and θz angles are then computed. this calculation is performed immediately following each equilibration step above. note that in the course of the equilibration the angular position of m will change. we refer to this as drift. re-docking at the end of each equilibration the ligand is removed from the structure file and then docked again to the protein using autodock. the docking code reports multiple poses with corresponding free energies of binding; we record the lowest of the latter as the gd, as we discuss below. the corresponding ligand coordinates are then appended to the protein structure file. terms of f: gd, rg, ef the fitness function f we devised to discriminate the holo structure from the many conformations generated includes four quantities which separately have strengths and weaknesses, but which together form a predictor. in this section we describe the terms, which are combined in a weighted sum to create f. gd, as mentioned, was computed using autodock at the end of each equilibration. this quantity was found early on to have significant ability to discriminate between near-holo structures and decoys. however in some cases an unnatural cleft was generated by the mentioned rotations, to which the ligand docked with high affinity. such unnatural conformations often have large radii of gyration and can be filtered out on that basis. we also found that this measure did not vary smoothly with srmsd, thus no gradient was available along which to minimize. for these reasons we added terms to f as follows. rg typically decreases during the cleft closure often associated with ligand binding, as has been observed by small angle x-ray scattering, and as we will further demonstrate here. this quantity is the distance rg from the center of mass at which all the mass of the protein's cα atoms can be concentrated to result in the same moment of inertia: rg=∑imri2m where i counts over all atoms, ri is the distance from the atom to the center of mass, and m is the total mass. it also varies smoothly with srmsd, and so helps deal with the noise issues mentioned above and leads to improved convergence. however rg alone is not a good predictive measure because it is trivially possible to minimize it with an unstable protein structure consisting of distorted interpenetrating domains. this problem led to the introduction of the next quantity into f. energy of folding, or ef, can be used to discriminate energetically feasible conformers from those with excessive interpenetration or unnatural orientation of domains. this quantity is estimated using foldx. the latter is a force field empirically fitted by guerois et al. using a database of mutationally induced changes in protein stability. however since a wide variety of conformers, including both the holo and the apo, are stable, this measure alone cannot find the holo structure in an ensemble. the point of computing this quantity is again to exclude unphysical structures. figure of merit: srmsd once the alternate conformers have been generated by rotation, equilibration, and re-docking, we are faced with the problem of determining how distant these are from the known target conformer, in our case the protein from the co-crystallized complex. for this purpose maiorov et al. used a "static" root-mean-square deviation, defined as the rmsd of the cα atoms in m, given that the cα atoms of s are optimally superimposed. we use the same measure in the current work. five-fold leave-one-out cross-validation of f we computed gd, ef, and rg for every structure in each of the five generated ensembles. we then divided the proteins into a single test protein and a training set consisting of the remaining four proteins. we did this five times, choosing a different test protein each time. each time, f was fitted to the training set as follows. we first define f as: f=xλ where x is a matrix with one row for each structure in the training set and one column for each of the quantities gd, ef, and rg the element in a given row and column is the physical quantity corresponding to that column, computed on the structure corresponding to that row. λ is a three-element column vector of weights to be fitted as follows. we assert that: xλ≈y where y is a column vector with each row containing the quantity srmsd, computed for the structure corresponding to that row. rows in y correspond to rows in x. λ is then given by: λ=-1xty five different λ's were computed in five-fold leave-one-out bootstrapping, with one λ used for each protein, except where noted below. holo, predicted, apo, and infeasible structures for each protein studied, we discuss three protein-ligand configurations. the holo structure is the crystallographically determined structure of the enzyme bound with ligand, which we use as a gold standard, and which would presumably be unknown in a practical situation. the apo structure is the protein structure which would be known in practice and which is used as a starting point in our analysis. the generated ensemble is the set of all structures generated by rotating m, equilibrating and docking. the predicted structure is the conformer chosen from the generated ensemble as the best approximation of the holo structure. in the ligand docking example explained here, it is selected based on the described f. an infeasible structure is one that could not be generated for one of two reasons: the rotation was effected but the steric strains generated were too large and so the equilibration step failed to converge, or the equilibration converged but the drift left the structure not in the targeted cell but in a neighboring one, even after two successive attempts at rotation and equilibration. in both cases the rotated but unequilibrated structure in the target cell was marked as infeasible. this marked the corresponding cell as being on the boundary of the accessible conformation space and prevented further attempts to generate conformers in it. minimization of f by line search algorithm we minimized f using the line search algorithm, which successively varies each of the angles over its full range. after each angle is varied, the conformer that minimized f over the full range of that angle is used as a starting point for the next minimization. the sampling of θy, θx, θz characteristic of this method can be appreciated graphically in figure. the line search minimization can be considered converged when no ry, rx, or rz rotation from the lowest-scoring conformer can further reduce the value of f. the algorithm is summarized in additional file in the following sections we describe the proteins we used for training, testing, and applying f. glutamine binding protein the motion of glnbp as it binds glutamine involves large-displacement domain hinge bending, experimentally determined to take ~5ns ns molecular dynamics simulations of the apo structure failed to result in domain closure, possibly in part because the dynamics of the apo structure were computed with no ligand present. the apo or starting structure used, and the holo structure both originate in e. coli. for glnbp flexoracle found a hinge at residues 86- and 182- we selected residues, and as the hinge location for the purpose of generating the rotations, but since the flexure occurs at the boundary between l and m, this adjustment makes little difference. biotin carboxylase acetyl-coa carboxylase, found in all animals, plants, and bacteria, catalyzes the carboxylation of acetyl-coa to malonyl-coa, the first committed step of fatty acid synthesis. the first half-reaction is the formation of carboxybiotin which is catalyzed by the biotin carboxylase subunit. pyruvate carboxylase is found in many eukaryotes and some prokaryotes. it plays a role in gluconeogenesis, mediating the carboxylation of pyruvate to oxaloacetate. it has three functional domains, of which biotin carboxylase is one. the half-reaction catalyzed by bc is common to acc and pc, although the second half-reaction catalyzed by a different subunit differs from enzyme to enzyme depending on the substrate. upon binding atp, m rotates approximately 45° with respect to s. the large scale of this change and the existence of structural domains makes bc a good test case for our method. we used the bc subunit of aquifex aeolicus pc as the starting or apo structure. the ligand-free bc subunit of e. coli acc has been crystallized, but many residues are disordered and so the structure of the m domain is not entirely clear. no atp-bound structure of a. aeolicus bc is available, and so we used atp-bound bc from e. coli acc as our holo structure. the sequence differences between the starting and holo structures had an impact on the results as we will discuss. flexoracle predicts a hinge at residues 86- and 182-185; this can be used to choose the location of l. ribose binding protein rbp is a member of a large family of bacterial periplasmic binding proteins, and displays clear domain hinge bending motion. it has a triple stranded hinge, with a small fragment of the n-terminus crossing over onto the m-domain. it binds ribose, a small, oxygen-rich ligand which we will discuss in more depth later. adenylate kinase adenylate kinase is a popular example of a domain hinge bending protein. even though the holo and starting structures come from the same organism, they are from different compartments and have significant structural differences. the holo structure was extracted from the mitochondrial intermembrane space, while the starting was taken from the mitochondrial matrix. it naturally binds mg-atp and amp, but it has been crys-tallized with so. apo holo the five proteins were each submitted to ce. for each protein, the parameters of f were fitted using the results of the other four proteins. in a preliminary round we computed closure of rbp binding its natural ligand, ribose, and of adk binding so but srmsd of the predicted structures was high. in the final round of calculations these two proteins were submitted with atp as a ligand, resulting in much lower srmsd. glnbp and bc were submitted with their natural ligands, and mura with an antibiotic. the motion of adk involves not two domains but three. here we predicted the motion of the lid domain with respect to the core domain. since the ampbd domain moves separately, for computing srmsd we aligned based only on the core domain. mura mura, which has been cyrstallized bound to the antibiotic t. the peculiarity of this ligand is that it binds to the open conformation of mura, rather than the closed, and thus interferes with the mechanism of motion. in this case if the predicted bound conformation is correct, it should not differ significantly from the open structure by the measure of srmsd. the starting structure used was taken from the complex of mura with an analogue of its natural ligand. abc transporter periplasmic ligand-binding protein abc-pbp is a very sparsely characterized protein. its natural ligand is unknown, and few other details are available about its function. it has been crystallized in the apo form. these characteristics make the protein a typical case in which a motion predictor like ce would be useful. in this work we make a blind prediction of its holo conformation.  natural vs. non-natural ligands in a first fitting, we docked each receptor with the ligand it was crystallized with in holo. we observed that for two of the proteins, the ligand penetrated deeply into one of the two domains rather than binding at the cleft between the two. observing that in both cases the ligand was small and oxygen-rich, we hypothesized: that the docking is not effective for any small and oxygen-rich ligand. that the precise choice of ligand is unimportant for the objective of predicting closed structures. we ran ce on glnbp with non-natural ligands. we found that the three highest srmsd's were obtained for so ribose, and oxalic acid, the three small oxygen-rich ligands. in the three cases the ligand penetrated one of the two domains rather than binding in the cleft.. finding that atp yielded the lowest srmsd's, we ran ce on all five proteins using atp. we found that in the cases of rbp and adk, the srmsd was lower than with their natural ligands. these two results support hypothesis. we also found that for glnbp and mura, the results were slightly worse using atp than the natural ligand, though only by an average of Å. thus the specific ligand matters, but only slightly. in our final result we use atp in place of so and ribose, but make no such substitution when the natural ligand is not small and oxygen-rich. the three small, oxygen-rich ligands led to the worst results, because the ligands bound deep inside one of the two domains rather than at the cleft. other ligands led to results comparable to those of glutamine. glutamine binding protein we ran ce on this protein to demonstrate the online results browser, and to rationalize the three terms of the fitness function. the structure with lowest gd is shown in inset e, and is also qualitatively similar to the holo. the structure with lowest gyration radius is shown in inset f. the holo tends to have a smaller rg than the apo structure. rg cannot be used alone in f since it is trivially possible to minimize this quantity with a compact structure which is unstable, has significantly distorted domains, and does not resemble the holo. the gyration radius appears squared for the physical reasons mentioned earlier. the structure with most favorable ef, inset g, is quite close to the starting structure. this should not be surprising, since such structures have gone through a minimum of structural manipulation and so have low steric strain. biotin carboxylase the significant sequence and structural differences between the starting and holo structures, which come from different organisms, made comparison on the basis of srmsd somewhat difficult. the structure that minimizes f bears clear similarities to the holo in relative domain orientation. ribose binding protein in figure the structure that minimizes f with atp as a ligand shows good agreement with the holo. the natural ligand led to poorer performance as mentioned. adenylate kinase the starting structure had 17Å, while the predicted structure had Å srmsd - almost Å lower. the qualitative agreement was excellent, as can be appreciated in figure on the other hand docking so a ligand it has been crystallized with, led to poor results. mura we have now shown that the conformation explorer can predict the bound conformation for four cases in which a large scale closing motion is required. but what if the ligand binds to the open conformation? one such case is that of mura. the predicted structure is more closed than the holo, indicating that f is biased towards closed structures. abc-pbp the above five proteins are well characterized, with a known holo complex, unlike cases of scientific interest for which ce is useful. in order to demonstrate a practical application, we submitted abc transporter perplasmic ligand binding protein to ce. since the natural ligand is unknown, we used atp. the predicted structure is shown in figure we await future experimental validation of the predicted coordinates, which can be downloaded from molmovdb. fitness function as mentioned in generating the predicted structure for each of glnbp, bc, adk, rbp, and mura, we trained on the remaining four proteins. the final fitness function, trained on all five proteins, is: f= ⋅gd+ ⋅ef+ ⋅rg computer time computational cost ranged from about hours to hours on a single processor. for comparison, the method of seeliger et al. requires about four days on a 50-node cluster. web server the described tool can be used through our online server at http://molmovdb.org. the user must first provide the structural coordinates of the receptor and ligand, as described in figure he/she must then select the hinge points. for this we make use of the hinge annotation tool on the database of macromolecular motions http://molmovdb.org morph page, as described in separate work. briefly, it consists of arrow buttons which control a moving window of two residues, highlighted in the jmol viewer. the selected residue numbers can also be queried using the "?" button. once the desired hinge location has been highlighted, the "submit" button causes it to be recorded. the user's job is queued upon submission. each time an angle is varied over its full range by the line search minimizer, the user is emailed a progress notice. this process can be tracked and inspected using a viewer similar to that described in figure a few iterations later the job will converge and the user will receive a link to a viewer that animates the ligand binding motion. links to structural coordinate files of the trajectory and the predicted holo structure are also provided. discussion prediction of protein conformational change for rbp, adk, glnbp, and bc, large scale motions of the m domain are required for binding. for three of those, we successfully predicted the closed state within 6Å srmsd. for bc, the predicted srmsd was higher, but still about 4Å lower than for the starting structure, and with qualitative similarities to the holo. ligands binding to the open conformation mura is a special case in which the ligand binds the open form. our predicted structure was closed, with srmsd higher than for its starting structure, by about 4Å. this is mostly because f, as a result of being trained on rbp, adk, glnbp, and bc for this protein, is explicitly biased towards closed structures through its gyration_radius term. if the docked_energy term were used alone in place of f, srmsd would be much better, though still Å higher than for the starting structure. small, oxygen-rich ligands the ligands so and ribose docked deep inside one of the two domains of adk and rbp, leading to poor srmsd. so ribose, and oxalic acid behaved similarly with glnbp, leading to the highest srmsd compared with other ligands. this may be due to an insufficient desolvation penalty in the docking force field, or to the lack of explicit ions in the docking, at least for this class of ligands. the above issue goes to the root of the utility of docking in predicting closed structures. we observed that the docked ligands in the predicted structures had relatively few receptor contacts in common with those in experimentally observed holo structures. and yet each of our five-fold leave-one-out studies verified the utility of the docked_energy term in f, for the purpose of detecting closed structures. it may be that the docking force field is biased to favor closed structures, in which the ligand can maximize the number of contacts in a non-specific manner. the cleft may also have non-specific ligand-binding properties. this would explain the success of using various non-natural ligands to predict the closed form of glnbp. it would also explain why ce performed well on rbp when atp was used instead of its natural ligand. ce's performance on glnbp, bc, and mura was somewhat worse when atp was substituted for the natural ligand, indicating at least some predicted interactions are specific.  we conclude that for many hinge bending proteins it is possible to generate conformers similar to the closed structure given the apo. we note that due to shortcomings in the docking method, any small, oxygen-rich ligands should be replaced, and atp is a good alternative. the hinge location was predicted using an existing hinge prediction server, and iterative rotations, equilibrations, and docking runs resulted in prediction of closed structures similar to those known crystallographically, in four of the five cases. the computational cost was moderate, permitting practical implementation on a single processor. authors' contributions sf conceived the idea, did the work, and wrote the manuscript. mg refined the idea, directed the work and provided critical reading and guidance on the manuscript. both authors have read and approved the work. supplementary material additional file 1 implementation of the line search algorithm. starting from the apo structure, we generate conformations in the pitch, yaw, and roll rotational directions. after exploring in each direction, the conformer that minimizes f is the starting point for exploration in the next direction. the algorithm is converged when no rotation is possible in any direction that further minimizes f. * a particular direction is exhaustively explored when one conformer has been generated or attempted every 15° in that direction, holding the other two direction angles constant. click here for file acknowledgements we greatly appreciate help from robert bjornson in running an early version of conformation explorer on yale's bulldog cluster. william grenawitzke wrote a script to dock ligands to generated protein conformers using autodock. qingping xu and johan Åqvist provided useful advice. trilok shahi provided editorial help. sf gratefully acknowledges support from the swedish research council through the essence project grant.
fourth, we score the generated structures using a novel fitness function which favors closed or holo structures. by iterating the second through fourth steps we systematically minimize the fitness function, thus predicting the conformational change required for small ligand binding for five well studied proteins. second, we apply an euler rotation to one of the domains about the hinge point. often the 3d structure of the protein is available in some conformation, but binding the ligand of interest may involve a large scale conformational change which is difficult to predict with existing methods. we demonstrate that the method in most cases successfully predicts the holo conformation given only an apo structure.
 knowledge of the structure of proteins bound to known or potential ligands is crucial for biological understanding and drug design. often the 3d structure of the protein is available in some conformation, but binding the ligand of interest may involve a large scale conformational change which is difficult to predict with existing methods.  we describe how to generate ligand binding conformations of proteins that move by hinge bending, the largest class of motions. first, we predict the location of the hinge between domains. second, we apply an euler rotation to one of the domains about the hinge point. third, we compute a short-time dynamical trajectory using molecular dynamics to equilibrate the protein and ligand and correct unnatural atomic positions. fourth, we score the generated structures using a novel fitness function which favors closed or holo structures. by iterating the second through fourth steps we systematically minimize the fitness function, thus predicting the conformational change required for small ligand binding for five well studied proteins.  we demonstrate that the method in most cases successfully predicts the holo conformation given only an apo structure. 