BACKGROUND
the analysis of protein folding by computational methods faces the enormous confor-mational space to be explored, and usually requires simulation lengths of at least several hundreds of ns.
simplified models  <cit> , higher simulation temperature  <cit>  or implicit solvation  <cit>  can accelerate the process, but induce a deformation of the free energy surface.
here, we propose a method based on the analysis of the nmr experimental restraints, to enhance the structure instability observed during explicit solvent molecular dynamics  simulations with an all-atom force field.
each nmr restraint is ranked according to its redundancy with respect to the other restraints.
calculation of nmr protein con-formations is then performed, based on a reduced set of restraints, from which the least redundant restraints are removed.
these conformations are modified starting points for md simulations having more chances to display structure instability.
md simulations are recorded in the absence of any distance restraints, as the different restraint sets are only used for the generation of the starting points.
sh <dig> domain of nephrocystin , a β barrel protein associated with juvenile nephronophthisis  <cit> , was used to illustrate this approach.
nph sh <dig> mutations associated to nephronophthisis, are located at the residue l <dig> in the strand β <dig> .
the structure of wt nph sh <dig> was determined by nmr  <cit>  and the nmr analysis of the wt protein and of point mutants revealed that l28p is unfolded, whereas l28a is only partially folded.
the impact of these mutations on the protein structure is quite intriguing, as residue l <dig> is pointing toward the solvent .
.
<cit>
drawn in cartoon with β strands colored in magenta and the l <dig> colored in green and drawn in sticks.
topology of the structure of the sh <dig> domain of the nephrocystine extracted from pdbsum  <cit> http://www.ebi.ac.uk/pdbsum/, with l <dig> written in green.
the structure of the nph sh <dig> domain  <cit>  reveals a classical β-barrel fold, surrounded by three loops: the rt loop , the n-src loop  and the distal hairpin .
this figure was produced with pymol  <dig>   <cit>
the effect of removing the least redundant restraints was tested on wt nph sh <dig> by using three nmr restraints sets:  the set full including all measured nmr restraints,  the set reduced where few least redundant restraints were removed from the set full,  the sets random, where restraints were randomly removed from the set full.
three series of nmr conformations were calculated with these restraint sets, and display only small differences.
notably, the md trajectories initiated from these different conformation series, display different trends.
similar β barrel shapes were observed during the trajectories starting from nmr conformations calculated with the sets full and random, whereas β barrel instability is found during trajectories started from the conformers calculated using the set reduced.
the effect of removing the least redundant restraints was also investigated on the nph sh <dig> mutants.
the leucine  <dig> sidechain was substituted by proline or alanine sidechains in the nmr conformers calculated using the restraint sets full and reduced.
all md trajectories started from these mutated conformations are less stable than the corresponding md trajectories recorded on wt protein, but the instability is enhanced if starting points obtained with the set reduced are used.
the redundancy of each nmr restraint was evaluated through the calculation of the information it brings to the total nmr restraint set.
the more information is brought by a restraint, the less redundant it is from the others.
conversely, low-information restraints depend a lot on the other restraints and are thus redundant.
the information brought by each restraint was calculated using the approach queen, introduced some years ago by nabuurs et al  <cit> .
the instability features observed on nph sh <dig> domain are in agreement with previous experimental and simulation information on the unfolding of sh <dig> domains, as the rt loop and the n and c terminal parts are the most destabilized regions.
the changes in β barrel architecture are revealed by an original geometric analysis and are related to a packing defect of the hydrophobic core residues induced by the mutation of l <dig>
<dig> methods
<dig>  nmr conformers and restraints
the nmr restraints measured on the nph sh <dig> domain  <cit>  were obtained from the pdb entry 1s1n.
the approach queen  <cit>  was used to calculate the information brought by each restraint.
from this information measure, the redundancy of each restraint with respect to the others was estimated.
three sets of restraints were then prepared:  the set full, including all restraints from 1s1n.mr,  the set reduced, determined by discarding the nine least redundant restraints displaying the largest queen information values,  the set random, determined by discarding nine randomly chosen restraints, displaying low queen information in the range 0- <dig> .
these three sets of restraints were used to produce the corresponding sets ,  and  of nmr conformations for the wt nph sh <dig> domain.
nmr conformers were generated using cns  <dig>   <cit> .
the simulated annealing procedure of cns was embedded into aria  <dig>   <cit> .
a single iteration in geometric force field produced  <dig> conformers, from which the  <dig> lowest-energy conformers were submitted to a water refinement step.
the nmr conformers were evaluated through: the coordinates rmsd between the backbone atoms, the number of violated distance restraints, the rms of restraints violations  <cit> , and the procheck v <dig> .
<dig>  <cit>  and whatif  <dig>   <cit>  quality scores.
<dig>  molecular dynamics simulations
nmr conformers of the nph sh <dig> were starting points of the molecular dynamics  simulations, performed with explicit solvent, and periodic boundary conditions at constant temperature  and pressure .
a cutoff distance of  <dig> Å was used to determine the water box size.
trajectories were recorded using the package amber  <dig>   <cit>  and the force field ff99sb  <cit> .
a cutoff of  <dig> Å was used for lennard-jones interactions, and long-range electrostatic interactions were calculated with the particle mesh ewald  protocol  <cit> .
the systems were neutralized using sodium counter-ions.
pressure was regulated with isotropic position scaling and a relaxation time of  <dig> ps, and temperature was regulated using a langevin thermostat  <cit>  with a collision frequency of  <dig> ps- <dig> or a berendsen thermostat  <cit>  with a coupling time of  <dig> ps.
the shake algorithm  <cit>  was used to keep rigid covalent bonds involving hydrogens, enabling a time step of  <dig> fs.
atom coordinates were saved at every ps.
simulations were initiated by rounds of semi-restrained and unrestrained system minimization.
the systems were thermalized to 298k for  <dig> ps at constant volume, while restraining the positions of the solute with a force constant of  <dig> kcal/.
the following equilibration protocol was then performed:  <dig> md round of  <dig> ps at constant volume and  <dig> md rounds of  <dig>  ps were run while reducing the position restraints from  <dig> kcal/ down to  <dig> kcal/; eventually a last md round of  <dig> ps was per-formed with a restraint of  <dig>  kcal/.
nmr conformers were selected from the conformers sets , ,  described in subsection "nmr conformers and restraints", as starting points for md simulations.
for each conformer set,  <dig> runs of  <dig> ns md trajectories were recorded varying: the thermo-stat controlling the system temperature , the seed for the initial velocities generation  and the starting point chosen among the nmr conformer set.
the mutations l28p and l28a were introduced in the wt nmr conformers calculated using the restraint sets full and reduced.
the generalized born surface area  approach  <cit>  was used to determine energy profiles per residues on the protein sequence.
the generalized born energy was calculated using the onufriev method  <cit>  and the following parameters: no counter-ions, a solvent dielectric constant of  <dig>  and a protein dielectric constant of  <dig>  the surface area was computed by recursively approximating a sphere around an atom, starting from an icosahedra  <cit>  .
<dig>  the β barrel architecture
an original geometric approach was used to analyze the sh <dig> β-barrel architecture .
first, the cylinder of maximum radius such that no protein backbone atom can be located inside  was determined to give the approximate main axis of the barrel.
second, a series of slabs  was defined along the main axis, and, within each slab, the maximum radius of the cylinder that can be fitted inside the volume defined by the protein backbone and the slab limits was calculated.
as the slab width is  <dig> Å and the slab center is moved by steps of  <dig>  Å, the cylinder radius obtained for each slab gives a local measure of the size of the inner β-barrel cavity, and the profile of the cylinder radius along the slabs thus provides an accurate and objective picture of the β barrel shape.
as only the backbone atoms of the protein are included in the optimization and fitting procedures, the analysis focuses on the geometry of the β barrel, and considers the protein interior as a cavity.
the maximization of cylinder radius was implemented as a constrained smooth op-timization problem, and solved with the algencan solver  <cit> .
analytic derivatives were computed.
multiple initial points were considered to guarantee convergence to global solutions.
this aspect is particularly relevant for the determination of the main axis, since, within each slab, the radius maximization is straightforward.
the algen-can solver is freely available at http://www.ime.usp.br/~egbirgin/tango and the package used in the present work is available at http://lm-utils.googlecode.com.
<dig> theory
<dig>  the restraint information
in a protein structure determination by nmr, the relationship between protein atomic coordinates and distance restraints can be described in the frame of the distance geometry approach  <cit> .
the ensemble of protein atomic coordinates is considered as a 3d euclidean object.
consequently, the upper and lower bounds of distance restraints verify the triangle inequalities, and the addition of a distance restraint potentially modifies the upper and lower bounds of all other restraints.
as each restraint keeps the distance value into the lower/upper bound interval, it brings information improving the definition of the protein structure.
not all restraints bring the same quantity of information, as they do not have the same influence on the structure definition.
in other words, low-information restraints are quite redundant and thus very dependent on the other restraints, whereas high-information restraints are much less redundant with respect to the other restraints defining the structure.
these properties of the nmr structure determination were known for a long time, but they were recently expressed in a quantitative way in the queen approach  <cit> , by proposing that the distance restraints defining a structure could be analyzed in the frame of information theory.
the uncertainty hij of the distance between spin nuclei i and j, can be expressed as:
hij=-∫ lijuij1uij-lij log1uij-lijddij= log
where lij and uij are respectively the lower and the upper bounds of the restraint.
the uncertainty of a given set r of restraints can be then expressed as:
hstructure|r=1n∑i,j=1nhij
where n is the number of spin nuclei.
therefore, the unique information iuni,r of a restraint r, applied between spins nuclei i and j, is defined as the difference of uncertainties between the sets r and r-1:
iuni,r=hstructure|r-1-hstructure|r
where r- <dig> represent the set of restraints others than r in the data-set, and r is the total set of restraints including r. the value iuni,r is different from the uncertainty hij of the distance between spin nuclei i and j, as the addition of a restraint may induce a variation of all upper and lower bounds by application of the triangle inequalities.
the present work intends to use the queen framework quantifying the redundancy of each restraint, in order to destabilize protein structure along its unfolding path.
<dig>  unique information and contact order
the unique information iuni,r of the queen approach  <cit>  can be related to the protein folding through the contact order parameter.
the contact order co  <cit>  is defined from the ranges of the inter-residues contacts:
co=1ncontactnseq ∑k,l|k-l|
where k and l are the residue numbers forming a native contact in the folded protein, ncontact is the number of native contacts and nseq is the sequence size.
the parameter co is positive and smaller than  <dig>  the parameter co was shown  <cit>  to be correlated to the folding rates, because of the decrease of conformational entropy associated with the formation of contacts  <cit> .
to see how co and iuni,r can be related, consider as a test case a folded protein, on which only the proximity information along the protein sequence was used to determine the upper and lower bound distances between residues.
the constant upper limit for se-quential proximity being constant  for all residues i, the upper uij and lower lij bounds between a given pair of residues i and j, would be given by:
uij=u|i-j|
lij= <dig>
as the only sequential information is considered for the calculation of lower and upper bounds, no overall variation of the bounds is introduced when a restraint is added.
thus, the quantity of information iuni,r can be reduced to the uncertainty term hij  of the distance between residues i and j. using equations  <dig>   <dig> and  <dig>  iuni,r can be written as:
iuni,r= log
the information iuni,r being related to the absolute difference value |i - j| also found in the co definition, a new definition of the contact order could then be written as:
co*=1ncontactnseq ∑rexp∕u
in co*, each term is proportional to the quantity of information brought by the corresponding restraint.
by analogy to the definition of co  <cit> , one can thus hypothesize that iuni,r values are correlated with the protein folding rates and folding energy removing the restraints of highest iuni,r during the nmr structure calculation will reduce co* and thus increase the probability of unfolding of the conformations calculated with this reduced set of restraints.
in the following, we will evaluate this tendency to instability by simulating md trajectory starting from nmr conformations calculated using reduced restraints sets.
<dig>  unique information and folding cooperativity
as said previously, the information brought by a nmr restraint depends on the redundancy of this restraint with respect to the other restraints defining the protein structure.
redundant restraints correspond to highly cooperative contacts, which may be involved in the formation of the transition state ensemble   <cit> .
on the other hand, the least redundant restraints correspond to contacts with low cooperativity, which are formed after the appearance of tse during the folding process.
reciprocally, when starting from the folded state, the removal of restraints displaying the largest iuni,r values should drive the protein conformation toward the tse.
in the frame of the funnel model of protein folding  <cit> , the reduction of conformational entropy from the funnel top to the bottom can be put in parallel with the increase of queen information brought by each distance restraint, the contribution of information being equivalent to a reduction of entropy in information theory  <cit> .
the removal of the most informative distance restraints induces an increase of conformational entropy, and thus a widening of the free-energy funnel, destabilizing the protein.
the nmr restraint set brings structural information, and, at a first sight, it seems paradoxical that such information can be related to protein unfolding, which depends on thermodynamics and kinetics parameters.
nevertheless, the lack of redundancy arises from the lack of noe restraints in some protein regions, due to larger amplitude of the protein internal dynamics, and thus to the conformational entropy of the protein.
nmr dynamical information was also proved in the past  <cit>  to give information about the thermodynamic aspects of protein stability.
a relationship was already put in evidence some years ago  <cit>  between the free energy changes related to hydrogen exchange, and the residue fluctuations predicted in the frame of the gaussian network model   <cit> .
the previous observation of this effect gives a strong support to the model proposed here.
indeed, the elastic network describing the mechanics of a protein in gnm can be put in parallel with the set of nmr distance restraints.
high redundancy of the nmr restraints involving the residue i corresponds in the frame of gnm to a large free energy contribution Δgi of the residue i upon distortion of the protein coordinates.
following the observations of  <cit> , the large free energy contributions Δgi were shown to be correlated to the unfolding penalties revealed by hydrogen-exchange nmr experiments.
the most redundant nmr restraints are thus the restraints whose disruption leads to the largest penalties for protein unfolding.
conversely, the disruption of the least redundant nmr restraints leads to lower unfolding penalties and is thus more prone to be disrupted early during the unfolding process.
the previous arguments are confirmed by the recent use  <cit>  of the gnm to give information on protein unfolding pathways.
in the work of su et al  <cit> , the contacts between protein residues, corresponding to the largest fluctuations calculated in the frame of gnm, are successively removed, and the corresponding variation of the protein contact map gives information about the possible protein unfolding pathways.
the large similarity between the method proposed here and the protocol applied by  <cit> , supports the use of nmr restraint redundancy to study protein unfolding.
<dig>
RESULTS
<dig>  effect of restraint redundancy on the nmr structure definition
the iuni,r values were calculated using queen  <cit>  on the 1s1n restraint list.
the distribution of these values is dominated by nine long-range restraints , distributed on different sh <dig> regions : three restraints involve the β <dig> or β <dig> strand, two restraints involve the rt loop, the β <dig>  β <dig> strand or the c terminal part, one restraint involves the distal hairpin or the n-src loop.
for each restraint partner, the protein region, the residue name and number and the hydrogen name are given.
the hydrogens without stereospecific assignment are marked with a sign #.
residues l <dig>  l <dig> and v <dig> as well as residues a <dig>  k <dig>  d <dig> and e <dig>  y <dig> are in-volved in these nine restraints.
residue  <dig> is the mutated residue and l <dig>  and v <dig> are its neighbors in the sequence, e <dig>  y <dig> and a <dig>  k <dig>  d <dig> being the neighboring residues in 3d space.
the mutated residue l <dig>  as well as its neighbors are involved in the least redundant restraints which were described in the section "unique information and folding cooperativity" to define the last folding steps.
it is thus not surprising that the perturbation of these residues by mutation, was shown by nmr  <cit>  to have a large influence on the protein stability.
the relationship between the iuni,r values and the distribution of restraints in the protein structure was explored.
two restraint sets were used:  the 1s1n nmr restraints with lower and upper bounds set to  <dig>  and  <dig>  Å,  a synthetic set of restraints between all hydrogen pairs closer than  <dig> Å in the 1s1n pdb structure, with lower and upper bounds of  <dig> and  <dig> Å. in set , the nine restraints displaying the best iuni,r values includes five restraints already detected during the iuni,r calculation on the 1s1n restraint list.
the redundancy of nmr restraints thus depends only partially on the distance bound values.
on the other hand, in the restraint set , all values iuni,r are null.
indeed, in the set  of restraints, for each couple of spins k, j for which a restraint r is considered, there are three other spins i and l and m which are connected to each other and connected to i and j by a restraint.
this set of spins i, j, k, l and m form a rigid pentahedron.
.
if the connection between k and j is removed, the pentahedron is still rigid , and this implies that the upper and lower bounds of other connections are not changing.
thus, the information brought by the connection between k and j is equal to zero.
from the comparison of the two calculations run on sets  and , queen thus seems to be the most sensitive to the restraints missing with respect to the full triangulation of the structure.
three series of aria conformers were calculated using three sets of nmr restraints:  the set full including all 1s1n nmr restraints,  the set reduced where the nine least redundant restraints described above were removed,  the set random where nine randomly picked-up restraints were removed.
the comparison of the aria conformers obtained from restraint sets full and reduced reveals that the removal of the nine least redundant restraints reduces slightly the convergence , but the quality scores  are similar for both conformers sets.
besides, the backbone mean rmsds cross-calculated between the conformers obtained from both restraints sets, are of  <dig>  ±  <dig>  Å, only slightly larger than the rmsd values found among the two sets of con-formers .
thus, the removal of least redundant restraints changes only slightly the obtained nmr structure.
unsurprisingly, the aria conformations obtained with the set random display features  quite similar to the ones observed on the conformations calculated with the set full.
procheck: percentage of residues located in core region, allowed region, and dis-allowed region of the ramachandran diagram.
the following whatif z scores are given: 2nd generation packing quality , ramachandran plot appearance , the χ1/χ <dig> rotamer normality , backbone conformation  and the number of bumps .
the convergence is analyzed using the coordinate rmsd value between the backbone and the heavy atoms, the fitting to the data is estimated from the number of restraint violations larger than  <dig>  Å and the rms of restraint violations.
the analysis of the redundancy of nmr restraints measured on the nph sh <dig> domains shows that the least redundant restraints involve the mutated residue as well as its neigh-bors.
the least redundant restraints are also related to the lack of structure triangulation.
nevertheless, the removal of the least redundant restraints does not have large influence on the definition of the calculated protein structure.
<dig>  the removal of the least redundant restraints induces β barrel instability during md trajectories
the conformational drift of the wt nph sh <dig> domain was compared during the md simulations starting from nmr conformers generated using three nmr restraint sets:  full ,  reduced  and  random .
the md simulations based on the set full display coordinates rmsd in the range 1- <dig> Å , whereas some of the runs based on the sets reduced  and random  display larger drifts up to  <dig> Å.
the distances between centers of mass of the secondary structure elements of the nph sh <dig> domain were computed through md simulations based on the sets full, reduced and random .
for the sets reduced  and random , the distances increase with respect to the observations made for the set full .
these variations of protein architecture agree with the observations made on the nph sh <dig> conformational drift.
the removal of nmr restraints, whatever they are randomly chosen or displaying the least redundancy, have a destabilizing effect on the sh <dig> structure.
nevertheless, to complete this picture, structural parameters of the sh <dig> β barrel ge-ometry were analyzed.
in the neighborhood of the residue  <dig>  several distances were monitored between the residues sidechains involved in the sh <dig> hydrophobic core: y <dig>  l <dig>  v <dig>  w <dig>  a <dig>  v <dig> and l <dig>  the mean distances between residues of the hy-drophobic core are similar whether md simulations are based on the sets of restraints full or random .
thus, if randomly chosen restraints are removed, no change of the hydrophobic packing is observed near residue  <dig>  on the other hand, if the set reduced is used , the distances v50/v <dig>  l55/v <dig>  y5/v <dig>  a40/l <dig> and v29/l <dig>  increase.
moreover, an original method was developed , to monitor the profile of the inner cavity of the β barrel  .
the lengths of the β barrel are increased in the set reduced  compared to the set full .
similar profiles of the cylinder radii are observed  for all simulations, in the range 3- <dig> Å, for the cylinder axis interval  Å and larger values in the range 5- <dig> Å, outside this interval.
the barrel hoop oscillates: a minimum of  <dig>  Å is observed at the origin, two maxima of about  <dig> Å are located around - <dig> and  <dig> Å, two minima smaller than  <dig> Å are located around - <dig> and  <dig> Å, the maxima located outside the cylinder axis interval  Å being less well defined.
the maxima observed for the radii profile arise from residues of different protein regions: n-src loop and c terminal part , n-src loop and α helix , rt loop and distal hairpin .
the relative positions of the valves formed by n-src loop/c terminal part, n-src loop/α helix and rt loop/distal hairpin are thus essential to define the profile of the β barrel inner cavity.
the inner cavity profiles are plotted as black curves for the md simulations of the wt nph sh <dig> domain .
one curve is obtained for each simulation by interlacing the values computed on each run.
the profiles observed for the simulations based on the restraint sets full  and random , have smaller amplitudes than the profile computed from the simulations based on the set reduced .
the β barrel is thus more unstable if the least redundant restraints are removed than in the case of full or randomly-reduced restraint sets.
the simulations starting from nmr conformations obtained with reduced restraint set, display a similar global conformational drift and architecture variation, either the least redundant or randomly picked up restraints are removed.
nevertheless, the removal of the least redundant restraints produces a specific alteration of the β barrel shape and of the packing of hydrophobic residues close to the mutated residue  <dig>
<dig>  the instability of nph sh <dig> mutants is enhanced by using a reduced nmr restraints set
the conformational drift increases for the mutated sequences l28a and l28p with respect to the wt sequence, for the sets full  and reduced .
the increase of conformational drift for mutants is in qualitative agreement with the experimental observations on sh <dig> stability  <cit> .
the comparison of regression lines, drawn as black lines, proves that similar conformational drifts are observed for a given mutated sequence whatever the set of restraints used.
but, the coordinate rmsd gives only a global description of the structural variations.
the mutations mostly increase distances between centers of mass  what-ever the restraint set full or reduced is used.
however, the ranges of distances variations  are enlarged when the set reduced is used, and the largest variation ranges  are observed for the distances  involving the n and c terminal regions or the rt loop.
as among the  <dig> restraints involving residues from the rt loop, only two are removed to produce the reduced set of restraints, the motion of the rt-loop apart from the core is not induced by a lack of restraints.
in the mutant proteins, the sidechains of residues neighboring the residue  <dig> and involved in the sh <dig> hydrophobic core   move apart from each other for the set reduced .
on the other hand, for the set full , the distances between sidechains display various be-haviors.
the homogeneous displacement observed for the reduced restraint set  reveals a larger packing disruption around position  <dig>  the inter-residue distances increase from wt to l28a and l28p sequences.
the β barrel shape was analyzed along the β axis  by quantifying the barrel inner cavity .
the β barrel length  increases by about  <dig> Å if the least redundant nmr restraints are removed.
moreover, the β barrel profiles vary more for the simulations based on the set reduced  than for the simulations based on the set full .
the removal of the least redundant restraints enhances the profile variations observed for the protein mutants.
the profile variability is concentrated around the local maxima and minima described in the previous section, and corresponding to the barrel valves: rt loop, n-src loop, distal hairpin.
the β barrel is thus destabilized by the changes of the barrel hoops inducing the disjunction of barrel staves.
the trajectories recorded on the nph sh <dig> mutants reveal tendencies of the rt-loop to move apart from the structure core.
the n and c terminal regions tend to move away from each other, and the packing of the core residues located in the vicinity of position  <dig> tends to be disrupted.
the mutant destabilizing effect is enhanced in trajectories based on the reduced set of nmr restraints.
<dig>  the variations of the nph sh <dig> conformation agrees with previous information on sh <dig> unfolding
the folding of several sh <dig> domains were studied in the literature by extensive molecular dynamics simulations  <cit> , nmr relaxation measurements  <cit> , and mutational analysis  <cit> .
a transition state ensemble   <cit>  similar for the protein unfolding  <cit>  and folding  <cit>  transitions is attained by the drifts of the n and c terminal parts, followed by the opening of the rt loop.
the local conformational drifts and the distances between the centers of mass observed here depict conformational features in agreement with the information previously available for the sh <dig> tse.
the results obtained on nph sh <dig> will be analyzed in more details with respect to the studies performed by borreguero et al  <cit> , and ding et al  <cit> .
borreguero et al  <cit>  showed that a cluster of contacts between the rt loop and the distal hairpin may be stabilized by the hydrogen bonds e16/m <dig> and l18/m <dig>  in nph sh <dig>  the mean distances between amide hydrogens and carbonyls of l <dig> and l <dig> are smaller than  <dig>  Å for the wt and l28a simulations based on the set full.
on the contrary, the l28p mutation or the use of the set reduced disrupt these hydrogen bonds, with corresponding mean distances larger than  <dig>  Å.
ding et al  <cit>  monitored on src sh <dig> the probability of contacts formation among the conformations observed along protein refolding, at the tse, before the tse formation  and after the tse formation .
the probabilities of corresponding contacts formation in nph sh <dig>  were derived from the mean and the standard deviation values of the cβ-cβ distances, a contact being formed if the distance is smaller than  <dig>  Å. these probabilities were compared to the probabilities values given in tables  <dig> and  <dig> of ref.
<cit> , and, for  <dig> cβ-cβ pairs analyzed, the number of contacts post-ts, pre-ts or tse, displays a two-fold increase in the simulations based on the restraint set reduced.
the removal of least redundant restraints thus induces the nph sh <dig> conformation to drift toward the vicinity of the tse described for src sh <dig>
<dig>  sh <dig> energetic stability and water bridges
the energetic stability of the different nph sh <dig> sequences were compared by calculating the gbsa profiles on the sequence residues, as described in materials and methods.
the profiles of total gbsa energies are superimposed whatever is the set of restraints used , except for the mutated residue  <dig> and the residue l <dig>  which display less favorable gbsa energy in the mutants.
as the residue l <dig> is pointing toward the protein interior, its energetic instability supports the distortion of the protein core, in agreement with the previous observations on the β-barrel shape and on the hydrophobic core residues.
the instability of the sh <dig> core, in the vicinity of the mutation, is furthermore supported by the van der waals energy profiles in the gbsa energy.
indeed, these profiles show a variability among the  <dig> runs of each simulation , and this variability is increased for residues located in the vicinity of l <dig>  other protein residues in the region 7- <dig> corresponding to a part of the rt loop, and in the region 45- <dig>  corresponding to the last β strand, show also variations of van der waals energies.
an increase of variability and a broader distribution along the sequence are observed if the set reduced  is used: this is in agreement with the increased distortion of the β barrel profile observed for the same restraints set .
conversely, the observations made on the van der waals profiles prove that the origin of the β-barrel distortion is effectively the disruption of hydrophobic interactions.
the md simulations in explicit solvent make possible the observation of the behavior of water molecules.
the water bridges corresponding to water molecules in which two atoms are closer than  <dig>  Å of a protein hydrogen bond donor or acceptor group, and present in more than 40% of the md simulation time, are located around two groups of hub residues :  d <dig>  d <dig>  d <dig>  colored in magenta in the rt loop, and  e <dig>  e <dig>  e <dig>  r <dig>  colored in cyan on the opposite side of the protein, involving residues from the n-src loop, the α helix and the strand β <dig>  and close to the mutated residue l <dig>
the detected water bridges are plotted as edges of a network relating the interacting protein residues.
for the sets full  and reduced , the most branched network is observed on the wt protein .
the less connected network, observed on protein mutants, points out the reduced cooperativity between water molecules, in agreement with the lower stability of the mutants.
<dig> discussion
the hierarchy of nmr restraints defined by the redundancy was used to produce three sets of starting points for md simulations, these starting points being nmr conformers calculated with the sets full, reduced and random.
all sets of nmr conformations display similar structural and convergence features, and thus equivalent enthalpic properties.
the protein fold appears to be stable during the md runs recorded starting from wt sh <dig> conformers generated with the full set  of restraints.
this proves that the force field ff99sb  <cit>  used for md simulations, is accurate enough to keep stable the structure of the wt nph sh <dig> domain, in agreement with the nmr experimental observations  <cit> .
the removal of the least redundant nmr restraints ) or of randomly chosen restraints ) induces similar conformational drifts of the protein.
however, a closer analysis of the simulations reveals that the β barrel architecture and the packing of hy-drophobic residues close to the mutated residue l <dig>  are more destabilized if the reduced restraint set is used.
the hierarchy of protein contacts, derived from the redundancy of the nmr restraints, is thus related to the funnel  <cit>  shape of free energy around the native state  <cit> .
the mutated residue l <dig> and its neighbors are involved in the least redundant nmr restraints.
the destabilizing effect of the l <dig> mutations can be thus predicted from the mutual influence of protein topology and nmr restraint geometry, and this could be a general trend of protein unfolding.
more physically, the less branched networks of water bridges found for the mutated sequences indicate that the mutations disturb the protein hydration.
this hydration instability may arise from the expansion of the β barrel as well as from the distortion of the local hydrophobic interactions, shown in the gbsa profiles, and in the packing defects of hydrophobic residues close to residue  <dig>
the use of nmr restraint redundancy to destabilize nph sh <dig> underscores the importance of the protein topology  <cit>  and of the sequence  <cit> .
indeed, the effect of restraint removal depends mainly on topology, as the wt sh <dig> displays β barrel instability when the least redundant restraints are removed.
on the other hand, mutational effects on the protein stability add to topological effect regardless of whether the full  or the reduced  restraint set is used.
the removal of the least redundant nmr restraints allows one to enhance protein instability in very short  <dig> ns md trajectories, recorded in explicit solvent and using an all-atom force-field.
the work described here is certainly a preliminary one, and the results should be confirmed on other proteins.
nevertheless, the previous use of nmr data-sets to obtain information on thermodynamics of protein stability  <cit>  and the relationship of the proposed method with the gaussian network model  <cit> , gives strong additional support to the protocol proposed here.
the simplified approaches, as the gō model  <cit> , were widely proved to be essential in the exploration of protein unfolding pathways  <cit> .
the approach proposed here does not allow the exploration of the full protein unfolding pathways, but rather concentrates on the initial steps of this pathway.
on the other hand, as the force field used here is not a simplified one, it allows higher precision in the description of the protein instability.
in particular, the use of explicit solvent allows one to analyze the role of water molecules in the protein destabilization.
this approach can be applied to any protein structure for which experimental nmr distance restraints are available.
restraint redundancy depends mainly on their distribution in the structure, and is thus indirectly related to the internal protein mobility.
the use of the nmr restraint redundancy to explore protein stability is thus similar to sample rare fluctuations in native protein states from the use of nmr hydrogen exchange information  <cit> .
<dig>
CONCLUSIONS
the instability of wt and mutants of nph sh <dig> domain was analyzed using md simulations in explicit solvent.
the profile of the β barrel inner cavity was found to be more sensitive to the md starting points than the global conformational drift.
the present work also permitted to validate a method for destabilizing the conformations of nmr protein structures without biasing the free energy surface.
indeed, the removal of the least redundant nmr restraints  <cit>  allows one to destabilize the protein over short md trajectories.
besides, the stability difference between wt protein and mutants, is enhanced by the removal of the least redundant restraints, showing that the mutant instability arises from the β barrel shape alteration coupled to an increase of hydrophobic packing defects and to a less connected network of water bridges.
<dig> abbreviations
md: molecular dynamics; nmr: nuclear magnetic resonance; wt: wild-type; nph: nephrocystin; noe: nuclear overhauser effect
<dig> authors' contributions
nds calculated the nmr conformers, performed the molecular dynamics simulations and most of the analyses.
lm implemented the analysis of the β-barrel profile, as a constrained smooth optimization problem.
mn gave conceptual and technical advice and implemented the structure calculation protocols.
tem designed the study and coordinated the project.
tem, nds and lm wrote the paper.
all authors read and approved the final manuscript.
supplementary material
additional file 1
analysis of the β barrel geometry.
the analyses was conducted in two steps :  determination of the main cylinder axis,  determination of local cylinders along the main axis, describing the profile of the β barrel inner cavity.
click here for file
additional file 2
distribution fo the iuni,r information obtained with queen as a function of the restraint index r. the iuni,r values of the nine least redundant restraints are colored in red.
click here for file
additional file 3
case study of the pentahedron.
scheme of the pentahedron defined by all possible distance nmr restraints between the spin nuclei i, j, k, l and m. the distance restraints are drawn as lines.
variation of the distance restraints in the pentahedron if the restraint between k and j is removed.
click here for file
additional file 4
mean distances  observed between residues of the hydrophobic core of nph sh <dig>  in the md trajectories.
they are based on the  the full set of nmr restraints of the pdb entry 1s1n,  the reduced set of restraints where the least redundant restraints were removed,  randomly-reduced sets of restraints, where randomly picked-up restraints were re-moved.
the distance values are plotted as bars, colored in black, red and green, for the wt, l28a and l28p sequences.
the distances were calculated between sidechain carbons of the residues quoted in abscissa.
v50-cγ2/v29-cγ <dig>  l55-cδ1/v29-cγ <dig>  y5-cδ1/v29-cγ <dig>  a40-cβ/l27-cδ <dig>  v29-cγ2/l27-cδ <dig>
click here for file
<dig> acknowledgements
the authors thank the cnrs and the institut pasteur for funding.
lm thanks the brazilian institution capes for a postdoctoral fellowship.
prof. frédéric dardel is acknowledged for fruitful discussions.
dr. ruth nussinov is deeply acknowledged for wise remarks about the manuscript conceptual organization.

 the analysis of protein folding by computational methods faces the enormous confor-mational space to be explored, and usually requires simulation lengths of at least several hundreds of ns. simplified models, higher simulation temperature or implicit solvation can accelerate the process, but induce a deformation of the free energy surface. here, we propose a method based on the analysis of the nmr experimental restraints, to enhance the structure instability observed during explicit solvent molecular dynamics simulations with an all-atom force field. each nmr restraint is ranked according to its redundancy with respect to the other restraints. calculation of nmr protein con-formations is then performed, based on a reduced set of restraints, from which the least redundant restraints are removed. these conformations are modified starting points for md simulations having more chances to display structure instability. md simulations are recorded in the absence of any distance restraints, as the different restraint sets are only used for the generation of the starting points. sh domain of nephrocystin, a β barrel protein associated with juvenile nephronophthisis, was used to illustrate this approach. nph sh mutations associated to nephronophthisis, are located at the residue l in the strand β. the structure of wt nph sh was determined by nmr and the nmr analysis of the wt protein and of point mutants revealed that l28p is unfolded, whereas l28a is only partially folded. the impact of these mutations on the protein structure is quite intriguing, as residue l is pointing toward the solvent.. drawn in cartoon with β strands colored in magenta and the l colored in green and drawn in sticks. topology of the structure of the sh domain of the nephrocystine extracted from pdbsum http://www.ebi.ac.uk/pdbsum/, with l written in green. the structure of the nph sh domain reveals a classical β-barrel fold, surrounded by three loops: the rt loop, the n-src loop and the distal hairpin. this figure was produced with pymol the effect of removing the least redundant restraints was tested on wt nph sh by using three nmr restraints sets: the set full including all measured nmr restraints, the set reduced where few least redundant restraints were removed from the set full, the sets random, where restraints were randomly removed from the set full. three series of nmr conformations were calculated with these restraint sets, and display only small differences. notably, the md trajectories initiated from these different conformation series, display different trends. similar β barrel shapes were observed during the trajectories starting from nmr conformations calculated with the sets full and random, whereas β barrel instability is found during trajectories started from the conformers calculated using the set reduced. the effect of removing the least redundant restraints was also investigated on the nph sh mutants. the leucine sidechain was substituted by proline or alanine sidechains in the nmr conformers calculated using the restraint sets full and reduced. all md trajectories started from these mutated conformations are less stable than the corresponding md trajectories recorded on wt protein, but the instability is enhanced if starting points obtained with the set reduced are used. the redundancy of each nmr restraint was evaluated through the calculation of the information it brings to the total nmr restraint set. the more information is brought by a restraint, the less redundant it is from the others. conversely, low-information restraints depend a lot on the other restraints and are thus redundant. the information brought by each restraint was calculated using the approach queen, introduced some years ago by nabuurs et al. the instability features observed on nph sh domain are in agreement with previous experimental and simulation information on the unfolding of sh domains, as the rt loop and the n and c terminal parts are the most destabilized regions. the changes in β barrel architecture are revealed by an original geometric analysis and are related to a packing defect of the hydrophobic core residues induced by the mutation of l methods nmr conformers and restraints the nmr restraints measured on the nph sh domain were obtained from the pdb entry 1s1n. the approach queen was used to calculate the information brought by each restraint. from this information measure, the redundancy of each restraint with respect to the others was estimated. three sets of restraints were then prepared: the set full, including all restraints from 1s1n.mr, the set reduced, determined by discarding the nine least redundant restraints displaying the largest queen information values, the set random, determined by discarding nine randomly chosen restraints, displaying low queen information in the range 0-. these three sets of restraints were used to produce the corresponding sets, and of nmr conformations for the wt nph sh domain. nmr conformers were generated using cns. the simulated annealing procedure of cns was embedded into aria. a single iteration in geometric force field produced conformers, from which the lowest-energy conformers were submitted to a water refinement step. the nmr conformers were evaluated through: the coordinates rmsd between the backbone atoms, the number of violated distance restraints, the rms of restraints violations, and the procheck v. and whatif quality scores. molecular dynamics simulations nmr conformers of the nph sh were starting points of the molecular dynamics simulations, performed with explicit solvent, and periodic boundary conditions at constant temperature and pressure. a cutoff distance of Å was used to determine the water box size. trajectories were recorded using the package amber and the force field ff99sb. a cutoff of Å was used for lennard-jones interactions, and long-range electrostatic interactions were calculated with the particle mesh ewald protocol. the systems were neutralized using sodium counter-ions. pressure was regulated with isotropic position scaling and a relaxation time of ps, and temperature was regulated using a langevin thermostat with a collision frequency of ps- or a berendsen thermostat with a coupling time of ps. the shake algorithm was used to keep rigid covalent bonds involving hydrogens, enabling a time step of fs. atom coordinates were saved at every ps. simulations were initiated by rounds of semi-restrained and unrestrained system minimization. the systems were thermalized to 298k for ps at constant volume, while restraining the positions of the solute with a force constant of kcal/. the following equilibration protocol was then performed: md round of ps at constant volume and md rounds of ps were run while reducing the position restraints from kcal/ down to kcal/; eventually a last md round of ps was per-formed with a restraint of kcal/. nmr conformers were selected from the conformers sets,, described in subsection "nmr conformers and restraints", as starting points for md simulations. for each conformer set, runs of ns md trajectories were recorded varying: the thermo-stat controlling the system temperature, the seed for the initial velocities generation and the starting point chosen among the nmr conformer set. the mutations l28p and l28a were introduced in the wt nmr conformers calculated using the restraint sets full and reduced. the generalized born surface area approach was used to determine energy profiles per residues on the protein sequence. the generalized born energy was calculated using the onufriev method and the following parameters: no counter-ions, a solvent dielectric constant of and a protein dielectric constant of the surface area was computed by recursively approximating a sphere around an atom, starting from an icosahedra. the β barrel architecture an original geometric approach was used to analyze the sh β-barrel architecture. first, the cylinder of maximum radius such that no protein backbone atom can be located inside was determined to give the approximate main axis of the barrel. second, a series of slabs was defined along the main axis, and, within each slab, the maximum radius of the cylinder that can be fitted inside the volume defined by the protein backbone and the slab limits was calculated. as the slab width is Å and the slab center is moved by steps of Å, the cylinder radius obtained for each slab gives a local measure of the size of the inner β-barrel cavity, and the profile of the cylinder radius along the slabs thus provides an accurate and objective picture of the β barrel shape. as only the backbone atoms of the protein are included in the optimization and fitting procedures, the analysis focuses on the geometry of the β barrel, and considers the protein interior as a cavity. the maximization of cylinder radius was implemented as a constrained smooth op-timization problem, and solved with the algencan solver. analytic derivatives were computed. multiple initial points were considered to guarantee convergence to global solutions. this aspect is particularly relevant for the determination of the main axis, since, within each slab, the radius maximization is straightforward. the algen-can solver is freely available at http://www.ime.usp.br/~egbirgin/tango and the package used in the present work is available at http://lm-utils.googlecode.com. theory the restraint information in a protein structure determination by nmr, the relationship between protein atomic coordinates and distance restraints can be described in the frame of the distance geometry approach. the ensemble of protein atomic coordinates is considered as a 3d euclidean object. consequently, the upper and lower bounds of distance restraints verify the triangle inequalities, and the addition of a distance restraint potentially modifies the upper and lower bounds of all other restraints. as each restraint keeps the distance value into the lower/upper bound interval, it brings information improving the definition of the protein structure. not all restraints bring the same quantity of information, as they do not have the same influence on the structure definition. in other words, low-information restraints are quite redundant and thus very dependent on the other restraints, whereas high-information restraints are much less redundant with respect to the other restraints defining the structure. these properties of the nmr structure determination were known for a long time, but they were recently expressed in a quantitative way in the queen approach, by proposing that the distance restraints defining a structure could be analyzed in the frame of information theory. the uncertainty hij of the distance between spin nuclei i and j, can be expressed as: hij=-∫ lijuij1uij-lij log1uij-lijddij= log where lij and uij are respectively the lower and the upper bounds of the restraint. the uncertainty of a given set r of restraints can be then expressed as: hstructure|r=1n∑i,j=1nhij where n is the number of spin nuclei. therefore, the unique information iuni,r of a restraint r, applied between spins nuclei i and j, is defined as the difference of uncertainties between the sets r and r-1: iuni,r=hstructure|r-1-hstructure|r where r- represent the set of restraints others than r in the data-set, and r is the total set of restraints including r. the value iuni,r is different from the uncertainty hij of the distance between spin nuclei i and j, as the addition of a restraint may induce a variation of all upper and lower bounds by application of the triangle inequalities. the present work intends to use the queen framework quantifying the redundancy of each restraint, in order to destabilize protein structure along its unfolding path. unique information and contact order the unique information iuni,r of the queen approach can be related to the protein folding through the contact order parameter. the contact order co is defined from the ranges of the inter-residues contacts: co=1ncontactnseq ∑k,l|k-l| where k and l are the residue numbers forming a native contact in the folded protein, ncontact is the number of native contacts and nseq is the sequence size. the parameter co is positive and smaller than the parameter co was shown to be correlated to the folding rates, because of the decrease of conformational entropy associated with the formation of contacts. to see how co and iuni,r can be related, consider as a test case a folded protein, on which only the proximity information along the protein sequence was used to determine the upper and lower bound distances between residues. the constant upper limit for se-quential proximity being constant for all residues i, the upper uij and lower lij bounds between a given pair of residues i and j, would be given by: uij=u|i-j| lij= as the only sequential information is considered for the calculation of lower and upper bounds, no overall variation of the bounds is introduced when a restraint is added. thus, the quantity of information iuni,r can be reduced to the uncertainty term hij of the distance between residues i and j. using equations and iuni,r can be written as: iuni,r= log the information iuni,r being related to the absolute difference value |i - j| also found in the co definition, a new definition of the contact order could then be written as: co*=1ncontactnseq ∑rexp∕u in co*, each term is proportional to the quantity of information brought by the corresponding restraint. by analogy to the definition of co, one can thus hypothesize that iuni,r values are correlated with the protein folding rates and folding energy removing the restraints of highest iuni,r during the nmr structure calculation will reduce co* and thus increase the probability of unfolding of the conformations calculated with this reduced set of restraints. in the following, we will evaluate this tendency to instability by simulating md trajectory starting from nmr conformations calculated using reduced restraints sets. unique information and folding cooperativity as said previously, the information brought by a nmr restraint depends on the redundancy of this restraint with respect to the other restraints defining the protein structure. redundant restraints correspond to highly cooperative contacts, which may be involved in the formation of the transition state ensemble. on the other hand, the least redundant restraints correspond to contacts with low cooperativity, which are formed after the appearance of tse during the folding process. reciprocally, when starting from the folded state, the removal of restraints displaying the largest iuni,r values should drive the protein conformation toward the tse. in the frame of the funnel model of protein folding, the reduction of conformational entropy from the funnel top to the bottom can be put in parallel with the increase of queen information brought by each distance restraint, the contribution of information being equivalent to a reduction of entropy in information theory. the removal of the most informative distance restraints induces an increase of conformational entropy, and thus a widening of the free-energy funnel, destabilizing the protein. the nmr restraint set brings structural information, and, at a first sight, it seems paradoxical that such information can be related to protein unfolding, which depends on thermodynamics and kinetics parameters. nevertheless, the lack of redundancy arises from the lack of noe restraints in some protein regions, due to larger amplitude of the protein internal dynamics, and thus to the conformational entropy of the protein. nmr dynamical information was also proved in the past to give information about the thermodynamic aspects of protein stability. a relationship was already put in evidence some years ago between the free energy changes related to hydrogen exchange, and the residue fluctuations predicted in the frame of the gaussian network model. the previous observation of this effect gives a strong support to the model proposed here. indeed, the elastic network describing the mechanics of a protein in gnm can be put in parallel with the set of nmr distance restraints. high redundancy of the nmr restraints involving the residue i corresponds in the frame of gnm to a large free energy contribution Δgi of the residue i upon distortion of the protein coordinates. following the observations of, the large free energy contributions Δgi were shown to be correlated to the unfolding penalties revealed by hydrogen-exchange nmr experiments. the most redundant nmr restraints are thus the restraints whose disruption leads to the largest penalties for protein unfolding. conversely, the disruption of the least redundant nmr restraints leads to lower unfolding penalties and is thus more prone to be disrupted early during the unfolding process. the previous arguments are confirmed by the recent use of the gnm to give information on protein unfolding pathways. in the work of su et al, the contacts between protein residues, corresponding to the largest fluctuations calculated in the frame of gnm, are successively removed, and the corresponding variation of the protein contact map gives information about the possible protein unfolding pathways. the large similarity between the method proposed here and the protocol applied by, supports the use of nmr restraint redundancy to study protein unfolding.  effect of restraint redundancy on the nmr structure definition the iuni,r values were calculated using queen on the 1s1n restraint list. the distribution of these values is dominated by nine long-range restraints, distributed on different sh regions: three restraints involve the β or β strand, two restraints involve the rt loop, the β β strand or the c terminal part, one restraint involves the distal hairpin or the n-src loop. for each restraint partner, the protein region, the residue name and number and the hydrogen name are given. the hydrogens without stereospecific assignment are marked with a sign #. residues l l and v as well as residues a k d and e y are in-volved in these nine restraints. residue is the mutated residue and l and v are its neighbors in the sequence, e y and a k d being the neighboring residues in 3d space. the mutated residue l as well as its neighbors are involved in the least redundant restraints which were described in the section "unique information and folding cooperativity" to define the last folding steps. it is thus not surprising that the perturbation of these residues by mutation, was shown by nmr to have a large influence on the protein stability. the relationship between the iuni,r values and the distribution of restraints in the protein structure was explored. two restraint sets were used: the 1s1n nmr restraints with lower and upper bounds set to and Å, a synthetic set of restraints between all hydrogen pairs closer than Å in the 1s1n pdb structure, with lower and upper bounds of and Å. in set, the nine restraints displaying the best iuni,r values includes five restraints already detected during the iuni,r calculation on the 1s1n restraint list. the redundancy of nmr restraints thus depends only partially on the distance bound values. on the other hand, in the restraint set, all values iuni,r are null. indeed, in the set of restraints, for each couple of spins k, j for which a restraint r is considered, there are three other spins i and l and m which are connected to each other and connected to i and j by a restraint. this set of spins i, j, k, l and m form a rigid pentahedron.. if the connection between k and j is removed, the pentahedron is still rigid, and this implies that the upper and lower bounds of other connections are not changing. thus, the information brought by the connection between k and j is equal to zero. from the comparison of the two calculations run on sets and, queen thus seems to be the most sensitive to the restraints missing with respect to the full triangulation of the structure. three series of aria conformers were calculated using three sets of nmr restraints: the set full including all 1s1n nmr restraints, the set reduced where the nine least redundant restraints described above were removed, the set random where nine randomly picked-up restraints were removed. the comparison of the aria conformers obtained from restraint sets full and reduced reveals that the removal of the nine least redundant restraints reduces slightly the convergence, but the quality scores are similar for both conformers sets. besides, the backbone mean rmsds cross-calculated between the conformers obtained from both restraints sets, are of ± Å, only slightly larger than the rmsd values found among the two sets of con-formers. thus, the removal of least redundant restraints changes only slightly the obtained nmr structure. unsurprisingly, the aria conformations obtained with the set random display features quite similar to the ones observed on the conformations calculated with the set full. procheck: percentage of residues located in core region, allowed region, and dis-allowed region of the ramachandran diagram. the following whatif z scores are given: 2nd generation packing quality, ramachandran plot appearance, the χ1/χ rotamer normality, backbone conformation and the number of bumps. the convergence is analyzed using the coordinate rmsd value between the backbone and the heavy atoms, the fitting to the data is estimated from the number of restraint violations larger than Å and the rms of restraint violations. the analysis of the redundancy of nmr restraints measured on the nph sh domains shows that the least redundant restraints involve the mutated residue as well as its neigh-bors. the least redundant restraints are also related to the lack of structure triangulation. nevertheless, the removal of the least redundant restraints does not have large influence on the definition of the calculated protein structure. the removal of the least redundant restraints induces β barrel instability during md trajectories the conformational drift of the wt nph sh domain was compared during the md simulations starting from nmr conformers generated using three nmr restraint sets: full, reduced and random. the md simulations based on the set full display coordinates rmsd in the range 1- Å, whereas some of the runs based on the sets reduced and random display larger drifts up to Å. the distances between centers of mass of the secondary structure elements of the nph sh domain were computed through md simulations based on the sets full, reduced and random. for the sets reduced and random, the distances increase with respect to the observations made for the set full. these variations of protein architecture agree with the observations made on the nph sh conformational drift. the removal of nmr restraints, whatever they are randomly chosen or displaying the least redundancy, have a destabilizing effect on the sh structure. nevertheless, to complete this picture, structural parameters of the sh β barrel ge-ometry were analyzed. in the neighborhood of the residue several distances were monitored between the residues sidechains involved in the sh hydrophobic core: y l v w a v and l the mean distances between residues of the hy-drophobic core are similar whether md simulations are based on the sets of restraints full or random. thus, if randomly chosen restraints are removed, no change of the hydrophobic packing is observed near residue on the other hand, if the set reduced is used, the distances v50/v l55/v y5/v a40/l and v29/l increase. moreover, an original method was developed, to monitor the profile of the inner cavity of the β barrel. the lengths of the β barrel are increased in the set reduced compared to the set full. similar profiles of the cylinder radii are observed for all simulations, in the range 3- Å, for the cylinder axis interval Å and larger values in the range 5- Å, outside this interval. the barrel hoop oscillates: a minimum of Å is observed at the origin, two maxima of about Å are located around - and Å, two minima smaller than Å are located around - and Å, the maxima located outside the cylinder axis interval Å being less well defined. the maxima observed for the radii profile arise from residues of different protein regions: n-src loop and c terminal part, n-src loop and α helix, rt loop and distal hairpin. the relative positions of the valves formed by n-src loop/c terminal part, n-src loop/α helix and rt loop/distal hairpin are thus essential to define the profile of the β barrel inner cavity. the inner cavity profiles are plotted as black curves for the md simulations of the wt nph sh domain. one curve is obtained for each simulation by interlacing the values computed on each run. the profiles observed for the simulations based on the restraint sets full and random, have smaller amplitudes than the profile computed from the simulations based on the set reduced. the β barrel is thus more unstable if the least redundant restraints are removed than in the case of full or randomly-reduced restraint sets. the simulations starting from nmr conformations obtained with reduced restraint set, display a similar global conformational drift and architecture variation, either the least redundant or randomly picked up restraints are removed. nevertheless, the removal of the least redundant restraints produces a specific alteration of the β barrel shape and of the packing of hydrophobic residues close to the mutated residue the instability of nph sh mutants is enhanced by using a reduced nmr restraints set the conformational drift increases for the mutated sequences l28a and l28p with respect to the wt sequence, for the sets full and reduced. the increase of conformational drift for mutants is in qualitative agreement with the experimental observations on sh stability. the comparison of regression lines, drawn as black lines, proves that similar conformational drifts are observed for a given mutated sequence whatever the set of restraints used. but, the coordinate rmsd gives only a global description of the structural variations. the mutations mostly increase distances between centers of mass what-ever the restraint set full or reduced is used. however, the ranges of distances variations are enlarged when the set reduced is used, and the largest variation ranges are observed for the distances involving the n and c terminal regions or the rt loop. as among the restraints involving residues from the rt loop, only two are removed to produce the reduced set of restraints, the motion of the rt-loop apart from the core is not induced by a lack of restraints. in the mutant proteins, the sidechains of residues neighboring the residue and involved in the sh hydrophobic core move apart from each other for the set reduced. on the other hand, for the set full, the distances between sidechains display various be-haviors. the homogeneous displacement observed for the reduced restraint set reveals a larger packing disruption around position the inter-residue distances increase from wt to l28a and l28p sequences. the β barrel shape was analyzed along the β axis by quantifying the barrel inner cavity. the β barrel length increases by about Å if the least redundant nmr restraints are removed. moreover, the β barrel profiles vary more for the simulations based on the set reduced than for the simulations based on the set full. the removal of the least redundant restraints enhances the profile variations observed for the protein mutants. the profile variability is concentrated around the local maxima and minima described in the previous section, and corresponding to the barrel valves: rt loop, n-src loop, distal hairpin. the β barrel is thus destabilized by the changes of the barrel hoops inducing the disjunction of barrel staves. the trajectories recorded on the nph sh mutants reveal tendencies of the rt-loop to move apart from the structure core. the n and c terminal regions tend to move away from each other, and the packing of the core residues located in the vicinity of position tends to be disrupted. the mutant destabilizing effect is enhanced in trajectories based on the reduced set of nmr restraints. the variations of the nph sh conformation agrees with previous information on sh unfolding the folding of several sh domains were studied in the literature by extensive molecular dynamics simulations, nmr relaxation measurements, and mutational analysis. a transition state ensemble similar for the protein unfolding and folding transitions is attained by the drifts of the n and c terminal parts, followed by the opening of the rt loop. the local conformational drifts and the distances between the centers of mass observed here depict conformational features in agreement with the information previously available for the sh tse. the results obtained on nph sh will be analyzed in more details with respect to the studies performed by borreguero et al, and ding et al. borreguero et al showed that a cluster of contacts between the rt loop and the distal hairpin may be stabilized by the hydrogen bonds e16/m and l18/m in nph sh the mean distances between amide hydrogens and carbonyls of l and l are smaller than Å for the wt and l28a simulations based on the set full. on the contrary, the l28p mutation or the use of the set reduced disrupt these hydrogen bonds, with corresponding mean distances larger than Å. ding et al monitored on src sh the probability of contacts formation among the conformations observed along protein refolding, at the tse, before the tse formation and after the tse formation. the probabilities of corresponding contacts formation in nph sh were derived from the mean and the standard deviation values of the cβ-cβ distances, a contact being formed if the distance is smaller than Å. these probabilities were compared to the probabilities values given in tables and of ref., and, for cβ-cβ pairs analyzed, the number of contacts post-ts, pre-ts or tse, displays a two-fold increase in the simulations based on the restraint set reduced. the removal of least redundant restraints thus induces the nph sh conformation to drift toward the vicinity of the tse described for src sh sh energetic stability and water bridges the energetic stability of the different nph sh sequences were compared by calculating the gbsa profiles on the sequence residues, as described in materials and methods. the profiles of total gbsa energies are superimposed whatever is the set of restraints used, except for the mutated residue and the residue l which display less favorable gbsa energy in the mutants. as the residue l is pointing toward the protein interior, its energetic instability supports the distortion of the protein core, in agreement with the previous observations on the β-barrel shape and on the hydrophobic core residues. the instability of the sh core, in the vicinity of the mutation, is furthermore supported by the van der waals energy profiles in the gbsa energy. indeed, these profiles show a variability among the runs of each simulation, and this variability is increased for residues located in the vicinity of l other protein residues in the region 7- corresponding to a part of the rt loop, and in the region 45- corresponding to the last β strand, show also variations of van der waals energies. an increase of variability and a broader distribution along the sequence are observed if the set reduced is used: this is in agreement with the increased distortion of the β barrel profile observed for the same restraints set. conversely, the observations made on the van der waals profiles prove that the origin of the β-barrel distortion is effectively the disruption of hydrophobic interactions. the md simulations in explicit solvent make possible the observation of the behavior of water molecules. the water bridges corresponding to water molecules in which two atoms are closer than Å of a protein hydrogen bond donor or acceptor group, and present in more than 40% of the md simulation time, are located around two groups of hub residues: d d d colored in magenta in the rt loop, and e e e r colored in cyan on the opposite side of the protein, involving residues from the n-src loop, the α helix and the strand β and close to the mutated residue l the detected water bridges are plotted as edges of a network relating the interacting protein residues. for the sets full and reduced, the most branched network is observed on the wt protein. the less connected network, observed on protein mutants, points out the reduced cooperativity between water molecules, in agreement with the lower stability of the mutants. discussion the hierarchy of nmr restraints defined by the redundancy was used to produce three sets of starting points for md simulations, these starting points being nmr conformers calculated with the sets full, reduced and random. all sets of nmr conformations display similar structural and convergence features, and thus equivalent enthalpic properties. the protein fold appears to be stable during the md runs recorded starting from wt sh conformers generated with the full set of restraints. this proves that the force field ff99sb used for md simulations, is accurate enough to keep stable the structure of the wt nph sh domain, in agreement with the nmr experimental observations. the removal of the least redundant nmr restraints ) or of randomly chosen restraints ) induces similar conformational drifts of the protein. however, a closer analysis of the simulations reveals that the β barrel architecture and the packing of hy-drophobic residues close to the mutated residue l are more destabilized if the reduced restraint set is used. the hierarchy of protein contacts, derived from the redundancy of the nmr restraints, is thus related to the funnel shape of free energy around the native state. the mutated residue l and its neighbors are involved in the least redundant nmr restraints. the destabilizing effect of the l mutations can be thus predicted from the mutual influence of protein topology and nmr restraint geometry, and this could be a general trend of protein unfolding. more physically, the less branched networks of water bridges found for the mutated sequences indicate that the mutations disturb the protein hydration. this hydration instability may arise from the expansion of the β barrel as well as from the distortion of the local hydrophobic interactions, shown in the gbsa profiles, and in the packing defects of hydrophobic residues close to residue the use of nmr restraint redundancy to destabilize nph sh underscores the importance of the protein topology and of the sequence. indeed, the effect of restraint removal depends mainly on topology, as the wt sh displays β barrel instability when the least redundant restraints are removed. on the other hand, mutational effects on the protein stability add to topological effect regardless of whether the full or the reduced restraint set is used. the removal of the least redundant nmr restraints allows one to enhance protein instability in very short ns md trajectories, recorded in explicit solvent and using an all-atom force-field. the work described here is certainly a preliminary one, and the results should be confirmed on other proteins. nevertheless, the previous use of nmr data-sets to obtain information on thermodynamics of protein stability and the relationship of the proposed method with the gaussian network model, gives strong additional support to the protocol proposed here. the simplified approaches, as the gō model, were widely proved to be essential in the exploration of protein unfolding pathways. the approach proposed here does not allow the exploration of the full protein unfolding pathways, but rather concentrates on the initial steps of this pathway. on the other hand, as the force field used here is not a simplified one, it allows higher precision in the description of the protein instability. in particular, the use of explicit solvent allows one to analyze the role of water molecules in the protein destabilization. this approach can be applied to any protein structure for which experimental nmr distance restraints are available. restraint redundancy depends mainly on their distribution in the structure, and is thus indirectly related to the internal protein mobility. the use of the nmr restraint redundancy to explore protein stability is thus similar to sample rare fluctuations in native protein states from the use of nmr hydrogen exchange information.  the instability of wt and mutants of nph sh domain was analyzed using md simulations in explicit solvent. the profile of the β barrel inner cavity was found to be more sensitive to the md starting points than the global conformational drift. the present work also permitted to validate a method for destabilizing the conformations of nmr protein structures without biasing the free energy surface. indeed, the removal of the least redundant nmr restraints allows one to destabilize the protein over short md trajectories. besides, the stability difference between wt protein and mutants, is enhanced by the removal of the least redundant restraints, showing that the mutant instability arises from the β barrel shape alteration coupled to an increase of hydrophobic packing defects and to a less connected network of water bridges. abbreviations md: molecular dynamics; nmr: nuclear magnetic resonance; wt: wild-type; nph: nephrocystin; noe: nuclear overhauser effect authors' contributions nds calculated the nmr conformers, performed the molecular dynamics simulations and most of the analyses. lm implemented the analysis of the β-barrel profile, as a constrained smooth optimization problem. mn gave conceptual and technical advice and implemented the structure calculation protocols. tem designed the study and coordinated the project. tem, nds and lm wrote the paper. all authors read and approved the final manuscript. supplementary material additional file 1 analysis of the β barrel geometry. the analyses was conducted in two steps: determination of the main cylinder axis, determination of local cylinders along the main axis, describing the profile of the β barrel inner cavity. click here for file additional file 2 distribution fo the iuni,r information obtained with queen as a function of the restraint index r. the iuni,r values of the nine least redundant restraints are colored in red. click here for file additional file 3 case study of the pentahedron. scheme of the pentahedron defined by all possible distance nmr restraints between the spin nuclei i, j, k, l and m. the distance restraints are drawn as lines. variation of the distance restraints in the pentahedron if the restraint between k and j is removed. click here for file additional file 4 mean distances observed between residues of the hydrophobic core of nph sh in the md trajectories. they are based on the the full set of nmr restraints of the pdb entry 1s1n, the reduced set of restraints where the least redundant restraints were removed, randomly-reduced sets of restraints, where randomly picked-up restraints were re-moved. the distance values are plotted as bars, colored in black, red and green, for the wt, l28a and l28p sequences. the distances were calculated between sidechain carbons of the residues quoted in abscissa. v50-cγ2/v29-cγ l55-cδ1/v29-cγ y5-cδ1/v29-cγ a40-cβ/l27-cδ v29-cγ2/l27-cδ click here for file acknowledgements the authors thank the cnrs and the institut pasteur for funding. lm thanks the brazilian institution capes for a postdoctoral fellowship. prof. frédéric dardel is acknowledged for fruitful discussions. dr. ruth nussinov is deeply acknowledged for wise remarks about the manuscript conceptual organization.
the present work aims to figure out whether nmr restraints data can be used to probe protein conformations in order to accelerate the unfolding simulation. the nmr restraints observed on the wt nph sh domain were sorted from the least redundant to the most redundant ones. protein nmr conformations were then calculated with: the set full including all nmr restraints measured on nph sh the set reduced where the least redundant restraints with respect to the set full were removed, the sets random where randomly picked-up restraints were removed. the β barrel architecture of nph sh in the trajectories starting from sets and appears to be stable. on the contrary, on trajectories based on the set, a displacement of the hydrophobic core residues and a variation of the β barrel inner cavity profile were observed. we conclude that the nmr restraint redundancy is connected to the instability of the sh nph domain. the relationship between the nmr restraint redundancy and the protein folding is also reminiscent of the previous use of the gaussian network model to predict protein folding parameters. this restraint redundancy generalizes the contact order parameter, which is calculated from the contact map of a folded protein and was shown in the literature to be correlated to the protein folding rate. the destabilizing effect of mutations was also found to be enhanced by the removal of the least redundant restraints.
 the simulation of protein unfolding usually requires recording long molecular dynamics trajectories. the present work aims to figure out whether nmr restraints data can be used to probe protein conformations in order to accelerate the unfolding simulation. the sh domain of nephrocystine was shown by nmr to be destabilized by point mutations, and was thus chosen to illustrate the proposed method.  the nmr restraints observed on the wt nph sh domain were sorted from the least redundant to the most redundant ones. protein nmr conformations were then calculated with: the set full including all nmr restraints measured on nph sh the set reduced where the least redundant restraints with respect to the set full were removed, the sets random where randomly picked-up restraints were removed. from each set of conformations, we recorded series of 5-ns md trajectories. the β barrel architecture of nph sh in the trajectories starting from sets and appears to be stable. on the contrary, on trajectories based on the set, a displacement of the hydrophobic core residues and a variation of the β barrel inner cavity profile were observed. the overall nph sh destabilization agrees with previous experimental and simulation observations made on other sh domains. the destabilizing effect of mutations was also found to be enhanced by the removal of the least redundant restraints.  we conclude that the nmr restraint redundancy is connected to the instability of the sh nph domain. this restraint redundancy generalizes the contact order parameter, which is calculated from the contact map of a folded protein and was shown in the literature to be correlated to the protein folding rate. the relationship between the nmr restraint redundancy and the protein folding is also reminiscent of the previous use of the gaussian network model to predict protein folding parameters. nmrprotein foldingsh domainmolecular dynamics simulationqueencontact order: gaussian network model 