We introduced a mechanistic Stress Model of protein sequence evolution. Mutations are modelled as random perturbations of the protein’s potential energy landscape, represented using Elastic Network Models. To model natural selection, we used basic statistical physics to derive the expected probability that a mutant samples a specific functional structure. From this, we deduced a linear relationship between a site’s mean evolutionary rate and the mean local mutational stress (MLmS) of the functional conformation. We compared this model with an empirical Flexibility Model that postulates that a site’s evolutionary rate is linearly dependent on its flexibility (measured by its MSF). We compared both models and found strong support for the Stress Model. Moreover, the independent contribution of flexibility is negligible once stress is controlled for.

The MLmS is proportional to Local Packing Density and, therefore, the Stress Model provides a mechanism for the connection between a site’s LPD and its evolutionary rate. Regarding the sequence-flexibility relationship, previous empirical correlation studies had already found that the sequence-flexibility relationship is nonlinear and either dismissed the nonlinear parts or attempted an interpretation in terms of different selection regimes [14, 32, 13]. We found the nonlinearity follows naturally from the Stress Model: evolutionary rates depend nonlinearly on *MSF* because they depend (approximately) linearly on *MLmS*, and *MSF* ≈ 1/*MLmS*. To summarize, the Stress Model accounts for the observed site-dependency of evolutionary rates and its relationship with packing density and flexibility.

A note of caution is in order here. For the Stress Model mutational stress was not postulated to be the determinant factor *a priori* but, rather, it was derived from the assumptions of the model that are essentially two (1) there is a single active conformation and (2) mutants are flexible and therefore can sample the active conformation so that they are at least partly functional. Therefore even though Stress Model was chosen to designate this mechanistic model, it should be kept in mind that it demonstrates the importance of protein flexibility.

It is worthwhile to mention some of the possible caveats and further developments of the Stress Model. First, we assume a single active conformation. In principle, it would be reasonable to assume that only changes of the active-site conformation should affect fitness. However, we note that if protein sites are strongly coupled, which is often the case, any conformational change will affect the active site conformation. For a strongly coupled elastic network forcing the active site to adopt a given conformation makes the rest of the protein move accordingly. Therefore, assuming that the whole protein conformation must be in the “active conformation” for the protein to function is not necessarily an important limitation. However, for cases where the coupling is not very strong, if the active site is known, this could be easily tackled using a modified version of the selection function that integrates away all coordinates except for those of the active site (i.e. uses marginal conformational distributions rather than the full ones in the definition of selection function).

Second, in Eq. (11) we performed a linear approximation of the exponential function. This is reasonable *a priori* only for weak selection, and *a posteriori* by the good fit of the resulting model to the data. We should note, however, this approximation can be easily removed, and the actual mean of the exponential can be calculated via simulation. Further work is needed to explore this possibility.

Third, we note that the z-normalized MLmS values, on which the Stress Model is based, are identical to the z-normalized LPD measures WCN (for the pfANM potential) and CN (for the ANM potential). For other potentials this need not be the case and it is for that reason that we chose to keep the notation MLmS in the present tables and figures, to make them comparable with further research based on estimating MLmS using different, perhaps better, potential energy functions.

To close, we note that the mutational part of the Stress Model accounts for observed patterns of evolutionary divergence of protein structure and dynamics [21–23]. Regarding structural divergence, unselected random mutations reproduce very well the evolutionary conservation of a “structural core” and account for the observation that structures diverge mainly within the space spanned by a few low-energy collective normal modes [21, 22]. Regarding protein motions, unselected random mutations explain the higher conservation of the low-energy normal modes in terms of their mutational robustness [31, 23]. In general, those studies could found no evidence of natural selection at the levels of structural or dynamical divergence. Clearly, without natural selection, all sites would evolve at the same rate, which is not the case. The Stress Model proposed here accounts rather well for the variation of rates of evolution among sites. It would be interesting to study the effect of the selection function introduced here on structural and dynamical divergence and compare the observed patterns with those that result from unselected mutations. This could advance our understanding of the effect of selection at the levels of structure and dynamics. In general, we think the Stress Model provides a possible unifying framework to study evolutionary protein divergence at the levels of sequence, structure, and dynamics.