As a first analysis, we investigated the composition of the dataset used, focusing in particular on the docking poses classified as well-predicted (i.e., structurally similar to the experimentally resolved native conformation of the complex) and on those incorrectly predicted (i.e., with structures that are considerably different from the corresponding experimental native conformation). To this end, we employed the DockQ score (see Methods for a more detailed description), which is able to classify native-like poses and decoy poses according to a threshold value.

In Figure 1a, we report the DockQ distribution for all docking poses generated by the HDOCK method. In particular, for each antibody–antigen docking simulation, we considered the top 10 docking poses ranked by score. The distribution shows that the two main peaks of the probability density function (PDF) are centered at low DockQ values, which are less then 0.24, and at high DockQ values, which are higher than 0.81, indicating that only a small fraction of poses are predicted as native-like (DockQ > 0.81 ), while the majority correspond to decoy poses (DockQ < 0.24 ). A cartoon representation of the structural alignment between the docking pose and the native structure, for different ranges of DockQ values, is shown in Figure 1c. This clearly highlights the difficulty of docking algorithms in accurately predicting the native conformation of interacting proteins. Very high DockQ values (close to 1) typically correspond to very small Root Mean Square Deviation (RMSD) values, which can be interpreted as structural fluctuations within thermal noise of experimentally determined native conformations [36, 37]. Therefore, the ability of the approaches proposed in the following sections must rely on identifying, based on specific interfacial structural properties, the decoy docking poses. In this way, the algorithm can be trained to discriminate between decoy and native-like conformations in a fully general manner, even when the predicted antibody–antigen complex exhibits an interface significantly different from those included in the training set.

In particular, according to the DockQ values calculated for each docking pose, the overall dataset is composed as follows: 19,406 decoy docking poses, 790 intermediate docking poses, and 1,684 native-like docking poses (see Figure 1a; Table 1). As shown by the bimodal trend in the distribution in Figure 1a, in most cases the docking method returns a pose that is structurally distant from the reference structure (i.e., the experimentally resolved complex). However, for 59% of the complexes in the dataset, the top-ranked pose generated by the algorithm is classified as native-like, in some cases with a very high DockQ score, making the docking model and the experimentally determined native structure nearly indistinguishable. This behavior may be due to the algorithm’s prior knowledge of the native structure (or its homologs), as well as particularly easy cases for the algorithm to predict. Nevertheless, this does not hinder the aims of the present work, which first seeks to classify docking poses according to their DockQ value and subsequently to predict the descriptor. In light of this, the development of computational methods capable of identifying decoy docking poses is crucial, as it helps reduce the space of possible binding conformations (by removing these from the candidate solutions) that require further investigation. A deeper insight on the overall dataset is presented in Table 1.

In this study, we defined two different subsets. The first, referred to as the Decoy–Native-like dataset (DNL dataset), includes only decoy and native-like docking poses and is used for the classification approach. The second, which comprises all three classes (decoy, intermediate, and native-like) and is referred to as the DINL dataset, is used to predict the DockQ value of a generic docking pose. The DNL dataset consists of 1,587 decoy docking poses and an equal number of native-like docking poses. Conversely, the DINL dataset, in which the DockQ value of each docking pose is taken into account, is composed of 1,000 decoy poses, 790 intermediate poses, and 1,000 native-like docking poses (see Figure 1a).

Each docking pose is characterized by a docking score, ITScore-PP [38], which is a numerical value used to rank the predicted binding modes of molecules—more negative scores indicate more stable and likely interactions. The distribution of ITScore-PP values is shown in Figure 1b for the native-like and decoy groups separately. The difference between the two distributions is evident, and the classification based on the ITscore-PP descriptor provided by the HDOCK docking method yields an Area Under the Receiver Operating Characteristic (ROC) Curve (AUC) of 0.78 (Equation 24).

The aim of this work is to investigate how a minimal neural network–based approach can improve the classification of docking poses into native-like and decoy categories when an appropriate selection of binding properties is adopted, and how it can directly predict the DockQ value on which this classification is based.

To evaluate the docking poses, an initial set of 21 features was defined and is listed and described in Table 2 (see Equations 4–19 for further details). The selected features primarily describe geometric properties at the interface between antibody–antigen docking poses, some of which are based on graph theory to better capture the complexity of the geometric organization of the residues involved in intermolecular interactions.

An initial Pearson correlation analysis was performed to remove pairs of features showing high correlation (absolute Pearson correlation coefficient > 0.75 , for both positive and negative correlations). The correlation matrices for all pairs of features, both for the initial 21 features and for the 15 features remaining after filtering, are shown in Figure 2a. In particular, in order to remove highly correlated feature pairs while minimizing feature pruning, the absolute Pearson correlation coefficients were mapped onto a graph in which pairs of strongly correlated features were connected. The resulting problem is equivalent to a minimum vertex cover problem, which was solved exactly using integer linear programming (ILP), given the small number of highly correlated features (see Equation 20). In addition, for both matrices, the corresponding graphs are displayed, where each node represents a feature and each edge between two features is weighted (using a red-to-blue color scale) according to their Pearson correlation.

The selection of 15 largely independent features, after appropriate normalization (see Methods), allowed us to perform a Principal Component Analysis (PCA). For each docking pose in the DNL dataset, which is used to classify docking poses as decoy or native-like, a vector of 15 normalized features was associated. The PCA results are shown in Figure 2b, and the proportion of variance explained by each eigenvector is reported in Figure 2c, showing that the first two principal components account for 32% of the total variance. The unsupervised PCA approach was employed here to explore a potential blind classification of the two docking pose classes (decoy and native-like). Each point in Figure 2b represents the projection of the 15-dimensional feature vector (associated with a single docking pose) onto the essential plane defined by the first two principal components (PC1–PC2). Points are colored red and blue according to their membership in the decoy or native-like class, respectively.

The analysis of the two distributions (decoy and native-like) along PC1 does not reveal a clear separation between the two classes, as evidenced by the strong overlap between distributions. This is also confirmed by the ROC curve shown in Figure 2b, with an AUC of 0.52, which is effectively close to random classification. By contrast, the distributions of decoy and native-like poses along PC2 are noticeably more separated, yielding an ROC AUC of 0.68 (see Figure 2e). Therefore, the use of PC2 alone, in a fully unsupervised manner, provides a moderate but non-negligible discriminative power between decoy and native-like classes.

The loading analysis in this context reveals the contribution of each feature to the definition of each principal component. In particular, the interface properties most relevant for the separation between decoy and native-like poses are pca_flatten_ratio, pca_alignment_score, pca_stretch_ratio_bs and pca_flatten_ratio_bs (features 2, 3, 9 and 10, see Table 2), which show a more pronounced difference compared to the corresponding loadings of PC1. The first two features are related to the geometry of the antibody–antigen interaction. Specifically, the first feature reflects the globularity of the complex, which increases when the interface lies in proximity to the CDR, while the second describes the relative orientation of the two molecules. Instead, the last two features capture the circularity and concavity of the binding interface, with higher PC2 values corresponding to a flatter interaction surface. A comparison between the ROC curves of PC1 and PC2 (with ROC AUCs of 0.52 and 0.68, respectively) and that of the HDOCK docking score (ROC AUC of 0.80, which was calculated using the DNL dataset) is reported in Figure 2e, highlighting the need to develop supervised methods to better evaluate each docking pose based on interfacial geometric properties.

In this section, a minimal neural network (see Methods and Equations 21–23) is employed with the aim of improving the classification between native-like and decoy docking poses as provided by the docking score. The goal is to investigate the contribution of neural network–based approaches to enhancing docking pose evaluation. As a first step, the training and test sets were randomly selected. Subsequently, in order to make the procedure as general as possible, multiple training and test datasets were generated so as to be maximally distinct with respect to the features selected for this study.

The first approach was therefore performed by considering one training set and multiple test groups, both drawn from the DNL dataset (see Methods for further details). The predictive performance, in terms of the Area Under the ROC Curve (AUC), was analyzed as a function of the number of network parameters, varying both the total amount of available data and the ratio between training and test samples. The results of this preliminary analysis are shown in Figure 3a, which highlights that, for a number of parameters equal to 70, a plateau in the test AUC curve (in red) is observed in almost all cases considered. By selecting the total number of available complexes—after verifying that this amount is sufficient to capture the information required to discriminate between the two classes—together with a 0.5 ratio between training and test data and a total of 70 network parameters, we obtained an average test set ROC AUC value of 0.90. This value exceeds the ROC AUC calculated using only the docking score. The comparison between the corresponding ROC curves is also reported in Figure 3c. An analogous analysis performed on the discrimination between native-like and intermediate docking poses in the DINL dataset yields a ROC AUC of 0.77, indicating a promising classification performance on this substantially more challenging and subtle task, where the structural differences between the pose classes are markedly smaller than in the native-like vs. decoy scenario.

The importance of the selected descriptors, as indicated by the PC2 loadings, can be assessed by training the NN after removing the descriptors with loadings greater than 0.20 (i.e., those contributing most to the separation between decoy and native-like docking poses according to PC2). This procedure results in a classification performance on the test set that is 14% lower than the performance obtained when retaining all selected features.

In addition, we propose training the NN on training and test sets that are as different as possible in terms of the selected features, in order to make the NN-based classification procedure as generalizable as possible. To this end, the entire dataset was split into two parts (training and test sets, and then swapped) according to the value of PC1, i.e., the projection of the feature vector onto the first principal component of each docking pose. Docking poses with PC1 values below the mean were assigned to one group, while those with higher PC1 values were assigned to the other.

The choice of PC1 as the reference distribution for defining the two groups was motivated by two considerations: (i) PC1, by definition, is the eigenvector associated with the largest proportion of explained variance, thus carrying the highest amount of information; and (ii) the distributions of the PC1 values for the decoy and native-like docking poses show no clear separation (and therefore no intrinsic discriminative power between the two classes), unlike PC2 (see Figure 2). This ensures that, before and after splitting by the PC1 mean, the relative proportion of decoy and native-like poses within each subset remains approximately the same, see Figure 3b.

The results are shown in Figure 3c, which illustrates a neural network discriminative capability between decoy and native-like docking poses that is intermediate between the ITScore-PP docking score provided by HDOCK (ROC AUC of 0.80) and the NN previously trained on randomly selected training and test subsets. A ROC AUC of 0.90 was measured for the NN trained and tested on randomly selected sets, decreasing to values that span between 0.81 and 0.82, when the training and test sets are separated based on the PC1 values associated with each docking pose. In particular, the improved classification capability of the proposed NN-based approach is further confirmed by the steep initial rise of the ROC curves corresponding to the NN-based methods, observed in the early phase (at low true positive rate and false positive rate values). This result highlights the ability of an NN-based approach, when coupled with properly selected features, to improve docking classification performance even when the training and test sets are deliberately constructed to have different underlying properties.

In the previous section, we demonstrated the importance of employing simple NN models for the classification of docking poses into decoy and native-like categories in antibody–antigen complexes. Furthermore, we emphasized the crucial role of accurate feature selection, which, combined with supervised machine learning methods, can significantly improve predictive performance. Here, a minimal feedforward neural network is trained to directly correlate (rather than classify) with the DockQ value, which is one of the standard metrics used to assess the quality of a docking pose. For this purpose, the DINL dataset was taken into account (see Methods for further details).

In this case as well, the scatterplot of the first two principal components obtained from the PCA of the feature vectors of all docking poses is shown in Figure 4a, where each point (docking pose) is colored according to its corresponding DockQ value. Given the inherent difficulty of capturing, through unsupervised approaches such as PCA, the relationship between the interface descriptors of predicted antibody–antigen complexes and their structural deviation from the corresponding experimentally resolved native structures (quantified by DockQ), we developed a neural network (NN) model trained on the DINL dataset.

By randomly selecting the training and test sets (see Methods for further details), we statistically analyzed how the Pearson correlation (Equation 25) in the test set between the experimental DockQ and the predicted DockQ (pDockQ) varies as a function of the number of NN parameters. The results for a NN trained on 80% of the available poses, reported in Figure 4a, show that a substantial performance gain is achieved by increasing the number of parameters up to approximately 70. Beyond this point, the correlation between DockQ and pDockQ increases much more slowly, while the mean square error (MSE) on the training set reaches a plateau for ≈ 86,400 parameters (Figure 4a).

Furthermore, the variation in correlation between pDockQ (computed with an 86,400 parameters NN) and DockQ has been studied as a function of the training set size, spanning from 50% of the DINL dataset to 93%, alongside the difference in MSE ( Δ MSE) between training and test set. For both the measures, the results, reported in Figure 4b, show an optimal average value for the 80% training set proportion.

For the NN trained with 86,400 parameters, on a set composed by 80% of the docking poses, the correlation between DockQ and pDockQ in the DINL set is 0.59.

The comparison between neural network predictions (pDockQ) and the docking score (ITScore-PP) was performed by evaluating the correlation with the DockQ score. In addition, we used projections onto the first two principal components (PC1 and PC2) of each docking score as potential predictors of DockQ. The scatter plots showing the relationship between DockQ and each proposed predictor (supervised and unsupervised) are reported in Figure 4c.

In particular, the correlations between DockQ and PC1, PC2, ITScore-PP (the docking score), and pDockQ are −0.04, −0.27, −0.41, 0.59 respectively (Table 3). This indicates that the NN-based approach, which uses as input the 15 selected features, substantially improves the quantitative evaluation of docking poses compared to the original docking score. Of particular note is the correlation between DockQ and the second principal component (PC2) of the PCA performed on the features. As a fully unsupervised descriptor, PC2 provides insight into the features that contribute most to the definition of the component (loadings), thereby offering the opportunity to further refine NN-based models through preliminary feature selection procedures.

To better illustrate the ability of the neural network–based predictive method to estimate DockQ values even in intermediate cases (0.24 < DockQ < 0.81), a probability density function (PDF) was computed for each DockQ range (analogous to the boxplot analysis shown in Figure 4c). The distributions of the pDockQ descriptor are progressively shifted with increasing DockQ ranges (see Figure 4d), thereby confirming the method’s ability not only to classify docking poses as native-like or decoy—as also supported by this DockQ estimation procedure—but also to correlate with intermediate DockQ values, with slightly lower yet satisfactory accuracy compared to the docking score.

Furthermore, to assess the overall quality of the pDockQ descriptor, it has been benchmarked in terms of Pearson correlation coefficient against both ITScore-PP and the predicted binding free energy ( Δ G ) , obtained using an MM/GBSA-based predictor [39] via the HawkDock server [40]. This comparison has been performed on a randomly selected small subset of the DINL dataset composed by 84 docking poses (30 decoy, 24 intermediate and 30 native-like poses, in order to maintain the proportions of the DINL dataset). For this analysis, the pDockQ values have been computed by a 86,400 parameters NN trained on all the DINL docking poses that do not share the reference native complex with any test set element. While Δ G and pDockQ show a comparable performance (respectively −0.62 and 0.67), both have a significantly larger Pearson correlation coefficient in magnitude than ITScore-PP (−0.43). Although the small size of the test set of this assessment does not allow a definitive statement, the pDockQ approach results are promising, even when compared with one of the state-of-the-art methods reported in the literature.

The initial dataset consisted of 9,780 experimentally resolved antibody–antigen complexes retrieved from the SAbDab database [41]. A first filtering step was applied to retain only complexes in which the antigen was classified as a protein or peptide and consisted of a single chain (thus preserving only monomeric antigens), resulting in 9,486 structures. Structures containing missing residues were either repaired or removed from the dataset, yielding 9,463 structures.

A multiple sequence alignment among all antibody–antigen complexes in the dataset was performed to remove redundancy. Specifically, for each complex we considered a single sequence obtained by concatenating the antibody sequence (heavy and light chains) with the sequence of the corresponding antigen. These sequences were then processed with CD-HIT [42–44] using a sequence identity cutoff of 0.9, resulting in 2,517 centroids, which represent the most representative sequences in the entire dataset. Since the study focuses on the calculation of interface properties, it was crucial to ensure that the interfaces were complete, i.e., without missing residues in the binding region. Therefore, complexes with incomplete interfaces were excluded, reducing the dataset from 2,517 to 2,244 structures.

Finally, energy minimization was performed on all structures, resulting in a final set of 2,188 properly minimized complexes.

Each antibody–antigen complex with a known experimental structure was split into two separate structures, antibody and antigen, which were then subjected to molecular docking simulations using the HDOCK method (thus considering the interacting structures in their bound conformations). For each antibody–antigen docking simulation, the top ten poses proposed by the method were retained. Each docking pose was evaluated using the DockQ metric, which is defined according to the following formula: D o c k Q = F nat + l R M S scaled + i R M S scaled 3 (1) with l R M S scaled = 1 1 + l R M S 8.5 A ̊ 2 (2) and i R M S scaled = 1 1 + i R M S 1.5 A ̊ 2 , (3) where F nat , l R M S and i R M S are the CAPRI-standard classification metrics [30, 31].

In particular, for the DNL dataset, we selected for each experimental complex the “decoy” docking pose as the one associated with the lowest DockQ value among the ten poses considered (ensuring in all cases that DockQ < 0.24 ), and identified the “native-like” pose as the one with the highest DockQ value among the ten poses generated by HDOCK (with DockQ > 0.81 ). Furthermore, 1,000 docking poses classified as decoy (based on their DockQ values) and 1,000 docking poses classified as native-like (also based on DockQ) were randomly selected from the docking poses obtained after the previous filtering steps. These poses were used to define the DINL dataset, which served as the training set for the neural network designed to predict DockQ values.

The features used throughout the whole paper can be divided into three groups: complex geometry features, interface features, complex graph features. The first group comprises all the measures related to the geometrical arrangement of the antibody-antigen complex α -carbon atoms. The first group is composed by

pca_stretch_ratio:

The second group is composed by features accounting for several properties of the antibody-antigen binding site (BS). For this group, we defined as BS residues those residues whose α -carbon is within 12 A ̊ to an α -carbon atom from a different molecule (i.e. the antibody residues closer than 12 A ̊ to an antigen residue and vice versa). The second group features are

bs_sasa_ratio:

The SASA values are measured using the Shrake-Rupley “rolling ball” algorithm with probe radius of 1.40 A ̊ and definition of 100 points/ A ̊ 2 via the Biophyton library [46].

The third group features are common graph theory descriptors measured on two networks: an unweighted network, where all nodes corresponding to residues whose α -carbons are within 12 A ̊ are linked, and a weighted network, where to any edge ( i , j ) of the unweighted network is assigned a weight W i j = 1 / d i s t ( i , j ) . The following features belong to the third group:

edge_density: edge density of the unweighted network

The network-related features are measured via the NetworkX library [47].

The set of 21 features has been reduced in order to avoid redundancy due to the presence of highly correlated features. In this instance, the least amount of features such that any remaining couple has absolute Pearson correlation < 0.75 was removed. Linking the highly correlated features in an undirected unweighted network, this problem results to be equivalent to a minimum vertex cover problem (pruning the least amount of nodes such that each remaining node is isolated), therefore exactly solvable via integer linear programming (ILP). The corresponding ILP formulation, with E the set of edges, Given:  x i = 1 , if node  i  is removed 0 , otherwise Minimize:  ∑ i = 1 N x i Subject to:  x i + x j ≥ 1 , ∀ i , j ∈ E (20) was solved via the PuLP modeler [48]. The remaining features were 15: pca_stretch_ratio, pca_flatten_ratio, pca_alignment_score, bs_sasa_ratio, ab_bs_size, ag_bs_size, pca_stretch_ratio_bs, pca_flatten_ratio_bs, bs_mean_hydrophobicity, mean_strength, network_radius, mean_assortativity, weighted_mean_clustering, network_transitivity.

For processing, each feature was normalized via the scikit-learn library [49], such that all the features share the same weight. The PCA was performed via scikit-learn, as well.

Every NN in this work has been defined via Tensorflow [50] and has the same structure: two-hidden-layers feed forward NN. Each hidden layer has reLu activaction function, furthermore, the output layer of the NNs used in Section 2.5 are provided with a sigmoid activation function, in order to retrieve p D o c k Q ∈ [ 0,1 ] . While varying the number of parameters the proportion of nodes in the first and second hidden layer is kept fixed at 1 : 2 . Therefore, naming M the number of first layer nodes, one can retrieve the number of parameters N : N = F + 1 M + 1 2 M 2 , (21) where F = 15 is the number of input features. The NN weights are fitted via AdamW algorithm with learning rate 0.001, through 300 epochs for the classifiers (Section 2.4) and 400 epochs for the predictor NNs (Section 2.5). In Section 2.4 binary cross-entropy was used as loss function: C p pred ‖ p true = − ∑ i ∈ 0,1 p pred i l n p true i (22) where { 0,1 } is the set of the possible classifications, i.e. “Decoy” or “Native-like”. In Section 2.5 mean square error (MSE) was the loss function, instead: MSE = 1 N ∑ i = 1 N p D o c k Q i − D o c k Q i 2 , (23) where N is the number of docking predictions in the dataset and p D o c k Q i and D o c k Q i are the values of pDockQ and DockQ associated to the i -th prediction. In Section 2.5, in order to obtain the pDockQ values for the whole DINL dataset, it has been split into several complementary subsets, according to the proportion of the training set. The pDockQ values of each subset have been computed using the others as training set.

The area under the receiver operating characteristic curve (ROC AUC) was used to assess the quality of the classifications throughout Section 2.4. Given two classes (Positive and Negative) and the distribution of a measure for each of the classes, the ROC curve is the parametric curve ROC     curve = ( f p r ( t ) , t p r ( t ) ) representing the variation of the false positive rate f p r ( t ) and the true positive rate t p r ( t ) in function of the measure threshold t used to split the classes, where t p r t = # true positives # positives , f p r t = # false positives # negatives . (24)

The ROC AUC equates the probability that, given a random negative element and a random positive element, the negative element correspond to a measure larger than the positive. In this instance the Negative and the Positive classes were “Decoy” and “Native-like”. The ROC curves and the ROC AUCs were computed via the scikit-learn library [49].

Regarding the regression tasks (Section 2.5), the assessment was done via Pearson correlation coefficient ( ρ ) between any measure x and the D o c k Q score of the docking prediction: ρ x = 〈 x − ⟨ x ⟩ D o c k Q − ⟨ D o c k Q ⟩ σ x σ D o c k Q . (25) The validity of the Pearson correlation was assessed performing a p-value test of the null hypothesis that the distributions underlying the samples are uncorrelated and normally distributed. Both the Pearson correlation coefficient and the p-value were computed via the SciPy library [51].