Computational Analysis of Mutations in the Receptor-Binding Domain of SARS-CoV-2 Spike and Their Effects on Antibody Binding

Currently, SARS-CoV-2 causing coronavirus disease 2019 (COVID-19) is responsible for one of the most deleterious pandemics of our time. The interaction between the ACE2 receptors at the surface of human cells and the viral Spike (S) protein triggers the infection, making the receptor-binding domain (RBD) of the SARS-CoV-2 S-protein a focal target for the neutralizing antibodies (Abs). Despite the recent progress in the development and deployment of vaccines, the emergence of novel variants of SARS-CoV-2 insensitive to Abs produced in response to the vaccine administration and/or monoclonal ones represent a potential danger. Here, we analyzed the diversity of neutralizing Ab epitopes and assessed the possible effects of single and multiple mutations in the RBD of SARS-CoV-2 S-protein on its binding affinity to various antibodies and the human ACE2 receptor using bioinformatics approaches. The RBD-Ab complexes with experimentally resolved structures were grouped into four clusters with distinct features at sequence and structure level. The performed computational analysis indicates that while single amino acid replacements in RBD may only cause partial impairment of the Abs binding, moreover, limited to specific epitopes, the variants of SARS-CoV-2 with multiple mutations, including some which were already detected in the population, may potentially result in a much broader antigenic escape. Further analysis of the existing RBD variants pointed to the trade-off between ACE2 binding and antigenic escape as a key limiting factor for the emergence of novel SAR-CoV-2 strains, as the naturally occurring mutations in RBD tend to reduce its binding affinity to Abs but not to ACE2. The results provide guidelines for further experimental studies aiming to identify high-risk RBD mutations that allow for an antigenic escape.


Introduction
The recent launch of vaccination campaigns in many countries allowed by the rapid development of several effective vaccines [1] gives hope for a forthcoming amelioration of the world pandemic of SARS-CoV-2. Vaccination will remain the main measure for antiviral protection against COVID19 for a long time, since the development of other types of antiviral drugs is much more time consuming [2].
In [3], it was firstly shown that ACE2 acted as the receptor for SARS-CoV, later in 2020 in [4], researchers provided one of two early experimental structures of the SARS-CoV-2 RBD-ACE2 complex, showing how the Spike protein recognizes its receptor. In [5,6], it was demonstrated that the Spike protein of SARS-CoV-2 and especially its receptor-binding domain (RBD) is one of the major targets of neutralizing antibodies elicited by natural infection or vaccination [7]. The titers of IgM and IgG antibodies against the receptorbinding domain (RBD) of the Spike protein of SARS-CoV-2 decrease significantly over 6 months, with IgA being less affected [8]. At the same time, a number of recent studies have identified viral mutations that escape neutralizing antibodies targeting the SARS-CoV-2 Spike protein [9,10]. Some of these mutations are already present in the human population [11], but many more may be present in natural reservoirs of coronaviruses and represent a potential threat [12,13]. These observations raise concerns about the potency of monoclonal antibodies as well as the protective efficacy of the existing vaccines [14].
Major efforts have been undertaken by the scientific community in order to classify existing data with bioinformatic resources [15][16][17][18] and map potentially hazardous mutations [19]. Particularly, several sites at the SARS-CoV-2 Spike protein, which reduce the neutralizing activity of monoclonal antibodies, and/or their cocktails/human sera, were identified, including E484K, K417N [11], N439K [20], E406W [19], N501Y [21], and others [10]. Many of these mutations occur in the receptor-binding domain (RBD) of Spike, which mediates binding to the angiotensin-converting enzyme 2 (ACE2) receptor, resulting in the virus entry into the cells. As mentioned above, the majority of leading anti-SARS-CoV-2 antibodies also target this domain [6,22,23] rendering these mutations especially risky.
Due to the central role of RBD domain as a key Ab target, in the present study, we carried out a comprehensive computational investigation of the effects which may be induced by its mutations on the affinity to various neutralizing antibodies and hACE2 exploiting structural data available to date. We classified RBD-targeting antibodies and revealed that their epitopes demonstrate remarkable variance in predicted binding energies, the degree of occlusion by the glycan chains, and conservativity. We further assessed the impact of naturally occurring residue replacements in RBD to possible antibody resistance and the ACE2 binding. Furthermore, analyzed the potential outcomes of all possible RBD mutations. We believe that the thorough virtual mutagenesis analysis reported here will guide further experimental studies of Ab resistance and identification of natural SARS-CoV-2 variants capable of antigenic escape.

Analysis of Atomic Structures and Clustering
The structures of all RBD-Ab complexes were retrieved from the PDB database (see Table 1). In order to classify Ab epitopes, each Ab-RBD interface was encoded as a binary vector with the length equal to the number of residues in the reference RBD structure (SARS-CoV-2 Spike protein in complex with the ACE2 receptor, the PDB code 6M17). Positions corresponding to residues in contact (distance between any pair of heavy atoms less than 6 Å) with an Ab were set to 1, while those not forming contacts were set to 0. The encoded epitopes were further clustered by means of the hierarchical clustering algorithm and split into four clusters based on the inspection of the inter-cluster distances. The PDB structure 6M17 of the Spike protein in complex with the ACE2 receptor resolved by Cryo-EM to 2.90 Å [24] was used to estimate the effects of mutations on the RBD binding to ACE2. Table 1. Overview of the analyzed antibody-RBD complexes. Complexes, which are the centers of the clusters in Figure 1, are highlighted in bold. Some structures represent complexes with multiple Abs; therefore, for each complex the specific protein chains corresponding to RBD and Ab heavy/light chains are indicated (only one chain is shown for nanobodies).

Estimation of Binding Energies
For the calculation of binding energies between Abs, ACE2 and RBD, w PRODIGY [42,43]. The contributions of amino acid mutations to the binding energ

Estimation of Binding Energies
For the calculation of binding energies between Abs, ACE2 and RBD, we used PRODIGY [42,43]. The contributions of amino acid mutations to the binding energies were estimated using the BeAtMuSiC server based on application of a statistical model to coarse-grained models of protein-protein complexes [44].
Briefly, the approach predicts a change in protein-protein binding affinity upon a mutation by combining the output of two binding models, one considering that both partners of the interaction are able to fold independently of each other with the change of affinity given by: and another assuming that the partners are unable to fold independently, and the change in binding-free energy on mutation is given just by: where ∆∆G P1 and ∆∆G P2 are the respective folding free energies of the two partners, and ∆∆G C is the folding free energy of the protein-protein complex.
In turn, the effects of mutations on folding free energies can be found using the following expression: where the first term is a sum of energetic changes, ∆∆W i , induced by the mutation according to a set of statistical potentials extracted from a data set of known protein structures and describing the correlations between amino acid types, pairwise inter-residue distances, backbone torsion angles and solvent accessibilities. The terms ∆V + and ∆V − account for the possible creation of packing defects and depend on the change of the side chain volume. The weights α i are sigmoid functions of the solvent accessibility A of the mutated residue, and they were fit on the basis of a data set of experimentally measured changes in folding free energy. We gave preference to the BeAtMuSiC approach showing satisfactory results [45], rather than potentially more accurate methods based on molecular dynamics calculations, in order to accelerate analysis and to be able to extend it to a large set of all possible RBD variants, which would otherwise become computationally infeasible. The effects of multiple mutations were estimated as the sum of individual contributions. Although multiple mutations may have a synergistic effect on binding, we took this approximation for the sake of speeding up the calculations.

Sequence Analysis
The sequences of the variants of the SARS-CoV-2 Spike protein were obtained from a database supported by the GISAID initiative [46] accessed on 15 November 2021. This database served as a source of information about novel RBD variants and the dates of their emergence. A multiple sequence alignment (MSA) was obtained using MAFFT version 7 [47] using the default settings and truncated to the region corresponding to the RBD (C336-L518) present in the reference structure (PDB code 6M17) for further analysis. The identical RBD sequences were sorted, resulting in 3370 unique sequences and the variability of MSA positions was estimated in terms of Shannon entropy using ProDy [48]: where P i is a probability of i-th amino acid at the given position in MSA. The lower and upper limits of Shannon entropy correspond to fully conserved and fully random (equal probability of all twenty amino acid types) amino acids at the given position. In turn, the interface variability was estimated as the average entropy of amino acid positions forming the interface.
Mutations in variable positions along the protein chain can occur because they are either accompanied or preceded by compensatory changes in other variable positions. Such compensation would result in a coupling between changes in the two positions, or coevolution. We used Mutual Information (MI) to identify such potentially co-evolving (i.e., dependent on each other) amino acid residues in RBD as implemented in ProDy [48]. MI corresponds to the reduction of uncertainty (as measured by Shannon entropy defined by Equation (4) of random variable X given random variable Y (and vice versa): Since the inclusion of related sequences in the MSA may lead to an overestimation of the conservation and coevolution propensities by including those amino acids that retain their identity due to common ancestry, we applied a correction the MI matrix as suggested in [49].
The positions were separated by less than 20 Å in the reference RBD structure (6M17) and those regarded to be in a direct contact were set to zero in the MI matrix in order to remove such direct effects from consideration.

Analysis of Epitope Accessibility in the Glycolisalted S Protein
To assess the alterations in epitope accessibility induced by the presence of posttranslational modifications (PTMs), we calculated the relative change of solvent-accessible surface area (SASA) for each residue in RBD in the glycosylated model of Spike protein and in the absence of PTMs. The analysis was performed for a molecular dynamics trajectory of the fully glycosylated Spike protein trimer in a viral membrane, based on the experimental structure of 6VSB reported in [50] and freely available at the CHARMM-GUI platform [51] (system 6VSB_1_1_1). SASA was estimated using the standard Gromacs [52] tool, gmx sasa. The resulting value corresponds to the relative decrease of a residue SASA in the presence of glycan chains averaged over MD trajectory and 3 individual monomers.

Analysis of Epitopes Reveals Distinct Clusters of Ab Binding Poses with Unique Features
The majority of human neutralizing antibodies (Ab) against SARS-CoV2 S-protein bind to diverse epitopes at the surface of its receptor-binding domain (RBD) preventing its attachment to the host cell and thus neutralizing the virus [53]. RBD accounts for 90% of serum neutralizing activity [31]. In order to explore the plasticity of their binding modes and classify them on the basis of residues involved in the complex formation with RBD, we applied hierarchical clustering to the binding interfaces as described in the Methods. The analysis was performed for 35 of the Ab-RBD complexes available in the PDB database (Table 1). We further considered four major clusters of epitopes based on the distances observed in the clustering dendrogram ( Figure 1A). Three of these clusters approximately correspond to three groups of RBD-targeting Abs suggested before [7]. However, our analysis revealed an additional distinct group of epitopes, cluster 3, including Abs such as synthetic nanobody MR17-K99Y and monoclonal Ab S2H14 (see Table 1).
A representative complex for each cluster was chosen based on the estimated bindingfree energy. In each cluster, the Ab-RBD complex with the lowest binding free energy was selected for further analysis (hereafter referred by their names or PDB codes with the prefix standing for the Ab chains in the corresponding structure: 52 Ab/7K9Z_HL (cluster 1, cyan), MR17-K99Y/7CAN_A (cluster 2, yellow), CR3022/6YLA_HL (cluster 3, blue), and CC12.1/6XC2_HL (cluster 4, red)). It may also be noted that complexes with lower binding energy Abs tend to form more contacts to RBD (see Figure S1) implying that focusing our analysis specifically on them should also allow us to examine the corresponding clusters of epitopes in a more complete way in addition to being concentrated at potentially the most affine Abs.
Detailed inspection of the interfaces in each cluster ( Figure 1B) suggests that three extend over the same region of RBD encompassing the receptor-binding motif (RBM) β-hairpin of the RBD and partially overlap with each other and with the ACE2 binding site (clusters 1, 2, and 4; the average overlap with the ACE2 site is 11.2, 53.5, and 74.5%, respectively). However, the fourth cluster (cluster 3) occupies a distinct area on the RBD surface including residues 368-388 forming two α-helices and an intervening β-strand. This cryptic epitope is remote from the ACE2 site and not overlapping with it or epitopes of other clusters of Abs at all. While the Abs binding to this site do not prevent the RBD interaction with ACE2 in the direct way they sterically clash with ACE2 interacting with the same protomer within an S trimer [31].
Interestingly, the identified groups of Abs appear distinct in terms of their predicted binding free energy and number of protein-protein contacts. While Abs belonging to clusters 1 and 2 form a similar or smaller number of contacts to RBD compared to ACE2, and they have comparable or lower predicted affinity ( Figure 1C,D), the Abs of clusters 3 and 4 have larger numbers of contacts resulting in lower binding free energies even exceeding the corresponding values for ACE2. The analysis of epitope shielding by glycans, based on the assessment of solvent-accessible solvent area in the molecular dynamics trajectory of the fully glycosylated Spike trimer [50], indicates that binding interfaces of Abs in cluster 4 and partially in cluster 3 are less occluded by the glycan chains as compared to epitopes belonging to three other clusters ( Figures 1F and 2D) suggesting that they can remain more accessible allowing bound Abs to make more contacts to RBD [54].
Further analysis of amino acid variability of the Spike RBD sequences available to date reveals that clusters 1, 2, and 4, which are spatially adjacent and overlap to some degree with the ACE2 interface, exhibit a significant level of sequence variability measured by means of Shannon entropy normalized by a number of amino acids comprising the epitope (see Figure 1E). Their average entropy is comparable with that of the ACE2 interface, suggesting that both the receptor interface and the epitopes of these Abs are present in multiple variants circulating in the population. This also implies a significant evolutionary adaptation undergone by these RBD epitopes due to the selection of variants, which are capable of escaping the immune response and/or enhance virulence [55]. It is also worth mentioning that the most variable amino acid positions in RBD demonstrate low occlusion by the glycans. Contrarily, the positions with the highest degree of relative glycan screening are quite conservative ( Figures 2C and S3). Altogether, it implies that the unshielded regions targeted by Abs experience the strongest evolutionary pressure, which, in its turn, results in their higher sequence variability.
In contrast to the Abs of clusters 1,2, and 4, the epitope of Abs belonging to cluster 3, which does not overlap with the ACE2 interface at all, demonstrates very low sequence variability, implying that the majority of Spike sequences lack diversity in this region, which is thus unlikely to experience strong selection. Given quite a high predicted affinity of Abs targeting this epitope ( Figure 1C), this observation seems surprising as it implies that these Abs do not significantly contribute to the natural human immune response despite their promising characteristics. One of the possible explanations comes from the structural analysis of CR3022 Ab belonging to this group in complex with the Spike. It shows that this epitope is only accessible when the RBD is in the up conformation in at least 2 Spike proteins within the trimer while the steric clashes with the bound antibodies appear otherwise [56]. Overall, the lack of strong evolutionary pressure on this epitope confirms previously proposed assumption that the group of Abs targeting it are likely immunosubdominant [7].
We also assessed how the variability of different epitopes/ACE2 interfaces has changed during the course of pandemic as well as how the novel missense mutations were accumulated in these regions of RBD. The results clearly show that novel mutations appear in the all analyzed epitopes regardless of their proximity to the ACE2 interface at about the same rate, which, as expected, decreases during the course of pandemics ( Figure S2) as the pool of novel unique mutations, which are beneficial or non-deleterious, is limited and rapidly depletes over time.

Single RBD Mutations Have Limited Effect on its Binding Affinity to Abs and ACE2
A substantial number of mutations in the RBD can be well tolerated witho noticeable consequences for ACE2 binding [62]. This raises the concern of antige escape, which may severely limit the therapeutic potential of neutralizing Abs [63]. W employed virtual mutagenesis analysis in order to assess the potential effects of po mutations in RBD on its binding affinity to antibodies belonging to four identified clust and to the SARS-CoV-2 natural receptor, ACE2. Firstly, we estimated the effects on binding energy for those single amino acid substitutions that are the most comm mutations in the reported viral genomes (Table 2). While almost all of them (except A475V, which might slightly increase affinity towards Ab CC12.1, see Figure 3) appear decrease the binding affinity to either all studied Abs (e.g., N439K) or some of them (e E484K), this effect is low. Notable exceptions are the G446V and K417N RBD variants, which the binding energy is increased by 1.3 kcal/mol and 0.9 kcal/mol towards MR At the same time, the variability of three out of the four investigated epitopes, i.e., those of 52, MR17-K99Y, and CC12.1 Abs, underwent a significant increase in 2021 compared to its rather low values observed before, which did not noticeably change till the end of 2020 (Figure 2A). Similar trend is observed for the ACE2 interface. Contrary to them, the variability of the CR3022 epitope does not experience any changes during the whole period of pandemics and it remains very low until now suggesting its high conservativity among the all of variants ever spread in the population.
The observed increase in variability of three epitopes overlapping with the ACE2 interface may result from the epistatic effects [57] of some RBD variants [58], which became widespread in the population during late 2020 and allowed for the fixation of otherwise compromised mutations. Intriguingly, this increase coincides in time with the period from ca. June 2020 to January 2021, when the rate at which the novel unique mutations emerged in RBD seems steady in contrast to the exponential decrease in the accumulation rate of novel mutations during both the previous and the consequent months (see Figure S2). One may hypothesize that it reflects the emergence of novel mutations during this time, which were epistatically dependent on previously emergent mutations and which could not be otherwise tolerated. Indeed, the analysis of coevolved (i.e., synchronously mutating) amino acid residues in RBD by means of estimation of mutual information (MI) reveals coupling between spatially remote positions separated by more than 20 Å (see Figure S4, the residue pairs with the highest coupling are N501-E484, N501-T478, N501-K417, E484-K417, and T478-L452). Such coupling may result from the indirect influence of one residue on other positions through some intermediate residues; however, it can also signify allosteric coupling in some cases [59] in line with the idea of possible epistatic interactions between certain positions in RBD. This quite intriguing hypothesis is worth further scrutiny.
Epistasis has previously been studied intensely in the influenza surface protein hemagglutinin, and positive epistasis was detected in several regions of the hemagglutinin receptor-binding domain [60]. In relation to the Spike protein of SARS-CoV-2, epistatic mutations could, for instance, allow the Spike protein to adopt a specific conformation enabled by the presence of a distinct set of mutations resulting in unique protein structural and dynamical features. For instance, a recent study demonstrated the impact of antigenicity and infectivity of the Spike D614G SARS-CoV-2 variant in combination with different other mutations [61]. The study shows that D614G alone only slightly affects the infectivity, but in combination with some mutations in Spike, they can both increase and decrease viral infectivity to a larger extent. Similar findings have been reported regarding potential escape from single-domain antibodies. Similarly, D614G alone cannot escape nanobodies. However, the combination of D614G with other Spike mutations can enable the antibody escape [61]. This data together with our coevolutionary analysis suggests that the emergence of epistatic mutations in Spike protein will likely be involved in further adaptability of SARS-CoV-2, including the enhanced transmissibility, stability, and antibody resistance.

Single RBD Mutations Have Limited Effect on its Binding Affinity to Abs and ACE2
A substantial number of mutations in the RBD can be well tolerated without noticeable consequences for ACE2 binding [62]. This raises the concern of antigenic escape, which may severely limit the therapeutic potential of neutralizing Abs [63]. We employed virtual mutagenesis analysis in order to assess the potential effects of point mutations in RBD on its binding affinity to antibodies belonging to four identified clusters and to the SARS-CoV-2 natural receptor, ACE2. Firstly, we estimated the effects on the binding energy for those single amino acid substitutions that are the most common mutations in the reported viral genomes (Table 2). While almost all of them (except for A475V, which might slightly increase affinity towards Ab CC12.1, see Figure 3) appear to decrease the binding affinity to either all studied Abs (e.g., N439K) or some of them (e.g., E484K), this effect is low. Notable exceptions are the G446V and K417N RBD variants, for which the binding energy is increased by 1.3 kcal/mol and 0.9 kcal/mol towards MR17-K99Y (PDB: 7CAN, cluster 2) and CC12.1 (PDB: 6XC2, cluster 4) Abs, respectively. These ∆∆G values should lead to~8.2 and~4.4 fold-change drops in K D of MR17-K99Y (PDB: 7CAN) and CC12.1 (PDB: 6XC2) Abs, respectively, at 310 K assuming the predicted binding-free energies for these Ab-RBD complexes (see Table 1). At the same time, most mutations also result in lower affinity to ACE2. Remarkably, the most widespread RBD mutation in the population to date, N501Y, which is particularly specific to α, β, γ, θ, µ, and o variants, either does not cause any changes in affinity (to 52, MR17-K99Y, and CR3022 Abs) or even results in a slight increase in the affinity to CC12.1 Ab and the native receptor, ACE2, which might account for the fixation of this mutation in the population. It is also worth mentioning that the affinity of CR3022 (PDB: 6YLA), the representative Ab of cluster 3, seems to be the least affected by the analyzed mutations compared to Abs representing the other three clusters. This is apparently due to the remoteness of its epitope from the most frequent mutations, which mainly occur in the variable RBM region of RBD. Table 2. SARS-CoV-2 variants of concern and variants of interest considered by WHO along with the corresponding mutations in RBD. demonstrated to reduce the affinity to sera/monoclonal Ab [64]. Table 2. SARS-CoV-2 variants of concern and variants of interest considered by WHO along with the corresponding mutations in RBD.  Overall, these results suggest that none of the currently widespread single mutations of RBD are able to completely break the binding between RBD and the broad spectrum of Abs targeting it. However, two mutations may lead to partial loss of affinity to some of the analyzed Abs, such as K417N and G446V. The former replacement, which is specific for the B.1.351 lineage (known as the South African variant) has already been demonstrated to reduce the affinity to sera/monoclonal Ab [64].

Variant(s) RBD Mutations
Since CoVs undergo continuous and extensive antigenic evolution [65], which can potentially lead to amino acid variations that are not currently observed but can appear in future, we did not limit our analysis to the existing variants but performed the complete virtual mutation scanning of RBD (notably, the least conserved domain of Spike [20]) by systematically changing all amino acid positions to 19 alternatives. The latter task was allowed here by exploiting the fast computational approach for virtual mutagenesis.
The effects of the majority of mutations on the binding affinity is very low, as indicated by corresponding distributions of the predicted ∆∆G values (Figure 4, middle panel). However, all of these distributions are not symmetrical but have a longer right tail corresponding to positive ∆∆G values and, thus, the decrease in affinity. In other words, the number of mutations that lower affinity and the extent of this lowering effect is greater than for mutations that strengthen the binding. This observation is also true for the RBD-ACE2 interaction (see Figure 4N).    Despite the very low effects induced by the vast majority of mutations on the RBD-Ab/ACE2 binding energy, some mutations may potentially lead to a change by up to +3-4 kcal/mol, which might already reduce the affinity between proteins drastically, given the typical Ab-RBD binding energies (see Table 1, ∆G average = −11.7 kcal/mol; experimental values fall into the same range, e.g., for CR3022 K D =~6.6 nM-~115 nM [36,56,66] corresponding to ∆G = −11.6-9.8 kcal/mol; for ACE2-K D =~15nM [22], i.e., ∆G = −11.1 kcal/mol at 310 K).
Importantly, mutations with strong effects on the binding energy usually occur directly or in a close vicinity of an Ab epitope, often also resulting in the reduced affinity of RBD towards ACE2, since for the majority of investigated Abs their epitopes overlap with the ACE2 interface. Indeed, we could observe a positive correlational trend between ∆∆Gs estimated for mutations in the RBD-ACE2 complex and ∆∆Gs for the same mutations in all of the studied Ab-RBD complexes (see Figure 5, Pearson's r = 0.59÷0.81). This fact may account for the lack of a major fraction of such deleterious mutations in real SARS-CoV-2 sequences, due to their interference with the RBD-ACE2 binding. Indeed, the comparison of the effects of real RBD mutations observed in deposited SARS-CoV-2 sequences onto Abs and ACE2 binding implies that they tend to reduce RBD-ACE2 binding to a lesser degree as compared to their effect on RBD-Ab affinities (see black dots in Figure 5). In other words, it appears that the mutations that promote antigenic escape but at the same time do not affect ACE2 binding are preferably selected during the SARS-CoV-2 evolution. Quite surprisingly, this trend is pronounced for 52 (PDB: 7K9Z), MR17-K99Y (PDB: 7CAN), and CR3022 (PDB: 6YLA) Abs but not for CC12.1 (PDB: 6XC2). A possible explanation for this discrepancy might be that the epitope of the latter Ab has the largest overlap with the interface of ACE2, leaving very little room for such adaptations in RBD, which lead to the antigenic escape from these Abs without affecting the ACE2 binding.

Antigenic Escape Assessment for Existing SARS-CoV2 Variants with Multipl
Finally, we estimated potential changes of the binding energy between R studied Abs caused by multiple mutations co-occurring in the real viral sequ in the GISAID database to date. The total affinity change was estimated as binding energy changes induced by all the individual mutations in a given s compared with the cognate change of the RBD affinity towards ACE2.
In general, the results indicate that a higher number of mutations in RB leads to larger ΔΔG (see Figure 6A, Pearson's r = 0.57 ÷ 0.84). At the same ti two substitutions are already sufficient to reduce the RBD affinity to indiv up to 2.5-3.0 kcal/mol, while six substitutions may result in ΔΔG up to 3.3-Some SARS-CoV2 variants, including the variants with the largest estimat RBD-Ab/ACE2 affinity and widely spread variants, are listed in Table 3.
None of the widely spread variants of concern except the recently emer (ο) seem to be able to escape the studied Abs efficiently. The RBD of the carries five times more mutations than any previously concerned variant of and may indeed be potentially resistant against some groups of Abs, inclu CC12.1. On the other hand, this developed resistance likely comes at a cos affinity to ACE2, as indicated by the corresponding ΔΔG values. Both obse support in the preliminary experimental reports, thereby indicating both ability of the omicron variant to avoid Abs and its decreased morbidity [67,6 The most drastic decrease in affinity was predicted for MR17-K99Y CC12.1 (6XC2) Abs towards the RBD variants found in animals, i.e., bats an These RBD sequences comprise as many as 20-25 mutations. At the same tim decrease in the affinity towards 52 (7K9Z) and CR3022 (6YLA) Abs was estim two sequences obtained from human samples in Iran and accommodat substitutions, respectively. In all of these cases, the predicted ΔΔG (~4.3-~ may potentially lead to the complete evasion of these variants from specific worth mentioning that all of the variants with the largest affinity decre particular Ab also , to some extent (at least by ~2.1 kcal/mol), demonstra affinity towards other examined Abs. The same is true for their affinity to A dropped by 2.7-6.9 kcal/mol. Furthermore, a positive correlation between ΔΔGs and the RBD-ACE2 ΔΔGs can be noted ( Figure 6B, Pearson's r = 0.63 to the correlation observed for single mutations and discussed in the previo

Antigenic Escape Assessment for Existing SARS-CoV2 Variants with Multiple Mutations
Finally, we estimated potential changes of the binding energy between RBD and four studied Abs caused by multiple mutations co-occurring in the real viral sequences stored in the GISAID database to date. The total affinity change was estimated as a sum of the binding energy changes induced by all the individual mutations in a given sequence, and compared with the cognate change of the RBD affinity towards ACE2.
In general, the results indicate that a higher number of mutations in RBD expectedly leads to larger ∆∆G (see Figure 6A, Pearson's r = 0.57 ÷ 0.84). At the same time, as few as two substitutions are already sufficient to reduce the RBD affinity to individual Abs by up to 2.5-3.0 kcal/mol, while six substitutions may result in ∆∆G up to 3.3-4.0 kcal/mol. Some SARS-CoV2 variants, including the variants with the largest estimated change in RBD-Ab/ACE2 affinity and widely spread variants, are listed in Table 3.

R REVIEW
14 of 18 towards specific Abs is decreased by ~2-4 kcal/mol while there is no notable change of affinity to ACE2 ( Figure 6B), although for CC12.1 (6XC2) Ab the effect is lower and does not exceed 2 kcal/mol, apparently for the same reason discussed with respect to the single mutations in the previous section, i.e., it is due to the largest overlap of the CC12.1 Ab epitope with the ACE2 interface, making those RBD mutations, which are unfavorable for the binding of this particular Ab, also adverse for the RBD interaction with ACE2. On the other hand, one cannot exclude a possibility that variants escaping even Abs targeting this epitope, which seem the most resistant to the development of resistance, and, at the same time, interacting in the normal way with ACE2 can appear in future given that over 2021 a significant number of SARS-CoV-2 variants with multiple mutations in RBD emerged comparing to 2020 (see Figure S5). The analysis of single mutations performed above indicates that as many as six mutations (P337G, R403P, T415C, N422P, N481P, and G502R) may potentially result in the substantial decrease in the RBD affinity to the CC12.1 Ab by ~5.6 kcal/mol while the affinity of this variant to ACE2 would even increase by ~0.5 kcal/mol. The full list of single mutations, which lower the binding affinities between RBD and Abs (ΔΔG > 0.6 kcal/mol) but do not lower the binding affinity to ACE2 (ΔΔG ≤ 0 kcal/mol), is provided in Table S1.
Overall, the analysis of RBD variants with multiple mutations suggests the emergence of novel SARS-CoV-2 variants with multiple co-occurring mutations as the major risk factor in terms of resistance development against Abs.   None of the widely spread variants of concern except the recently emerged B.1.1.529 (o) seem to be able to escape the studied Abs efficiently. The RBD of the latter variant carries five times more mutations than any previously concerned variant of SARS-CoV-2 and may indeed be potentially resistant against some groups of Abs, including 52 and CC12.1. On the other hand, this developed resistance likely comes at a cost of reduced affinity to ACE2, as indicated by the corresponding ∆∆G values. Both observations find support in the preliminary experimental reports, thereby indicating both an increased ability of the omicron variant to avoid Abs and its decreased morbidity [67,68].
The most drastic decrease in affinity was predicted for MR17-K99Y (7CAN) and CC12.1 (6XC2) Abs towards the RBD variants found in animals, i.e., bats and pangolins. These RBD sequences comprise as many as 20-25 mutations. At the same time, the largest decrease in the affinity towards 52 (7K9Z) and CR3022 (6YLA) Abs was estimated for the two sequences obtained from human samples in Iran and accommodated 8 and 11 substitutions, respectively. In all of these cases, the predicted ∆∆G (~4.3-~9.5 kcal/mol) may potentially lead to the complete evasion of these variants from specific Abs. It is also worth mentioning that all of the variants with the largest affinity decrease towards particular Ab also, to some extent (at least by~2.1 kcal/mol), demonstrate decreased affinity towards other examined Abs. The same is true for their affinity to ACE2 which is dropped by 2.7-6.9 kcal/mol. Furthermore, a positive correlation between the RBD-Ab ∆∆Gs and the RBD-ACE2 ∆∆Gs can be noted ( Figure 6B, Pearson's r = 0.63-0.90) similar to the correlation observed for single mutations and discussed in the previous section. It does not exclude, however, the presence of some RBD variants for which the affinity towards specific Abs is decreased bỹ 2-4 kcal/mol while there is no notable change of affinity to ACE2 ( Figure 6B), although for CC12.1 (6XC2) Ab the effect is lower and does not exceed 2 kcal/mol, apparently for the same reason discussed with respect to the single mutations in the previous section, i.e., it is due to the largest overlap of the CC12.1 Ab epitope with the ACE2 interface, making those RBD mutations, which are unfavorable for the binding of this particular Ab, also adverse for the RBD interaction with ACE2. On the other hand, one cannot exclude a possibility that variants escaping even Abs targeting this epitope, which seem the most resistant to the development of resistance, and, at the same time, interacting in the normal way with ACE2 can appear in future given that over 2021 a significant number of SARS-CoV-2 variants with multiple mutations in RBD emerged comparing to 2020 (see Figure S5). The analysis of single mutations performed above indicates that as many as six mutations (P337G, R403P, T415C, N422P, N481P, and G502R) may potentially result in the substantial decrease in the RBD affinity to the CC12.1 Ab by~5.6 kcal/mol while the affinity of this variant to ACE2 would even increase by~0.5 kcal/mol. The full list of single mutations, which lower the binding affinities between RBD and Abs (∆∆G > 0.6 kcal/mol) but do not lower the binding affinity to ACE2 (∆∆G ≤ 0 kcal/mol), is provided in Table S1.
Overall, the analysis of RBD variants with multiple mutations suggests the emergence of novel SARS-CoV-2 variants with multiple co-occurring mutations as the major risk factor in terms of resistance development against Abs.

Conclusions
The great effort with which the scientific community is studying SARS-CoV-2 has generated a huge amount of data about this novel coronavirus since its emergence about two years ago. The open databases of sequences and atomic structures provide an opportunity to perform the high-throughput analysis of interactions between various monoclonal neutralizing antibodies (Abs), hACE2, serving as the main entry point for coronavirus, and their principal viral binding partner, the receptor-binding domain (RBD) of Spike. Here, we employed a theoretical approach to predict the effects of the known and conceivable RBD mutations on the binding energy between RBD, ACE2 and a set of selected representative Abs targeting structurally distinct epitopes, and augmented this analysis with the sequencelevel data about the existing SARS-CoV-2 variants.
An analysis of epitopes revealed four major clusters of Ab binding sites differing in both of the involved RBD residues, their conservativity, and degree of the glycan shielding. The assessment of epitope conservativity suggested that Abs belonging to one of the identified clusters (including, among others, CR3022) are likely immunosubdominant despite their seemingly high affinity to RBD. Although most mutations that are currently widely spread appear in the least conservative part of RBD, the receptor-binding motif (RBM), we show that there also exist potentially deleterious mutations in more conservative regions of RBD, which can significantly impair its binding to some Abs. While certain point replacements in RBD, including those previously characterized experimentally and/or observed in the population, can reduce the binding affinity for specific Abs or even several types of Abs to a certain degree, the cumulative effects of multiple mutations may be much more hazardous, thereby leading to resistance against a broader spectrum of Abs as demonstrated by the recently emerged omicron variant. By analyzing sequences of the existing SARS-CoV-2 variants, we identified a number of such naturally occurring mutants potentially capable of extensive antigenic escape.
We expect that our theoretical study will promote future experimental work, complement the eventual prognosis of the properties of the novel SARS-CoV-2 variants, and help to design more efficient treatments for COVID-19. Particularly, the knowledge-based selection of antibody mixtures with non-overlapping escape mutations should reduce the emergence of resistance and prolong the utility of antibody therapies.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/v14020295/s1, Figure S1: Theoretical estimation of the free energy of the Ab-RBD complexes plotted as a function of the contact count. The linear regression trend is shown.; Figure S2: Number of novel unique mutations appeared at the RBD interfaces of selected Abs and ACE2 during the course of pandemics; Figure S3: Shannon entropy and per-residue change of solvent accessible surface area (SASA) for RBD residues; Figure S4: Mutual Information (MI) matrix estimated from the MSA of RBD; Figure S5: Number of RBD sequences stored in the GISAID database during 2020 and during 2020-2021 periods with the given mutation counts with respect to the original Wuhan variant. Number of multiple mutants increases significantly in 2021 compared to 2020; Table S1: List of mutations, which lower the binding affinities between RBD and antibodies (∆∆G > 0.6 kcal/mol) but do not lower the binding affinity to ACE2 (∆∆G ≤ 0 kcal/mol).