In the present section, we discuss the results we obtained applying our recently developed computational protocol (Milanetti et al., 2021a) on a structural dataset composed of 229 antibody–antigen complexes (see Methods).

In particular, we firstly identified for each complex the paratope (epitope) as the set of residues with at least one atom closer than 4 A ̊ to an antigen (antibody) atom. Therefore, after separately computing the molecular surface (Richards, 1977) for both the proteins in interaction, we extracted the portions of the molecular surfaces relating to the binding site residues to properly characterize the shape of the interacting regions of antibodies and antigens (see Figure 1A).

Once we identified the interacting regions, we characterized them through the 2D Zernike polynomials, summarizing their geometrical properties in an ordered set of numerical descriptors (see Materials and Methods). By definition, two perfectly fitting surfaces are characterized by the same shape, meaning that, in principle, the difference between their Zernike descriptors is zero. Therefore, the shape complementarity between two molecular regions is compactly evaluated by such formalism. The lower is the distance between the Zernike descriptors and the higher is the shape complementarity between the corresponding protein regions (Milanetti et al., 2021a; Di Rienzo et al., 2021b).

We described all the paratopes and epitopes in the dataset with the Zernike formalism. In summary, we deal with 229 (the number of structures in our dataset) sets of 121 (the number of invariant descriptors when the order of expansion is set to 20) numerical descriptors for the paratopes and 229 for the epitopes. We thus defined the specific complementarities as the Euclidean distance between the descriptors of the 229 pairs of interacting (extracted from the same structure) paratope and epitope. Diversely, Non-Specific complementarity is the Euclidean distance between all the pairs of unrelated paratopes and epitopes (i.e. a paratope extracted from one structure and an epitope extracted from another one). In the end, we deal with 229 specific complementarities (one for each complex) and a high number (25,000) of non-specific complementarities (all the possible paratope–epitope pairs given a dataset of 229 items). In other words, in dealing with N antibody–antigen complexes, the specific complementarity, C _s , is defined as D p i , e j | i = j = ∑ k = 1 121 p i k − e j k 2 | i = j , (1) where D is the Euclidean distance between the vectors of the paratope (p _i ) and epitope (e _i ) Zernike descriptors. Since we expanded all the paratopes and epitopes to the order 20, we dealt with 121 descriptors for each binding region. On the other hand, the non-specific complementarities, C _ns , can be computed as D p i , e j | i ≠ j = ∑ k = 1 121 p i k − e j k 2 | i ≠ j . (2)

In Figure 1B, we reported a boxplot highlighting the differences between C _s and C _ns distributions. As expected, the distribution of C _s is statistically lower than the distribution of C _ns (Kolmogorov–Smirnov test p value < 2.2 e − 16 ), testifying the sensitivity of Zernike in recognizing regions actually in contact from the non-interacting ones. In the second step, we normalized the Zernike complementarities with the Z-score. In particular, for each paratope we have 1 C _s and 228 C _ns complementarity values, each of which is related to a specific patch pair. Normalizing over this set of numbers and looking at the Z-score regarding the specific one, we assessed the propensity that characterizes each antibody toward its specific antigen. In Figure 1C, we report the distribution of such specific Z-scores. As evident, they are mostly negatives (86% of the cases have a Z-score lower than 0, and 28% have a Z-score lower than -1), providing evidence that specific interactions are characterized by lower distances (higher complementarity) than non-specific ones.

Taken together, these results confirm the ability of our method to correctly capture the main determinant of molecular shape complementarity.

The Zernike formalism enjoys several advantageous features in representing the molecular surface: mainly, the invariance under rotation that makes such descriptors an “absolute” characterization of local protein morphology and the low computational cost of its calculation. Indeed, in this section, we present our algorithm that exploiting these advantages aims to optimize the shape complementarity of an antibody toward a given molecular target region. A similar procedure has already been presented and tested in our previous work (Di Rienzo et al., 2020b), and here for the first time, it is applied to antibody–antigen interaction systems.

Figure 2A illustrates the main steps of the algorithm. We defined the target region as the portion of the antigen surface toward which an antibody will be optimized. It is thus necessary to identify the antibody chosen as a starting point of the algorithm. Summarized with the Zernike descriptors of the target region we compute the shape complementarity with all the paratopes of our dataset; here, the template, i.e., the antibody selected as the starting point for the mutagenesis process, can be chosen among the paratopes characterized by a high initial complementarity.

After establishing the template, we perform a Monte Carlo simulation employing computational mutagenesis on the paratope residues. In each step, we randomly selected a residue mutating it into another one of the 19 possible ones. The mutation generates a different paratope, characterized by a different shape of its molecular surface. Consequently, re-computing the Zernike descriptors, we can evaluate the effect of the mutation on the complementarity with the target region; indeed we can define the complementarity balance as Δ C = C m u t − C w t = D p m u t , e t a r − D p o r i g , e t a r , (3) where p _orig and p _mut are the Zernike descriptors of the original and the mutated paratopes, respectively, while the e _tar represent the Zernike descriptors of the target epitope and D represent the distance between 2 sets of descriptors. Since, as said, a high complementarity is reached when D is low, ΔC < 0 means a higher complementary surface, and ΔC > 0 is obtained when the mutation is deleterious since it causes a worsening of the shape complementarity.

The number of combinations of possible mutations in an interacting region, composed usually of tens of residues, is huge. Therefore, to effectively sample the space of the possible mutants, we perform a Monte Carlo Metropolis simulation, iterating the procedure described previously, where the mutation in each step is accepted according to the following rules: P = 1 i f Δ C < 0 e − β Δ C i f Δ C ≥ 0 , (4) where β ∼ 1 T is the temperature factor that determines the probability of acceptance of a step worsening complementarity. Note that this aspect is crucial to properly exploring a large number of different mutants: β is thus progressively increased during the Monte Carlo simulation, progressively confining the system to minimum energy in a simulated annealing process (Kirkpatrick et al., 1983).

To observe how many steps are necessary to reach the equilibrium for each β, we preliminary ran several fixed-temperature Monte Carlo simulations. We selected the epitope of an antigen structure in our dataset (PDB id: 1AR1) and, excluding its one, we choose as the starting template, the most complementary paratope in the structural dataset. We thus performed six different Monte Carlo simulations, each for a different β, where the acceptation probability in each mutagenesis step is given by Eq. 4. Performing 10 independent simulations of M = 2,000 steps for each β, the averaged results we obtained are summarized in Figure 2B. In the top panel, we reported the energy (i.e., the complementarity) as a function of the Monte Carlo steps. In the low panel, we reported the standard deviation of the energy of the remaining part of the Monte Carlo simulations as a function of the number of steps (i.e., σ(E) for m = 1,000 means standard deviations of the energy obtained in steps 1,001–2,000). As expected, for low values of β (i.e., high temperature) the system lives in a condition of indifferent equilibrium, where whatever mutation has an equal likelihood to be accepted, independently from its effect on complementarity. When, on the contrary, β is high (low temperature), the energy of the system rapidly decreases to a stationary local minimum. This trend is confirmed by looking at the stationary value of energies or, equivalently, noting that standard deviation tends progressively to zero. In the light of these results, in our protocol, we set N = 700 for each temperature. In this way, we preliminary allow the system to move away from the starting local conformation, thus freezing it in a new energy minimum, characterized by an increased shape complementarity with the target region.

The approach described here is general and can be applied to whatever protein. This notwithstanding that we applied it to the SARS-CoV-2 spike protein is a very relevant case of macromolecular interaction. Indeed, the severe acute respiratory syndrome Corona virus 2 infection is still a very serious danger for public health (Huang et al., 2020; Zhu et al., 2020).

Many therapeutic strategies are devoted to the SARS-CoV-2 spike protein, protruding from the viral envelope and responsible for the cell entry mechanism (Walls et al., 2020; Wan et al., 2020; and Zhou et al., 2020). Thus, we obtained, using the dedicated section of Cov-AbDab (Raybould et al., 2021), a structural dataset of 145 spike–antibody complexes (we will call it the “Spike dataset”). We thus characterized the paratopes (antibody binding residues region) and the epitopes (various regions on the Spike) with the Zernike formalism. This allows us to compute, also for the Spike dataset, the specific complementarities and the non-specific complementarities, defined in Eqs 1, 2. The result of this analysis is shown in the following Figure 3A. These complementarities are reported as boxplots (light green for specific complementarities, light red for non-specific ones), even showing the boxplots regarding the general dataset (dark green and dark red, already shown in Figure 1B). Its results are evident that, also in the Spike dataset, specific interactions are characterized by a complementarity much higher than non-specific interactions (k.s. test p-value < 2.2 e-16). As expected, the non-specific interactions are represented by very similar distributions, since in both datasets they would represent an ensemble of non-interactions. Looking at these results, we noticed that the population of spike–antibody complexes showed the same behavior as the general population of the protein–antibody complexes.

In this framework, we selected on spike molecular surfaces, three different regions as targets for the optimization protocol. On one hand, we targeted two different molecular regions involved in the interaction between spike and angiotensin-converting enzyme 2 (ACE-2), the well-known cellular receptor responsible for viral cell invasion. Moreover, we also optimized an antibody toward a very exposed region in the N-terminal spike domain, responsible for contacting sialic acid molecules. Indeed, such an interaction can confer to the virus, as occurred for the Middle East respiratory syndrome coronavirus (MERS-CoV) (Li et al., 2017), an additional molecular mechanism for cell intrusion. The responsible spike region represents a promising therapeutic target (Baker et al., 2020; Milanetti et al., 2021b).

We selected the residues constituting such epitopes and we characterized their molecular surfaces through Zernike formalism. Thus, we calculated the complementarity between these regions and all the antibody binding sites in our original dataset. To begin the optimization from a favorite starting point, we selected as templates antibody-binding sites characterized by the highest complementarity with each identified target.

We applied the procedure described in the previous section obtaining optimized paratopes whose molecular images are shown in Figure 3B, where we reported both the optimized antibodies and antigen interacting surfaces. In Figure 3C, we reported, for each of the Monte Carlo simulations performed, the shape complementarity as a function of the steps of the simulation, where the dashed lines enclose ranges with different β values. Each simulation significantly optimizes shape complementarity, obtaining a Zernike distance decrease of 43% on average. Significantly, all the designed binding sites are characterized by a very high final shape complementarity, in terms of Zernike descriptors. Indeed, it is worth noting that the values obtained by all the three designed binding sites are lower than all the specific complementarities obtained in our structural dataset.

Moreover, in Figure 3D, we reported the probability density function of specific interactions in the Spike dataset. The vertical dashed lines represent the final complementarity values we get after the optimization procedures. It has to be noted that our protocol can effectively optimize the shape complementarity, obtaining final shape complementarity similar to the best cases observed in the Spike dataset.

The computational protocol we developed does not take into account several properties, known to be important in molecular recognition, such as electrostatics or hydrophobicity. In particular, our working hypothesis is that the shape complementarity plays a primary role as a perfect match between molecular surfaces due to an optimal structural rearrangement, which is probably caused by the compatibility of amino acid compositions of the interacting patches. However, the relationship between shape complementarity and chemical–physical properties is not always trivial, requiring a further test for the patches proposed as interacting, to also analyze the compatibility of a chemical nature. This means that a residue substitution can, in principle, worsen the chemical compatibility between molecules, even if the shape complementarity is enhanced. For this reason, as a further step of our optimization protocol, we performed a molecular docking analysis using HDOCK (Yan et al., 2020). More specifically, we docked the spike protein and the antibodies, both in the original and in the optimized versions, to study the effects our computational protocol has produced. We constrained docking to interact with the residues composing the spike target epitopes and the antibodies’ optimized regions. We summarized in Figure 4 the results we obtained.

Thus, we selected the 10 best docking poses regarding both the original and the optimized antibodies. To assess whether the optimization protocol has been effective, some estimators of binding compatibility have been calculated. In particular, the number of residue–residue inter-molecular contacts, the surface buried in the complex, the mean distance of the closest atoms between the two interfaces, the Coulomb inter-molecular energy, the Lennard–Jones inter-molecular energy, and the HDock binding score. For each of these observables, we computed the relative percentage of gaining after optimization so that the positive values indicate an increased binding tightness (See Materials and Methods). As shown in Figure 4A, even if in two applications we note a worsening, in one case, the optimization procedure has produced an antibody with better values of all such estimators, indicating the importance of including the molecular docking approach as a filter of selected patches based on geometric compatibility.

We focused therefore on this case and we analyzed how the introduced residue substitutions were responsible for this better compatibility. Analyzing the best docking poses, in Figure 4B, we reported the gaining in terms of the number of inter-molecular contacts and inter-molecular energy each residue registered before and after the optimization. Each residue is represented by a blue bar, while the residues mutated in the protocol are depicted in orange. As evident from the main effect regarding the residues “H 31” and “H 32” indeed, to increase the shape complementarity, the optimization protocol preferred to switch the exposition of these residues. It can be noted that the “H 31” residue, characterized by a high increase in the number of contacts, gains a very high amount of favorable (negative) Coulomb energy. Moreover, even if the number of inter-molecular contacts gained by H 54 is negligible, such a residue (and its neighborhood) acquired in the docking poses an increment of favorable Lennard–Jones energies.

Lastly, we assessed the difference in residue–residue inter-molecular interaction networks. In Figure 4C we reported the contacts between the main couples of residues, where a higher number of occurrences in the docking poses is testified by the thickness and the color of the edge. The spike residues are shown in green, and the antibody ones are in red. As further proof of the goodness of the proposed mutants, it can be noted that the interface of the optimized antibody (lower figure) is much more interconnected than the one of the original antibody (upper figure), indicating a possible effect on binding stability.

We selected 229 protein-binding antibodies with a sequence identity lower than 90% and resolution < 3.0 Å using the SabDab database (Dunbar et al., 2013). The Spike dataset, i.e., the structural dataset of the spike–antibody complexes was built using CoV-AbDab (Raybould et al., 2021). It results in 145 complex structures with a sequence identity lower than 90%, as calculated using CD-HIT (Huang et al., 2010).

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.

The structure of the SARS-CoV-2 spike protein used for the identification of the ACE-2 interacting region has the PDB code 6vw1. When we investigated the N-terminal domain, we used the structure 7jji.

We identified on the spike protein two epitopes in the ACE-2 binding region: spike A and spike B. spike A epitope is constituted by the residues: “TYR 453, LEU 455, PHE 456, ALA 475, GLY 476, PHE 486, ASN 487, TYR 489, and GLN 493.” Spike B epitope is constituted by the residues: “TYR 449, GLY 496, GLN 498, THR 500, ASN 501, GLY 502, and TYR 505.” We identified another epitope in the N-terminal domain, in the region involved in sialic acid binding. That region is defined as the set of residues whose CA atoms are closer than 8 A ̊ to the TRP 258 CA (Milanetti et al., 2021b). Such an epitope is constituted by these residues: “LEU 244, HIS 245, ARG 246, SER 247, TYR 248, LEU 249, THR 250, PRO 251, GLY 252, ASP 253, SER 254, SER 255, SER 256, GLY 257, TRP 258, THR 259, and ALA 260”.

Using as reference the experimental structures, computational mutagenesis has been performed using SCWRL4 (Krivov et al., 2009).

For each protein structure, solvent accessible surface is computed using DMS software with the standard option (Richards, 1977). The interacting surface is constituted by the surface points belonging to interacting residues, defined as the set of residues having at least one atom closer than 4 A ̊ to any atoms of the molecular partner.

Given a 2D function f (r, ϕ) in the unitary circle (region r < 1), it can be expanded in the Zernike polynomials basis. Therefore, f r , ϕ = ∑ n = 0 ∞ ∑ m = 0 m = n c n m Z n m , (5) where c n m = n + 1 π Z n m | f = = n + 1 π ∫ 0 1 d r r ∫ 0 2 π d ϕ Z n m * r , ϕ f r , ϕ (6) are the expansion coefficients (Zernike moments). The complex functions Z _nm (r, ϕ) are the Zernike polynomials, each composed of a radial and an angular part: Z n m = R n m r e i m ϕ . (7)

The radial dependence, given n and m, can be written as follows: R n m r = ∑ k = 0 n − m 2 − 1 k n − k ! k ! n + m 2 − k ! n − m 2 − k ! r n − 2 k . (8)

For each couple of polynomials, the following rule holds: Z n m | Z n ′ m ′ = π n + 1 δ n n ′ δ m m ′ . (9)

Therefore, the set of polynomials forms a basis. Knowing the coefficients, {c _nm } allows the reconstruction of the original function. The level of the detail can be modified by modulating the order of expansion, N = max (n).

The norm of each coefficient (z _nm = |c _nm |) does not depend on the phase; therefore, it is invariant under rotations around the origin.

The shape complementarity between two regions can be evaluated by comparing their Zernike invariants. In particular, we measured the complementarity between regions i and j as the Euclidean distance between the invariant vectors, i.e., d i j = ∑ k = 1 M = 121 z i k − z j k 2 . (10)

We adopted N = 20, therefore, dealing with 121 invariant descriptors for each patch.

In Table 1, we report the mutations proposed as a result of the three Monte Carlo simulations.

We docked the three original and the three optimized antibody structures with spike using HDOCK (Yan et al., 2020), indicating as interacting residues the one written in the Dataset section.

We selected, for all the 6 docking simulations, the best 10 poses according to the Hdock binding score, an iterative knowledge-based scoring function. For each pose we get the following:

• The number of inter-molecular residue–residue contacts. Two residues are in contact if the minimum distance between their atoms is less than 4 A ̊ .

The comparisons between the results regarding original and optimized antibodies are performed so as a positive value means an increase in binding compatibility after optimization. Therefore, the relative percentage of gaining is defined as follows:

• Number of contacts: < C o n t > o p t − < C o n t > o r i g < C o n t > o r i g