Too packed to change: side-chain packing and site-specific substitution rates in protein evolution

Theory

In this section, we show that the mechanistic stress model of protein evolution predicts that the substitution rate of a protein site is determined by the packing density of its side chain. This prediction and its empirical assessment are the point of this paper.

The stress model was proposed by Huang et al. (2014) to explain the observed correlation between site-specific substitution rates and packing density. The model is based on the idea that a mutant is viable to the extent that it spends time in the active conformation. In turn, this time will depend on mutational changes of the stability of the active conformation. The fixation probability of a mutant is modeled as (1)${\displaystyle p_{\text{fix}}\propto \frac{C_{\text{mut}}^F{\rho }_{\text{mut}}\left({\mathbf{r}}_{\text{active}}\right)}{C_{\text{wt}}^F{\rho }_{\text{wt}}\left({\mathbf{r}}_{\text{active}}\right)}}$ where wt stands for wild-type, mut for mutant, C^F is the concentration of folded protein and ρ(r_active) its probability of adopting the active conformation. Assuming that C_mut/C_wt is equal to the ratio of partition functions, from basic statistical physics it follows that: (2)${\displaystyle p_{\text{fix}}\propto e^{-\beta \delta V^{\ast }},}$ where β represents the selection pressure and (3)${\displaystyle \delta V^{\ast }=V_{\text{mut}}\left({\mathbf{r}}_{\text{active}}\right)-V_{\text{wt}}\left({\mathbf{r}}_{\text{active}}\right)}$ is the energy difference between mutant and wild-type in the active conformation.

Assuming that βδV^∗ ≪ 1 (weak selection), from (2) we find: (4)${\displaystyle K^i\propto -\langle \delta {V{}^{\ast }\rangle }^i,}$ i.e., the rate of substitution of site i, Kⁱ, is proportional to (minus) the change in stability of the active conformation averaged over mutations at i, 〈δV^∗〉ⁱ. This is the basic equation of the stress theory.

In Huang et al. (2014), mutational stability changes were calculated using an elastic network model in which each residue is represented by a single node. Within such a one-node-per-residue representation, there is no differentiation between main chain and side chain. Therefore, we cannot predict whether evolutionary rates will be determined by main chain packing or side chain packing. To address this issue, here we represent each residue using two nodes: a main-chain node α, placed at the residue’s C_α, and a side-chain node ρ, placed at the geometric center of the residue’s side chain (Gly’s are represented using only one node at C_α). The energy function is: (5)${\displaystyle V\left(\mathbf{r}\right)=\frac{1}{2}\sum _i\sum _{j>i}{k}_{{\alpha }_i{\alpha }_j}{\left(r_{{\alpha }_i{\alpha }_j}-d_{{\alpha }_i{\alpha }_j}\right)}^2+\frac{1}{2}\sum _i\sum _{j>i}{k}_{{\alpha }_i{\rho }_j}{\left(r_{{\alpha }_i{\rho }_j}-d_{{\alpha }_i{\rho }_j}\right)}^2+\frac{1}{2}\sum _i\sum _{j>i}{k}_{{\rho }_i{\alpha }_j}{\left(r_{{\rho }_i{\alpha }_j}-d_{{\rho }_i{\alpha }_j}\right)}^2+\frac{1}{2}\sum _i\sum _{j>i}{k}_{{\rho }_i{\rho }_j}{\left(r_{{\rho }_i{\rho }_j}-d_{{\rho }_i{\rho }_j}\right)}^2,}$ where r_ninj is the distance between nodes n_i and n_j (n is α or ρ), k_ninj is the force constant of the spring connecting these nodes, and d_ninj the equilibrium spring length.

A mutation at site i will replace ρ_i, affecting only the parameters of the energy function related to this node. We emphasize: while the mutation may well induce global structural changes involving the backbone and other side chains, the only parameters that will change are those of the mutated side chain. Following Echave (2008) and Echave & Fernández (2010), we model a mutation at i by adding random perturbations to the lengths of the springs connected to ρ_i: d_ρiρj → d_ρiρj + δ_ρiρj and d_ρiαj → d_ρiαj + δ_ρiαj, to find, using (3) and (5): (6)${\displaystyle \delta V^{\ast }=\frac{1}{2}\sum _{j\ne i}\left(k_{{\rho }_i{\alpha }_j}{\delta }_{{\rho }_i{\alpha }_j}^2+k_{{\rho }_i{\rho }_j}{\delta }_{{\rho }_i{\rho }_j}^2\right).}$ Assuming that perturbations are drawn independently from the same distribution, averaging (6) over mutations at i we find: (7)${\displaystyle \langle \delta {V{}^{\ast }\rangle }^i\propto \sum _{j\ne i}\left(k_{{\rho }_i{\alpha }_j}+k_{{\rho }_i{\rho }_j}\right).}$

To finish, we assume, as in the parameter-free Anisotropic Network Model (pfANM) of Yang, Song & Jernigan (2009), that $k_{n_i{n}_j}=\frac{1}{d_{n_i{n}_j}^2}$. Then, from (4) and (7) we obtain: (8)${\displaystyle K^i\propto -{\text{WCN}}_{\rho }^{\alpha \rho },}$ where (9)${\displaystyle {\text{WCN}}_{\rho }^{\alpha \rho }=\sum _{j\ne i}\left(\frac{1}{d_{{\rho }_i{\alpha }_j}^2}+\frac{1}{d_{{\rho }_i{\rho }_j}^2}\right).}$ ${\text{WCN}}_{\rho }^{\alpha \rho }$, derived here, is the side-chain weighted contact number. It depends on contacts between node ρ of the site considered (subscript) and nodes α and ρ of the other sites (superscript). Therefore, the stress model, combined with a two-nodes-per-site pfANM energy function, predicts that site-specific rates will depend on the contact density of the side chain ${\text{WCN}}_{\rho }^{\alpha \rho }$.

By analogy with (9) we can calculate the main-chain weighted contact number: (10)${\displaystyle {\text{WCN}}_{\alpha }^{\alpha \rho }=\sum _{j\ne i}\left(\frac{1}{d_{{\alpha }_i{\alpha }_j}^2}+\frac{1}{d_{{\alpha }_i{\rho }_j}^2}\right).}$ We expect ${\text{WCN}}_{\alpha }^{\alpha \rho }$ to correlate with ${\text{WCN}}_{\rho }^{\alpha \rho }$, which may result in indirect correlations with substitution rates. However, if the stress model is correct, rates will be determined only by ${\text{WCN}}_{\rho }^{\alpha \rho }$ and there should not be any independent effect of ${\text{WCN}}_{\alpha }^{\alpha \rho }$.

Other structural predictors

To assess the prediction of the previous section, we also consider the following structural properties. First, the Weighted Contact Number WCN, which was introduced by Lin et al. (2008) and found to be among the best structural predictors of site-dependent evolutionary rates (Yeh et al., 2014a; Yeh et al., 2014b). It is defined as: (11)${\displaystyle \text{WCN}={\text{WCN}}_{\alpha }^{\alpha }=\sum _{j\ne i}\frac{1}{d_{{\alpha }_i{\alpha }_j}^2}}$ where d_αiαj is the distance between the the alpha carbons of sites i and j. For the sake of clarity, wherever it is convenient we will use the notation ${\text{WCN}}_{\alpha }^{\alpha }$ to make explicit that the distances between the C_α of a site (subscript) and the C_α of the other sites (superscript) are considered. Therefore, ${\text{WCN}}_{\alpha }^{\alpha }$ can be considered a measure of main-chain packing density (based only on C_α–C_α interactions).

Second, by analogy with (11) we can use side-chain centers of mass ρ rather than C_α to define: (12)${\displaystyle {\text{WCN}}_{\rho }^{\rho }=\sum _{j\ne i}\frac{1}{d_{{\rho }_i{\rho }_j}^2}}$ ${\text{WCN}}_{\rho }^{\rho }$ quantifies the packing density of the side chain including only ρ–ρ interactions.

Finally, we also consider the Relative Solvent Accessibility, RSA, which is the most studied structural determinant of evolutionary rates. The RSA of a residue is obtained by dividing its area accessible to the solvent (SA) by the maximum SA for the given amino acid type (Tien et al., 2013).

Dataset and empirical substitution rates

To test our theory, we used the data set of Echave, Jackson & Wilke (2015). The set consists of 209 monomeric enzymes of known structure covering diverse structural and functional classes. Each structure is accompanied by up to 300 homologous sequences.

We used the empirical site-specific rates of evolution of Echave, Jackson & Wilke (2015). They were calculated as follows. First, the homologous sequences for each structure were aligned using MAFFT (Multiple Alignment using Fast Fourier Transform) (Katoh et al., 2005; Katoh & Standley, 2013). Second, using the resulting alignments as input, Maximum Likelihood phylogenetic trees were inferred with RAxML (Randomized Axelerated Maximum Likelihood), using the LG substitution matrix (named after Le and Gascuel) and the CAT model of rate heterogeneity (Stamatakis, 2014). Third, the alignment and phylogenetic tree for each structure was used as input of Rate4Site to obtain the site-specific rates of substitution using the empirical Bayesian method and the amino-acid Jukes-Cantor mutational model (aaJC) (Mayrose et al., 2004). Finally, site-specific relative rates were obtained by dividing site-specific rates by their average over all sites of the protein. We denote the empirical rates by K_R4S.

Comparison of empirical rates with structural properties

For each protein of the dataset, we used the pdb structure to calculate the five site-dependent structural properties defined above: ${\text{WCN}}_{\rho }^{\alpha \rho }$, ${\text{WCN}}_{\alpha }^{\alpha \rho }$, ${\text{WCN}}_{\rho }^{\rho }$, ${\text{WCN}}_{\alpha }^{\alpha }$ (=WCN), and RSA. For a given predictor x, we quantified its predictive power using the squared Pearson correlation coefficient R²(K_R4S, x). According to the theoretical predictions, ${\text{WCN}}_{\rho }^{\alpha \rho }$ should be the sole determinant of site-specific rates. We quantified the independent contribution of each of the other structural descriptors by partialing out the effect of ${\text{WCN}}_{\rho }^{\alpha \rho }$ using semipartial correlations. The squared semipartial correlation ${\rho }^2\left(K_{\text{R4S}},x|{\text{WCN}}_{\rho }^{\alpha \rho }\right)$ represents the unique contribution of predictor x. Also, it is the amount by which the explained variation of K_R4S (R²) would increase when going from the single-variable linear fit $K\sim {\text{WCN}}_{\rho }^{\alpha \rho }$ to the two-variable fit $K\sim {\text{WCN}}_{\rho }^{\alpha \rho }+x$. Expected values, standard deviations, and p-values were obtained by averaging protein correlations and semipartial correlations for 10,000 bootstrapped replicas of the dataset of 209 proteins.

For statistical analysis we used R (R Core Team, 2014). Correlation coefficients and their p-values were calculated using cor.test(). Semipartial correlation coefficients and p-values were calculated using spcor.test(). For bootstrapping with used boot() with default options.

Side-chain contact density (${\text{WCN}}_{\rho }^{\alpha \rho }$) vs. main-chain contact density (${\text{WCN}}_{\alpha }^{\alpha \rho }$)

According to the stress model, site-specific substitution rates depend only on side-chain packing, so that main-chain packing should not be directly related to substitution rates. To test this prediction, we compared empirical substitution rates K_R4S with ${\text{WCN}}_{\rho }^{\alpha \rho }$ and ${\text{WCN}}_{\alpha }^{\alpha \rho }$.

Consider, for example, Human Carbonic Anhidrase II (pdb code 1CA2). As we mentioned in the Introduction, main chain environments and side-chain environments are different (Fig. 1). Accordingly, ${\text{WCN}}_{\rho }^{\alpha \rho }$ and ${\text{WCN}}_{\alpha }^{\alpha \rho }$ result in different predicted rates (Fig. 2). The two site-dependent profiles of predicted rates are similar to the empirical K_R4S profile. However, ${\text{WCN}}_{\rho }^{\alpha \rho }$-based predictions look somewhat better (Fig. 2) and are better (Fig. 3): the R² values are 0.56 for ${\text{WCN}}_{\rho }^{\alpha \rho }$ and 0.41 for ${\text{WCN}}_{\alpha }^{\alpha \rho }$. Thus, for 1CA2 ${\text{WCN}}_{\rho }^{\alpha \rho }$ outperforms ${\text{WCN}}_{\alpha }^{\alpha \rho }$ as predictor of evolutionary rates.

We repeated the previous assessment for each of the 209 enzymes of the data set (Fig. 4). Empirical rates correlate with ${\text{WCN}}_{\rho }^{\alpha \rho }$ better than with ${\text{WCN}}_{\alpha }^{\alpha \rho }$ for 204 of the 209 proteins studied. Aggregating the data over all proteins, we obtained expected R² values of 0.392 ± 0.008 and 0.319 ± 0.008 for ${\text{WCN}}_{\rho }^{\alpha \rho }$ and ${\text{WCN}}_{\alpha }^{\alpha \rho }$, respectively. The difference ΔR² = 0.073 ± 0.003 is significantly positive (p < 10⁻³, bootstrapping). Therefore, ${\text{WCN}}_{\rho }^{\alpha \rho }$ outperforms ${\text{WCN}}_{\alpha }^{\alpha \rho }$ as predictor of site-specific substitution rates.

The stress model predicts ${\text{WCN}}_{\rho }^{\alpha \rho }$ to be the sole predictor of substitution rates. Any correlation between rates and ${\text{WCN}}_{\alpha }^{\alpha \rho }$ should be indirect. We measured the direct association between empirical rates and ${\text{WCN}}_{\alpha }^{\alpha \rho }$ using the squared semipartial correlation ${\rho }^2\left(K_{\text{R4S}},{\text{WCN}}_{\alpha }^{\alpha \rho }|{\text{WCN}}_{\rho }^{\alpha \rho }\right)$, where the variation of rates due to ${\text{WCN}}_{\rho }^{\alpha \rho }$ is partialed out. This measure is the unique contribution of ${\text{WCN}}_{\alpha }^{\alpha \rho }$ and it represents how much R² would increase when going from the one variable model $K\sim {\text{WCN}}_{\rho }^{\alpha \rho }$ to the two-variable model $K\sim {\text{WCN}}_{\rho }^{\alpha \rho }+{\text{WCN}}_{\alpha }^{\alpha \rho }$. Averaging over the 209 proteins studied, we found ${\rho }^2\left(K_{\text{R4S}},{\text{WCN}}_{\alpha }^{\alpha \rho }|{\text{WCN}}_{\rho }^{\alpha \rho }\right)=0.0024\pm 0.0005$. This value is statistically significant (p < 10⁻³, bootstrapping), but very small: ${\text{WCN}}_{\alpha }^{\alpha \rho }$’s unique contribution to rate variation among sites is just 0.2%. As predicted by the stress model, the independent contribution of ${\text{WCN}}_{\alpha }^{\alpha \rho }$ is negligible.

${\text{WCN}}_{\rho }^{\alpha \rho }$ vs. ${\text{WCN}}_{\rho }^{\rho }$

${\text{WCN}}_{\rho }^{\alpha \rho }$, Eq. (9), is based on a two-nodes-per-site network representation of the protein. It considers the contacts between the node ρ that represents the side chain of a site with all other nodes, ρ and α of the network. ${\text{WCN}}_{\rho }^{\rho }$, Eq. (12), is an alternative alternative measure of side-chain packing based only on ρ − ρ contacts. ${\text{WCN}}_{\rho }^{\alpha \rho }$ is a better rate predictor than ${\text{WCN}}_{\rho }^{\rho }$ for 122 of the 209 proteins. The expected correlations are $R^2\left(K_{\text{R4S}},{\text{WCN}}_{\rho }^{\alpha \rho }\right)=0.392\pm 0.008$ and $R^2\left(K_{\text{R4S}},{\text{WCN}}_{\rho }^{\rho }\right)=0.389\pm 0.008$. The average difference ΔR² = 0.0024 ± 0.007 is significant (p < 10⁻³, bootstrapping), but very small (just 0.24% of explained variation). Thus, ${\text{WCN}}_{\rho }^{\rho }$-based predictions are (almost) as good as ${\text{WCN}}_{\rho }^{\alpha \rho }$ predictions. However, while ${\text{WCN}}_{\rho }^{\rho }$ was posed ad hoc, ${\text{WCN}}_{\rho }^{\alpha \rho }$ was theoretically derived.

${\text{WCN}}_{\rho }^{\alpha \rho }$ vs. WCN

Currently, WCN ($={\text{WCN}}_{\alpha }^{\alpha }$), the original weighted contact number (Lin et al., 2008), is one of the two main structural predictors of site-dependent evolutionary rates (Yeh et al., 2014a; Yeh et al., 2014b). It is worthwhile to consider whether the new measure presented here, ${\text{WCN}}_{\rho }^{\alpha \rho }$ provides an improvement over WCN.

We found that ${\text{WCN}}_{\rho }^{\alpha \rho }$ outperforms WCN for 206 out of the 209 proteins studied. The expected correlations are $R^2\left(K_{\text{R4S}},{\text{WCN}}_{\rho }^{\alpha \rho }\right)=0.392\pm 0.008$ and R²(K_R4S, WCN) = 0.314 ± 0.007. The difference is ΔR² = 0.078 ± 0.003, which is statistically significant (p ≪ 10⁻³, bootstrapping). Thus, not only does ${\text{WCN}}_{\rho }^{\alpha \rho }$ outperform WCN for almost all proteins, but by a rather large amount: while WCN explains 31.4 % of the variation of evolutionary rates, ${\text{WCN}}_{\rho }^{\alpha \rho }$ explains 39.2%, an increase by a factor of 1.25. Moreover, ${\rho }^2\left(K_{\text{R4S}},\text{WCN}|{\text{WCN}}_{\rho }^{\alpha \rho }\right)=0.0024\pm 0.005$ (p < 10⁻³, bootstrapping). Despite statistical significance, the unique contribution of WCN is just 0.2%, which is negligible. Thus, ${\text{WCN}}_{\rho }^{\alpha \rho }$ is a better predictor and the independent contribution of WCN is negligible.

${\text{WCN}}_{\rho }^{\alpha \rho }$ vs. RSA

The most studied structural predictor of site-dependent evolutionary rates is the relative solvent accessibility RSA (Bustamante, Townsend & Hartl, 2000; Conant & Stadler, 2009; Franzosa & Xia, 2009; Ramsey et al., 2011; Shahmoradi et al., 2014). Therefore, we compare the new measure ${\text{WCN}}_{\rho }^{\alpha \rho }$ with RSA.

According to protein-by-protein results, $R^2\left(K_{\text{R4S}},{\text{WCN}}_{\rho }^{\alpha \rho }\right)>R^2\left(K_{\text{R4S}},\text{RSA}\right)$ for 175 of the 209 proteins. The expected square correlations are $R^2\left(K_{\text{R4S}},{\text{WCN}}_{\rho }^{\alpha \rho }\right)=0.392\pm 0.008$ and R²(K_R4S, RSA) = 0.327 ± 0.007. The difference is ΔR² = 0.065 ± 0.004, which is statistically significant (p < 10⁻³, bootstrapping). Thus, ${\text{WCN}}_{\rho }^{\alpha \rho }$ outperforms RSA as rate predictor for 84% of the proteins and ${\text{WCN}}_{\alpha }^{\alpha }$ explains 6.5% more of the rate variation among sites, an improvement by a factor of 1.2 over the explaining power of RSA. Moreover, the expected value of the independent contribution of RSA is $\rho \left(K_{\text{R4S}},\text{RSA}|{\text{WCN}}_{\rho }^{\alpha \rho }\right)=0.005\pm 0.001$). This is statistically significant (p < 10⁻³, bootstrapping), but very small. Therefore, ${\text{WCN}}_{\rho }^{\alpha \rho }$ is a better predictor and the independent contribution of RSA is minor.