### Elastic network models

Let the conformation of an N-sites protein be represented by the column vector of the 3 N Cartesian coordinates of its N *C*_{
α
} atoms: *r* = (*x*_{1} *y*_{1} *z*_{1} *x*_{2} *y*_{2} *z*_{2} … *x*_{
N
} *y*_{
N
} *z*_{
N
} )^{T}. *r*_{
i
} = (*x*_{
i
} *y*_{
i
} *z*_{
i
})^{T} is the position vector of the *i*th *C*_{
α
}. The vector joining sites *i* and *j* is *d*_{
ij
} = *r*_{
j
} - *r*_{
i
} with length *d*_{
ij
} = *d*_{
ij
}. We use *r*^{0} for the protein’s equilibrium conformation in which the *i*th site is at {\mathit{r}}_{\mathit{i}}^{0}.

An Elastic Network Model (ENM) represents the folded protein as a network of sites connected by springs. They have proved accurate and useful in a variety of applications [17, 19]. The potential energy landscape is given by:

\mathit{V}\left(\mathit{r}\right)=\frac{1}{2}{\displaystyle \sum _{\mathit{i}=1}^{\mathit{N}-1}}{\displaystyle \sum _{\mathit{j}=\mathit{i}+1}^{\mathit{N}}}{\mathit{k}}_{\mathit{ij}}{\left({\mathit{d}}_{\mathit{ij}}-{\mathit{d}}_{\mathit{ij}}^{0}\right)}^{2}

(1)

where {\mathit{d}}_{\mathit{ij}}^{0} and *k*_{
ij
} are, respectively, the equilibrium length and force constant of spring *ij*. As far as we know, all models proposed so far assume that {\mathit{d}}_{\mathit{ij}}^{0}={\mathit{d}}_{\mathit{ij}}\left({\mathit{r}}^{0}\right)={\mathit{r}}_{\mathit{j}}^{0}-{\mathit{r}}_{\mathit{i}}^{0}, i.e. that at the equilibrium conformation *r*^{0}, all springs are relaxed.

### Fluctuations and flexibility

No protein is frozen at its equilibrium conformation. At non-zero absolute temperature, the folded protein fluctuates around *r*^{0} sampling conformational space with equilibrium Boltzmann’s probability density function:

\mathit{\rho}\left(\mathit{r}\right)=\frac{{\mathit{e}}^{-\mathit{\beta V}\left(\mathit{r}\right)}}{{\mathit{Z}}_{\mathit{F}}}

(2)

where \mathit{\beta}=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{${\mathit{k}}_{\mathit{B}}\mathit{T}$}\right., with *T* the absolute temperature and *k*_{
B
} Boltzmann’s constant. The denominator of Eq. (2) is the partition function of the folded protein:

{\mathit{Z}}_{\mathit{F}}={\displaystyle \int}{\mathit{e}}^{-\mathit{\beta V}\left(\mathit{r}\right)}\mathit{d\tau}

(3)

where ∫ … *dτ* stands for integration over the whole of conformational space.

The dynamical flexibility (mobility) of a site is ordinarily quantified using its Mean Square Fluctuation:

\mathit{MS}{\mathit{F}}_{\mathit{i}}\equiv \u3008{\u2551{\mathit{r}}_{\mathit{i}}-{\mathit{r}}_{\mathit{i}}^{0}\u2551}^{2}\u3009={\displaystyle \int}{\u2551{\mathit{r}}_{\mathit{i}}-{\mathit{r}}_{\mathit{i}}^{0}\u2551}^{2}\mathit{\rho}\left(\mathit{r}\right)\mathit{d\tau}

(4)

To calculate *MSF*_{
i
} using Eq. (4), the potential energy function Eq. (1) is approximated using a second-order Taylor expansion around its equilibrium conformation. First, the Hessian matrix *H* of second derivatives of the potential Eq. (1) with respect to the atoms’ Cartesian coordinates is calculated. Then, *H* is inverted to obtain the 3*N* × 3*N* variance-covariance matrix *C*, which is composed of a 3 × 3 *C*_{
ij
} block for each pair of sites. Finally, a site’s MSF is given by [20]:

\mathit{MS}{\mathit{F}}_{\mathit{i}}=\mathit{Tr}\left({\mathit{C}}_{\mathit{ii}}\right)

(5)

### An empirical flexibility model

Several studies have investigated the correlation between site-specific rates of evolution or other sequence-variability measures and the corresponding flexibility. Since such studies use Pearson’s correlation coefficients as measure of association, the underlying assumption is that there is a linear relationship between rate of evolution and flexibility. To make such assumption explicit, here we postulate the following Flexibility Model:

{\tilde{{\mathit{\omega}}_{\mathit{i}}}}^{\mathit{FLEX}}={\mathit{a}}_{\mathit{P}}^{\mathit{FLEX}}+{\mathit{b}}_{\mathit{P}}^{\mathit{FLEX}}\tilde{\mathit{MS}{\mathit{F}}_{\mathit{i}}}

(6)

where \tilde{{\mathit{\omega}}_{\mathit{i}}}\phantom{\rule{0.25em}{0ex}} is the relative rate of substitution of the *i*th site. In general, for site-specific scalar properties we will use relative values obtained by z-score normalization. For any given site-specific property *x*_{
i
}, we the z-score normalized values are \tilde{{\mathit{x}}_{\mathit{i}}}=\raisebox{1ex}{$\left({\mathit{x}}_{\mathit{i}}-\u3008\mathit{x}\u3009\right)$}\!\left/ \!\raisebox{-1ex}{$\sqrt{\u3008{\mathit{x}}^{2}\u3009-\u3008{\mathit{x}}^{2}\u3009}$}\right., where the averages are calculated over all sites of the same protein. The subscript *P* is used to note that *a priori* the coefficients may depend on the protein considered. We emphasize that the Flexibility Model is empirical: rather than derived from first principles, it is postulated, based on the intuitive notion that flexible sites should accommodate mutations more easily.

### A mechanistic stress model

We introduce here a mechanistic model that includes explicitly the effects of mutations and natural selection. We consider mutations as random perturbations of the wild-type ENM potential [21–23]. A random mutation at site *i* results in a mutant whose potential *V*_{
mut
} is obtained from Eq. (1) by adding perturbations to the equilibrium length of each of its springs: {\mathit{d}}_{\mathit{ij}}^{0}\to {\mathit{d}}_{\mathit{ij}}^{0}+{\mathit{\delta}}_{\mathit{ij}}. We further assume that the springs are independently perturbed and that perturbations are spring-independent, randomly drawn from a distribution with zero mean and constant variance *α*^{2}:

\u3008{\mathit{\delta}}_{\mathit{ij}}\u3009=0;\u3008{\mathit{\delta}}_{\mathit{ij}}^{2}\u3009={\mathit{\alpha}}^{2}.

(7)

As we mentioned above, when the wild type is at its equilibrium conformation {\mathit{r}}_{\mathit{wt}}^{0}, all springs are relaxed by construction. In contrast, when the mutant is at {\mathit{r}}_{\mathit{wt}}^{0}, the mutated site’s springs will be stressed (stretched or compressed). For further reference, we define the Mean Local mutational Stress (MLmS) as follows:

\mathit{MLm}{\mathit{S}}_{\mathit{i}}\equiv {\u3008{\mathit{V}}_{\mathit{mut}}\left({\mathit{r}}_{\mathit{wt}}^{0}\right)-{\mathit{V}}_{\mathit{wt}}\left({\mathit{r}}_{\mathit{wt}}^{0}\right)\u3009}_{\mathit{mut}@\mathit{i}}

(8)

where 〈 … 〉_{mut @ i} stands for averaging over random mutations at the *i*th site.

To complete the model, we derive a simple selection function. First, we assume that there is a single specific active conformation *r*_{
active
}. Next, we acknowledge fluctuations and assume that the protein’s activity (either the wild-type’s or a mutant’s) is proportional to the concentration of the active conformation *r*_{
active
}. Finally, we assume that {\mathit{r}}_{\mathit{active}}={\mathit{r}}_{\mathit{wt}}^{0} and, accordingly, we model the acceptance probability of a mutant as:

{\mathit{p}}^{\mathit{accept}}\equiv \frac{{\mathit{C}}_{\mathit{mut}}^{\mathit{F}}{\mathit{\rho}}_{\mathit{mut}}\left({\mathit{r}}_{\mathit{wt}}^{0}\right)}{{\mathit{C}}_{\mathit{wt}}^{\mathit{F}}{\mathit{\rho}}_{\mathit{wt}}\left({\mathit{r}}_{\mathit{wt}}^{0}\right)}

(9)

Where {\mathit{C}}_{\mathit{mut}}^{\mathit{F}} and {\mathit{C}}_{\mathit{wt}}^{\mathit{F}} are the concentrations of folded protein for the mutant and wild type, respectively. From statistical mechanics, the Folded-Unfolded equilibrium constants for the wild-type and mutant proteins are, respectively, \raisebox{1ex}{${\mathit{C}}_{\mathit{wt}}^{\mathit{F}}$}\!\left/ \!\raisebox{-1ex}{${\mathit{C}}_{\mathit{wt}}^{\mathit{U}}$}\right.=\raisebox{1ex}{${\mathit{Z}}_{\mathit{wt}}^{\mathit{F}}$}\!\left/ \!\raisebox{-1ex}{${\mathit{Z}}_{\mathit{wt}}^{\mathit{U}}$}\right. and \raisebox{1ex}{${\mathit{C}}_{\mathit{mut}}^{\mathit{F}}$}\!\left/ \!\raisebox{-1ex}{${\mathit{C}}_{\mathit{mut}}^{\mathit{U}}$}\right.=\raisebox{1ex}{${\mathit{Z}}_{\mathit{mut}}^{\mathit{F}}$}\!\left/ \!\raisebox{-1ex}{${\mathit{Z}}_{\mathit{mut}}^{\mathit{U}}$}\right.. We further assume that the partition function and concentration of unfolded protein is the same for the mutant and wild type. Therefore \raisebox{1ex}{${\mathit{C}}_{\mathit{mut}}^{\mathit{F}}$}\!\left/ \!\raisebox{-1ex}{${\mathit{C}}_{\mathit{wt}}^{\mathit{F}}$}\right.=\raisebox{1ex}{${\mathit{Z}}_{\mathit{mut}}$}\!\left/ \!\raisebox{-1ex}{${\mathit{Z}}_{\mathit{wt}}$}\right.. Replacing this relationship and Eq. (2) into Eq. (9) we find:

{\mathit{p}}^{\mathit{accept}}={\mathit{e}}^{-\mathit{\beta}\phantom{\rule{0.25em}{0ex}}\left[{\mathit{V}}_{\mathit{mut}}\left({\mathit{r}}_{\mathit{wt}}^{0}\right)-{\mathit{V}}_{\mathit{wt}}\left({\mathit{r}}_{\mathit{wt}}^{0}\right)\right]}

(10)

Finally, averaging over random mutations at site *i* and using Eq. (8) we obtain the *acceptance rate*:

{\mathit{\omega}}_{\mathit{i}}\equiv {\u3008{\mathit{p}}_{\mathit{i}}^{\mathit{accept}}\u3009}_{\mathit{mut}@\mathit{i}}={\u3008{\mathit{e}}^{-\mathit{\beta}\phantom{\rule{0.25em}{0ex}}\left[{\mathit{V}}_{\mathit{mut}}\left({\mathit{r}}_{\mathit{wt}}^{0}\right)-{\mathit{V}}_{\mathit{wt}}\left({\mathit{r}}_{\mathit{wt}}^{0}\right)\right]}\u3009}_{\mathit{mut}@\mathit{i}}\cong 1-\mathit{\beta}\phantom{\rule{0.25em}{0ex}}\mathit{MLm}{\mathit{S}}_{\mathit{i}}

(11)

Where *β* may be thought of as representing not just temperature but also selection pressure, and we have assumed that *β*Δ*V* << 1 (mild selection) to approximate the exponential to first order. To finish, we z-normalize the variables of Eq. (11) to get the relative rates of evolution:

{\tilde{{\mathit{\omega}}_{\mathit{i}}}}^{\mathit{STRESS}}={\mathit{a}}_{\mathit{P}}^{\mathit{STRESS}}+{\mathit{b}}_{\mathit{P}}^{\mathit{STRESS}}\tilde{\mathit{MLm}{\mathit{S}}_{\mathit{i}}}.

(12)

This equation specifies the stress model.

### Relationship of flexibility and stress with packing density

The purpose of this work is to investigate why LPD correlates with rate of evolution at site level. The previous models relate rates of evolution with MSF (Eq. 6) and MLmS (Eq. 12). Here we derive the relationship between these properties and LPD measures.

First, we relate flexibility and stress with the potential energy parameters of Eq. (1). Let us define:

{\mathit{k}}_{\mathit{i}}\equiv {\displaystyle \sum _{\mathit{j}\ne \mathit{i}}}{\mathit{k}}_{\mathit{ij}}

(13)

Regarding flexibility, replacing Eqs. (1), (2), and (3) into Eq. (4), following [12], and using Eq. (13), it can be found that:

\mathit{MS}{\mathit{F}}_{\mathit{i}}\cong \frac{3}{2\mathit{\beta}{\mathit{k}}_{\mathit{i}}}

(14)

Regarding stress, from Eqs. (1), (7), and (8), after some algebra, we get:

\mathit{MLm}{\mathit{S}}_{\mathit{i}}=\frac{1}{2}{\mathit{\alpha}}^{2}{\mathit{k}}_{\mathit{i}}

(15)

Note that Eq. (14) is an approximation while Eq. (15) is an identity.

Second, to relate the previous models to LPD we need to specify the ENMs spring constants. A variety of ENMs have been developed (see [24] for a recent comparison). Here, we consider two models. First, the “parameter-free Anisotropic Network Model” (pfANM) [25], which uses:

{\mathit{k}}_{\mathit{ij}}=\frac{1}{{\left({\mathit{d}}_{\mathit{ij}}^{0}\right)}^{2}}

(16)

Second, the “Anisotropic Network Model” (ANM) [20], for which:

{\mathit{k}}_{\mathit{ij}}=\left\{\begin{array}{cc}\hfill 1\hfill & \hfill {\mathit{d}}_{\mathit{ij}}^{0}\le {\mathit{R}}_{\mathit{cut}}\hfill \\ \hfill 0\hfill & \hfill {\mathit{d}}_{\mathit{ij}}^{0}>{\mathit{R}}_{\mathit{cut}}\hfill \end{array}\right.

(17)

where *R*_{
cut
} is typically between 10 Å and 18 Å.

From Eqs. (13), (16), and (17) and z-normalizing we find:

\tilde{{\mathit{k}}_{\mathit{i}}}=\tilde{\mathit{LP}{\mathit{D}}_{\mathit{i}}}

(18)

where for the pfANM, LPD is the Weighted Contact Number (WCN) of [26], and for the ANM, it is the Contact Number (CN): the number of sites closer than *R*_{
cut
}. Finally, from Eqs. (14) and (18) it follows:

\tilde{\mathit{MS}{\mathit{F}}_{\mathit{i}}}\cong {\tilde{{\mathit{k}}_{\mathit{i}}}}^{-1}={\tilde{\mathit{LP}{\mathit{D}}_{\mathit{i}}}}^{-1}

(19)

Similarly, from Eqs. (15) and (18) we get:

\tilde{\mathit{MLm}{\mathit{S}}_{\mathit{i}}}=\tilde{{\mathit{k}}_{\mathit{i}}}=\tilde{\mathit{LP}{\mathit{D}}_{\mathit{i}}}

(20)

Note that while MSF is approximately equal to 1/LPD, MLmS is exactly equal to LPD (for relative z-normalized values).

### Calculation details

We used the dataset of 213 monomeric enzymes of Yeh et al. [11]. The dataset includes proteins of diverse sizes, functional, and structural classes (Additional file 1: Table S1).

We used the evolutionary rates of [11]. They were inferred from multiple alignments of homologous sequences using Rate4Site, which builds the phylogenetic tree using a neighbour-joining algorithm and estimates rates with an empirical Bayesian approach and the JTT model of sequence evolution [27, 28]. To keep in mind that we are not dealing with the (unknown) “true rates”, but with Rate4Site-inferred rates, we use the notation {\tilde{\mathit{\omega}}}_{\mathit{i}}^{\mathit{R}4\mathit{S}}.

From the pdb equilibrium structure of each protein we calculated the spring constants of pfANM (Eq. 16) and ANM (Eq. 17), for which we used a cut-off distance of 13 Å [11]. Given a protein and ENM model, we calculated the Hessian matrix, inverted it to obtain the variance-covariance matrix, and calculated the site-specific flexibility values \tilde{\mathit{MS}{\mathit{F}}_{\mathit{i}}} using Eq. (5) and z-normalizing. Regarding stress, we obtained the relative site-specific values \tilde{\mathit{MLm}{\mathit{S}}_{\mathit{i}}} using Eq. (15) and z-normalizing.

Since we always use z-normalized relative values, for the sake of notational simplicity, we shall use *ω*^{R4S}, *MSF*, and *MLmS* to refer to z-normalized values from now on.

We performed two analyses. In a protein-by-protein analysis, we performed linear fits of *ω*^{R4S} with either *MSF* (Flexibility Model) or *MLmS* (Stress Model) using the lm() function of the base package of R for each protein. In a global analysis we pooled together all sites of all proteins and performed similar global fits.

To assess the goodness-of-fit of a model to the data, we used the Akaike Information Criterion *AIC* = 2*k* - 2 ln *L*, where *k* is the number of parameters and *L* is the model’s likelihood given the data. When comparing models, the AIC weight of evidence for each model is given by \mathit{w}\left(\mathit{AIC}\right)\propto {\mathit{e}}^{\mathit{\u2012}\frac{1}{2}\mathrm{\Delta}\left(\mathit{AIC}\right)}, where Δ(*AIC*) = *AIC* ‒ min(*AIC*) [29, 30].

We also calculated Pearson’s correlation coefficients between evolutionary rates and the independent variable that defines each model. When comparing two models, we calculated partial correlation coefficients of evolutionary rates with the independent variable of each model controlling that of the other.