This work describes a computational method to evaluate the effect of different SARS-CoV-2 mutations on the binding affinity of available NAbs and on the stability of the complex. As a working pattern, we developed a reference complex model between the NAb S309 and the original SARS-CoV-2 Wuhan strain. We evaluated the other variants and subvariants based on the differences of their specific molecular interactions and Gibbs free energy (ΔG) with S309. Figure 1 outlines the methods used to determine the anti–SARS-CoV-2 antibody neutralization potential of S309.

Several reports have discussed the neutralizing effect of NAbs and possible antibody escape of some new Omicrons subvariants [34-39]. Here we used our developed method to evaluate the binding affinity of several of these new subvariants including AY.1, XBB.1.5, BF.7, BQ.1.1, BA.1.1, BA.3, BA.2.3.20, BM.1.1.1, BA.5.6.2, BA.2.75.2, and CH.1.1 (Orthrus), with the NAb S309. Additionally, the effect of several amino acid substitutions in the NAb epitope have been tested experimentally using the enzyme-linked immunosorbent assay and/or pseudovirus neutralization assays. Several mutations are reportedly resistant to inhibition by S309 leading to an antibody escape. These key residues include R346, P337, G339, N440, and S371 [40,41]. Therefore, we applied our method to computationally test the effect of some mutations on these residues. As we already generated parent RBD sequences, newly emerged mutations were introduced, new models and complexes were built, and the mutation’s effect on binding energy with the NAb was predicted by recalculating complex’s ΔG in reference to that of the parent complex and binding affinity with the original SARS-CoV-2 Wuhan strain.

This study was exempt from ethical review since it was conducted in silico and no human subjects were involved.

Figure 1 outlines the methods for assessing the anti–SARS-CoV-2 neutralization potential of S309. The blueprint of the method we developed using monoclonal antibody S309—an experimentally proven neutralizing monoclonal antibody for SARS-CoV-2 and its variants—is described in Figure 1A. We proceeded by modifying the available model 7YAD to generate a reference model that can be used to measure neutralization potential in terms of binding affinity ΔG (Figure 1B). Several in silico 3D models representing spike monomer chain of each variant were generated. The quality of the generated 3D model was evaluated based on the homology modeling report and SWISS-MODEL structural assessment. The generated models showed a QMEAN z score between –1.0 and –3.2 indicating a good-quality model where z scores of around 0.0 are ideal and any value below –4.0 indicates a low-quality model [42]. The QMEANDisCo global score represents the combined scoring of global (for the entire structure) and local (per residue) absolute quality estimates of a single model [43]. Our models’ QMEANDisCo global scores ranged from 0.64 to 0.76 (SD 0.05). These values reflect a good-quality model (any value below 0.6 represents a low-quality model). Each complex was built by molecular replacement of chain M of the reference model with the extracted RBD, followed by binding affinity and interaction analyses.

The generated complexes were energy-minimized and polar and hydrophobic interactions were analyzed. Several interactions were identified between the S309 Fv domain and spike RBD with more interactions toward the heavy chain. Interacting residues of the spike protein include residue 321-428 in SARS-CoV and 334-441 in SARS-CoV-2 and its variants. SARS-CoV showed 4 polar interactions compared to the original SARS-CoV-2 Wuhan strain that shares a total of 3 polar interactions with S309. Interestingly, variant Kappa showed the highest number of polar interactions (n=6), while variant Delta-21J showed the lowest (n=1) number of polar interactions. Variant Kappa showed 2 unique salt bridges between residues R346 and K356 with the S309 heavy chain residue E108. All the variants share the same polar interaction between E340 and S309 heavy chain A104 except for variant Delta-21J. All Omicron subvariants showed the same interaction pattern except for BA.2.75 with 1 missing polar interaction between T343 and S309 heavy chain S109. Variant Gamma showed more hydrophobic interactions with the light chain of S309. All polar interactions are represented in Figure 3 and detailed interactions are listed in Multimedia Appendix 1.

The thermodynamic stability of the generated complexes was measured via computational prediction of ΔG using the CSM-AB tool. ΔG reflects energy differences between coupled and decoupled antibody-antigen complexes. This difference in energy indicates complex stability where a negative normalized energy (ΔG<0) indicates spontaneous and exergonic reactions and hence more stable complexes and more efficient protein–ligand interactions. Thus, the lower the value of ΔG, the more stable the (antibody-antigen) complex. In our model, we found that the NAb S309 has a binding affinity of –8.26 kcal/mol with SARS-CoV and –7.13.26 kcal/mol with SARS-CoV-2, indicating a loss in binding affinity. However, comparing SARS-CoV-2 variants to the binding affinity of the first Wuhan strain showed an improvement in the binding affinity of S309 with variants Alpha, Beta, Gamma, and Kappa. This improvement in affinity, when compared to the interaction profile, can be related to the increased number of polar and hydrophobic interactions and more similar interaction profiles with SARS-CoV than with the original SARS-CoV-2 Wuhan strain. In contrast, variant Delta showed a substantial decrease in binding affinity as it exhibited only 1 polar interaction. All Omicron subvariants shared similar interaction profiles; however, they exhibited different binding affinities. Although they showed a significant decrease in binding affinity compared to original SARS-CoV-2 Wuhan strain, they can be clustered in 2 groups: those with a G339H mutation (BA.2.75, XBB, and XBB.1) and those with a G339D mutation (BA.1, BA.2, BA.4/5, BQ.1, and BA.2.12.1; Figure 4 and Multimedia Appendix 1). The data show that the H339 residue slightly enhanced binding affinity compared to the D residue substitution. This residue is located in the middle of the interaction loop and hence plays a marked role in maintaining the complex’s stability and binding affinity. In addition, our results are in line with the reported effect of the G339D mutation and its role in escaping antibody neutralization [41,44,45].

Furthermore, to test the impact of a mutation in residue G339, we analyzed the effect of reverse mutagenesis. We used the generated models and in silico tools to test the effect of reverse mutation at residue G339 on complex stability in subvariants BA2.75, XBB, and XBB.1. They have an aspartic acid residue at position 339. By reversing this residue to either glycine or histidine (G339 or H339), we calculated the effect in the form of the ΔG value. Our results showed an increase in the stability of the SARS-CoV-2/S309 complex and hence enhanced binding affinity with the glycine residue. However, reverse mutagenesis to histidine has no to a very low effect, except for subvariant BA.2.12.1 where there was a slight increase in binding affinity (Table 2).

The effect of several amino acid substitutions in the NAb S309 epitope have been tested experimentally using the enzyme-linked immunosorbent assay and/or pseudovirus neutralization assays. These mutations resulted in resistance to neutralization by S309, leading to antibody escape. These key substitutions include R346S and P337L, G339D, N440K, and S371L [40,41]. Here we used our developed method to evaluate this effect computationally. By generating models with the newly reported mutations and CSM-AB tool, we predicted the effect of the reported mutations on the binding affinity of the complex and hence on neutralizing effect of S309. Interestingly, our computational results are comparable with the experimentally reported effect of these mutations on the S309 evasion from monoclonal antibodies. Additionally, we predicted a possible effect of hypothetical mutations on some of the proteoglycan epitopes (Table 3).

Antibody-based therapies have proven effective against SARS-CoV-2 infection and appear to be the most promising approach to control the COVID-19 pandemic. A number of neutralizing monoclonal antibodies used in the clinical setting have shown highly favorable results, particularly in stopping disease progression [46,47]. However, the constant emergence of new virus variants has hindered the potency of available anti–SARS-CoV-2 antibodies and urged the continuous development of improved, more effective NAbs. In this study, we describe an in silico rapid method that we developed to predict a possible effect of newly emerged mutations on the efficacy of available neutralizing anti–SARS-CoV-2 antibodies. We used the monoclonal antibody S309 as an example. S309 recognizes a proteoglycan epitope embedded in a structural loop located on the outer side the SARS-CoV-2 spike protein and encompasses residues 334-441 (Multimedia Appendix 1). This specific epitope location permits the binding to RBD in both the up and down configurations without affecting binding to the ACE2 receptor. Indeed, this epitope does not overlap with the ACE2 binding site. However, several newly emerged RBD mutations were reported to have an impact on the neutralizing effect of S309. To further explore this, we developed this computational method to evaluate and compare the neutralization potential of S309 against different SARS-CoV-2 variants and possible new emerging mutations (Figure 1).

Using bioinformatics tools, we developed spike models for several new SARS-CoV-2 variants and evaluated the effect of several emerged mutations on the interaction with the neutralizing monoclonal antibody S309 used for the treatment of mild-to-moderate COVID-19. In addition, by applying this method, we foresee the effect of some predicted or not yet observed mutations. Interestingly, the predicted significantly decreased computational neutralization values of the monoclonal antibody S309 (from 10% to 50%) for some new Omicron subvariants are confirmed by the newly published clinical results indicating a reduction in its effectiveness against these same new Omicron subvariants and possible immune evasion [39,48-51]. Early on, S309 was clinically considered one of the most effective monoclonal antibodies against all SARS-CoV-2 variants [7]. However, this statement has been proven wrong as recent convergent evolution of Omicron and its subvariants has led to a new set of spike mutations within the S309 epitope, and, consequently, the new subvariants became increasingly resistant [52]. Several mutations were identified to be critical, and others are yet to be investigated. For example, a substitution in the nonpolar G339 residue located at the center of the antibody epitope to the acidic charged aspartic acid residue (G339D) has been shown to have a remarkable impact on the binding affinity of Omicron’s subvariants [44,53], with a predicted reduction in neutralization power of 30% for BA.1; 45% for BA.4, BA.5, and BQ.1; 50% for BA.2.12.1 and BA.2; and 60% for BF.7 and BQ.1.1. We reported a similar effect in our proposed computational method and we found that the impact was less intense with the G339H mutation (Table 2 and Multimedia Appendix 1). However, the combination of multiple mutations in Omicron subvariants has a more profound effect on binding affinity, indicating increased antibody resistance. This effect was clearly detected in the subsequent, potentially dominant new subvariants BM.1.1.1, BA.2.3.20, and CH.1.1 (Orthrus) [54] (Table 3). Furthermore, to test our method, we examined some experimentally evaluated mutations in residues P337, R346, G339, and S371 that are located in the S309 epitope, and once more, our computational method was compatible with the experimental results (Table 3). This reduced susceptibility of S309 with mutations in residues P337, R346, and other residues has been experimentally recognized [13,40,41]. Considering the clinical observations of the efficiency of Sotrovimab in neutralizing SARS-CoV, SARS-CoV-2 variants, and Omicron subvariants, a 50% reduction in binding affinity, compared to that in the reference model, may be considered the cutoff for determining whether a monoclonal antibody will neutralize a new variant, using the method described in this paper. Comparison of the predicted values of the evaluation of neutralizing power with a larger number of clinical observations about the efficiency of a neutralizing monoclonal antibody would help refine this theoretical cutoff value and further validates the method. Ultimately, molecular dynamics simulations can be performed to more accurately define the most stable conformation of monoclonal antibody/spike protein–RBD complexes.

This in silico method provides significant insights into possible antibody escape following the emergence of new SARS-CoV-2 mutants and helps evaluate the usefulness of existing NAbs in combating new emerging variants and subvariants. This method is straightforward, rapid, and applicable ahead of obtaining statistically significant clinical observations. In addition, this method highlights the advantages of computational approaches in viral the rapid surveillance and for the development of novel monoclonal antibody therapies.