We selected 326 antibodies with redundancy lower than 90% and resolution <3.0 Å using the SabDab database (Dunbar et al., 2013b). 229 antibodies were solved in complex with protein antigens, 71 with haptens, 19 with carbohydrates and 7 with nucleic acids. The sequence of each antibody was renumbered according to the Chothia numbering scheme (Chothia and Lesk, 1987; Chothia et al., 1989) using an in-house python script.

For each antibody and protein antigen 3D structure, atomic partial charges and radii were assigned using PDB2PQR with default parameters (Dolinsky et al., 2004). Solvent Accessible Surface (SAS) was computed using GROMACS (Abraham et al., 2015). Electrostatic surface (ES) potential was computed using the Bluues software (options -srf and -srfpot) (Fogolari et al., 2012). Each molecular surface point was assigned to the electrostatic potential of the corresponding residue. The “geometry” (Habel et al., 2019) and “Bio3D” (Grant et al., 2006) packages available in R were used for PDB structure processing and analysis.

The set of selected molecular surface points was scaled to the unit sphere and placed into a 3D grid of dimension 128³. To avoid boundary effects, the size of the bounding box of the point cloud was set so as to be contained within 80% of the unit sphere (Grandison et al., 2009). Voxelization was performed separately for SAS and ES. In SAS voxelization, each voxel was assigned a value of 1 if the center of the voxel was closer than 1.7 to any SAS point, 0 otherwise. In ES voxelization, each voxel was assigned the mean ES value of the enclosed points, 0 otherwise.

Since the Zernike formalism does not differentiate positive and negative values (Chikhi et al., 2010; Daberdaku and Ferrari, 2018), but only patterns of non-zero values in the 3D space, voxels were initialized for positive and negative patterns separately using a similar approach as done in (Chikhi et al., 2010), as follows: f e l e c + = 0 i f f e l e c < 0 f e l e c + = f e l e c i f f e l e c > 0 (1) f e l e c − = f e l e c i f f e l e c < 0 f e l e c − = 0 i f f e l e c > 0 (2)

In summary, voxels with positive electrostatics values were initialized to 1 and all other voxels with negative electrostatics values were set to zero, and vice versa. The resulting voxels, one for SAS values, and two for positive and negative ES values, respectively, were considered as three different 3D functions, f(x), each expanded into the 3DZD as described in the next section.

For the quantitative description of the binding sites, we rely on a representation based on the Zernike polynomials and their corresponding moments. Moment-based representations are a class of mathematical descriptors of shape, originally developed for pattern recognition and subsequently generalized to three-dimensions (Ming-Kuei Hu, 1962; Canterakis, 1999; Novotni and Klein, 2004).

A surface described by a function f (r, θ, ϕ) in polar coordinates can be represented by a series expansion in an orthonormal sequence of polynomials (Canterakis, 1999): f ( r , θ , ϕ ) = ∑ n = 0 ∞ ∑ l = 0 n ∑ m = − l l C n l m Z n l m ( r , θ , ϕ ) (3) where the indices n, m and l are the order, degree and repetition, respectively.

The Zernike polynomials can be written as: Z n l m ( r , θ , ϕ ) = R n l ( r ) Y l m ( θ , ϕ ) (4) where the Y functions are complex spherical harmonics depending on both θ and ϕ while R only depends on the radius r, which is given by R n l ( r ) = ∑ k = 0 ( n − l ) 2 N n l k r n − 2 k (5) where N is a normalization factor.

The 3D Zernike moments of a surface described by a function f (r, θ, ϕ) are defined as the coefficients of the expansion of f(r) in the Zernike polynomial basis, i.e.: C n l m = ∫ ∣ r ∣ ≤ 1 f ( r ) Z n l m ( r , θ , ϕ ) ̄ d r (6) where Z ̄ is the polynomial complex conjugate.

Their rotation invariant norms, i.e. the 3DZD, are defined as: D n l = ‖ C n l m ‖ = ∑ m = − l l ( C n l m ) 2 . (7)

The Zernike formalism can be as detailed as desired by modulating the order of the expansion n. In our implementation, the function f represents the geometric or the (positive or negative) electrostatic potential of the molecular surface, and the maximum order of expansion was set to 20, giving a total of 121 invariants.

Given the dataset of Antibody-Antigen complexes containing protein antigens, the native geometric epitope was defined as the set of residues of the antigen having a distance lower than 6 Å to any residue of the antibody. The pivot residue was defined as the residue with the lowest mean distance to any residue of the native geometric epitope. The native electrostatic epitope was defined as the set of residues of the antigen having a distance shorter than 15 Å to any residue of the antibody. For the set of native geometric epitope residues, the Solvent Accessible Surface Area (SASA) was computed using GROMACS. The mean and standard deviation values of the computed global and residue-based SASA were used to generate an alternative set of surface patches, i.e. decoy epitopes. The algorithm first selects a decoy pivot residue, i.e. by randomly selecting any solvent exposed residue having a value of SASA within half standard deviation of the mean SASA value measured over all pivot residues of the native epitopes, i.e. SASA = 0.48 ± 0.33 nm² (Supplementary Figure S1). The algorithm proceeds by adding neighboring solvent accessible residues, i.e. having relative SASA >0.2 (Tien et al., 2013), until the decoy geometric epitope reaches a similar global SASA to that of the native epitope. To ensures continuous coverage of the antigen protein surface (Supplementary Figure S2) and diversity of the generated patches, a maximum 50% surface patch overlap was allowed between native and decoy epitopes. Electrostatic decoy epitopes were defined by calculating the electrostatic potential over the region defined by a geometric decoy epitope considering all the charged residues within 15 A ̊ to the pivot residue.

Given a pair of ordered set of 3DZD, x and y, their cosine distance is measured as: D ( x , y ) = 1 − S c ( x , y ) = 1 − x y ‖ x ‖ ‖ y ‖ (8) where S _c (x, y) is the cosine similarity as measured by the “proxy” R package (Meyer and Buchta, 2019).

Given two patches A and B, the similarity between their 3DZD is computed as: [ A − B ] s h a p e = D ( X s h a p e A , X s h a p e B ) (9) [ A − B ] e l e c = ( D ( X e l e c + , A , X e l e c + , B ) + D ( X e l e c − , A , X e l e c − , B ) ) 2 (10) where X _shape , X e l e c + and X e l e c − are, respectively, the shape, the electrostatic positive potential, and the electrostatic negative potential 3DZD.

The surface complementarity between A and B is defined as follows: [ A − B ] s h a p e = D ( X s h a p e A , X s h a p e B ) (11) [ A − B ] e l e c = ( D ( X e l e c + , A , X e l e c − , B ) + D ( X e l e c − , A , X e l e c + , B ) ) 2 (12)

We have previously shown that a 3DZD-based description of the surface of the antibody CDRs provides an effective metric for antibody classification according to their specificity towards protein and non-protein antigens (Di Rienzo et al., 2017). Here we extend this approach to the analysis of both the shape and electrostatic properties of the CDRs and analyze the classification performance of both descriptors at different orders of the Zernike expansion. For each CDR we generated two sets of 121-dimensional vectors, representing the 3DZD of the shape and the electrostatic surface, similar to what done in (Chikhi et al., 2010; Di Rienzo et al., 2020). The similarity between each set of descriptors is then computed to perform an all-against-all comparison of CDRs, according to Eq. 9, 10 in Methods section. For each CDR, we then selected the nearest neighbors set as the 5% most similar CDRs in terms of shape and electrostatic surface and analyse the number of protein binding antibodies (N _pb ) in the neighbours set. As it is shown in Figures 2A,B, protein-binding antibodies (green curve) are typically characterized by an higher number of N _pb ( m e a n ( N p b s h a p e ) = 13.37 ± 2.61 , m e a n ( N p b e l e c ) = 13.54 ± 3.24 ) in the neighbors set as compared to non protein-binding antibodies (orange curve) ( m e a n ( N p b s h a p e ) = 10.31 ± 2.99 , m e a n ( N p b e l e c ) = 9.93 ± 3.13 ) and to random expectation (i.e., Ex [N _pb ] = N _Prot /N _tot , where Ex [N _pb ] is the expected number of protein-binding antibodies if they were distributed uniformly, N _prot represents the number of protein-binding in the dataset and N _tot is the total number of antibodies in the dataset.).

We next analyzed the performance of each descriptor in classifying the CDRs as a function of the antigen type, using a leave-one-out approach. In summary, for each CDR, if the N _pb was greater than Ex (N _pb ) the CDR was labeled as protein-binding, non protein-binding otherwise. The obtained classification accuracy for the shape and electrostatic descriptors at order n = 20 is 75 and 73%, respectively. Using a Receiver Operating Curve (ROC) analysis, both descriptors achieved an Area Under the Curve (AUC) of 0.78. We next analyzed the classification performance when assigning the class label based on the weighted contribution of shape and electrostatics, as follows: N p b ̄ = A N p b e l e c + ( 1 − A ) N p b s h a p e A ∈ [ 0,1 ] (13) where N p b s h a p e and N p b e l e c correspond to the N _pb computed based on shape and electrostatic descriptors, respectively, and A is the weight ranging from 0 to 1. The results are shown in Figure 2C, where the ROC AUC is reported as a function of the weight A and the order n of the Zernike expansion. As it can be noticed, overall performance increases with increasing values of n. Higher AUC values are achieved when both descriptors contribute with similar weight in the classification. Top classification performance indeed is obtained with A = 0.4 and n = 17, leading to an AUC = 0.85 and accuracy of 81%. A very similar performance is obtained with n = 20 and A = 0.4 (AUC = 0.83).

The sequence and structure analysis of antibodies, as well as antibody engineering experiments, crucially rely on the precise identification of the CDRs from the antibody sequence (Chothia and Lesk, 1987; Chothia et al., 1989; Kabat et al., 1992; MacCallum et al., 1996; Lefranc, 2011). On the other hand, it is well known that the CDRs only provide a proxy of the actual antigen-binding site, i.e. the antibody paratope (Kunik et al., 2012; Olimpieri et al., 2013). Indeed, early studies showed that only 20–30% of residues within the CDRs are directly involved in the interaction with the antigen (Padlan, 1994; Sela-Culang et al., 2013). To quantify to what extent this approximation affects our predictions, we analyzed the classification performance as a function of distance from the center of the antibody-antigen interface. For each Antibody-Antigen complex, we defined a centerpoint, b, as the centroid of the 10 interface atoms of the antibody closer to the antigen and computed the 3DZD for increasing concentric shells around b.

We then applied the same classification procedure as described previously, by fixing the order n = 20 for both shape and electrostatic 3DZD. The results are shown in Figure 3 where the ROC AUC of the individual classifiers are reported as a function of the percentage of the CDR surface included in the analysis.

As it can be noticed in Figure 3B, the performance of the shape-based classifier shows a maximum when the selected surface region around b extends up to including 20% of the CDRs (ROC AUC = 0.88) followed by a linear decrease when larger surfaces are considered. These results are consistent with the previous notion that shape recognition of the antigen is largely mediated by smaller interacting surfaces contained within the CDR, i.e. the antibody paratope. In summary, while the overall CDR surface can inform about the function of the antibody, this analysis highlights that the information of the paratope can significantly increase our ability to predict antibody specificity. On the other hand, in Figure 3C, the classification performance based on the electrostatic descriptor shows a different trend. Indeed, while the classifier shows an overall lower performance compared to the shape-based classifier, performance increases when larger CDR surfaces are considered, reaching a maximum when almost the entire CDR surface is included in the analysis.

A key feature of the 3DZD description is that it is invariant under rotation and translation of the represented surface. This implies that two interacting protein regions with perfect surface complementarity yield identical sets of 3DZD descriptors (Venkatraman et al., 2009a). In line with this principle, here we focus on the application of 3DZD to the analysis of surface complementarity between antibody CDRs and their cognate protein antigens (Details in Methods). The results are shown in Figure 4, where the average surface shape and electrostatic complementarity computed on 229 antibody-antigen complexes are reported as a function of the interaction cutoff distance between the antibody and the antigen, and the order n of the series expansion. As expected, shape complementarity decreases at higher values of the cutoff distance, i.e. as regions of the antibody/antigen that are distant from the interaction interface are progressively included in the analysis. On the other hand, electrostatic complementarity increases at higher distances, reaching a maximum when the distance cutoff is 15 A ̊ . Notably, in both cases, results are consistent at different orders n of the series expansion. These results indicate that the two descriptors are competent in capturing both short- and long-range effects occurring during antibody-antigen recognition. As further validation of our approach, we measured the surface complementarity at the paratope-epitope interface and compared it with that measured between the paratope and a set of non-epitope, solvent-exposed regions of the antigen, i.e. surface decoys. The results are reported in Figure 5, where both shape and electrostatic complementarity are reported for each paratope as normalized Z-score distances to native epitopes and surface decoys, respectively. Notably, while in terms of amino acid composition, surface size, and solvent accessibility the antibody epitopes are essentially not distinguishable from the decoys (Supplementary Figure S3), they display significantly higher surface shape and electrostatics complementarity to the paratope. In summary, the metric is able to distinguish the correct paratope-epitope pair among the set of decoys with a classification performance of AUC = 0.75 based on the shape descriptor, and AUC = 0.61 based on the electrostatic 3DZD. Additionally, we compared the 3DZD complementarity observed between specific paratope-epitope pairs and that between the antibody paratopes and non-native epitopes. The results (Supplementary Figure S4) show that only a relatively low number, i.e. 68% (72%) of the antibodies in our dataset show a higher shape (electrostatic) complementarity to their cognate epitope compared to non-native epitopes, highlighting the limitation of this metric in the very elusive task of predicting which antibody recognises specifically a given antigen.