BACKGROUND
protein-coding sequences play a central role in cellular functions necessary to produce the vast variation in organismal phenotypes observed in nature.
to function, most proteins must fold into a unique and stable structure  <cit> .
the structure is responsible for orienting residues necessary for proper function, including binding and catalysis.
to maintain fitness, protein function and the underlying structure must be preserved.
protein structure and binding  are determined by the interactions of individual amino acid residues.
therefore, to mechanistically model the evolution of protein sequences, these residue-residue interactions must be considered, relaxing the common assumption of independent evolution of each amino acid position  <cit> .
structural models for sequence evolution where function is protein-protein intermolecular interaction will now be considered.
to evaluate if a sequence will fold into a unique and stable structure, it must be demonstrated that the sequence prefers the folded state to both the unfolded state  as well as to alternative folded states  <cit> .
from a practical perspective, it is impossible to enumerate the enormous space of all possible alternative conformations  <cit> .
therefore, sets of representative "decoy conformations" must be used to approximate folding intermediates and/or other stable or meta-stable states  <cit> .
alternatively, without explicitly considering decoy structures, the inter-residue contacts from different conformations can be randomly sampled  <cit> .
additionally, a scoring function is needed to quantify differences between the possible states.
strategies developed for generating such scoring functions include those derived statistically from existing structures  <cit>  and those derived from physical first principles  <cit> .
both types of scoring function assume that the native sequence lies close to  a thermodynamic optimum, and that there is a gap in energy between the native state and the closest non-native state  <cit> .
formally, the protein sequence evolution problem relates to the inverse folding problem  rather than the folding problem.
underlying both problems are similar physical assumptions and similar models  <cit> .
protein biological function depends not only on proper folding, but also on the ability to bind the target ligands.
therefore, to study sequence evolution, binding must also be evaluated.
from a physical perspective this is easier, as the only states to consider are the bound and the unbound states.
however, it has recently been suggested that selective pressures on proteins to avoid non-specific binding are also an important aspect of protein fitness, necessitating consideration of binding "decoys" as well  <cit> .
identifying the actual selective pressures against non-specific binding in a cell is a daunting task, but the number and nature of binding decoys is a major determinant of the ease of evolving new binding interactions.
the binding decoy characteristics also affect the level of selective constraint on the binding interface.
site and function independent models clearly do not capture critical elements of protein evolution.
for structure-based models, a good scoring function must produce sequences with similar properties to real proteins.
this includes a hydrophobic core that evolves slower than the hydrophilic surface  <cit> .
a dn/ds ratio  much smaller than unity, particularly for functional residues, should also emerge  <cit> .
more particularly, amino acid substitution rates should show heterogeneity, expressed by a gamma distribution across positions  <cit> , reflecting their relative importance to folding and function.
while the gamma function for rate heterogeneity was adopted for model fit rather than mechanistic purposes, it is one of the most common parameters in standard evolutionary models and accounts for heterogeneity in the substitution rate driven by structural signals as well as other signals in evolutionary sequence data  <cit> .
in addition to rate heterogeneity, sequences must also have an energy gap between the native and alternative conformations to ensure rapid and stable folding  <cit> , although structurally disordered proteins represent a distinct class of proteins that do not have this property  <cit> .
the model should place the native sequence near an optimum so as not to provide a signal of directional selection when function is not changing.
consequently, most mutations should be deleterious or nearly neutral rather than adaptive  <cit> .
lastly, proteins must retain their binding function.
population size dictates what fitness changes are neutral as well as what mutations become substitutions  <cit>  and must therefore be taken into account.
to enable use in forward  and backward  studies of evolution, especially in a population genetic context, in a complete genome context, or in studies of interactome evolution, the model needs to be coarse grained at a level that allows for sufficient computational speed.
in this context the computational cost makes state of the art models and methods from the protein structure community  <cit>  intractable.
here a novel physics-based scoring function is developed and compared with a commonly used informational approach on the criteria described above.
the physics-based model is more specific and predicts less drastic stability changes.
both types of models violate the assumption of high native sequence stability and produce many adaptive mutations through directional selection towards a scoring function optimum, even though the model works with truncation selection.
these properties are discussed in light of their effect on modeling accuracy and improvements to the models are suggested.
RESULTS
protein folding and function cause interdependence between sites in the protein sequence.
for instance, a mutation that removes a cysteine involved in a disulfide bridge  is likely to be deleterious, whereas mutation of other cysteines  may be more neutral.
of course, disulfide bridges can be lost and interconverted to other types of stabilizing interactions over evolutionary time, albeit with a large transitional barrier due to the strong site interdependence.
such selective pressures can be modeled by scoring the thermodynamic effects of mutations.
two approaches to calculating such scores are the statistically motivated informational methods  and the first principles physics-based methods.
informational models score the likelihood of observing specific types of contacts in the folded protein based on those seen in previously known structures  <cit> .
physics-based approaches evaluate structures by measuring the fit of residues to geometrical properties such as backbone torsion and residue-residue distances  <cit> .
in both cases the fit to both the native and the many possible non-native conformations must be measured to ensure specificity.
this manuscript sets out to evaluate the ability to replicate biological properties of proteins using both approaches described above through evolutionary simulation.
a new physics-based approach is described.
native sequences were evaluated in both approaches to evaluate if they fell within the neutral network of an optimum or if they underwent directional selection initially.
similarly, sequences sampled from within the optimum were evolved to examine the evolutionary properties of each model.
with the physics-based model, the relative importance of each term in the folding and the binding functions was evaluated.
together, these analyses give a picture of the evolutionary performance of each model and the appropriateness for use in addressing various questions in molecular evolution.
the native state in physics-based and informational scoring functions
in this study a novel physics-based scoring function for a coarse-grained representation of protein structure  was developed.
the model  includes terms for backbone bond angles, helix and sheet propensities, a new disulfide bridge term, a lennard-jones potential to account for van der waals interactions, a coloumbic potential, and a solvation potential.
the last three terms were also included as part of the binding score function.
the newly developed physics-based potential energy function  for a coarse grained protein representation , was evaluated on structure fit to sequence  and sequence fit to structure .
a previously developed informational scoring function  <cit> , which is widely used for evolutionary studies was similarly tested.
both models produce an energy gap zfold  between native and random conformations for the native sequence for  <dig> proteins in a structurally diverse test set.
also, both models produce larger energy gaps for the native sequence than random sequences in the native conformation from the same set of proteins: an average native zfold of  <dig>  for the physics-based model and  <dig>  for the informational model.
hence, both model types solve the fold recognition and inverse folding problems to some degree.
binding function can be scored similarly , and a gap zbind  is produced for both models with an average of  <dig>  for the informational model and  <dig>  for the physics-based model.
not surprisingly, both models score folding  more specifically than binding.
natural protein sequences are not near the modeled thermodynamic minimum
however, evaluation of folding stability  of near-native sequences for a representative protein sap  <cit>  shows that neither model is highly specific for the native sequence .
within 15% divergence from the native sequence, 42% of sequences are more stable than the native under the physics-based model, and for the informational model this rises to 73%.
this is not consistent with biological protein sequences, for which the vast majority of mutations are destabilizing  <cit> .
the sampling procedure also uncovers sequences up to  <dig> × more stable than the native within 30% divergence .
the landscapes for both models are qualitatively similar: the native sequence sits on a steep slope and there are local minima elsewhere .
the physics-based function yields a somewhat funnel-like landscape, whereas the informational model produces a surface with multiple similar minima.
in both cases the overall slope and large fraction of adaptive changes would be expected to cause directional selection away from the native state.
simulations corroborate this by showing rapid divergence from the native state for both models .
sequences near the modeled thermodynamic minimum have protein-like properties
as the native sequence of sap is not near a stable equilibrium for folding stability on the landscape, stability-biased markov chain monte carlo  sampling  was used to find such minima.
as seen in figure  <dig> , these highly stable sequences are very dissimilar to the native state, with only 5% and 11% average identity for the physics-based and informational models, respectively.
however, both models show high conservation at many positions and a marked signal for hydrophobicity.
the sequences show a rugged energy landscape across sequence space and expectedly complex patterns of change relative to the native sequence .
in the informational model the sequences are dominated by stretches of hydrophobic residues, interspersed with charged residues and a few tyrosines .
the remaining positions are highly randomized and do not maintain a specific biochemical character.
when mapping the most informative positions onto the structure  a fairly simple pattern emerges.
the β-sheet core is dominated by small hydrophobic residues, surrounded by a shell of larger such residues with the occasional tyrosine embedded.
clusters of charged residues maintain the stability of some loops, although most of the highly exposed positions are randomized.
the physics-based model produces qualitatively similar results, but differs in some aspects.
core residues are prolines rather than leucines, and the larger hydrophobic residues surrounding the core tend to be tryptophans instead of phenylalanines .
tyrosines are again favored in semi-exposed positions.
however, in place of ionic interactions the physics-based model prefers glycines.
mapping the informative positions onto the structure  reveals that these conserved glycines are located in or near portions of the backbone that are β-sheet-like and highly exposed.
another notable difference is that the physics-based model seems more residue-specific as it conserves fewer positions overall , but those residues are typically more strongly conserved with larger co-evolutionary barriers to change.
both models impart roughly the same selective pressure overall.
evolutionary simulations near minima reproduce biological sequence properties
using the stable sequences as starting points for population-based simulations allows comparison of the models with respect to other properties than folding stability.
as can be seen in figure  <dig>  evolutionary rate ratios in both cases are comparable to those found in real orthologous proteins  and indicative of the negative selective pressure applied  <cit> .
the rates also vary as expected depending on the position of each residue in the 3d structure.
exposed surface positions evolve up to three times faster than the buried and highly constrained core residues.
residues in the binding site, which are exposed but under functional selective pressure, show an intermediate rate.
interestingly, under both models the rates in the selectively important regions  become very similar over time.
the physics-based model consistently shows lower rates throughout the structure, with the exception of initial changes in the binding site.
surface mutations are noticeably more selected against compared to the information-based approach, which is consistent with the explicit consideration of solvation in this model.
the physics-based model also restricts movement on the sequence-stability landscape to one or a few well-defined  basins, whereas the informational model appears to diffuse more neutrally across a mostly flat surface .
real protein sequences also show rate heterogeneity across positions , such that most positions evolve slowly but some evolve quite rapidly  <cit> .
the fit of such a model to the simulated sequences was found to be reasonable .
the relative support for including gamma-distributed rates is generally high using the informational model, and somewhat lower in the physics-based model.
the latter instead shows higher relative support for a model with some invariant sites among common parameters in site-independent substitution models.
the drop in support for rate heterogeneity toward the end of the simulations may indicate emergence of rate heterotachy  <cit> , caused by migration to a non-native stable basin through intramolecular co-evolution.
this is consistent with a view of proteins where an initial equal rates across sites model with a small number of substitutions gives rise to a gamma distribution as substiutions accumulate, driven by the three dimensional structure.
with further substitution, intramolecular co-evolution gives rise to support for a covarion process rather than for a gamma distribution  <cit> .
when gamma is supported, the shape parameter of the gamma distribution  is consistent with those found in biological sequences in general  <cit>  and with the shape parameter  <dig>  obtained from the pfam sh <dig> domain seed alignment in particular.
simulation time is listed in generations.
core and surface hydrophobicities shown reflect the population average percentage of hydrophobic residues in buried and exposed positions.
gamma support indicates the fractional support for a model including a gamma distribution of rates, and the shape parameter shows the α value of such a distribution.
the sequences and structures produced by each simulation  further highlight the differences between the models.
the sequence collection from the informational model with truncation selection on folding and binding  shows the same basic features as those found by sampling from stable structures  but with added emphasis on tyrosines rather than ionic interactions.
there is some indication of increased conservation of the binding site .
during the simulations core hydrophobicity increased by ~25%, and surface hydrophobicity almost tripled .
this resulted in a similar ratio of core and surface hydrophobicities as in the native sap sequence .
by mapping conserved positions onto the structure  we recover the starting pattern  of conserved β-sheet core, a shell of larger hydrophobic residues around this, and tyrosines at medium exposure.
as noted above, the physics-based model sampled with truncation selection for folding and binding is substantially more restrictive in how many non-synonymous substitutions it allows, and the simulated sequences reflect this .
nearly 1/ <dig> of the positions are effectively invariant over shorter evolutionary timeframes, compared to only a few residues in the informational model.
structural mapping of highly conserved positions  again recovers the starting pattern , although a few β-sheet-associated glycines have been replaced by charges on the most exposed loops.
the binding site is very highly conserved.
only very exposed positions  have started to randomize, mostly containing polar and charged residues.
core and surface hydrophobicities started close to those of the native sequence , and did not change appreciably during the simulation .
comparing the simulated sequences to the hidden markov model of the pfam  <cit>  sh <dig> family , some similarities are noticeable.
both models reproduce the overall conservation of the core residues, particularly the very buried ones.
the informational model is more accurate in predicting the residue type within the core , whereas the physics-based model correctly conserves some glycines.
neither model captures the defining functional feature of the family , but this is due to the sap sh <dig> domain binding a non-phosphorylated ligand.
instead the models more generally preserve residues within the binding interface.
this type of conservation is absent in the hmm since sh <dig> domains have unique specificity-determining binding site geometries.
similarly, since the hmm reflects several dozen distinct structures and sequences it is less conserved overall than simulated sequences derived from a single structure evolved over relatively short evolutionary distances.
components of the physics-based model have variable importance for scoring
the parameters of both models would be expected to have an influence on the sequences produced.
the properties of the informational model have been described elsewhere  <cit> .
here, the nature of the physics-based model is examined in more detail.
the term weights are re-adjusted for each individual protein in the test set, and the resulting distributions are shown in figure 8a.
solvation, side chain angles  and pair-wise interactions  are quite important, consistent with a well-packed hydrophobic core with biological side chain rotamers and no steric clashes.
buried salt bridges  and disulfide bonds  seem less influential, but are also comparatively rare in the investigated protein structures.
of the secondary structure terms, helix propensity is somewhat more specific than β-sheet tendency.
for binding, shape complementarity  is by far the most effective term.
however, these parameter values are found in the context of the complete scoring function and may have strong synergistic effects in combination with other terms.
an examination of the contribution of each term in isolation to the total gap between native and random states  reveals some additional trends.
the lj term is nearly uninformative on its own, indicating that it needs to be combined with at least one other term to be specific.
the difference between representations of the non-native state is visible in the helix and beta terms.
the helix term is more informative when considering the gap to decoy structures, whereas the beta term is more effective in discriminating against random sequences.
repeating the parameterization step many times on the sap protein reveals further differences between the terms.
although the algorithm converges on roughly the same optimality score , there is substantial variation in parameter values.
the co-dependence of roughly equally likely solution sets for term weights is shown in figure  <dig>  some terms  have pronounced optima around which there is little variation, but others  produce a relatively flat optimality landscape.
particularly the less important terms have weights that vary over several orders of magnitude.
this can cause difficulty with finding realistic stable sequences for a particular fold, as small changes in some terms  can have large effects on the sequence composition.
a representative sample of parameter sets was used to generate collections of stable sequences, and the most protein-like set was selected .
the parameter values fall close to the center of each distribution, although charges are de-emphasized and solvation has increased importance.
discussion
in evaluating our ability to model the evolution of sequences using structural and functional  constraints, both informational and physical models show protein-like properties, but also need further improvement.
despite a stability gap between random sequences  from the native for both methods, the native sequence is not near enough to a thermodynamic stability optimum to prevent directional selection.
sequences near such an optimum are more homopolymeric in composition.
the evolution of such sequences, however, does show many attributes of real protein evolution.
the effect of negative  selection, a dn/ds ratio smaller than  <dig>  is well modeled by both approaches, even reproducing the known correlation between exposure and evolutionary rate  <cit> .
the physics-based model produces the stronger selective pressure overall, particularly so on the protein surface .
variation in this rate across positions is mimicked by the informational model, which more consistently supports a gamma parameter throughout the simulation.
the physics-based model instead conserves some positions completely over shorter evolutionary distances, but ultimately supports a covarion process rather than simply rate heterogeneity.
in particular, hydrophobic residue content is much closer to that of the native sequence.
when it comes to the exact identity of such hydrophobic residues  the informational model has higher correspondence to native residues, presumably due to being derived from known sequence/structure combinations.
overall the physics-based model appears more specific for local but crucial contributions to free energy of folding, whereas the informational model produces sequences with protein-like properties but without fold specificity.
in the context of fold specificity, a major difference between models is consideration of secondary structure.
the physics-based model is enriched in prolines in the β-sheet core, and glycines in other sheet-like portions of the backbone.
in principle any small hydrophobic residue should suffice in the tightly packed and highly buried core positions, but prolines in particular score well in the helical term  <cit>  and are mostly acceptable in the sheet term.
the glycines in half-exposed positions with sheet-like geometry provide an excellent compromise between the helix, sheet and solvation terms, and have little effect on the remainder of the scoring function.
although somewhat unrealistic in the context of the native sequence, this does demonstrate the strong impact of the particular backbone conformation under selection.
the physics-based model is also different in that its parameters are adjusted for each protein.
this is beneficial in adding specificity to the model  <cit> , but also presents issues with parameter selection  <cit> .
the lennard-jones term in particular is surprisingly variable and offers little benefit in isolation.
when testing against random sequences, this is mostly due to the side chain replacement protocol.
it minimizes a lennard-jones function  <cit> , thereby removing much of the signal a priori.
side chain optimization with a complete scoring function including other terms might remedy this problem.
the vastly larger utility in a binding context , where no such pre-minimization occurs, further illustrates this point.
the test against random compact structures with constant sequence  is more illuminating.
it shows that the l-j term needs the restrictions imposed by other factors  to efficiently identify a protein-like heteropolymer.
in other words, it adds information to the total function, but only when the solution space is already somewhat restricted.
generally, the parameterization of both folding and binding scoring functions is fairly difficult.
for instance, the choice of the reference state  has a substantial impact on the calculations.
the ruggedness of the optimality landscape in parameter space makes it quite difficult to find the global optimum, and parameterizing folding and binding functions separately may lead to inconsistencies between the two.
future work will explore more sophisticated approaches to solving these problems.
the specificity of either model is ultimately reflected by the shape of its stability landscape.
as seen both near the native sequence  and non-native optima , the physics-based model produces a landscape that is more rugged and has fewer equivalent minima.
this is consistent with theoretical expectations of the real free energy landscape in both conformational  <cit>  and sequence  <cit>  space.
the result is that the simulated sequences tend to be more conserved, stay near the starting  point during simulation, and show mostly deleterious or neutral mutations.
in contrast, under the informational model populations simply diffuse across a relatively flat sequence landscape with less co-evolutionary pressure between interacting sites, rapidly moving away from the starting point.
in this aspect the physics-based model shows more protein-like and fold-specific behavior overall, despite selecting sequences with some compositional "quirks", as discussed above.
informational methods, especially those built only on pairwise contacts, necessarily lack specificity.
for example, all lys-phe interactions are not equivalent.
those oriented in a plane will be repulsive, while those that are orthogonal will generate a cation-pi effect and be attractive.
this is but one example.
physical models, at the right level of granularity, will not suffer from this problem and are thereby inherently attractive for this interaction specificity.
further, they enable study of the forces driving evolutionary processes at a physical level.
why then are the native states so far from sequence optima?
the answer is found in the approximations made by the models.
packing is perhaps the most severely affected, as can be seen from the accumulation of large hydrophobic residues in both models.
these substitutions fill up the empty space that results from the coarse-graining procedure.
spherical cβ beads, for instance, are a particularly poor approximation for elongated  or flattened  residues.
a structural model with more than one bead per side chain, such as that of hills and co-workers  <cit> , could remedy this issue.
the representation of ions as spherical point charges is also problematic, and salt bridge formation could be better addressed by using a model of induced dipoles  <cit> .
many entropic contributions are not considered, although the physics-based model includes some of them implicitly in the solvation term.
in addition to entropic considerations, explicit hydrogen bonding  <cit> , cation-pi interactions  <cit> , and other interactions are not considered, but would increase model complexity.
finally, an overarching problem may be found in the effect of individual mutations.
a typical disadvantageous mutation under the informational model has a very small effect on the stability score e  , and the physics-based model generally predicts even smaller stability changes.
in real proteins the average mutation is strongly deleterious  <cit> , indicating that the models are much too lenient when it comes to the cost of mutation.
appropriate weighting of contributing terms or more precise modeling of packing would improve this aspect, and make the energy landscape more realistically rugged.
early attempts at considering alternative states in the above models of folding used explicit decoys.
because of the large number of potential alternative states for a sequence, both kinetic traps and thermodynamic optima, evaluating even a large finite number of such states proved inadequate.
further, identifying the states that were most informative near the native confirmation was not clearly driven by fold sequence or geometry and it may be that the relevant structures originate from smaller regions of diverse structures.
the random contacts model  <cit>  proved powerful and efficient in sampling many contacts from diverse structures as an implicit solution.
for binding, there is a greater conundrum in that while the real biological decoys are not known, a random contacts model has the potential to be too restrictive in selecting against binding interactions that are not biologically deleterious.
for an alternative binding interaction to be deleterious, the proteins need to have a deleterious interaction  and to have the potential to interact at the same time, in the same cellular location, at the right concentration  <cit> .
however, the proper decoys are frequently not known and this restricts the use of the alternative approach of explicit binding decoys.
further, when evaluating binding specificity as a function of affinity it is assumed that complex thermodynamic stability  must be maximized.
while this is a relevant constraint for cellular structures such as the nuclear pore complex  <cit> , protein complexes involved in signaling or even in forming the cytoskeleton have transient associations that are necessary for functional signaling  <cit> .
a requirement of transient binding means that dissociation kinetics are selected to be just fast enough, not as slow  as possible.
therefore, even in the limit of perfect predictions one should not expect to recover the native binding site sequence based on thermodynamics alone.
from the perspective of protein structural space  <cit> , the folding decoys all lie near the native conformation.
the total amount of secondary structure is conserved, as are the relative percentages of helices and sheets.
this means that the decoys can be thought of as being sampled from the same superfamily, or at least fold, in the scop hierarchy  <cit> .
conservation of residue exposures and sequence length also guarantees that the conformations have roughly the same radius of gyration as the native state.
this should cause selection against subtle misfolds and accumulating pathway intermediates, but may be a poor representation of the diversity of conformations a sequence may fold into and thus may miss important local alternative conformations that are close in energy.
for example, an increase in residues with high helical propensity could shift the lowest-energy conformation toward an all-α topology instead of α+β.
the current approach does not allow modeling this type of event, but evaluation of the partition function of the boltzmann probability of folding is a notoriously difficult  problem  <cit> .
a critical aspect of characterizing the energy landscape is the role of decoy structures in restricting the funnel shape as well as providing fold specificity, and this is a challenging task.
recently, the utility of site-interdependent structural models of sequence evolution was evaluated in the context of phylogenetics  <cit> .
despite use of a variety of scoring functions and probability measures, the utility of such models was found to be generally low.
minimizing a coarse-grained potential failed to recover much of the native sequence and did not outperform site-independent models  <cit> .
simulating the effects of substitutions on folding energy revealed a rapid divergence from the native state even when modeling solvation in addition to residue-residue contacts  <cit> , and including the contribution of folding energy to transition path probability can even decrease tree reconstruction accuracy  <cit> .
these results are all consistent with the nature of the coarse-grained energy landscape seen above.
there is no global minimum near the native sequence, and its position on a steep slope results in rapid divergence in terms of sequence, energy or both.
models that accurately characterize the structural effects of substitutions are urgently needed to make progress in this direction.
it is also worth noting that these approaches all relied on informational potentials, which we have shown above to be less structurally specific than physics-based scoring.
with models that properly characterize the evolution of proteins, many important evolutionary biological questions can be addressed.
what are the patterns of sequence evolution associated with different mechanisms of duplicate gene retention?
relatedly, how easily does orthologous neofunctionalization occur and how dependent upon protein fold and binding interface size is it?
what are the roles of positive and negative pleiotropy in restricting neofunctionalization?
how do structural transitions occur between neighboring folds?
at the contact level, how are different stabilizing interactions interconverted ?
what is the interplay between population size, fold distribution, and functional evolution?
to address these important evolutionary biology questions, continued progress on models such as described here is needed.
CONCLUSIONS
proteins evolve according to the laws of physical chemistry, biochemistry, and population genetics.
classic phenomenological site independent models, while computationally simple to use, do not adequately describe evolutionary processes for either phylogenetic or functional evolutionary analysis  purposes.
mechanistic models built upon the type described here will be necessary.
biological  pathways dictate the level of constraint on physical constants and such considerations will also need to be taken.
while the currently available models remain inadequate, an understanding of the problems in physical chemistry associated with evolutionary models will lead to future improvements, with many downstream applications.

 protein-coding sequences play a central role in cellular functions necessary to produce the vast variation in organismal phenotypes observed in nature. to function, most proteins must fold into a unique and stable structure. the structure is responsible for orienting residues necessary for proper function, including binding and catalysis. to maintain fitness, protein function and the underlying structure must be preserved. protein structure and binding are determined by the interactions of individual amino acid residues. therefore, to mechanistically model the evolution of protein sequences, these residue-residue interactions must be considered, relaxing the common assumption of independent evolution of each amino acid position. structural models for sequence evolution where function is protein-protein intermolecular interaction will now be considered. to evaluate if a sequence will fold into a unique and stable structure, it must be demonstrated that the sequence prefers the folded state to both the unfolded state as well as to alternative folded states. from a practical perspective, it is impossible to enumerate the enormous space of all possible alternative conformations. therefore, sets of representative "decoy conformations" must be used to approximate folding intermediates and/or other stable or meta-stable states. alternatively, without explicitly considering decoy structures, the inter-residue contacts from different conformations can be randomly sampled. additionally, a scoring function is needed to quantify differences between the possible states. strategies developed for generating such scoring functions include those derived statistically from existing structures and those derived from physical first principles. both types of scoring function assume that the native sequence lies close to a thermodynamic optimum, and that there is a gap in energy between the native state and the closest non-native state. formally, the protein sequence evolution problem relates to the inverse folding problem rather than the folding problem. underlying both problems are similar physical assumptions and similar models. protein biological function depends not only on proper folding, but also on the ability to bind the target ligands. therefore, to study sequence evolution, binding must also be evaluated. from a physical perspective this is easier, as the only states to consider are the bound and the unbound states. however, it has recently been suggested that selective pressures on proteins to avoid non-specific binding are also an important aspect of protein fitness, necessitating consideration of binding "decoys" as well. identifying the actual selective pressures against non-specific binding in a cell is a daunting task, but the number and nature of binding decoys is a major determinant of the ease of evolving new binding interactions. the binding decoy characteristics also affect the level of selective constraint on the binding interface. site and function independent models clearly do not capture critical elements of protein evolution. for structure-based models, a good scoring function must produce sequences with similar properties to real proteins. this includes a hydrophobic core that evolves slower than the hydrophilic surface. a dn/ds ratio much smaller than unity, particularly for functional residues, should also emerge. more particularly, amino acid substitution rates should show heterogeneity, expressed by a gamma distribution across positions, reflecting their relative importance to folding and function. while the gamma function for rate heterogeneity was adopted for model fit rather than mechanistic purposes, it is one of the most common parameters in standard evolutionary models and accounts for heterogeneity in the substitution rate driven by structural signals as well as other signals in evolutionary sequence data. in addition to rate heterogeneity, sequences must also have an energy gap between the native and alternative conformations to ensure rapid and stable folding, although structurally disordered proteins represent a distinct class of proteins that do not have this property. the model should place the native sequence near an optimum so as not to provide a signal of directional selection when function is not changing. consequently, most mutations should be deleterious or nearly neutral rather than adaptive. lastly, proteins must retain their binding function. population size dictates what fitness changes are neutral as well as what mutations become substitutions and must therefore be taken into account. to enable use in forward and backward studies of evolution, especially in a population genetic context, in a complete genome context, or in studies of interactome evolution, the model needs to be coarse grained at a level that allows for sufficient computational speed. in this context the computational cost makes state of the art models and methods from the protein structure community intractable. here a novel physics-based scoring function is developed and compared with a commonly used informational approach on the criteria described above. the physics-based model is more specific and predicts less drastic stability changes. both types of models violate the assumption of high native sequence stability and produce many adaptive mutations through directional selection towards a scoring function optimum, even though the model works with truncation selection. these properties are discussed in light of their effect on modeling accuracy and improvements to the models are suggested.  protein folding and function cause interdependence between sites in the protein sequence. for instance, a mutation that removes a cysteine involved in a disulfide bridge is likely to be deleterious, whereas mutation of other cysteines may be more neutral. of course, disulfide bridges can be lost and interconverted to other types of stabilizing interactions over evolutionary time, albeit with a large transitional barrier due to the strong site interdependence. such selective pressures can be modeled by scoring the thermodynamic effects of mutations. two approaches to calculating such scores are the statistically motivated informational methods and the first principles physics-based methods. informational models score the likelihood of observing specific types of contacts in the folded protein based on those seen in previously known structures. physics-based approaches evaluate structures by measuring the fit of residues to geometrical properties such as backbone torsion and residue-residue distances. in both cases the fit to both the native and the many possible non-native conformations must be measured to ensure specificity. this manuscript sets out to evaluate the ability to replicate biological properties of proteins using both approaches described above through evolutionary simulation. a new physics-based approach is described. native sequences were evaluated in both approaches to evaluate if they fell within the neutral network of an optimum or if they underwent directional selection initially. similarly, sequences sampled from within the optimum were evolved to examine the evolutionary properties of each model. with the physics-based model, the relative importance of each term in the folding and the binding functions was evaluated. together, these analyses give a picture of the evolutionary performance of each model and the appropriateness for use in addressing various questions in molecular evolution. the native state in physics-based and informational scoring functions in this study a novel physics-based scoring function for a coarse-grained representation of protein structure was developed. the model includes terms for backbone bond angles, helix and sheet propensities, a new disulfide bridge term, a lennard-jones potential to account for van der waals interactions, a coloumbic potential, and a solvation potential. the last three terms were also included as part of the binding score function. the newly developed physics-based potential energy function for a coarse grained protein representation, was evaluated on structure fit to sequence and sequence fit to structure. a previously developed informational scoring function, which is widely used for evolutionary studies was similarly tested. both models produce an energy gap zfold between native and random conformations for the native sequence for proteins in a structurally diverse test set. also, both models produce larger energy gaps for the native sequence than random sequences in the native conformation from the same set of proteins: an average native zfold of for the physics-based model and for the informational model. hence, both model types solve the fold recognition and inverse folding problems to some degree. binding function can be scored similarly, and a gap zbind is produced for both models with an average of for the informational model and for the physics-based model. not surprisingly, both models score folding more specifically than binding. natural protein sequences are not near the modeled thermodynamic minimum however, evaluation of folding stability of near-native sequences for a representative protein sap shows that neither model is highly specific for the native sequence. within 15% divergence from the native sequence, 42% of sequences are more stable than the native under the physics-based model, and for the informational model this rises to 73%. this is not consistent with biological protein sequences, for which the vast majority of mutations are destabilizing. the sampling procedure also uncovers sequences up to × more stable than the native within 30% divergence. the landscapes for both models are qualitatively similar: the native sequence sits on a steep slope and there are local minima elsewhere. the physics-based function yields a somewhat funnel-like landscape, whereas the informational model produces a surface with multiple similar minima. in both cases the overall slope and large fraction of adaptive changes would be expected to cause directional selection away from the native state. simulations corroborate this by showing rapid divergence from the native state for both models. sequences near the modeled thermodynamic minimum have protein-like properties as the native sequence of sap is not near a stable equilibrium for folding stability on the landscape, stability-biased markov chain monte carlo sampling was used to find such minima. as seen in figure, these highly stable sequences are very dissimilar to the native state, with only 5% and 11% average identity for the physics-based and informational models, respectively. however, both models show high conservation at many positions and a marked signal for hydrophobicity. the sequences show a rugged energy landscape across sequence space and expectedly complex patterns of change relative to the native sequence. in the informational model the sequences are dominated by stretches of hydrophobic residues, interspersed with charged residues and a few tyrosines. the remaining positions are highly randomized and do not maintain a specific biochemical character. when mapping the most informative positions onto the structure a fairly simple pattern emerges. the β-sheet core is dominated by small hydrophobic residues, surrounded by a shell of larger such residues with the occasional tyrosine embedded. clusters of charged residues maintain the stability of some loops, although most of the highly exposed positions are randomized. the physics-based model produces qualitatively similar results, but differs in some aspects. core residues are prolines rather than leucines, and the larger hydrophobic residues surrounding the core tend to be tryptophans instead of phenylalanines. tyrosines are again favored in semi-exposed positions. however, in place of ionic interactions the physics-based model prefers glycines. mapping the informative positions onto the structure reveals that these conserved glycines are located in or near portions of the backbone that are β-sheet-like and highly exposed. another notable difference is that the physics-based model seems more residue-specific as it conserves fewer positions overall, but those residues are typically more strongly conserved with larger co-evolutionary barriers to change. both models impart roughly the same selective pressure overall. evolutionary simulations near minima reproduce biological sequence properties using the stable sequences as starting points for population-based simulations allows comparison of the models with respect to other properties than folding stability. as can be seen in figure evolutionary rate ratios in both cases are comparable to those found in real orthologous proteins and indicative of the negative selective pressure applied. the rates also vary as expected depending on the position of each residue in the 3d structure. exposed surface positions evolve up to three times faster than the buried and highly constrained core residues. residues in the binding site, which are exposed but under functional selective pressure, show an intermediate rate. interestingly, under both models the rates in the selectively important regions become very similar over time. the physics-based model consistently shows lower rates throughout the structure, with the exception of initial changes in the binding site. surface mutations are noticeably more selected against compared to the information-based approach, which is consistent with the explicit consideration of solvation in this model. the physics-based model also restricts movement on the sequence-stability landscape to one or a few well-defined basins, whereas the informational model appears to diffuse more neutrally across a mostly flat surface. real protein sequences also show rate heterogeneity across positions, such that most positions evolve slowly but some evolve quite rapidly. the fit of such a model to the simulated sequences was found to be reasonable. the relative support for including gamma-distributed rates is generally high using the informational model, and somewhat lower in the physics-based model. the latter instead shows higher relative support for a model with some invariant sites among common parameters in site-independent substitution models. the drop in support for rate heterogeneity toward the end of the simulations may indicate emergence of rate heterotachy, caused by migration to a non-native stable basin through intramolecular co-evolution. this is consistent with a view of proteins where an initial equal rates across sites model with a small number of substitutions gives rise to a gamma distribution as substiutions accumulate, driven by the three dimensional structure. with further substitution, intramolecular co-evolution gives rise to support for a covarion process rather than for a gamma distribution. when gamma is supported, the shape parameter of the gamma distribution is consistent with those found in biological sequences in general and with the shape parameter obtained from the pfam sh domain seed alignment in particular. simulation time is listed in generations. core and surface hydrophobicities shown reflect the population average percentage of hydrophobic residues in buried and exposed positions. gamma support indicates the fractional support for a model including a gamma distribution of rates, and the shape parameter shows the α value of such a distribution. the sequences and structures produced by each simulation further highlight the differences between the models. the sequence collection from the informational model with truncation selection on folding and binding shows the same basic features as those found by sampling from stable structures but with added emphasis on tyrosines rather than ionic interactions. there is some indication of increased conservation of the binding site. during the simulations core hydrophobicity increased by ~25%, and surface hydrophobicity almost tripled. this resulted in a similar ratio of core and surface hydrophobicities as in the native sap sequence. by mapping conserved positions onto the structure we recover the starting pattern of conserved β-sheet core, a shell of larger hydrophobic residues around this, and tyrosines at medium exposure. as noted above, the physics-based model sampled with truncation selection for folding and binding is substantially more restrictive in how many non-synonymous substitutions it allows, and the simulated sequences reflect this. nearly 1/ of the positions are effectively invariant over shorter evolutionary timeframes, compared to only a few residues in the informational model. structural mapping of highly conserved positions again recovers the starting pattern, although a few β-sheet-associated glycines have been replaced by charges on the most exposed loops. the binding site is very highly conserved. only very exposed positions have started to randomize, mostly containing polar and charged residues. core and surface hydrophobicities started close to those of the native sequence, and did not change appreciably during the simulation. comparing the simulated sequences to the hidden markov model of the pfam sh family, some similarities are noticeable. both models reproduce the overall conservation of the core residues, particularly the very buried ones. the informational model is more accurate in predicting the residue type within the core, whereas the physics-based model correctly conserves some glycines. neither model captures the defining functional feature of the family, but this is due to the sap sh domain binding a non-phosphorylated ligand. instead the models more generally preserve residues within the binding interface. this type of conservation is absent in the hmm since sh domains have unique specificity-determining binding site geometries. similarly, since the hmm reflects several dozen distinct structures and sequences it is less conserved overall than simulated sequences derived from a single structure evolved over relatively short evolutionary distances. components of the physics-based model have variable importance for scoring the parameters of both models would be expected to have an influence on the sequences produced. the properties of the informational model have been described elsewhere. here, the nature of the physics-based model is examined in more detail. the term weights are re-adjusted for each individual protein in the test set, and the resulting distributions are shown in figure 8a. solvation, side chain angles and pair-wise interactions are quite important, consistent with a well-packed hydrophobic core with biological side chain rotamers and no steric clashes. buried salt bridges and disulfide bonds seem less influential, but are also comparatively rare in the investigated protein structures. of the secondary structure terms, helix propensity is somewhat more specific than β-sheet tendency. for binding, shape complementarity is by far the most effective term. however, these parameter values are found in the context of the complete scoring function and may have strong synergistic effects in combination with other terms. an examination of the contribution of each term in isolation to the total gap between native and random states reveals some additional trends. the lj term is nearly uninformative on its own, indicating that it needs to be combined with at least one other term to be specific. the difference between representations of the non-native state is visible in the helix and beta terms. the helix term is more informative when considering the gap to decoy structures, whereas the beta term is more effective in discriminating against random sequences. repeating the parameterization step many times on the sap protein reveals further differences between the terms. although the algorithm converges on roughly the same optimality score, there is substantial variation in parameter values. the co-dependence of roughly equally likely solution sets for term weights is shown in figure some terms have pronounced optima around which there is little variation, but others produce a relatively flat optimality landscape. particularly the less important terms have weights that vary over several orders of magnitude. this can cause difficulty with finding realistic stable sequences for a particular fold, as small changes in some terms can have large effects on the sequence composition. a representative sample of parameter sets was used to generate collections of stable sequences, and the most protein-like set was selected. the parameter values fall close to the center of each distribution, although charges are de-emphasized and solvation has increased importance. discussion in evaluating our ability to model the evolution of sequences using structural and functional constraints, both informational and physical models show protein-like properties, but also need further improvement. despite a stability gap between random sequences from the native for both methods, the native sequence is not near enough to a thermodynamic stability optimum to prevent directional selection. sequences near such an optimum are more homopolymeric in composition. the evolution of such sequences, however, does show many attributes of real protein evolution. the effect of negative selection, a dn/ds ratio smaller than is well modeled by both approaches, even reproducing the known correlation between exposure and evolutionary rate. the physics-based model produces the stronger selective pressure overall, particularly so on the protein surface. variation in this rate across positions is mimicked by the informational model, which more consistently supports a gamma parameter throughout the simulation. the physics-based model instead conserves some positions completely over shorter evolutionary distances, but ultimately supports a covarion process rather than simply rate heterogeneity. in particular, hydrophobic residue content is much closer to that of the native sequence. when it comes to the exact identity of such hydrophobic residues the informational model has higher correspondence to native residues, presumably due to being derived from known sequence/structure combinations. overall the physics-based model appears more specific for local but crucial contributions to free energy of folding, whereas the informational model produces sequences with protein-like properties but without fold specificity. in the context of fold specificity, a major difference between models is consideration of secondary structure. the physics-based model is enriched in prolines in the β-sheet core, and glycines in other sheet-like portions of the backbone. in principle any small hydrophobic residue should suffice in the tightly packed and highly buried core positions, but prolines in particular score well in the helical term and are mostly acceptable in the sheet term. the glycines in half-exposed positions with sheet-like geometry provide an excellent compromise between the helix, sheet and solvation terms, and have little effect on the remainder of the scoring function. although somewhat unrealistic in the context of the native sequence, this does demonstrate the strong impact of the particular backbone conformation under selection. the physics-based model is also different in that its parameters are adjusted for each protein. this is beneficial in adding specificity to the model, but also presents issues with parameter selection. the lennard-jones term in particular is surprisingly variable and offers little benefit in isolation. when testing against random sequences, this is mostly due to the side chain replacement protocol. it minimizes a lennard-jones function, thereby removing much of the signal a priori. side chain optimization with a complete scoring function including other terms might remedy this problem. the vastly larger utility in a binding context, where no such pre-minimization occurs, further illustrates this point. the test against random compact structures with constant sequence is more illuminating. it shows that the l-j term needs the restrictions imposed by other factors to efficiently identify a protein-like heteropolymer. in other words, it adds information to the total function, but only when the solution space is already somewhat restricted. generally, the parameterization of both folding and binding scoring functions is fairly difficult. for instance, the choice of the reference state has a substantial impact on the calculations. the ruggedness of the optimality landscape in parameter space makes it quite difficult to find the global optimum, and parameterizing folding and binding functions separately may lead to inconsistencies between the two. future work will explore more sophisticated approaches to solving these problems. the specificity of either model is ultimately reflected by the shape of its stability landscape. as seen both near the native sequence and non-native optima, the physics-based model produces a landscape that is more rugged and has fewer equivalent minima. this is consistent with theoretical expectations of the real free energy landscape in both conformational and sequence space. the result is that the simulated sequences tend to be more conserved, stay near the starting point during simulation, and show mostly deleterious or neutral mutations. in contrast, under the informational model populations simply diffuse across a relatively flat sequence landscape with less co-evolutionary pressure between interacting sites, rapidly moving away from the starting point. in this aspect the physics-based model shows more protein-like and fold-specific behavior overall, despite selecting sequences with some compositional "quirks", as discussed above. informational methods, especially those built only on pairwise contacts, necessarily lack specificity. for example, all lys-phe interactions are not equivalent. those oriented in a plane will be repulsive, while those that are orthogonal will generate a cation-pi effect and be attractive. this is but one example. physical models, at the right level of granularity, will not suffer from this problem and are thereby inherently attractive for this interaction specificity. further, they enable study of the forces driving evolutionary processes at a physical level. why then are the native states so far from sequence optima? the answer is found in the approximations made by the models. packing is perhaps the most severely affected, as can be seen from the accumulation of large hydrophobic residues in both models. these substitutions fill up the empty space that results from the coarse-graining procedure. spherical cβ beads, for instance, are a particularly poor approximation for elongated or flattened residues. a structural model with more than one bead per side chain, such as that of hills and co-workers, could remedy this issue. the representation of ions as spherical point charges is also problematic, and salt bridge formation could be better addressed by using a model of induced dipoles. many entropic contributions are not considered, although the physics-based model includes some of them implicitly in the solvation term. in addition to entropic considerations, explicit hydrogen bonding, cation-pi interactions, and other interactions are not considered, but would increase model complexity. finally, an overarching problem may be found in the effect of individual mutations. a typical disadvantageous mutation under the informational model has a very small effect on the stability score e, and the physics-based model generally predicts even smaller stability changes. in real proteins the average mutation is strongly deleterious, indicating that the models are much too lenient when it comes to the cost of mutation. appropriate weighting of contributing terms or more precise modeling of packing would improve this aspect, and make the energy landscape more realistically rugged. early attempts at considering alternative states in the above models of folding used explicit decoys. because of the large number of potential alternative states for a sequence, both kinetic traps and thermodynamic optima, evaluating even a large finite number of such states proved inadequate. further, identifying the states that were most informative near the native confirmation was not clearly driven by fold sequence or geometry and it may be that the relevant structures originate from smaller regions of diverse structures. the random contacts model proved powerful and efficient in sampling many contacts from diverse structures as an implicit solution. for binding, there is a greater conundrum in that while the real biological decoys are not known, a random contacts model has the potential to be too restrictive in selecting against binding interactions that are not biologically deleterious. for an alternative binding interaction to be deleterious, the proteins need to have a deleterious interaction and to have the potential to interact at the same time, in the same cellular location, at the right concentration. however, the proper decoys are frequently not known and this restricts the use of the alternative approach of explicit binding decoys. further, when evaluating binding specificity as a function of affinity it is assumed that complex thermodynamic stability must be maximized. while this is a relevant constraint for cellular structures such as the nuclear pore complex, protein complexes involved in signaling or even in forming the cytoskeleton have transient associations that are necessary for functional signaling. a requirement of transient binding means that dissociation kinetics are selected to be just fast enough, not as slow as possible. therefore, even in the limit of perfect predictions one should not expect to recover the native binding site sequence based on thermodynamics alone. from the perspective of protein structural space, the folding decoys all lie near the native conformation. the total amount of secondary structure is conserved, as are the relative percentages of helices and sheets. this means that the decoys can be thought of as being sampled from the same superfamily, or at least fold, in the scop hierarchy. conservation of residue exposures and sequence length also guarantees that the conformations have roughly the same radius of gyration as the native state. this should cause selection against subtle misfolds and accumulating pathway intermediates, but may be a poor representation of the diversity of conformations a sequence may fold into and thus may miss important local alternative conformations that are close in energy. for example, an increase in residues with high helical propensity could shift the lowest-energy conformation toward an all-α topology instead of α+β. the current approach does not allow modeling this type of event, but evaluation of the partition function of the boltzmann probability of folding is a notoriously difficult problem. a critical aspect of characterizing the energy landscape is the role of decoy structures in restricting the funnel shape as well as providing fold specificity, and this is a challenging task. recently, the utility of site-interdependent structural models of sequence evolution was evaluated in the context of phylogenetics. despite use of a variety of scoring functions and probability measures, the utility of such models was found to be generally low. minimizing a coarse-grained potential failed to recover much of the native sequence and did not outperform site-independent models. simulating the effects of substitutions on folding energy revealed a rapid divergence from the native state even when modeling solvation in addition to residue-residue contacts, and including the contribution of folding energy to transition path probability can even decrease tree reconstruction accuracy. these results are all consistent with the nature of the coarse-grained energy landscape seen above. there is no global minimum near the native sequence, and its position on a steep slope results in rapid divergence in terms of sequence, energy or both. models that accurately characterize the structural effects of substitutions are urgently needed to make progress in this direction. it is also worth noting that these approaches all relied on informational potentials, which we have shown above to be less structurally specific than physics-based scoring. with models that properly characterize the evolution of proteins, many important evolutionary biological questions can be addressed. what are the patterns of sequence evolution associated with different mechanisms of duplicate gene retention? relatedly, how easily does orthologous neofunctionalization occur and how dependent upon protein fold and binding interface size is it? what are the roles of positive and negative pleiotropy in restricting neofunctionalization? how do structural transitions occur between neighboring folds? at the contact level, how are different stabilizing interactions interconverted? what is the interplay between population size, fold distribution, and functional evolution? to address these important evolutionary biology questions, continued progress on models such as described here is needed.  proteins evolve according to the laws of physical chemistry, biochemistry, and population genetics. classic phenomenological site independent models, while computationally simple to use, do not adequately describe evolutionary processes for either phylogenetic or functional evolutionary analysis purposes. mechanistic models built upon the type described here will be necessary. biological pathways dictate the level of constraint on physical constants and such considerations will also need to be taken. while the currently available models remain inadequate, an understanding of the problems in physical chemistry associated with evolutionary models will lead to future improvements, with many downstream applications. 
sampling and simulation show that the physics-based model is more specific for fold-defining interactions that vary less among residue type. we find that both models violate the assumption of the native sequence being close to a thermodynamic optimum, causing directional selection away from the native state. the informational model diffuses further in sequence space with fewer barriers and tends to provide less support for an invariant sites model, although amino acid substitutions are generally conservative. this is partly due to improper packing of the hydrophobic core. we suggest possible improvements on the representation of structure, folding energy, and binding function, as regards both native and non-native conformations, and describe a large number of possible applications for such a model.
 protein sequence evolution is constrained by the biophysics of folding and function, causing interdependence between interacting sites in the sequence. however, current site-independent models of sequence evolutions do not take this into account. recent attempts to integrate the influence of structure and biophysics into phylogenetic models via statistical/informational approaches have not resulted in expected improvements in model performance. this suggests that further innovations are needed for progress in this field.  here we develop a coarse-grained physics-based model of protein folding and binding function, and compare it to a popular informational model. we find that both models violate the assumption of the native sequence being close to a thermodynamic optimum, causing directional selection away from the native state. sampling and simulation show that the physics-based model is more specific for fold-defining interactions that vary less among residue type. the informational model diffuses further in sequence space with fewer barriers and tends to provide less support for an invariant sites model, although amino acid substitutions are generally conservative. both approaches produce sequences with natural features like dn/ds < and gamma-distributed rates across sites.  simple coarse-grained models of protein folding can describe some natural features of evolving proteins but are currently not accurate enough to use in evolutionary inference. this is partly due to improper packing of the hydrophobic core. we suggest possible improvements on the representation of structure, folding energy, and binding function, as regards both native and non-native conformations, and describe a large number of possible applications for such a model. 