BACKGROUND
predicting the sequences compatible with a given structure defines what is traditionally called the inverse folding problem, or more often, protein design  <cit> .
as suggested by the terminology, this question is usually considered in an engineering perspective: the aim is then to determine a sequence, or a set of sequences, that stably fold into a pre-specified conformation.
in a thermodynamic perspective, this requirement translates into eliciting sequences that have lowest free energy under the target fold, compared to all possible alternative conformations.
in principle, such a criterion would imply a search through the joint structure-sequence space, which is not feasible but for small on-lattice model proteins  <cit> .
as an alternative to the engineering approach, a more evolutionary stance can be taken towards the inverse folding problem, in which case the aim would rather be to predict the sequences of natural proteins having the conformation of interest.
seen from this new point of view, the design problem raises new questions: natural proteins are the result of a complex evolutionary process, involving an intricate interplay between mutation and selection, and this probably entails many constraints directly related to the native conformation, but nevertheless not equivalent to the mere requirement of structural stability.
for instance, the requirement of fast and cooperative folding has an impact on the dispersion of contact energies  <cit> .
for this and many other potential reasons, among all sequences predicted by classical engineering-oriented protein design, probably only a subset will look like natural proteins.
the evolutionary approach to protein design is particularly relevant to phylogenetic studies, where one of the current motivations is to develop the so-called structurally constrained models of protein evolution, i.e.
models explicitly dependent on the protein's conformation, either for simulation purposes  <cit> , or in the context of phylogenetic inference  <cit> .
in this framework, each substitution undergone by a protein during evolution has to be tested for its compatibility with the structure, in the context of the sequence that the protein displays at all other sites when the substitution occurs.
such repeated evaluation of the structure-sequence compatibility along a phylogenetic tree requires relevant and computationally very efficient scoring schemes/functions.
it is interesting to compare the different methods proposed thus far for performing protein design in light of this engineering/evolutionary distinction.
a first direction of research has consisted in using all-atom semi-empirical force fields to evaluate the conformational free energy .
these empirical methods have been applied to many theoretical and experimental cases, reaching a high level of accuracy.
on the other hand, they are computationally heavy, mainly because of the side-chain positioning problem, and thus cannot be easily applied to structurally constrained phylogenetic models  <cit> .
concerns may also be expressed about their over-sensitivity to the native conformation, in particular in the core of the target structures and when the flexibility of the backbone is not accounted for  <cit> .
but more importantly, approaches based on physical force fields are, by definition, exclusively focussed on the conformational stability, and thereby, completely oversee other potential factors shaping the sequences of biological proteins.
as such, they are well suited for engineering synthetic proteins  <cit> , or for testing to what extent natural sequences are shaped by selection for protein stability  <cit> , but may not be sufficient for more general evolutionary purposes.
an alternative to the semi-empirical strategy consists in relying on knowledge-based, or statistical, potentials.
these scoring functions mimic physical boltzmann distributions, but merely encode statistical patterns present in the databases.
some of these potentials were obtained under the quasi-chemical approximation, whereby frequencies of patterns, such as contacts between each pair of amino-acids, are transformed into energies using the boltzmann law  <cit> .
alternatively, contact energies can be obtained by maximizing the potential's predictive accuracy in a threading test  <cit> .
in the present context, an advantage of these knowledge-based potentials, compared to semi-empirical force-fields, is that they should in principle capture all kinds of patterns that true biological sequences have, in relation to their conformation, and not only those directly related to thermodynamic stability.
furthermore, statistical potentials need not be defined at the atomic level, but can be based on a coarse-grained description of the protein's configuration, essentially by omitting the degrees of freedom associated to side chains.
this allows faster computations, by avoiding the problem of searching through the rugged landscape of side-chain conformations.
in addition, coarse-grained potentials could turn out to be an advantage, in that they will not recover the native sequence too faithfully.
most protein design procedures based on statistical potentials proposed until now have relied on coarse-grained, pairwise contact pseudo-energies  <cit> .
yet, irrespective of the level of description adopted, currently available statistical potentials may not be ideal for protein design, since they have generally been optimized in the context of the folding problem, i.e.
for maximizing the rate of correct structure prediction, given the sequence.
in contrast, we would like to optimize the reciprocal prediction, namely, the sequences given the conformation.
several approaches have been proposed in this direction, consisting in maximizing the z-score between the energy of the native sequence on the target conformation and its energy on a set of decoy sequences  <cit> , or, alternatively, in applying a mean-square criterion on the values taken by the scoring function on each structure-sequence pair of the database  <cit> .
however, these methods have thus far only been tested in cubic lattice protein models.
in addition, they lack a firm theoretical basis.
in particular, it would be interesting to guarantee optimal predictive power, and to have a robust methodology available to assess and compare the performance of alternative forms of statistical potentials.
standard statistical theory provides such theoretical guarantees  <cit> .
in the present case, the inverse folding problem can be formulated directly in terms of the probability of observing a sequence s given a conformation c, i.e.
p. this probability explicitly depends on the pre-specified model through a series of parameters, represented here by θ. these may be, for instance, the coefficients of a pairwise potential, parameters describing compositional effects, secondary structure environment, solvent accessibility, etc.
taking the product over a database of p independent sequence-conformation pairs, s = p =  <dig> .p and c = p =  <dig> .p, yields a joint probability
p=∏pp
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgwbacdaqadaqaaiabdofatjabcyha8jabdoeadjabcycasggaciab=h7axbgaayjkaiaawmcaaiabg2da9maarafabagaemicaa3aaewaaeaacqwgzbwcdaahaawcbeqaaiabdchawbaakiabcyha8jabdogajnaacaaaleqabagaemicaahaaogaeiilawiae8hudehacagloagaayzkaaaaleaacqwgwbacaeqaniabg+givdgccawljagaaczcamaabmaabagaegymaedacagloagaayzkaaaaaa@4b42@
which, as a function of θ, can be seen as a likelihood.
the parameter θ is then learnt by maximizing the likelihood with respect to θ. once this is done, sequences can be assessed, or sampled, under the optimal parameter value θ^
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadaqiaaqaaggaciab=h7axbgaaypadaaaaa@2f2b@, by direct numerical evaluation of their probability, or by monte carlo sampling methods.
reformulated in this way, the method maximizes the predictive power of the potential, now in the structure-seeks-sequence direction.
by construction, it yields the optimal parameter values that can be obtained for a given form of the potential.
in addition, the fit of the model can be directly evaluated, based on the value of the likelihood obtained on a test data set, distinct from the learning set , giving a means of rigorous model selection.
finally, the statistical framework proposed here allows one to explicitly combine together, in a model dependent manner, all kinds of factors that we surmise may induce correlations between the structure and the sequence of proteins.
we have implemented this maximum likelihood  procedure in a markov chain monte carlo framework, and applied it to a simple case, using a contact potential, supplemented with a solvent accessibility term.
using cross-validation, we show that the resulting potentials yield a better fit than currently available potentials of the same form, and that combining solvent-accessibility considerations with contact energies is better than either alone.
furthermore, we find that solvent accessibility requires a more complex description than what is currently used.
ultimately, the overall method proposed in this work can be extended to a large spectrum of alternative models and statistical potentials.
RESULTS
the probabilistic model
let us consider a sequence s = i =  <dig> .n, of length n, and of conformation c. in its most general form, the method introduced here can work with any model m specifying the conditional probability of s given c, in terms of an unnormalized non negative function q:
p=q∑sq.
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgwbacdaqadaqaaiabdohazjabcyha8jabdogajjabcycasiabd2eanbgaayjkaiaawmcaaiabg2da9maalaaabagaemycae3aaewaaeaacqwgzbwccqggsaalcqwgjbwyaiaawicacaglpaaaaeaadaaeqaqaaiabdghaxnaabmaabagaem4camnaeiilawiaem4yamgacagloagaayzkaaaaleaacqwgzbwcaeqaniabgghildaaaogaeiola4iaaczcaiaaxmaadaqadaqaaiabikdayagaayjkaiaawmcaaaaa@4c18@
to illustrate the method, we will apply it to a simple case, using a pairwise contact potential.
the argument is as follows.
first, by bayes' theorem:
p=pp∑spp.
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgwbacdaqadaqaaiabdohazjabcyha8jabdogajjabcycasiabd2eanbgaayjkaiaawmcaaiabg2da9maalaaabagaemicaa3aaewaaeaacqwgjbwycqgg8bafcqwgzbwccqggsaalcqwgnbqtaiaawicacaglpaaacqwgwbacdaqadaqaaiabdohazjabcyha8jabd2eanbgaayjkaiaawmcaaaqaamaaqababagaemicaa3aaewaaeaacqwgjbwycqgg8bafcqwgzbwccqggsaalcqwgnbqtaiaawicacaglpaaacqwgwbacdaqadaqaaiabdohazjabcyha8jabd2eanbgaayjkaiaawmcaaawcbagaem4camhabeqdcqghris5aaaakiabc6cauiaaxmaacawljawaaewaaeaacqaizawmaiaawicacaglpaaaaaa@5f64@
if, in addition, we assume a uniform prior on s, we can simply relate equations  <dig> and  <dig> by posing q = p. next, given a statistical potential e, the conformational probability p can be expressed as a boltzmann distribution:
p=e−e/ktzs=e−−f)/kt,
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaafaqadegacaaabagaemicaa3aaewaaeaacqwgjbwycqgg8bafcqwgzbwccqggsaalcqwgnbqtaiaawicacaglpaaacqgh9aqpdawcaaqaaiabdwgalnaacaaaleqabagaeyoei0iaemyrau0aaewaaeaacqwgzbwccqggsaalcqwgjbwyaiaawicacaglpaaacqggvawlcqwgrbwacqwgubavaaaakeaacqwgabgwdawgaawcbagaem4camhabeaaaaaakeaadaqadaqaaiabisda0agaayjkaiaawmcaaaqaaiabg2da9iabdwgalnaacaaaleqabagaeyoei0yaaewaaeaacqwgfbqrdaqadaqaaiabdohazjabcycasiabdogajbgaayjkaiaawmcaaiabgkhitiabdaeagnaabmaabagaem4camhacagloagaayzkaaaacagloagaayzkaagaei4la8iaem4aasmaemivaqfaaogaeiilawcabawaaewaaeaacqai1aqnaiaawicacaglpaaaaaaaaa@5fa0@
where
zs=∑ce−e/kt
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgabgwdawgaawcbagaem4camhabeaakiabg2da9maaqafabagaemyzau2aawbaasqabeaacqghsislcqwgfbqrdaqadaqaaiabdohazjabcycasiabdogajbgaayjkaiaawmcaaiabc+caviabdugarjabdsfaubaaaeaacqwgjbwyaeqaniabgghildgccawljagaaczcamaabmaabagaegonaydacagloagaayzkaaaaaa@43e9@
is a normalization constant, and
f = - ln zs.
t and k are the absolute temperature and the boltzmann constant, respectively.
without loss of generality, it is possible to rescale the potential so that kt =  <dig>  which we will do in the following.
then, by defining the inverse potential:
g = e - f,
the conditional probability of sequence s reads as
p=e−gy,
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgwbacdaqadaqaaiabdohazjabcyha8jabdogajjabcycasggaciab=h7axjabcycasiabd2eanbgaayjkaiaawmcaaiabg2da9maalaaabagaemyzau2aawbaasqabeaacqghsislcqwghbwrdaqadaqaaiabdohazjabcycasiabdogajbgaayjkaiaawmcaaaaaaoqaaiabdmfazbaacqggsaalcawljagaaczcamaabmaabagaegyoakdacagloagaayzkaaaaaa@482f@
where
y=∑s′e−g
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgzbqwcqgh9aqpdaaeqbqaaiabdwgalnaacaaaleqabagaeyoei0iaem4rac0aaewaaeaacuwgzbwcgaqbaiabcycasiabdogajbgaayjkaiaawmcaaaaaaeaacuwgzbwcgaqbaaqab0gaeyyeiuoakiaaxmaacawljawaaewaaeaacqaixaqmcqaiwaamaiaawicacaglpaaaaaa@3fec@
is the normalization factor.
note that, contrary to the zs factor of equation  <dig>  which was a sum over all conformations, the present factor y is a sum over sequence space .
statistical potentials
in the present work, we used a statistical potential made of two terms:
e=∑1≤i<j≤nΔijεsisj+∑1≤i≤nαsivi.
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgfbqrdaqadaqaaiabdohazjabcycasiabdogajbgaayjkaiaawmcaaiabg2da9maaqafabagaeuildq0aasbaasqaaiabdmgapjabdqgaqbqabaaccigccqwf1oqzdawgaawcbagaem4cam3aasbaawqaaiabdmgapbqabawccqwgzbwcdawgaaadbagaemoaaogabeaaasqabaaabagaegymaejaeyizimqaemyaakmaeyipawjaemoaaomaeyizimqaemota4eabeqdcqghris5aogaey4kasyaaabuaeaacqwfxoqydaqhaawcbagaem4cam3aasbaawqaaiabdmgapbqabaaaleaacqwg2bgddawgaaadbagaemyaakgabeaaaaaaleaacqaixaqmcqghkjyocqwgpbqacqghkjyocqwgobgtaeqaniabgghildgccqgguaglcawljagaaczcamaabmaabagaegymaejaegymaedacagloagaayzkaaaaaa@62bd@
the first term is a contact free energy: Δij =  <dig> if positions i and j are closer in space than a certain cut-off distance, and  <dig> otherwise, and εab defines the contact energy between amino acids a and b. the second term encodes a solvent-accessibility free energy: for each position, αad
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaaiigacqwfxoqydaqhaawcbagaemyyaegabagaemizaqgaaaaa@311b@ represents the free energy of amino acid a in the solvent accessibility class d, a =  <dig> .
<dig>  and d =  <dig> .d, where d is the total number of solvent accessibility classes considered.
deriving the inverse potential requires the calculation of f, which is already entirely specified by the potential e as a sum over all conformations.
however, this computation is difficult in practice.
as an alternative, we can give it a simple phenomenological form, inspired from the random energy model  <cit> :
f=−∑1≤i≤nμsi,
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwggbgrdaqadaqaaiabdohazbgaayjkaiaawmcaaiabg2da9iabgkhitmaaqafabaaccigae8hvd02aasbaasqaaiabdohaznaabaaameaacqwgpbqaaeqaaawcbeaakiabcycasawcbagaegymaejaeyizimqaemyaakmaeyizimqaemota4eabeqdcqghris5aogaaczcaiaaxmaadaqadaqaaiabigdaxiabikdayagaayjkaiaawmcaaaaa@4637@
where the a =  <dig> .
<dig> are unknown parameters, analogous to "chemical potentials" for the  <dig> amino acids.
altogether, our parameter vector is made of three components: θ = , and the inverse potential reads as:
note that the probability defined by equation  <dig> is invariant under the following transformation:
μ′a
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaaiigacuwf8oqbgaqbamaabaaaleaacqwghbqyaeqaaaaa@2fec@ = μa + j <dig>
ε′ab
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaaiigacuwf1oqzgaqbamaabaaaleaacqwghbqycqwgibgyaeqaaaaa@312a@ = εab + j <dig>
α′ad
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaaiigacuwfxoqygaqbamaadaaaleaacqwghbqyaeaacqwgkbazaaaaaa@3127@ = αad
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaaiigacqwfxoqydaqhaawcbagaemyyaegabagaemizaqgaaaaa@311b@ + j <dig>
where j <dig>  j <dig> and j <dig> are arbitrary real constants.
therefore, to ensure identifiability of our probabilistic model, we enforce the following constraints:
∑aμa= <dig>
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadaaeqbqaaggaciab=x7atnaabaaaleaacqwghbqyaeqaaaqaaiabdggahbqab0gaeyyeiuoakiabg2da9iabicdawiabcycasiaaxmaacawljawaaewaaeaacqaixaqmcqai3awnaiaawicacaglpaaaaaa@3ada@
∑abεab= <dig>
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadaaeqbqaaggaciab=v7alnaabaaaleaacqwghbqycqwgibgyaeqaaaqaaiabdggahjabdkgaibqab0gaeyyeiuoakiabg2da9iabicdawiabcycasiaaxmaacawljawaaewaaeaacqaixaqmcqai4aaoaiaawicacaglpaaaaaa@3d67@
∑aαad= <dig>  d= <dig> .d.
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadaaeqbqaaggaciab=f7ahnaadaaaleaacqwghbqyaeaacqwgkbazaaaabagaemyyaegabeqdcqghris5aogaeyypa0jaegimaajaeiilawiaeeiiaaiaemizaqmaeyypa0jaegymaejaeiola4iaeiola4iaemiraqkaeiola4iaaczcaiaaxmaadaqadaqaaiabigdaxiabimda5agaayjkaiaawmcaaaaa@43e4@
a series of alternative inverse potentials can be obtained by suppressing the first or the second of the components of equation  <dig>  in the present work, we tested the following combinations:
• μ,
• α + μ,
• ε + μ,
• ε + α + μ.
we also explored various numbers of accessibility classes, with d ranging from  <dig> to  <dig>  alternatively, the ε component can be fixed to values of a contact potential obtained by other authors   <cit> .
in this case, we must add a multiplicative scaling factor λ in front of the contact component to account for the fact that these potentials are normalized differently:
the scaling factor is optimized by ml, along with μ.
optimizing the potentials by gradient descent
if we now consider a database, made of p protein sequences s = p =  <dig> .p, of respective lengths np and their corresponding three dimensional structures c = p =  <dig> .p, the probability of observing the whole database, which we define as the likelihood l, is the product of the probabilities of observing each protein independently:
l=p=∏pp=e−gy
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaafaqadewacaaabagaemitaw0aaewaaeaaiigacqwf4oqcaiaawicacaglpaaacqgh9aqpcqwgwbacdaqadaqaaiabdofatjabcyha8jabdoeadjabcycasiab=h7axbgaayjkaiaawmcaaaqaamaabmaabagaegomaijaegymaedacagloagaayzkaaaabagaeyypa0zaaebuaeaacqwgwbacdaqadaqaaiabdohaznaacaaaleqabagaemicaahaaogaeiifawnaem4yam2aawbaasqabeaacqwgwbacaagccqggsaalcqwf4oqcaiaawicacaglpaaaasqaaiabdchawbqab0gaey4diunaaoqaamaabmaabagaegomaijaegomaidacagloagaayzkaaaabagaeyypa0zaasaaaeaacqwglbqzdaahaawcbeqaaiabgkhitiabdeeahnaabmaabagaem4uamlaeiilawiaem4qameacagloagaayzkaaaaaagcbagaemywakfaaaqaamaabmaabagaegomaijaeg4mamdacagloagaayzkaaaaaaaa@61c8@
where
g=∑pg
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwghbwrdaqadaqaaiabdofatjabcycasiabdoeadbgaayjkaiaawmcaaiabg2da9maaqafabagaem4rac0aaewaaeaacqwgzbwcdaahaawcbeqaaiabdchawbaakiabcycasiabdogajnaacaaaleqabagaemicaahaaagccagloagaayzkaaaaleaacqwgwbacaeqaniabgghildgccawljagaaczcamaabmaabagaegomaijaeginaqdacagloagaayzkaaaaaa@4539@
is the inverse potential summed over the database, and
y=∑s′e−g
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgzbqwcqgh9aqpdaaeqbqaaiabdwgalnaacaaaleqabagaeyoei0iaem4rac0aaewaaeaacuwgtbwugaqbaiabcycasiabdoeadbgaayjkaiaawmcaaaaaaeaacuwgtbwugaqbaaqab0gaeyyeiuoakiaaxmaacawljawaaewaaeaacqaiyagmcqai1aqnaiaawicacaglpaaaaaa@3f38@
is the corresponding normalization constant.
since it is more convenient to work on minus the logarithm of the probability, we define the score ω:
ω = - ln l
= g + ln y.
we wish to maximize the likelihood, or equivalently, minimize ω, with respect to θ. we do this by gradient descent, based on a numerical evaluation of the derivative of ω .
the overall method is akin to an expectation maximization algorithm  <cit> .
in fact, it can be seen as a differential version of dempster's method, and therefore, we call it differential em.
the derivative of ω reads as:
∂ω∂θ=∂g∂θ+∂ln⁡y∂θ.
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadawcaaqaaggaciab=jgi2kab=l8a3bqaaiab=jgi2kab=h7axbaacqgh9aqpdawcaaqaaiab=jgi2kabdeeahnaabmaabagaem4uamlaeiilawiaem4qameacagloagaayzkaaaabagae8naiyrae8hudehaaiabgucarmaalaaabagae8naiyragiibawmaeioba4maemywakfabagae8naiyrae8hudehaaiabc6cauiaaxmaacawljawaaewaaeaaieaacqgfyagmcqgfzawmaiaawicacaglpaaaaaa@4d30@
applying the partition function formalism to equation  <dig>  we can express the second term as an expectation over p:
∂ln⁡y∂θ=1y∂y∂θ=−1y∑s′∂g∂θe−g=−∑s′∂g∂θp=−〈∂g∂θ〉
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaafaqadeabcaaaaeaadawcaaqaaggaciab=jgi2kgbcygasjabc6gaujabdmfazbqaaiab=jgi2kab=h7axbaacqgh9aqpdawcaaqaaiabigdaxaqaaiabdmfazbaadawcaaqaaiab=jgi2kabdmfazbqaaiab=jgi2kab=h7axbaaaeaadaqadaqaaiabikdayiabimda5agaayjkaiaawmcaaaqaaiabg2da9iabgkhitmaalaaabagaegymaedabagaemywakfaamaaqafabawaasaaaeaacqwfciitcqwghbwrdaqadaqaaiqbdofatzaafagaeiilawiaem4qameacagloagaayzkaaaabagae8naiyrae8hudehaaawcbagafm4uamlbauaaaeqaniabgghildgccqwglbqzdaahaawcbeqaaiabgkhitiabdeeahnaabmaabagafm4uamlbauaacqggsaalcqwgdbwqaiaawicacaglpaaaaaaakeaadaqadaqaaiabiodaziabicdawagaayjkaiaawmcaaaqaaiabg2da9iabgkhitmaaqafabawaasaaaeaacqwfciitcqwghbwrdaqadaqaaiqbdofatzaafagaeiilawiaem4qameacagloagaayzkaaaabagae8naiyrae8hudehaaawcbagafm4uamlbauaaaeqaniabgghildgccqwgwbacdaqadaqaaiqbdofatzaafagaeiifawnaem4qamkaeiilawiae8hudehacagloagaayzkaaaabawaaewaaeaacqaizawmcqaixaqmaiaawicacaglpaaaaeaacqgh9aqpcqghsislcqghpms4dawcaaqaaiab=jgi2kabdeeahbqaaiab=jgi2kab=h7axbaacqghqms8aeaadaqadaqaaiabiodaziabikdayagaayjkaiaawmcaaaaaaaa@8a2a@
which leads us to the following expression for the derivative of ω:
∂ω∂θ=∂g∂θ−〈∂g∂θ〉.
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadawcaaqaaggaciab=jgi2kab=l8a3bqaaiab=jgi2kab=h7axbaacqgh9aqpdawcaaqaaiab=jgi2kabdeeahnaabmaabagaem4uamlaeiilawiaem4qameacagloagaayzkaaaabagae8naiyrae8hudehaaiabgkhitiabgmyihpaalaaabagae8naiyraem4raceabagae8naiyrae8hudehaaiabgqyixlabc6cauiaaxmaacawljawaaewaaeaacqaizawmcqaizawmaiaawicacaglpaaaaaa@4de1@
the computation of the first term in this equation is straightforward, while the second term must be estimated numerically.
in order to do so, we obtain a sample  drawn from p by a gibbs sampling algorithm similar to that of robinson et al.
<cit>  .
applying formula  <dig> on the inverse potential  <dig> yields the following expressions for the derivatives:
∂ω∂εab=−,
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadawcaaqaaggaciab=jgi2kab=l8a3bqaaiab=jgi2kab=v7alnaabaaaleaacqwghbqycqwgibgyaeqaaaaakiabg2da9iabgkhitmaadmaabagaemoba42aasbaasqaaiabdggahjabdkgaibqabagccqghsisldaaadaqaaiabd6gaunaabaaaleaacqwghbqycqwgibgyaeqaaagccaglpmiaaypkjaaacaglbbgaayzxaagaeiilawiaaczcaiaaxmaadaqadaqaaiabiodaziabisda0agaayjkaiaawmcaaaaa@4a5d@
where nab is the number of contacts between amino acids a and b observed in the database, and 〈nab〉 is its expectation over the probability distribution p. formula  <dig> thus leads to an intuitive characterization of the maximum likelihood estimate ε^
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadaqiaaqaaggaciab=v7albgaaypadaaaaa@2f1c@: it is the value of ε such that the average number of each type of contact predicted by the potential matches the number observed in the database.
following a similar derivation:
∂ω∂μa=−,
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadawcaaqaaggaciab=jgi2kab=l8a3bqaaiab=jgi2kab=x7atnaabaaaleaacqwghbqyaeqaaaaakiabg2da9iabgkhitmaadmaabagaemyba02aasbaasqaaiabdggahbqabagccqghsisldaaadaqaaiabd2gatnaabaaaleaacqwghbqyaeqaaagccaglpmiaaypkjaaacaglbbgaayzxaagaeiilawiaaczcaiaaxmaadaqadaqaaiabiodaziabiwda1agaayjkaiaawmcaaaaa@4683@
where ma is the total number of amino acids of type a, and
∂ω∂αad=−,
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadawcaaqaaggaciab=jgi2kab=l8a3bqaaiab=jgi2kab=f7ahnaadaaaleaacqwghbqyaeaacqwgkbazaaaaaogaeyypa0jaeyoei0yaamwaaeaacqwgsbabdaqhaawcbagaemyyaegabagaemizaqgaaogaeyoei0yaaawaaeaacqwgsbabdaqhaawcbagaemyyaegabagaemizaqgaaagccaglpmiaaypkjaaacaglbbgaayzxaagaeiilawiaaczcaiaaxmaadaqadaqaaiabiodaziabiada2agaayjkaiaawmcaaaaa@4a60@
where lad
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgsbabdaqhaawcbagaemyyaegabagaemizaqgaaaaa@30d6@ is the total number of amino acids of type a belonging to solvent-accessibility class d.
we first performed an optimization of the pure contact potential  on each data set.
figure  <dig> shows the evolution of the scoring function ω and of the contact potential during the gradient descent.
as can be seen from these traceplots, the differential em algorithm converges after a few hundred cycles.
the scoring function stabilizes at around  <dig>  natural units of logarithm , and then fluctuates by up to  <dig> nits around this value.
these fluctuations are mainly due to the finite size of the sample of sequences on which the derivative of ln y is evaluated and, to a lesser extent, to the error on the estimation of ln y by thermodynamic integration.
in any case, these errors are small compared to the differences between scores obtained with alternative models .
the evolution of the potential for some residue pairs is shown in figure 1b and 1c.
effects in the final values due to residue polarity are easily seen: known favorable interactions such as glutamate-lysine or the hydrophobic isoleucine-valine have a lower contact energy, while known unfavorable interactions, such as glutamate-glutamate, have higher energies, indicating that the potentials obtained are biologically reasonable.
the potentials obtained in two independent runs are virtually identical , indicating that the gradient descent does not get trapped into local minima.
we can also compare the values of the potential for two distinct data sets of equivalent size, ds <dig> and ds <dig> , which uncovers a greater discrepancy than for two independent runs on the same data set ds <dig>  the correlation is high, however, suggesting that data sets are large enough for the learning procedure to reach stability.
in addition, these differences are small compared to the discrepancy between the potential obtained by our method and that of miyazawa & jernigan .
model comparison
the same optimization procedure was applied to the potential consisting only of the solvent accessibility term , with an increasing number of accessibility classes, and to the combined  potential.
the resulting log likelihood scores cannot directly be compared, since the models do not have the same dimensionality.
we therefore applied a 2-fold cross-validation procedure , consisting in learning the potential on ds <dig>  and testing it on ds <dig>  and vice versa.
the evolution of the cv score as a function of the number of accessibility classes  is shown in figure  <dig>  when d increases, the fit of the model improves, until a point is reached where the penalization for model dimensionality starts to dominate the score.
the optimal number of classes obtained is  <dig> to  <dig>  depending on the form of the potential studied, although  <dig> to  <dig> classes is sufficient to attain 90% of the fit improvement.
the scores obtained for the different models tested are reported in figure  <dig>  we also included in the comparison the miyazawa and jernigan potential  <cit> .
the contact potential performs better than the pure solvent accessibility potential, and the combination of both terms is the most informative.
miyazawa and jernigan's potential results in a poorer fit improvement than any of the other models.
specificity of the designed sequences
once an optimal value of θ is obtained, properties of the sequences induced by the models can be investigated by sampling sequences from p, using this optimal value of θ. in particular, we tested to what extent the sequences proposed by our method met the requirement of specificity, i.e.
the condition that the sequences designed on a given conformation c indeed have c as their unique ground state.
more precisely, we generated  <dig> sequences by gibbs sampling for  <dig> randomly chosen structures , i.e.
<dig>  sequences for each potential, and performed a fold recognition experiment for the designed sequences, monitoring the score for the target fold using threader  <cit>  .
the solvent accessibility potential alone  is not sufficient to provide specificity to the designed sequences, and behaves almost as poorly as the flat potential .
a mild improvement is seen when using the contact potential : for 10% of the designed sequences the target fold is found among the best scoring folds , and the distribution of this ranking is skewed towards lower values.
however, it is only with the combined potential  that a significant improvement is observed: for more than half of the designed sequences the target fold is found among the best 1% scoring folds, even though the average sequence identity with the native sequence is less than 10% in all cases .
we also tested a subset of  <dig> randomly chosen designed sequences using another fold recognition program, loopp  <cit> .
loopp is based on a combination of several structure prediction methods, based on threading, secondary structure, sequence profile and exposed surface area prediction.
the results obtained with this program were similar to those of threader: for  <dig> % of the designed sequences using the combined  potential, the target fold was found as the first hit, and for  <dig> % the target fold was found among the first  <dig> hits.
in contrast, many of the current fold recognition programs based on sequence profile methods produced no significant hits , which is not surprising, given that our sampling algorithm produces highly divergent sequences, with no similarity to any natural protein.
discussion
the central idea of the present work is to reformulate the problem of devising statistical potentials for protein design as a statistical inference problem.
this reformulation, based on the maximum likelihood  principle, led us naturally to a gradient descent method, with the only additional aspect being that the gradient to follow is itself estimated by monte-carlo averaging.
the main advantage of this ml framework is that it guarantees an optimal predictive power of the resulting potential.
in addition, it is very general, and can in principle be applied to any form of statistical potential.
in particular, it is not restricted to coarse grained descriptions of proteins, and it could also be applied at the atomic level.
interestingly, our gradient descent method turns out to be similar in spirit to an iterative scheme proposed by thomas and dill  <cit> , although in that case the purpose was to optimize a potential in the context of the folding problem.
specifically, thomas and dill tune the potential so as to match the observed and expected number of contacts of each type, except that their expectation is taken on a set of alternative conformations, for a fixed sequence, whereas we take the expectation on a set of alternative sequences, on the conformation of interest.
note that thomas and dill derived their method from intuitive arguments, and not as a mathematical consequence of the ml principle.
these two alternative optimization schemes, obtained by normalizing either over the sequence or over the structure space, are quite distinct, at least conceptually.
how the resulting potentials would differ in practice is more difficult to evaluate.
among other things, it will depend on how the approximation of lnzs based on the random energy model works.
in the eventuality that it does not work well, it is likely that the contact term of our inverse potential will in fact combine two things: the information corresponding to the conformational energy of the sequence itself, which is also encoded in classical potentials optimized for threading, plus some information coming from the decoy term ln zs.
a way to settle this question would be to optimize a contact potential using, on the same learning set, both normalization schemes, and then compare the resulting values as well as their predictive powers.
model assessment and comparison
the methodological framework proposed here offers reliable criteria for comparing the empirical fit of alternative models on real data.
in this respect, it should be noted that the lack of a reliable objective criterion for evaluating different statistical potentials has often been invoked for justifying the use of on-lattice idealized models  <cit> .
however, on-lattice approaches are only moderately interesting, as they completely ignore the problem of the robustness of the learning method to model violation.
coarse-grained statistical potentials are by definition over-simplified models of proteins, and therefore, model violation is an intrinsic feature of the protein design problem.
in this respect, the statistical language is interesting, since it is still valid, even for fitting and assessing models that are known to be imperfect.
on the other hand, the intuitive idea underlying cross-validation, i.e.
measuring the rate of prediction of the native sequence, is quite simple, and has been invoked and used several times previously  <cit> .
what we propose here is a better formalization of this idea.
note that in contrast to previous methods, we do not measure the marginal native prediction rate at each site, but the joint probability of the native sequence.
this can be important, as it accounts for possible correlations in the predictive distribution.
for instance, two given positions may not display any particular pattern, when considered marginally, but may jointly follow charge or steric compensatory patterns.
these phenomena will not be taken into account in the overall fit of the potential when measuring the marginal prediction rate, as is usually done.
technically speaking, the joint probability of the native sequence on the corresponding structure is extremely small, and cannot be evaluated just by counting the frequency at which the native sequence appears in the sample obtained by gibbs sampling.
for this, more elaborate numerical methods, such as thermodynamic integration, are required.
in the present case, the comparison between alternative models has allowed us to measure the relative contribution of each term of the potential and to refine the protein representation.
the contact component turns out to be the most informative , although it should be complemented with other energetic forms.
here, we have tested the addition of a solvent accessibility component, which significantly improves the fit of the model.
contact information and solvent exposure are correlated, which is reflected in the fact that the fit improvement of each term is not additive.
our model comparison method also gives us a direct way of choosing the optimal number of solvent accessibility classes .
here, we found a number of  <dig> to  <dig> classes, which is higher than what one may have expected and than what is usually used.
note that this number depends on the way the classes are defined; here, the classes are based on quantiles, but as an alternative, we also tried a linear definition , which gave us an even higher optimal number of classes .
in general, the present methodology could be used to investigate different definitions of accessibility classes, to refine the pairwise contact definition, or any other elements of the structure representation included in the potential.
the fact that our potential has a significantly better predictive power than that of miyazawa and jernigan  is trivially expected, by construction of the ml potential.
what is more surprising is that the mj matrix is less fit than a simple solvent-accessibility profile.
a possible explanation would be that miyazawa and jernigan's potential is based on the quasi-chemical approximation, which is now known to be somewhat drastic  <cit> , as it neglects correlations between observed pairing frequencies, due to chain connectivity and multiple contacts.
alternatively, it could mean that potentials optimized for folding are really not suited for protein design purposes.
testing other pairwise contact potentials, in particular those that do not rely on the quasi-chemical approximation  <cit> , would be a way to address this issue.
sequence sampling
the method that we propose in this work is probabilistic in essence.
as such, it offers a very natural framework for investigating the patterns induced by the models on distributions of sequences.
specificity of the designed sequences
a sequence s designed for a target conformation c should not only be compatible with c, but also incompatible with competing folds.
a rigorous solution to this problem involves a simultaneous search over the sequence and conformation space.
it is possible, however, to achieve specificity without explicitly seeking to penalize competing states , if we rely on the approximation based on the random energy model, where the normalization constant of equation  <dig> can be considered as a function of the sequence composition only  <cit> .
in our case, the normalization of the likelihood will also play an important role: since the total probability over all possible sequences has to be  <dig>  maximizing the probability for a given sequence s <dig> on its native conformation c <dig> will lower the probability that another natural sequence s <dig>  with native conformation c <dig>  also gets a high probability on c <dig>  when many sequences are learnt in parallel, this phenomenon should ultimately favor specificity of s <dig> on c <dig>  compared to all other conformations of the data set.
on the other hand, the extent to which the specificity is achieved will depend on the actual form of the potential used, as well as on the data base used for learning.
to address this question, we produced a large number of sequences with four different potentials, and checked their ability to recognize the target fold, as measured by the z-score ratio or by the ranking of the target structure in a fold recognition experiment.
indeed, an improvement of specificity is observed when using better potentials, suggesting that the method is effectively capturing specific dependencies between the conformation and the sequence of the proteins in the learning set, even for the simple forms of potentials tested here.
for the combined  potential, the average z-score ratio of the designed sequences is similar to what has been reported for other protein design algorithms  <cit> .
conversely, this also suggests that a more sophisticated potential may further improve the specificity of the sequences designed using our algorithm.
conformation-dependent site-specific profiles
to compare natural protein sequences with those predicted by the optimized potentials, marginal, leave-one-out and empirical profiles  were generated for the  <dig> proteins used in the design specificity experiment described above; the profiles obtained for the best and the worst scoring structures are provided as supplementary materials .
overall, leave-one-out profiles  and marginal profiles  do not display significant differences in the discriminative power between sites: the mean shannon entropy per site is  <dig>  ±  <dig>  for marginal profiles, and  <dig>  ±  <dig>  for leave-one-out profiles.
it is worth noting that the mean entropy per site for each protein, and the corresponding standard deviation, i.e.
the average amount of information at each site and the variation between sites, are both correlated with the performance of the particular protein in the fold recognition experiment, and this, only for the combined  potential .
a detailed analysis of the leave-one-out profiles for a particular case, an alpha-aminotransferase, may be useful to understand which type of information is effectively captured by the potential, and which is not captured at all, thereby suggesting possible ways of improving the current form of potential.
first, regions of the protein that show little secondary structure  contain less information  than regions with local structure .
this is not surprising, since these regions typically have fewer contacts between residues, and thus the amount of information included in the protein representation is lower.
concerning regions with defined secondary structure, residue polarity is the information most easily captured.
charged residues are also distinctively inferred, as well as glycines, to a lesser extent .
in contrast, prolines are rarely correctly predicted, which is expected, since the properties most distinctive of prolines  are not included in this particular form of potential.
interestingly, some residues that have a crucial importance for the protein structure or function fail to be predicted, simply because the properties conferring their importance are not included in the protein description.
this is the case of the amino acids that are in close interaction with a ligand .
finally, the leave-one-out profiles display an interesting behavior with respect to positions where the amino-acid present in the reference sequence is not at all conserved in other members of the family.
in some cases, they simply do not predict anything , which suggests that their limited importance in structure stability or function is recognized by the inverse potential.
in other cases, the natural profile is even reproduced in the leave-one-out profile, instead of the amino acid of the reference sequence; such is the case for phenylalanine  <dig>
CONCLUSIONS
as illustrated by the sequence logos and the fold recognition experiments performed above, the predictive power of the models proposed here is encouraging, but nevertheless still weak.
it is not yet clear to what extent this is due to the specific choice made concerning the form of the statistical potential, to the approximation of ln zs as a function of the sole composition of the sequence, or to yet other reasons.
most probably, we are facing a combination of several factors.
the methods proposed here can now be used to address these difficult questions empirically.
in one direction, other approximations of ln zs, less drastic than the random energy model, but still accessible in practice, can be investigated.
for instance, following deutsch and kurozky , the conditional probability of a sequence could be defined as:
p ∝ e- p
where the expectation 〈·〉 is taken over a pre-defined set of decoy conformations.
more sophisticated monte carlo methods, jointly sampling the sequence and conformation spaces, can also be imagined, in order to get more precise evaluations of ln zs, while staying in the same global maximum likelihood formalism.
on the other hand, all the many statistical potentials that have been proposed over the last fifteen years may in principle be investigated in the same way as we have done here.
in particular, distance-dependent potentials  <cit>  and main-chain dihedral angle potentials  <cit> , which imply a richer representation of the protein structure, may result in models of greater predictive power.
other ways of implicitly considering side-chain conformation may also be easily incorporated into the model.
in a completely different perspective, it is possible to devise probabilistic models that are not exclusively defined in terms of a conformational free energy, even in a formal way.
for instance, additional terms, concerning secondary structure aspects, interactions between successive positions along the sequence, or terms related to the folding constraints, can all be combined in an additive manner in the inverse potential.
in fact, the model need not even be formulated in terms of a boltzmann distribution, as long as the parameters are fitted by ml, and the predictive power of the resulting models is evaluated in a systematic way.
altogether, this amounts to setting up a robust statistical framework helping us to understand how, and to what extent, the sequences of natural proteins are determined by protein structure.

 predicting the sequences compatible with a given structure defines what is traditionally called the inverse folding problem, or more often, protein design. as suggested by the terminology, this question is usually considered in an engineering perspective: the aim is then to determine a sequence, or a set of sequences, that stably fold into a pre-specified conformation. in a thermodynamic perspective, this requirement translates into eliciting sequences that have lowest free energy under the target fold, compared to all possible alternative conformations. in principle, such a criterion would imply a search through the joint structure-sequence space, which is not feasible but for small on-lattice model proteins. as an alternative to the engineering approach, a more evolutionary stance can be taken towards the inverse folding problem, in which case the aim would rather be to predict the sequences of natural proteins having the conformation of interest. seen from this new point of view, the design problem raises new questions: natural proteins are the result of a complex evolutionary process, involving an intricate interplay between mutation and selection, and this probably entails many constraints directly related to the native conformation, but nevertheless not equivalent to the mere requirement of structural stability. for instance, the requirement of fast and cooperative folding has an impact on the dispersion of contact energies. for this and many other potential reasons, among all sequences predicted by classical engineering-oriented protein design, probably only a subset will look like natural proteins. the evolutionary approach to protein design is particularly relevant to phylogenetic studies, where one of the current motivations is to develop the so-called structurally constrained models of protein evolution, i.e. models explicitly dependent on the protein's conformation, either for simulation purposes, or in the context of phylogenetic inference. in this framework, each substitution undergone by a protein during evolution has to be tested for its compatibility with the structure, in the context of the sequence that the protein displays at all other sites when the substitution occurs. such repeated evaluation of the structure-sequence compatibility along a phylogenetic tree requires relevant and computationally very efficient scoring schemes/functions. it is interesting to compare the different methods proposed thus far for performing protein design in light of this engineering/evolutionary distinction. a first direction of research has consisted in using all-atom semi-empirical force fields to evaluate the conformational free energy. these empirical methods have been applied to many theoretical and experimental cases, reaching a high level of accuracy. on the other hand, they are computationally heavy, mainly because of the side-chain positioning problem, and thus cannot be easily applied to structurally constrained phylogenetic models. concerns may also be expressed about their over-sensitivity to the native conformation, in particular in the core of the target structures and when the flexibility of the backbone is not accounted for. but more importantly, approaches based on physical force fields are, by definition, exclusively focussed on the conformational stability, and thereby, completely oversee other potential factors shaping the sequences of biological proteins. as such, they are well suited for engineering synthetic proteins, or for testing to what extent natural sequences are shaped by selection for protein stability, but may not be sufficient for more general evolutionary purposes. an alternative to the semi-empirical strategy consists in relying on knowledge-based, or statistical, potentials. these scoring functions mimic physical boltzmann distributions, but merely encode statistical patterns present in the databases. some of these potentials were obtained under the quasi-chemical approximation, whereby frequencies of patterns, such as contacts between each pair of amino-acids, are transformed into energies using the boltzmann law. alternatively, contact energies can be obtained by maximizing the potential's predictive accuracy in a threading test. in the present context, an advantage of these knowledge-based potentials, compared to semi-empirical force-fields, is that they should in principle capture all kinds of patterns that true biological sequences have, in relation to their conformation, and not only those directly related to thermodynamic stability. furthermore, statistical potentials need not be defined at the atomic level, but can be based on a coarse-grained description of the protein's configuration, essentially by omitting the degrees of freedom associated to side chains. this allows faster computations, by avoiding the problem of searching through the rugged landscape of side-chain conformations. in addition, coarse-grained potentials could turn out to be an advantage, in that they will not recover the native sequence too faithfully. most protein design procedures based on statistical potentials proposed until now have relied on coarse-grained, pairwise contact pseudo-energies. yet, irrespective of the level of description adopted, currently available statistical potentials may not be ideal for protein design, since they have generally been optimized in the context of the folding problem, i.e. for maximizing the rate of correct structure prediction, given the sequence. in contrast, we would like to optimize the reciprocal prediction, namely, the sequences given the conformation. several approaches have been proposed in this direction, consisting in maximizing the z-score between the energy of the native sequence on the target conformation and its energy on a set of decoy sequences, or, alternatively, in applying a mean-square criterion on the values taken by the scoring function on each structure-sequence pair of the database. however, these methods have thus far only been tested in cubic lattice protein models. in addition, they lack a firm theoretical basis. in particular, it would be interesting to guarantee optimal predictive power, and to have a robust methodology available to assess and compare the performance of alternative forms of statistical potentials. standard statistical theory provides such theoretical guarantees. in the present case, the inverse folding problem can be formulated directly in terms of the probability of observing a sequence s given a conformation c, i.e. p. this probability explicitly depends on the pre-specified model through a series of parameters, represented here by θ. these may be, for instance, the coefficients of a pairwise potential, parameters describing compositional effects, secondary structure environment, solvent accessibility, etc. taking the product over a database of p independent sequence-conformation pairs, s = p =.p and c = p =.p, yields a joint probability p=∏pp      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgwbacdaqadaqaaiabdofatjabcyha8jabdoeadjabcycasggaciab=h7axbgaayjkaiaawmcaaiabg2da9maarafabagaemicaa3aaewaaeaacqwgzbwcdaahaawcbeqaaiabdchawbaakiabcyha8jabdogajnaacaaaleqabagaemicaahaaogaeiilawiae8hudehacagloagaayzkaaaaleaacqwgwbacaeqaniabg+givdgccawljagaaczcamaabmaabagaegymaedacagloagaayzkaaaaaa@4b42@ which, as a function of θ, can be seen as a likelihood. the parameter θ is then learnt by maximizing the likelihood with respect to θ. once this is done, sequences can be assessed, or sampled, under the optimal parameter value θ^ mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadaqiaaqaaggaciab=h7axbgaaypadaaaaa@2f2b@, by direct numerical evaluation of their probability, or by monte carlo sampling methods. reformulated in this way, the method maximizes the predictive power of the potential, now in the structure-seeks-sequence direction. by construction, it yields the optimal parameter values that can be obtained for a given form of the potential. in addition, the fit of the model can be directly evaluated, based on the value of the likelihood obtained on a test data set, distinct from the learning set, giving a means of rigorous model selection. finally, the statistical framework proposed here allows one to explicitly combine together, in a model dependent manner, all kinds of factors that we surmise may induce correlations between the structure and the sequence of proteins. we have implemented this maximum likelihood procedure in a markov chain monte carlo framework, and applied it to a simple case, using a contact potential, supplemented with a solvent accessibility term. using cross-validation, we show that the resulting potentials yield a better fit than currently available potentials of the same form, and that combining solvent-accessibility considerations with contact energies is better than either alone. furthermore, we find that solvent accessibility requires a more complex description than what is currently used. ultimately, the overall method proposed in this work can be extended to a large spectrum of alternative models and statistical potentials.  the probabilistic model let us consider a sequence s = i =.n, of length n, and of conformation c. in its most general form, the method introduced here can work with any model m specifying the conditional probability of s given c, in terms of an unnormalized non negative function q: p=q∑sq.      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgwbacdaqadaqaaiabdohazjabcyha8jabdogajjabcycasiabd2eanbgaayjkaiaawmcaaiabg2da9maalaaabagaemycae3aaewaaeaacqwgzbwccqggsaalcqwgjbwyaiaawicacaglpaaaaeaadaaeqaqaaiabdghaxnaabmaabagaem4camnaeiilawiaem4yamgacagloagaayzkaaaaleaacqwgzbwcaeqaniabgghildaaaogaeiola4iaaczcaiaaxmaadaqadaqaaiabikdayagaayjkaiaawmcaaaaa@4c18@ to illustrate the method, we will apply it to a simple case, using a pairwise contact potential. the argument is as follows. first, by bayes' theorem: p=pp∑spp.      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgwbacdaqadaqaaiabdohazjabcyha8jabdogajjabcycasiabd2eanbgaayjkaiaawmcaaiabg2da9maalaaabagaemicaa3aaewaaeaacqwgjbwycqgg8bafcqwgzbwccqggsaalcqwgnbqtaiaawicacaglpaaacqwgwbacdaqadaqaaiabdohazjabcyha8jabd2eanbgaayjkaiaawmcaaaqaamaaqababagaemicaa3aaewaaeaacqwgjbwycqgg8bafcqwgzbwccqggsaalcqwgnbqtaiaawicacaglpaaacqwgwbacdaqadaqaaiabdohazjabcyha8jabd2eanbgaayjkaiaawmcaaawcbagaem4camhabeqdcqghris5aaaakiabc6cauiaaxmaacawljawaaewaaeaacqaizawmaiaawicacaglpaaaaaa@5f64@ if, in addition, we assume a uniform prior on s, we can simply relate equations and by posing q = p. next, given a statistical potential e, the conformational probability p can be expressed as a boltzmann distribution: p=e−e/ktzs=e−−f)/kt, mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaafaqadegacaaabagaemicaa3aaewaaeaacqwgjbwycqgg8bafcqwgzbwccqggsaalcqwgnbqtaiaawicacaglpaaacqgh9aqpdawcaaqaaiabdwgalnaacaaaleqabagaeyoei0iaemyrau0aaewaaeaacqwgzbwccqggsaalcqwgjbwyaiaawicacaglpaaacqggvawlcqwgrbwacqwgubavaaaakeaacqwgabgwdawgaawcbagaem4camhabeaaaaaakeaadaqadaqaaiabisda0agaayjkaiaawmcaaaqaaiabg2da9iabdwgalnaacaaaleqabagaeyoei0yaaewaaeaacqwgfbqrdaqadaqaaiabdohazjabcycasiabdogajbgaayjkaiaawmcaaiabgkhitiabdaeagnaabmaabagaem4camhacagloagaayzkaaaacagloagaayzkaagaei4la8iaem4aasmaemivaqfaaogaeiilawcabawaaewaaeaacqai1aqnaiaawicacaglpaaaaaaaaa@5fa0@ where zs=∑ce−e/kt      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgabgwdawgaawcbagaem4camhabeaakiabg2da9maaqafabagaemyzau2aawbaasqabeaacqghsislcqwgfbqrdaqadaqaaiabdohazjabcycasiabdogajbgaayjkaiaawmcaaiabc+caviabdugarjabdsfaubaaaeaacqwgjbwyaeqaniabgghildgccawljagaaczcamaabmaabagaegonaydacagloagaayzkaaaaaa@43e9@ is a normalization constant, and f = - ln zs.       t and k are the absolute temperature and the boltzmann constant, respectively. without loss of generality, it is possible to rescale the potential so that kt = which we will do in the following. then, by defining the inverse potential: g = e - f,       the conditional probability of sequence s reads as p=e−gy,      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgwbacdaqadaqaaiabdohazjabcyha8jabdogajjabcycasggaciab=h7axjabcycasiabd2eanbgaayjkaiaawmcaaiabg2da9maalaaabagaemyzau2aawbaasqabeaacqghsislcqwghbwrdaqadaqaaiabdohazjabcycasiabdogajbgaayjkaiaawmcaaaaaaoqaaiabdmfazbaacqggsaalcawljagaaczcamaabmaabagaegyoakdacagloagaayzkaaaaaa@482f@ where y=∑s′e−g      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgzbqwcqgh9aqpdaaeqbqaaiabdwgalnaacaaaleqabagaeyoei0iaem4rac0aaewaaeaacuwgzbwcgaqbaiabcycasiabdogajbgaayjkaiaawmcaaaaaaeaacuwgzbwcgaqbaaqab0gaeyyeiuoakiaaxmaacawljawaaewaaeaacqaixaqmcqaiwaamaiaawicacaglpaaaaaa@3fec@ is the normalization factor. note that, contrary to the zs factor of equation which was a sum over all conformations, the present factor y is a sum over sequence space. statistical potentials in the present work, we used a statistical potential made of two terms: e=∑1≤i<j≤nΔijεsisj+∑1≤i≤nαsivi.      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgfbqrdaqadaqaaiabdohazjabcycasiabdogajbgaayjkaiaawmcaaiabg2da9maaqafabagaeuildq0aasbaasqaaiabdmgapjabdqgaqbqabaaccigccqwf1oqzdawgaawcbagaem4cam3aasbaawqaaiabdmgapbqabawccqwgzbwcdawgaaadbagaemoaaogabeaaasqabaaabagaegymaejaeyizimqaemyaakmaeyipawjaemoaaomaeyizimqaemota4eabeqdcqghris5aogaey4kasyaaabuaeaacqwfxoqydaqhaawcbagaem4cam3aasbaawqaaiabdmgapbqabaaaleaacqwg2bgddawgaaadbagaemyaakgabeaaaaaaleaacqaixaqmcqghkjyocqwgpbqacqghkjyocqwgobgtaeqaniabgghildgccqgguaglcawljagaaczcamaabmaabagaegymaejaegymaedacagloagaayzkaaaaaa@62bd@ the first term is a contact free energy: Δij = if positions i and j are closer in space than a certain cut-off distance, and otherwise, and εab defines the contact energy between amino acids a and b. the second term encodes a solvent-accessibility free energy: for each position, αad mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaaiigacqwfxoqydaqhaawcbagaemyyaegabagaemizaqgaaaaa@311b@ represents the free energy of amino acid a in the solvent accessibility class d, a =. and d =.d, where d is the total number of solvent accessibility classes considered. deriving the inverse potential requires the calculation of f, which is already entirely specified by the potential e as a sum over all conformations. however, this computation is difficult in practice. as an alternative, we can give it a simple phenomenological form, inspired from the random energy model: f=−∑1≤i≤nμsi,      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwggbgrdaqadaqaaiabdohazbgaayjkaiaawmcaaiabg2da9iabgkhitmaaqafabaaccigae8hvd02aasbaasqaaiabdohaznaabaaameaacqwgpbqaaeqaaawcbeaakiabcycasawcbagaegymaejaeyizimqaemyaakmaeyizimqaemota4eabeqdcqghris5aogaaczcaiaaxmaadaqadaqaaiabigdaxiabikdayagaayjkaiaawmcaaaaa@4637@ where the a =. are unknown parameters, analogous to "chemical potentials" for the amino acids. altogether, our parameter vector is made of three components: θ =, and the inverse potential reads as: note that the probability defined by equation is invariant under the following transformation: μ′a mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaaiigacuwf8oqbgaqbamaabaaaleaacqwghbqyaeqaaaaa@2fec@ = μa + j     ε′ab mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaaiigacuwf1oqzgaqbamaabaaaleaacqwghbqycqwgibgyaeqaaaaa@312a@ = εab + j     α′ad mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaaiigacuwfxoqygaqbamaadaaaleaacqwghbqyaeaacqwgkbazaaaaaa@3127@ = αad mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaaiigacqwfxoqydaqhaawcbagaemyyaegabagaemizaqgaaaaa@311b@ + j     where j j and j are arbitrary real constants. therefore, to ensure identifiability of our probabilistic model, we enforce the following constraints: ∑aμa=       mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadaaeqbqaaggaciab=x7atnaabaaaleaacqwghbqyaeqaaaqaaiabdggahbqab0gaeyyeiuoakiabg2da9iabicdawiabcycasiaaxmaacawljawaaewaaeaacqaixaqmcqai3awnaiaawicacaglpaaaaaa@3ada@ ∑abεab=       mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadaaeqbqaaggaciab=v7alnaabaaaleaacqwghbqycqwgibgyaeqaaaqaaiabdggahjabdkgaibqab0gaeyyeiuoakiabg2da9iabicdawiabcycasiaaxmaacawljawaaewaaeaacqaixaqmcqai4aaoaiaawicacaglpaaaaaa@3d67@ ∑aαad=  d=.d.      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadaaeqbqaaggaciab=f7ahnaadaaaleaacqwghbqyaeaacqwgkbazaaaabagaemyyaegabeqdcqghris5aogaeyypa0jaegimaajaeiilawiaeeiiaaiaemizaqmaeyypa0jaegymaejaeiola4iaeiola4iaemiraqkaeiola4iaaczcaiaaxmaadaqadaqaaiabigdaxiabimda5agaayjkaiaawmcaaaaa@43e4@ a series of alternative inverse potentials can be obtained by suppressing the first or the second of the components of equation in the present work, we tested the following combinations: • μ, • α + μ, • ε + μ, • ε + α + μ. we also explored various numbers of accessibility classes, with d ranging from to alternatively, the ε component can be fixed to values of a contact potential obtained by other authors. in this case, we must add a multiplicative scaling factor λ in front of the contact component to account for the fact that these potentials are normalized differently: the scaling factor is optimized by ml, along with μ. optimizing the potentials by gradient descent if we now consider a database, made of p protein sequences s = p =.p, of respective lengths np and their corresponding three dimensional structures c = p =.p, the probability of observing the whole database, which we define as the likelihood l, is the product of the probabilities of observing each protein independently: l=p=∏pp=e−gy mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaafaqadewacaaabagaemitaw0aaewaaeaaiigacqwf4oqcaiaawicacaglpaaacqgh9aqpcqwgwbacdaqadaqaaiabdofatjabcyha8jabdoeadjabcycasiab=h7axbgaayjkaiaawmcaaaqaamaabmaabagaegomaijaegymaedacagloagaayzkaaaabagaeyypa0zaaebuaeaacqwgwbacdaqadaqaaiabdohaznaacaaaleqabagaemicaahaaogaeiifawnaem4yam2aawbaasqabeaacqwgwbacaagccqggsaalcqwf4oqcaiaawicacaglpaaaasqaaiabdchawbqab0gaey4diunaaoqaamaabmaabagaegomaijaegomaidacagloagaayzkaaaabagaeyypa0zaasaaaeaacqwglbqzdaahaawcbeqaaiabgkhitiabdeeahnaabmaabagaem4uamlaeiilawiaem4qameacagloagaayzkaaaaaagcbagaemywakfaaaqaamaabmaabagaegomaijaeg4mamdacagloagaayzkaaaaaaaa@61c8@ where g=∑pg      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwghbwrdaqadaqaaiabdofatjabcycasiabdoeadbgaayjkaiaawmcaaiabg2da9maaqafabagaem4rac0aaewaaeaacqwgzbwcdaahaawcbeqaaiabdchawbaakiabcycasiabdogajnaacaaaleqabagaemicaahaaagccagloagaayzkaaaaleaacqwgwbacaeqaniabgghildgccawljagaaczcamaabmaabagaegomaijaeginaqdacagloagaayzkaaaaaa@4539@ is the inverse potential summed over the database, and y=∑s′e−g      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgzbqwcqgh9aqpdaaeqbqaaiabdwgalnaacaaaleqabagaeyoei0iaem4rac0aaewaaeaacuwgtbwugaqbaiabcycasiabdoeadbgaayjkaiaawmcaaaaaaeaacuwgtbwugaqbaaqab0gaeyyeiuoakiaaxmaacawljawaaewaaeaacqaiyagmcqai1aqnaiaawicacaglpaaaaaa@3f38@ is the corresponding normalization constant. since it is more convenient to work on minus the logarithm of the probability, we define the score ω: ω = - ln l       = g + ln y.       we wish to maximize the likelihood, or equivalently, minimize ω, with respect to θ. we do this by gradient descent, based on a numerical evaluation of the derivative of ω. the overall method is akin to an expectation maximization algorithm. in fact, it can be seen as a differential version of dempster's method, and therefore, we call it differential em. the derivative of ω reads as: ∂ω∂θ=∂g∂θ+∂ln⁡y∂θ.      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadawcaaqaaggaciab=jgi2kab=l8a3bqaaiab=jgi2kab=h7axbaacqgh9aqpdawcaaqaaiab=jgi2kabdeeahnaabmaabagaem4uamlaeiilawiaem4qameacagloagaayzkaaaabagae8naiyrae8hudehaaiabgucarmaalaaabagae8naiyragiibawmaeioba4maemywakfabagae8naiyrae8hudehaaiabc6cauiaaxmaacawljawaaewaaeaaieaacqgfyagmcqgfzawmaiaawicacaglpaaaaaa@4d30@ applying the partition function formalism to equation we can express the second term as an expectation over p: ∂ln⁡y∂θ=1y∂y∂θ=−1y∑s′∂g∂θe−g=−∑s′∂g∂θp=−〈∂g∂θ〉 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaafaqadeabcaaaaeaadawcaaqaaggaciab=jgi2kgbcygasjabc6gaujabdmfazbqaaiab=jgi2kab=h7axbaacqgh9aqpdawcaaqaaiabigdaxaqaaiabdmfazbaadawcaaqaaiab=jgi2kabdmfazbqaaiab=jgi2kab=h7axbaaaeaadaqadaqaaiabikdayiabimda5agaayjkaiaawmcaaaqaaiabg2da9iabgkhitmaalaaabagaegymaedabagaemywakfaamaaqafabawaasaaaeaacqwfciitcqwghbwrdaqadaqaaiqbdofatzaafagaeiilawiaem4qameacagloagaayzkaaaabagae8naiyrae8hudehaaawcbagafm4uamlbauaaaeqaniabgghildgccqwglbqzdaahaawcbeqaaiabgkhitiabdeeahnaabmaabagafm4uamlbauaacqggsaalcqwgdbwqaiaawicacaglpaaaaaaakeaadaqadaqaaiabiodaziabicdawagaayjkaiaawmcaaaqaaiabg2da9iabgkhitmaaqafabawaasaaaeaacqwfciitcqwghbwrdaqadaqaaiqbdofatzaafagaeiilawiaem4qameacagloagaayzkaaaabagae8naiyrae8hudehaaawcbagafm4uamlbauaaaeqaniabgghildgccqwgwbacdaqadaqaaiqbdofatzaafagaeiifawnaem4qamkaeiilawiae8hudehacagloagaayzkaaaabawaaewaaeaacqaizawmcqaixaqmaiaawicacaglpaaaaeaacqgh9aqpcqghsislcqghpms4dawcaaqaaiab=jgi2kabdeeahbqaaiab=jgi2kab=h7axbaacqghqms8aeaadaqadaqaaiabiodaziabikdayagaayjkaiaawmcaaaaaaaa@8a2a@ which leads us to the following expression for the derivative of ω: ∂ω∂θ=∂g∂θ−〈∂g∂θ〉.      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadawcaaqaaggaciab=jgi2kab=l8a3bqaaiab=jgi2kab=h7axbaacqgh9aqpdawcaaqaaiab=jgi2kabdeeahnaabmaabagaem4uamlaeiilawiaem4qameacagloagaayzkaaaabagae8naiyrae8hudehaaiabgkhitiabgmyihpaalaaabagae8naiyraem4raceabagae8naiyrae8hudehaaiabgqyixlabc6cauiaaxmaacawljawaaewaaeaacqaizawmcqaizawmaiaawicacaglpaaaaaa@4de1@ the computation of the first term in this equation is straightforward, while the second term must be estimated numerically. in order to do so, we obtain a sample drawn from p by a gibbs sampling algorithm similar to that of robinson et al.. applying formula on the inverse potential yields the following expressions for the derivatives: ∂ω∂εab=−,      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadawcaaqaaggaciab=jgi2kab=l8a3bqaaiab=jgi2kab=v7alnaabaaaleaacqwghbqycqwgibgyaeqaaaaakiabg2da9iabgkhitmaadmaabagaemoba42aasbaasqaaiabdggahjabdkgaibqabagccqghsisldaaadaqaaiabd6gaunaabaaaleaacqwghbqycqwgibgyaeqaaagccaglpmiaaypkjaaacaglbbgaayzxaagaeiilawiaaczcaiaaxmaadaqadaqaaiabiodaziabisda0agaayjkaiaawmcaaaaa@4a5d@ where nab is the number of contacts between amino acids a and b observed in the database, and 〈nab〉 is its expectation over the probability distribution p. formula thus leads to an intuitive characterization of the maximum likelihood estimate ε^ mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadaqiaaqaaggaciab=v7albgaaypadaaaaa@2f1c@: it is the value of ε such that the average number of each type of contact predicted by the potential matches the number observed in the database. following a similar derivation: ∂ω∂μa=−,      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadawcaaqaaggaciab=jgi2kab=l8a3bqaaiab=jgi2kab=x7atnaabaaaleaacqwghbqyaeqaaaaakiabg2da9iabgkhitmaadmaabagaemyba02aasbaasqaaiabdggahbqabagccqghsisldaaadaqaaiabd2gatnaabaaaleaacqwghbqyaeqaaagccaglpmiaaypkjaaacaglbbgaayzxaagaeiilawiaaczcaiaaxmaadaqadaqaaiabiodaziabiwda1agaayjkaiaawmcaaaaa@4683@ where ma is the total number of amino acids of type a, and ∂ω∂αad=−,      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadawcaaqaaggaciab=jgi2kab=l8a3bqaaiab=jgi2kab=f7ahnaadaaaleaacqwghbqyaeaacqwgkbazaaaaaogaeyypa0jaeyoei0yaamwaaeaacqwgsbabdaqhaawcbagaemyyaegabagaemizaqgaaogaeyoei0yaaawaaeaacqwgsbabdaqhaawcbagaemyyaegabagaemizaqgaaagccaglpmiaaypkjaaacaglbbgaayzxaagaeiilawiaaczcaiaaxmaadaqadaqaaiabiodaziabiada2agaayjkaiaawmcaaaaa@4a60@ where lad mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgsbabdaqhaawcbagaemyyaegabagaemizaqgaaaaa@30d6@ is the total number of amino acids of type a belonging to solvent-accessibility class d. we first performed an optimization of the pure contact potential on each data set. figure shows the evolution of the scoring function ω and of the contact potential during the gradient descent. as can be seen from these traceplots, the differential em algorithm converges after a few hundred cycles. the scoring function stabilizes at around natural units of logarithm, and then fluctuates by up to nits around this value. these fluctuations are mainly due to the finite size of the sample of sequences on which the derivative of ln y is evaluated and, to a lesser extent, to the error on the estimation of ln y by thermodynamic integration. in any case, these errors are small compared to the differences between scores obtained with alternative models. the evolution of the potential for some residue pairs is shown in figure 1b and 1c. effects in the final values due to residue polarity are easily seen: known favorable interactions such as glutamate-lysine or the hydrophobic isoleucine-valine have a lower contact energy, while known unfavorable interactions, such as glutamate-glutamate, have higher energies, indicating that the potentials obtained are biologically reasonable. the potentials obtained in two independent runs are virtually identical, indicating that the gradient descent does not get trapped into local minima. we can also compare the values of the potential for two distinct data sets of equivalent size, ds and ds, which uncovers a greater discrepancy than for two independent runs on the same data set ds the correlation is high, however, suggesting that data sets are large enough for the learning procedure to reach stability. in addition, these differences are small compared to the discrepancy between the potential obtained by our method and that of miyazawa & jernigan. model comparison the same optimization procedure was applied to the potential consisting only of the solvent accessibility term, with an increasing number of accessibility classes, and to the combined potential. the resulting log likelihood scores cannot directly be compared, since the models do not have the same dimensionality. we therefore applied a 2-fold cross-validation procedure, consisting in learning the potential on ds and testing it on ds and vice versa. the evolution of the cv score as a function of the number of accessibility classes is shown in figure when d increases, the fit of the model improves, until a point is reached where the penalization for model dimensionality starts to dominate the score. the optimal number of classes obtained is to depending on the form of the potential studied, although to classes is sufficient to attain 90% of the fit improvement. the scores obtained for the different models tested are reported in figure we also included in the comparison the miyazawa and jernigan potential. the contact potential performs better than the pure solvent accessibility potential, and the combination of both terms is the most informative. miyazawa and jernigan's potential results in a poorer fit improvement than any of the other models. specificity of the designed sequences once an optimal value of θ is obtained, properties of the sequences induced by the models can be investigated by sampling sequences from p, using this optimal value of θ. in particular, we tested to what extent the sequences proposed by our method met the requirement of specificity, i.e. the condition that the sequences designed on a given conformation c indeed have c as their unique ground state. more precisely, we generated sequences by gibbs sampling for randomly chosen structures, i.e. sequences for each potential, and performed a fold recognition experiment for the designed sequences, monitoring the score for the target fold using threader. the solvent accessibility potential alone is not sufficient to provide specificity to the designed sequences, and behaves almost as poorly as the flat potential. a mild improvement is seen when using the contact potential: for 10% of the designed sequences the target fold is found among the best scoring folds, and the distribution of this ranking is skewed towards lower values. however, it is only with the combined potential that a significant improvement is observed: for more than half of the designed sequences the target fold is found among the best 1% scoring folds, even though the average sequence identity with the native sequence is less than 10% in all cases. we also tested a subset of randomly chosen designed sequences using another fold recognition program, loopp. loopp is based on a combination of several structure prediction methods, based on threading, secondary structure, sequence profile and exposed surface area prediction. the results obtained with this program were similar to those of threader: for % of the designed sequences using the combined potential, the target fold was found as the first hit, and for % the target fold was found among the first hits. in contrast, many of the current fold recognition programs based on sequence profile methods produced no significant hits, which is not surprising, given that our sampling algorithm produces highly divergent sequences, with no similarity to any natural protein. discussion the central idea of the present work is to reformulate the problem of devising statistical potentials for protein design as a statistical inference problem. this reformulation, based on the maximum likelihood principle, led us naturally to a gradient descent method, with the only additional aspect being that the gradient to follow is itself estimated by monte-carlo averaging. the main advantage of this ml framework is that it guarantees an optimal predictive power of the resulting potential. in addition, it is very general, and can in principle be applied to any form of statistical potential. in particular, it is not restricted to coarse grained descriptions of proteins, and it could also be applied at the atomic level. interestingly, our gradient descent method turns out to be similar in spirit to an iterative scheme proposed by thomas and dill, although in that case the purpose was to optimize a potential in the context of the folding problem. specifically, thomas and dill tune the potential so as to match the observed and expected number of contacts of each type, except that their expectation is taken on a set of alternative conformations, for a fixed sequence, whereas we take the expectation on a set of alternative sequences, on the conformation of interest. note that thomas and dill derived their method from intuitive arguments, and not as a mathematical consequence of the ml principle. these two alternative optimization schemes, obtained by normalizing either over the sequence or over the structure space, are quite distinct, at least conceptually. how the resulting potentials would differ in practice is more difficult to evaluate. among other things, it will depend on how the approximation of lnzs based on the random energy model works. in the eventuality that it does not work well, it is likely that the contact term of our inverse potential will in fact combine two things: the information corresponding to the conformational energy of the sequence itself, which is also encoded in classical potentials optimized for threading, plus some information coming from the decoy term ln zs. a way to settle this question would be to optimize a contact potential using, on the same learning set, both normalization schemes, and then compare the resulting values as well as their predictive powers. model assessment and comparison the methodological framework proposed here offers reliable criteria for comparing the empirical fit of alternative models on real data. in this respect, it should be noted that the lack of a reliable objective criterion for evaluating different statistical potentials has often been invoked for justifying the use of on-lattice idealized models. however, on-lattice approaches are only moderately interesting, as they completely ignore the problem of the robustness of the learning method to model violation. coarse-grained statistical potentials are by definition over-simplified models of proteins, and therefore, model violation is an intrinsic feature of the protein design problem. in this respect, the statistical language is interesting, since it is still valid, even for fitting and assessing models that are known to be imperfect. on the other hand, the intuitive idea underlying cross-validation, i.e. measuring the rate of prediction of the native sequence, is quite simple, and has been invoked and used several times previously. what we propose here is a better formalization of this idea. note that in contrast to previous methods, we do not measure the marginal native prediction rate at each site, but the joint probability of the native sequence. this can be important, as it accounts for possible correlations in the predictive distribution. for instance, two given positions may not display any particular pattern, when considered marginally, but may jointly follow charge or steric compensatory patterns. these phenomena will not be taken into account in the overall fit of the potential when measuring the marginal prediction rate, as is usually done. technically speaking, the joint probability of the native sequence on the corresponding structure is extremely small, and cannot be evaluated just by counting the frequency at which the native sequence appears in the sample obtained by gibbs sampling. for this, more elaborate numerical methods, such as thermodynamic integration, are required. in the present case, the comparison between alternative models has allowed us to measure the relative contribution of each term of the potential and to refine the protein representation. the contact component turns out to be the most informative, although it should be complemented with other energetic forms. here, we have tested the addition of a solvent accessibility component, which significantly improves the fit of the model. contact information and solvent exposure are correlated, which is reflected in the fact that the fit improvement of each term is not additive. our model comparison method also gives us a direct way of choosing the optimal number of solvent accessibility classes. here, we found a number of to classes, which is higher than what one may have expected and than what is usually used. note that this number depends on the way the classes are defined; here, the classes are based on quantiles, but as an alternative, we also tried a linear definition, which gave us an even higher optimal number of classes. in general, the present methodology could be used to investigate different definitions of accessibility classes, to refine the pairwise contact definition, or any other elements of the structure representation included in the potential. the fact that our potential has a significantly better predictive power than that of miyazawa and jernigan is trivially expected, by construction of the ml potential. what is more surprising is that the mj matrix is less fit than a simple solvent-accessibility profile. a possible explanation would be that miyazawa and jernigan's potential is based on the quasi-chemical approximation, which is now known to be somewhat drastic, as it neglects correlations between observed pairing frequencies, due to chain connectivity and multiple contacts. alternatively, it could mean that potentials optimized for folding are really not suited for protein design purposes. testing other pairwise contact potentials, in particular those that do not rely on the quasi-chemical approximation, would be a way to address this issue. sequence sampling the method that we propose in this work is probabilistic in essence. as such, it offers a very natural framework for investigating the patterns induced by the models on distributions of sequences. specificity of the designed sequences a sequence s designed for a target conformation c should not only be compatible with c, but also incompatible with competing folds. a rigorous solution to this problem involves a simultaneous search over the sequence and conformation space. it is possible, however, to achieve specificity without explicitly seeking to penalize competing states, if we rely on the approximation based on the random energy model, where the normalization constant of equation can be considered as a function of the sequence composition only. in our case, the normalization of the likelihood will also play an important role: since the total probability over all possible sequences has to be maximizing the probability for a given sequence s on its native conformation c will lower the probability that another natural sequence s with native conformation c also gets a high probability on c when many sequences are learnt in parallel, this phenomenon should ultimately favor specificity of s on c compared to all other conformations of the data set. on the other hand, the extent to which the specificity is achieved will depend on the actual form of the potential used, as well as on the data base used for learning. to address this question, we produced a large number of sequences with four different potentials, and checked their ability to recognize the target fold, as measured by the z-score ratio or by the ranking of the target structure in a fold recognition experiment. indeed, an improvement of specificity is observed when using better potentials, suggesting that the method is effectively capturing specific dependencies between the conformation and the sequence of the proteins in the learning set, even for the simple forms of potentials tested here. for the combined potential, the average z-score ratio of the designed sequences is similar to what has been reported for other protein design algorithms. conversely, this also suggests that a more sophisticated potential may further improve the specificity of the sequences designed using our algorithm. conformation-dependent site-specific profiles to compare natural protein sequences with those predicted by the optimized potentials, marginal, leave-one-out and empirical profiles were generated for the proteins used in the design specificity experiment described above; the profiles obtained for the best and the worst scoring structures are provided as supplementary materials. overall, leave-one-out profiles and marginal profiles do not display significant differences in the discriminative power between sites: the mean shannon entropy per site is ± for marginal profiles, and ± for leave-one-out profiles. it is worth noting that the mean entropy per site for each protein, and the corresponding standard deviation, i.e. the average amount of information at each site and the variation between sites, are both correlated with the performance of the particular protein in the fold recognition experiment, and this, only for the combined potential. a detailed analysis of the leave-one-out profiles for a particular case, an alpha-aminotransferase, may be useful to understand which type of information is effectively captured by the potential, and which is not captured at all, thereby suggesting possible ways of improving the current form of potential. first, regions of the protein that show little secondary structure contain less information than regions with local structure. this is not surprising, since these regions typically have fewer contacts between residues, and thus the amount of information included in the protein representation is lower. concerning regions with defined secondary structure, residue polarity is the information most easily captured. charged residues are also distinctively inferred, as well as glycines, to a lesser extent. in contrast, prolines are rarely correctly predicted, which is expected, since the properties most distinctive of prolines are not included in this particular form of potential. interestingly, some residues that have a crucial importance for the protein structure or function fail to be predicted, simply because the properties conferring their importance are not included in the protein description. this is the case of the amino acids that are in close interaction with a ligand. finally, the leave-one-out profiles display an interesting behavior with respect to positions where the amino-acid present in the reference sequence is not at all conserved in other members of the family. in some cases, they simply do not predict anything, which suggests that their limited importance in structure stability or function is recognized by the inverse potential. in other cases, the natural profile is even reproduced in the leave-one-out profile, instead of the amino acid of the reference sequence; such is the case for phenylalanine  as illustrated by the sequence logos and the fold recognition experiments performed above, the predictive power of the models proposed here is encouraging, but nevertheless still weak. it is not yet clear to what extent this is due to the specific choice made concerning the form of the statistical potential, to the approximation of ln zs as a function of the sole composition of the sequence, or to yet other reasons. most probably, we are facing a combination of several factors. the methods proposed here can now be used to address these difficult questions empirically. in one direction, other approximations of ln zs, less drastic than the random energy model, but still accessible in practice, can be investigated. for instance, following deutsch and kurozky, the conditional probability of a sequence could be defined as: p ∝ e- p       where the expectation 〈·〉 is taken over a pre-defined set of decoy conformations. more sophisticated monte carlo methods, jointly sampling the sequence and conformation spaces, can also be imagined, in order to get more precise evaluations of ln zs, while staying in the same global maximum likelihood formalism. on the other hand, all the many statistical potentials that have been proposed over the last fifteen years may in principle be investigated in the same way as we have done here. in particular, distance-dependent potentials and main-chain dihedral angle potentials, which imply a richer representation of the protein structure, may result in models of greater predictive power. other ways of implicitly considering side-chain conformation may also be easily incorporated into the model. in a completely different perspective, it is possible to devise probabilistic models that are not exclusively defined in terms of a conformational free energy, even in a formal way. for instance, additional terms, concerning secondary structure aspects, interactions between successive positions along the sequence, or terms related to the folding constraints, can all be combined in an additive manner in the inverse potential. in fact, the model need not even be formulated in terms of a boltzmann distribution, as long as the parameters are fitted by ml, and the predictive power of the resulting models is evaluated in a systematic way. altogether, this amounts to setting up a robust statistical framework helping us to understand how, and to what extent, the sequences of natural proteins are determined by protein structure. 
the aim of protein design is to predict amino-acid sequences compatible with a given target structure. furthermore, the model comparison method presented here allows one to measure the relative contribution of each component of the potential, and to choose the optimal number of accessibility classes, which turns out to be much higher than classically considered. altogether, this reformulation makes it possible to test a wide diversity of models, using different forms of potentials, or accounting for other factors than just the constraint of thermodynamic stability. we propose an implementation based on markov chain monte carlo, in which the likelihood is maximized by gradient descent and is numerically estimated by thermodynamic integration. our framework uses the maximum likelihood principle to optimize the unknown parameters of a statistical potential, which we call an inverse potential to contrast with classical potentials used for structure prediction. we propose a formulation of the protein design problem in terms of model-based statistical inference. we apply this to a simple pairwise contact potential, supplemented with a solvent-accessibility term, and show that the resulting models have a better predictive power than currently available pairwise potentials. traditionally envisioned as a purely thermodynamic question, this problem can also be understood in a wider context, where additional constraints are captured by learning the sequence patterns displayed by natural proteins of known conformation.
 the aim of protein design is to predict amino-acid sequences compatible with a given target structure. traditionally envisioned as a purely thermodynamic question, this problem can also be understood in a wider context, where additional constraints are captured by learning the sequence patterns displayed by natural proteins of known conformation. in this latter perspective, however, we still need a theoretical formalization of the question, leading to general and efficient learning methods, and allowing for the selection of fast and accurate objective functions quantifying sequence/structure compatibility.  we propose a formulation of the protein design problem in terms of model-based statistical inference. our framework uses the maximum likelihood principle to optimize the unknown parameters of a statistical potential, which we call an inverse potential to contrast with classical potentials used for structure prediction. we propose an implementation based on markov chain monte carlo, in which the likelihood is maximized by gradient descent and is numerically estimated by thermodynamic integration. the fit of the models is evaluated by cross-validation. we apply this to a simple pairwise contact potential, supplemented with a solvent-accessibility term, and show that the resulting models have a better predictive power than currently available pairwise potentials. furthermore, the model comparison method presented here allows one to measure the relative contribution of each component of the potential, and to choose the optimal number of accessibility classes, which turns out to be much higher than classically considered.  altogether, this reformulation makes it possible to test a wide diversity of models, using different forms of potentials, or accounting for other factors than just the constraint of thermodynamic stability. ultimately, such model-based statistical analyses may help to understand the forces shaping protein sequences, and driving their evolution. 