Fast and Interpretable Estimation of Amino AcidResidue Surface Accessibility Based on Protein Contact Graph

1. Introduction

The solvent accessible surface area (SASA) of amino acid residues is a key characteristic in the analysis of the spatial structure of a protein. It is widely used to identify surface residues, predict binding sites and epitopes, and assess the quality of protein packaging [1,2]. Traditional SASA calculation methods, such as the Shrake–Rupley algorithm [3] implemented in DSSP [4] or FreeSASA [5], require complete three-dimensional protein geometry at atomic resolution. This makes them computationally expensive for large-scale analyses (e.g., virtual screening or training machine learning models on thousands of PDB files).

An alternative is to use empirical functions based on simple characteristics such as node degree in the RIN. RIN is based on Cα–Cα distances; node degree reflects the number of spatial neighbors of the residue, which is closely related to its packing and therefore inversely proportional to solvent accessibility [6]. Despite the obvious correlation between node degrees in the RIN and SASA values, practical use requires quantification of this relationship and the selection of empirical parameters. In this paper, we explore two approaches to constructing a protein contact graph and evaluate their applicability for SASA prediction.

The first approach is a classical Cα-graph with a threshold of 8 Å. Each residue is represented by its Cα atom coordinates; edges are added when the Cα–Cα distance is < 8 Å. This approach has minimal computational overhead and is widely used in the literature [6].

The second approach is a graph of heavy atoms (HAG) with a threshold of 5.0 Å. Each residue is represented by a set of coordinates of all heavy atoms (C, N, O, S, etc., excluding hydrogens), edges are added at a minimum distance between any heavy atoms of two residues < 5.0 Å.

For each of the two types of graphs, we propose and validate three empirical features for estimating relative SASA: approx_sasa (linear saturation model), surface_score (piecewise linear function) and exp_sasa (exponential model). Formal parameter optimization was carried out on an extended sample of 509 protein structures obtained from PDB with a resolution of ≤ 2.0 Å and an R-factor of ≤ 0.25. This ensures reliability and enables assessment of model generalization across diverse proteins.

In addition, we explore the possibility of improving prediction accuracy by using an extended set of 11 RIN topological features, including basic node characteristics, centrality (betweenness, eigenvector, closeness), and hydrophobic subgraph features. Training ensemble models (Random Forest, XGBoost) with these features using the LOSO (Leave-One-Structure-Out) scheme enables precise quantification of potential improvement compared to empirical formulas.

A comparison with graph neural networks (GCN, GAT, GraphSAGE) was carried out to assess the applicability of deep learning on the graph structure. The models were trained using the train/val/test scheme (70/15/15) with early stopping and hyperparameter optimization.

Machine learning methods, including neural networks and ensemble models, have been increasingly applied to predict protein properties [7,8,9,10,11,12,13,14,15,16]. However, such methods often require significant computational resources and remain "black boxes", hindering interpretation. Accordingly, we focus on simple, interpretable empirical formulas that use only the degree of the node in the contact graph, and further compare them with more complex regression models on an extended set of topological features. We show that simple formulas provide high accuracy with minimal computational overhead, and their parameters can be universally calibrated on a diverse sample of proteins. Comparing ML methods and graph types enables systematic evaluation of trade-offs for practical tool selection in structural bioinformatics.

Main results of the work:

• For the Cα-graph, the optimal parameters of the empirical formulas provide the Pearson correlation r = 0.711, MAE = 0.117.
• For HAG, the empirical formulas give r = 0.813, MAE = 0.094.
• Random Forest on HAG achieves record values: r = 0.915, MAE = 0.057, which is 39% better than the best empirical formula on a Cα graph and 13% better than Random Forest on a Cα graph.
• GraphSAGE on HAG showed a result close to Random Forest: MAE = 0.0715, Pearson r = 0.8917, which demonstrates the potential of graph neural networks when using HAG.
• GCN and GAT performed significantly worse (MAE = 0.14–0.15, Pearson r = 0.51–0.61).
• Empirical formulas are calculated in 0.008 ms per structure (~26,000× faster than FreeSASA), Random Forest in prediction mode is calculated in 36.5 ms (~6× faster than FreeSASA).
• Constructing a HAG takes 21 times longer than a Cα graph.

2. Materials and Methods

2.1. Data

For validation, a sample of 509 structures obtained from the Protein Data Bank (PDB) with a resolution of ≤ 2.0 Å and an R-factor of ≤ 0.25 was used. The sample includes proteins of various origins: viral proteins (hemagglutinin H3N2, H1N1, SARS-CoV-2 spike, HIV-1 gp120), enzymes, immunoglobulins, and others. We processed 128,794 residues.

For each structure, contact graphs of two types (Cα-graph and HAG) were constructed, node degrees and an extended set of topological features were calculated. The real relative SASA was calculated using the FreeSASA library [5] and then normalized to maximum values for each amino acid [17].

2.2. Construction of Protein Contact Graphs (RINs)

Let protein consist of amino acid residues.$N$

2.2.1. Cα-Graph (Threshold 8 Å)

Each residue is represented by the coordinates of its Cα atom. The undirected contact graph comprises a set of vertices, each corresponding to a residue, and a set of edges defined by the condition:${\mathbf{r}}_i\in {\mathbb{R}}^{3}i=1,\dots ,N{G}_{\mathrm{C}\mathsf{\alpha }}=(V,E)V=\{1,\dots ,N\}$

$$E=\{(i,j)\mid \parallel {\mathbf{r}}_i-{\mathbf{r}}_j\parallel <d_{\mathrm{cutoff}},\hspace{0.25em} i\ne j\},$$

with a threshold distance $d_{\mathrm{cutoff}}=8$Å. The 8 Å threshold is based on typical Cα–Cα contact distances corresponding to secondary structure packing [6].

Node degree$i$:

$$\mathrm{d}\mathrm{e}\mathrm{g}(i)=\mid \{j\mid (i,j)\in E\}\mid .$$

The degree reflects the number of spatial neighbors of the residue, which is inversely proportional to its availability to the solvent.

2.2.2. HAG (Threshold 5.0 Å)

In contrast to the Cα-graph, each residue is represented by the set of coordinates of all its heavy atoms $A_i$ (atoms C, N, O, S, etc., excluding hydrogens). Atomic coordinates are obtained from PDB files. The edge between the residues is added if the minimum distance between any heavy atoms of these residues is less than 5.0 Å:$ij$

$$E=\left\{\left(i,j\right)\hspace{0.25em} ∣\hspace{0.25em} \underset{a\in A_i,\hspace{0.25em} b\in A_j}{\mathrm{m}\mathrm{i}\mathrm{n}}\parallel {\mathbf{r}}_a-{\mathbf{r}}_b\parallel <5.0\stackrel{\mathrm{˚}}{\mathrm{A}}\right\}.$$

The threshold of 5.0 Å corresponds to typical van der Waals interactions and real interatomic contacts between residues.

2.3. Empirical Features

The same three empirical formulas are used for both types of graphs. The parameters are calibrated independently for the Cα-graph and HAG.

2.3.1. Feature Approx Sasa

Linear model with saturation:

$${\widehat{s}}_i=1-\min \left({\displaystyle \frac{\mathrm{d}\mathrm{e}\mathrm{g}(i)}{M}},\hspace{0.17em}1\right),M>0.$$

The parameter $M$ (max_degree) is interpreted as the threshold at which the residue is considered fully packed.

2.3.2. Feature Surface_Score

Piecewise linear model:

$${\widehat{t}}_i=\max \hspace{0.17em}\left(0,\hspace{0.25em} a-{\displaystyle \frac{\mathrm{d}\mathrm{e}\mathrm{g}(i)}{b}}\right),a>0,\hspace{0.25em} b>0.$$

Here defines the maximum feature value at deg=0, and $b$ controls the decay rate.

2.3.3. Feature Exp_Sasa

Exponential Model:

$${\widehat{u}}_i=\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{\mathrm{d}\mathrm{e}\mathrm{g}(i)}{\lambda }}\right),\lambda >0.$$

This functional form produces smooth, non-linear decay characteristic of actual SASA dependence on the number of contacts.

2.4. Calculating the True Relative SASA

The true relative SASA for the residue was calculated as$i$

$$s_i={\displaystyle \frac{{\mathrm{SASA}}_{\mathrm{abs}}(i)}{S_{\max }({\mathrm{aa}}_i)}},$$

where $SAS{A}_{\mathrm{abs}}(i)$ is the absolute SASA of the residue, obtained by summing the contributions of all its atoms using the FreeSASA library [5], and $S_{\mathrm{m}\mathrm{a}\mathrm{x}}(a{a}_i)$ is the maximum SASA for this amino acid type (tabular values from [17]). Thus, $s_i\in \left[\mathrm{0,1}\right]$.

2.5. Extended set of RIN Topological Features

To improve the accuracy of the relative SASA prediction, an extended set of 11 topological features calculated directly from the protein contact graph was used. The features are divided into three categories.

2.5.1. Basic Topological Features (6 Elements)

• The degree of the node is the number of contacts of the residue:

$$\mathrm{d}\mathrm{e}\mathrm{g}(i)=\mid \{j\mid (i,j)\in E\}\mid .$$

3.

The clustering coefficient is a measure of the density of the local environment:

$$C(i)={\displaystyle \frac{2\cdot \mid \{(j,k)\in E:j,k\in N(i)\}\mid }{\mathrm{d}\mathrm{e}\mathrm{g}(i)(\mathrm{d}\mathrm{e}\mathrm{g}(i)-1)}},$$

where N(i) is the set of neighbors of node i.

4.

The average degree of neighbors characterizes the environment of the residue:

$${\overline{\mathrm{d}\mathrm{e}\mathrm{g}}}_{\mathrm{neigh}}\left(i\right)={\displaystyle \frac{1}{\mathrm{deg}\left(i\right)}}{\sum }_{j\in N\left(i\right)}\mathrm{deg}\left(j\right).$$

6.

The k-core number is the maximum value of k such that a node belongs to a k-core (a subgraph in which all nodes have degree ≥ k).

$$\mathrm{core}(i)=\mathrm{m}\mathrm{a}\mathrm{x}\{k\mid i\in k-\mathrm{core}(G)\}.$$

• Local density is the number of Cα atoms within the radius $R_{\mathrm{local}}=12\mathrm{\AA }$

:

$${\rho }_{\mathrm{local}}(i)=\mid \{j\ne i:\parallel {\mathbf{r}}_i-{\mathbf{r}}_j\parallel <R_{\mathrm{local}}\}\mid .$$

For HAG, this feature is calculated from the same Cα coordinates as for the Cα graph. This ensures comparability between the two graph types and avoids ambiguity associated with different heavy atom packing densities.

9.

Normalized distance to the protein center:

$$d_{\mathrm{center}}(i)={\displaystyle \frac{\parallel {\mathbf{r}}_i-\stackrel{\mathbf{ˉ}}{\mathbf{r}}\parallel }{\underset{j}{\mathrm{m}\mathrm{a}\mathrm{x}}\parallel {\mathbf{r}}_j-\stackrel{\mathbf{ˉ}}{\mathbf{r}}\parallel }},\stackrel{\mathbf{ˉ}}{\mathbf{r}}={\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^N}{\mathbf{r}}_i.$$

2.5.2. Centrality (3 Elements)

• Betweenness centrality is the proportion of the shortest paths that pass through a node:

$$B(i)=\sum _{s\ne i\ne t}{\displaystyle \frac{{\sigma }_{st}(i)}{{\sigma }_{st}}},$$

where ${\sigma }_{st}$ is the total number of shortest paths between s and t, and ${\sigma }_{st}(i)$ is the number of such paths that pass through i.

2.

Eigenvector centrality is the importance of a node as determined by the importance of its neighbors:

$$x_i={\displaystyle \frac{1}{\lambda }}\sum _{j\in N(i)}{x}_j,$$

where λ is the largest eigenvalue of the adjacency matrix.

3.

Closeness centrality is the inverse of the average shortest path distance to all other nodes:

$$C_{\mathrm{close}}(i)={\displaystyle \frac{N-1}{\sum _{j\ne i}d(i,j)}}.$$

2.5.3. Features of a Hydrophobic Subgraph (2 Elements)

Hydrophobic residues (A, V, L, I, M, F, W, Y, P) [17] form the protein core. To analyze their interactions, a hydrophobic subgraph is constructed $G_{\mathrm{hydro}}=(V_{\mathrm{hydro}},E_{\mathrm{hydro}})$, where is the set of hydrophobic residues, and the edges correspond to the same contacts as in the original graph.$V_{\mathrm{hydro}}$

1.

Degree in hydrophobic subgraph:

$${\mathrm{d}\mathrm{e}\mathrm{g}}_{\mathrm{hydro}}(i)=\mid \{j\in V_{\mathrm{hydro}}:(i,j)\in E\}\mid .$$

3.

Clustering coefficient in hydrophobic subgraph:

$$C_{\mathrm{hydro}}(i)=\left\{\begin{array}{cc}0,& \mathrm{if} {\mathrm{d}\mathrm{e}\mathrm{g}}_{\mathrm{hydro}}(i)<2,\\ {\displaystyle \frac{2\cdot \mid \{(j,k)\in E:j,k\in N_{\mathrm{hydro}}(i)\}\mid }{{\mathrm{d}\mathrm{e}\mathrm{g}}_{\mathrm{hydro}}(i)\hspace{0.17em}({\mathrm{d}\mathrm{e}\mathrm{g}}_{\mathrm{hydro}}(i)-1)}},& \mathrm{otherwise}.\end{array}\right.$$

where $N_{\mathrm{hydro}}\left(i\right)=N\left(i\right)\cap V_{\mathrm{hydro}}$.

2.6. Quality Metrics

To compare the predicted values ${\hat{S}}_i$ (or ${\hat{t}}_i,{\hat{u}}_i$) with the actual values S_i, the following metrics were used:

• Pearson correlation coefficient (r)
• Spearman's rank correlation coefficient $(\rho )$
• Mean Absolute Error (MAE)
• Root Mean Square Error (RMSE)
• Mean error (bias)

2.7. Optimization of the Parameters of Empirical Formulas

Empirical formula parameters were optimized using the Leave-One-Structure-Out (LOSO) scheme. In each fold, the parameters were calibrated on the training data, and then applied to the test structure. The optimal parameters obtained were stable (standard deviation = 0) and summarized in Table 1.

2.8. ML Models and Validation Scheme

The following regression models were trained to predict relative SASA for an extended set of 11 features:

• Linear Regression
• Ridge regression (Ridge, α = 1.0)
• Random Forest (100 trees)
• XGBoost (100 trees)

Feature standardization was used for linear models. To assess generalization ability, the Leave-One-Structure-Out (LOSO) scheme was used - each protein alternately acted as a test protein, the rest of the structures were used for training. Quality metrics were averaged across all folds. Training was carried out separately for features computed on the Cα-graph and on the HAG.

2.9. Graph Neural Networks

For comparison with ML models, three graph neural network architectures were trained:

• GCN (Graph Convolutional Network) — 2 layers, 64 hidden neurons.
• GAT (Graph Attention Network) – 2 layers, 64 hidden neurons, 4 attention heads.
• GraphSAGE (Graph Sample and Aggregator) - 2 layers, 64 hidden neurons.

All models used the same 11 features as ML models. The training was conducted according to the train/val/test scheme (70/15/15) with an early stopping (patience = 30), an Adam optimizer (lr = 0.001) and an MSE loss function. BatchNorm, dropout (0.2) and L2 regularization (5e-4) were used.

3. Results

3.1. Optimal Parameters of Empirical Formulas

Calibration on 509 structures showed that the optimal parameters for the Cα graph were: M = 12 (approx_sasa), λ = 6.0 (exp_sasa), a = 1.00, b = 12.00 (surface_score). For HAG: m = 13, λ = 6.6, a = 1.25, b = 9.50. The parameters for HAG have shifted towards higher values, reflecting the denser packing in the heavy atom graph.

3.2. Aggregated Metrics of Empirical Formulas

Table 2 shows the aggregated metrics for all three models, averaged over 509 structures.

Switching to HAG improved approx_sasa by 13% on MAE (from 0.117 to 0.094) and by 14% on correlation (from 0.711 to 0.813).

3.3. Results of ML Models (LOSO)

To estimate the potential gain in accuracy when using more complex models, we trained regression algorithms on an extended set of RIN topological features. Training and validation were carried out using the Leave-One-Structure-Out (LOSO) scheme on the same sample of 509 structures. Comparison results for the Cα-graph and HAG are shown in Table 3.

Random Forest on HAG showed the best result: MAE = 0.057, Pearson r = 0.915, which is 13% better than Random Forest on a Cα graph (MAE = 0.066) and 39% better than the best empirical formula on a Cα graph (MAE = 0.117).

3.4. Comparison with Graph Neural Networks

Table 4 shows the data from experiments with several GNN architectures.

GraphSAGE on HAG showed a result close to Random Forest: MAE = 0.0715, Pearson r = 0.8917. This indicates that graph neural networks can be effective when using HAG. GCN and GAT performed significantly worse (MAE > 0.14, Pearson r < 0.61).

3.5. Analysis of the Importance of Features

To interpret the contribution of each feature to SASA prediction, feature importance analysis in the Random Forest model was performed (Table 5). The Mean Decrease in Impurity at tree splitting was calculated.

The largest contributors in both graphs are the degree of the node, the local density, and the average degree of the neighbors.

Key observation: the betweenness centrality was the second most important feature on the HAG (15.8% vs. 11.4% on the Cα graph), which reflects the more complex topology of the heavy atom graph and its greater informative value for characterizing the global role of the residue in the protein structure.

Centralities (betweenness, eigenvector, closeness) add up to:

• ~18% of the contribution on the Cα-graph,
• ~22% on HAG.

Features of a hydrophobic subgraph give ~11% and ~13%, respectively.

3.6. Computational Efficiency

Table 6 compares computation times for different methods in a sample of 509 structures. Measurements were taken on a PC with an Intel Core i7-1165G7 (Python 3.8) processor, averaging over 3 repetitions.

Empirical formulas are calculated in 0.008 ms per structure (~26,000× faster than FreeSASA). Random Forest in prediction mode requires 36.5 ms (~6× faster than FreeSASA). GraphSAGE in prediction mode is faster than Random Forest (2.68 ms), but inferior to it in accuracy.

3.7. Final Comparison of Approaches

Table 7 summarizes the results of a comparison of all the studied approaches.

4. Discussion

4.1. Comparison of Cα Graph and HAG: Accuracy

Transitioning from a Cα-graph to HAG significantly improves accuracy for all methods. For empirical formulas, the correlation increase was +0.10 (from 0.711 to 0.813), and the MAE decreased by 19% (from 0.117 to 0.094). For Random Forest, the correlation increase was +0.025 (from 0.890 to 0.915), and the MAE decreased by 13% (from 0.066 to 0.057). For GraphSAGE, the correlation increase was +0.051 (from 0.841 to 0.892), and the MAE decreased by 20% (from 0.090 to 0.072). This is due to three factors: side-chain representation, a lower threshold (5.0 Å vs 8.0 Å), and a wider degree range.

4.2. Comparison of the Cα-Graph and HAG: Computational Efficiency

GraphSAGE on HAG showed a result close to Random Forest: MAE = 0.0715, Pearson r = 0.8917. This is an important result that demonstrates that graph neural networks can be effective when using HAG. GCN and GAT performed significantly worse (MAE > 0.14, Pearson r < 0.61), which may be attributed to their lower ability to aggregate information in dense graphs. Thus, GraphSAGE on HAG can be considered as an alternative to Random Forest when using deep learning. However, Random Forest remains the preferred choice due to its higher accuracy (MAE = 0.057 vs 0.072) and interpretability.

4.3. Selection of Optimal Parameters of Empirical Formulas

Optimization on 509 structures showed a shift in the parameters for HAG towards higher values: M = 13 (vs. 11 in the previous work), λ = 6.6 (vs. 5.4), a = 1.25, b = 9.50 (vs. a = 1.5, b = 5.5). This reflects the denser packing in the heavy atom graph.

4.4. Practical Recommendations for Choosing a Method

Based on this analysis, we formulate the following practical recommendations (Table 8).

4.5. Trade-Off Between Speed and Accuracy: GraphSAGE as an Alternative to Random Forest

The inference time comparison reveals a clear trade-off. Random Forest on HAG achieves top accuracy (MAE = 0.057, r = 0.915) but requires 36.5 ms per structure (~6× faster than FreeSASA). GraphSAGE on HAG shows only a slight accuracy drop (MAE = 0.072, r = 0.892) while being significantly faster — 2.68 ms per structure, which is 13.6× faster than Random Forest and ~77× faster than FreeSASA. Thus, GraphSAGE on HAG is preferable for large-scale screening where speed is critical, whereas Random Forest on HAG remains the best choice when maximum accuracy is required.

4.6. Restrictions

The main limitation of both approaches is the use of only contact information without explicitly accounting for the types of interactions (hydrogen bonds, ionic pairs, hydrophobic contacts). This can lead to errors for irregularly packed residues or in the presence of ligands/ions.

• For HAG, an additional limitation is the significant increase in graph construction time (21×), which makes it less suitable for large-scale screening tasks where tens of thousands of structures need to be processed.
• For ML approaches, the limitation is the loss of interpretability and the need for pre-training on a representative sample. Although feature importance analysis (Table 5) partially addresses this, the full interpretability of empirical formulas remains unattainable.

5. Conclusion

In this work, three empirical features for assessing amino acid residue surface accessibility based on the degree of the node in the protein contact graph are proposed and validated. A comparative study of two approaches to graph construction was carried out: the classical Cα-graph (threshold 8 Å) and HAG (threshold 5.0 Å). Calibration was performed on an extended sample of 509 structures (128,794 residues).

Main results:

• Empirical formulas on HAG are MAE = 0.094, Pearson r = 0.813 (19% improvement on MAE compared to Cα graph).
• Random Forest on HAG — MAE = 0.057, Pearson r = 0.915 (absolute best result, 13% better than RF on Cα graph).
• GraphSAGE on HAG — MAE = 0.072, Pearson r = 0.892 (close to Random Forest, the best result among GNN).
• GCN and GAT showed significantly worse results (MAE > 0.14, Pearson r < 0.61).
• Computational efficiency - empirical formulas: 0.008 ms/structure (~26,000× faster than FreeSASA), Random Forest: 36.5 ms/structure (~6× faster than FreeSASA).

Table 9 summarizes the practical recommendations.

The proposed approach can be recommended for large-scale pipelines of protein structure analysis.