Computational approach for binding prediction of SARS-CoV-2 with neutralizing antibodies.

Coronavirus disease 2019 (COVID-19) caused by a novel severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has spread worldwide as a severe pandemic and caused enormous global health and economical damage. Since December 2019, more than 197 million cases have been reported, causing 4.2 million deaths. In the settings of pandemic it is an urgent necessity for the development of an effective COVID-19 treatment. While in-vitro screening of hundreds of antibodies isolated from convalescent patients is challenging due to its high cost, use of computational methods may provide an attractive solution in selecting the top candidates. Here, we developed a computational approach (SARS-AB) for binding prediction of spike protein SARS-CoV-2 with monoclonal antibodies. We validated our approach using existing structures in the protein data bank (PDB), and demonstrated its prediction power in antibody-spike protein binding prediction. We further tested its performance using antibody sequences from the literature where crystal structure is not available, and observed a high prediction accuracy (AUC = 99.6%). Finally, we demonstrated that SARS-AB can be used to design effective antibodies against novel SARS-CoV-2 mutants that might escape the current antibody protections. We believe that SARS-AB can significantly accelerate the discovery of neutralizing antibodies against SARS-CoV-2 and its mutants.

Introduction

Three highly pathogenic human coronaviruses have emerged during the past 20 years, including the Middle East respiratory syndrome coronavirus (MERS-CoV), the severe acute respiratory syndrome coronavirus (SARS-CoV), and a 2019 novel coronavirus (SARS-CoV-2) [1], [2], [3]. Among them, SARS-CoV-2 has exceeded both SARS-CoV and MERS-CoV in its rate of transmission among humans [4], [5], and caused the global pandemic of COVID-19.

Recently, three COVID-19 vaccines have been authorized in the United States for emergency use by the US Food and Drug Administration (FDA). Pfizer and Moderna vaccines work by delivering mRNA into host cells to allow expression of the SARS-CoV-2 S antigen [6], [7], while Johnson&Johnson vaccine represents a viral vector. The vaccine elicits an immune response to the S antigen, which protects against COVID-19. In addition to the prophylactic approach, COVID-19 therapies based on monoclonal antibodies have also been developed to treat mild to moderate COVID-19 in adults and pediatric patients with positive results of direct SARS-CoV-2 viral testing and who are at high risk for progressing to severe COVID-19 [8]. Nonetheless, the development of antibody therapies against COVID-19 is still at the beginning phase, with the purpose to overcome multiple virus variants that might escape the current reagents [9].

Although antibody therapy can be extremely effective in treating COVID-19 and preventing severe disease symptoms, finding high-affinity antibodies usually requires the labor intensive screening assays, which cannot scale up to more than a few thousand of target sequences. Therefore, a computational method that reliably predicts RBD-binding antibodies accelerate antibody discovery and reduce the cost. In this work, we developed a novel computational approach (SARS-AB) based on molecular dynamics (MD) simulations and an energy estimation method of antibodies-RBD-SARS-CoV-2 contacts. MD simulations have been successfully used in drug discovery [10], [11], [12], [13], lead optimization, exploring structural changes in proteins, providing energetic information about protein and ligand interactions, and validation of protein–ligand models deposited to PDB [14]. The proposed SARS-AB performs de novo prediction of the binding energy of a given antibody to the RBD region directly from the full-length antibody DNA or protein sequences, which can be applied to most antibody screening studies to identify promising therapeutic candidates against SARS-CoV-2. We benchmarked SARS-AB using putative RBD neutralizing antibodies from the recent literature to evaluate its ability to distinguish SARS-CoV-2 binders from non-binders, where SARS-AB achieved 100% accuracy. Finally, SARS-AB was validated using independent datasets of experimentally acquired antibodies without available crystal structures, and consistently reached high prediction accuracies (AUC = 99.6%).

New variants of SARS-CoV-2 may cause rapid spread and antibody neutralization escape of the virus [15]. The greatest concern is caused by a new lineages of SARS-CoV-2 in South Africa, 501Y.V2, which contains multiple mutations in RBD (N501Y, K417N and E484K) and NTD domains (L18F, D80A, D215G, Δ242–244, and R246I) [16], Brazil, P.1, which harbors 21 lineage-defining mutations including three in RBD (N501Y, K417T, E484K) [17] India, B.1.617.2 [18], containing two mutations in RBD (L452R, T478K) and the latest severe acute SARS-CoV-2 Omicron variant B.1.1.529 containing 15 mutations in the RBD domain. The 501Y.V2, P.1, B.1.617.2 and B.1.1.529 variants have been associated with virus increased transmissibility and neutralization escape from some classes of SARS-CoV-2 monoclonal antibodies [18], [19], [20], [21]. In this work we made an attempt to predict the impact of 501Y.V2, P.1 and B.1.1.529 spike amino acid changes on a binding affinity of potent antibody, Regdanvimab, and proposed the possible modifications of this antibody to maintain a high neutralization activity against new variant of SARS-CoV-2 found in South Africa and Brazil.

Methods

Study design

The goal of this study was to develop a computational approach, which would allow to implement a virtual screening of antibodies against SARS-CoV-2 protein. This goal was achieved through the development of SARS-AB method that is based on MD simulations and ranking of antibodies based on the energy of antibody-S protein contacts. Two datasets were generated and analyzed in this study: (a) a dataset of 20 antibodies obtained from the protein data bank, (b) additional dataset of 30 experimentally acquired antibodies. The details of SARS-AB and dataset categories are described below.

Antibody selection

In order to test the ability of SARS-AB to distinguish SARS-CoV-2 binders from non-binders a total of 20 structural entries were obtained from the PDB. The accession numbers of all 20 structures are shown in the Table 1. The dataset contains 10 antibodies neutralizing SARS-CoV-2 protein and 10 anti-Flu antibodies, which were considered as a negative control.

For independent validation of SARS-AB we obtained 3 additional datasets:

First dataset contained 15 anti-SARS-CoV-2 antibodies, 11 of them were experimentally acquired antibodies without available crystal structures and 4 were recently deposited to the PDB (Table 2). To represent non-binders we randomly selected antibodies with different neutralization activities from the PDB. They included 2 antibodies effective against cytomegalovirus (pdb id 5C6T, 5VOB), 3 anti-meningitides antibodies (pdb id 5T5F, 2YPV, 2MPA), 4 anti-hepatitis C virus antibodies (pdb id 6meh, 4WHT, 5KZP, 5VXR) and 4 anti-HIV-1 antibodies (pdb id 1E6J, 5F9W, 1RZ8, 1RZG) (Table 2). Due to the highly specific interactions of these antibodies with their targets, selected antibodies are not expected to be cross-reactive with SARS-CoV-2.

A second dataset was extracted from a work of Li et al. [80]. The authors did single cell RNA and immune repertoire profiling of 16 COVID-19 patients at the time of hospital discharge. As a result of bioinformatics analysis authors identified 347 potential antigen-specific B cell receptors, from which 100 antibodies were chosen for the screening assay. Here, we selected top 70 of these 100 antibodies for our computational analysis.

A third dataset contained 37 humanized mouse antibodies from our in-house analysis.

Obtaining 3D model of antibody.

When available, the 3D model of antibody was downloaded from the PDB. In other cases antibody was built using SWISS-MODEL software by providing the full-length antibody sequence. SWISS-MODEL is a bioinformatics web-server dedicated to homology modeling of 3D protein structures [72]. In order to provide objective assessments of modelling performance, SWISS-MODEL participates in the CAMEO project (Continuous Automated Model Evaluation, https://cameo3d.org) [81]. Based on the CAMEO results in the ‘3D Structure Prediction’ category, SWISS-MODEL is consistently ranked among the top-modelling servers [72].

SARS-AB pipeline

Molecular dynamics simulations

The crystal structure of SARS-CoV-2 protein was obtained from Protein Data Bank (PDB code 7C01). In order to speed up MD only RBD domain of SARS-CoV-2 protein was used during simulations. The initial binding mode of SARS-CoV-2 with antibody was modeled similar to the interaction mode of SARS-CoV-2 with neutralizing antibody CB6 (PDB code 7C01). For that each built antibody was aligned with CB6-SARS-CoV-2 model from 7C01 in Chimera software [73]. The aligned structure was saved and minimized for 3000 steps using NAMD [82]. Next, structure was solvated in a rectangular box such that the minimum distance to the edge of the box was 10 Å under periodic boundary conditions in VMD [83]. Na and Cl ions were added to neutralize the protein charge, then further ions were added to mimic a salt solution concentration of 0.15 M. The system was minimized for 200 steps and gradually heated to 310 K. Each complex was equilibrated with NVT ensemble at 310 K for 1 ns. Further production run was performed for 10 ns with NPT ensemble using NAMD software [82]. A cutoff distance of 12 Å for Coulomb and van der Waals interactions was used. Long-range electrostatics were evaluated through the Particle Mesh Ewald method. The integration time step was 2 fs with all bonds involving hydrogen atoms constrained using the SHAKE algorithm [84]. Trajectory snapshots were saved every 50 000 steps for further analysis. UCSF Chimera [73] were employed to analyze and visualize the MD trajectories and to render the molecular graphics.

Ranking of antibodies using energy of antibody-SARS-CoV-2 contacts

A number of software tools have been proposed for predicting binding affinities of complexes [85], [86], [87]. These methods are based on the evaluation of energies calculated from protein structures and utilize different scoring functions. While several scoring functions exist, semi-empirical force field function provides a fast estimation of the interaction energy in biological systems. The most common protein force-field approaches include AMBER, CHARMM, GROMOS and OPLS [88]. All these approaches use similar form of the potential energy function for estimation of the interaction energy of non-bonded atoms.

Here, we estimated energy of antibody-SARS-CoV-2 contacts using the semi-empirical force field as implemented in AutoDockTools and LigEnergy [14], [74], and described in Huey et al. [75]. This scoring function is based on the AMBER force field and was calculated as following:

The first term represents the 12-6 Lennard-Jones potential for dispersion/repulsion interactions. The second term is a hydrogen-bond energy estimated by 10-12 potential, where E(t) is directionality of the hydrogen-bond interactions, which estimated according to Boobbyer et al. and Huey et al. Energy of electrostatic interactions is calculated based on the Coulomb potential. The last term is desolvation potential, which includes the volume (V) surrounding a given atom, weighted by the solvation parameter (S) and an exponential term based on the distance [75].

The values of energies of antibody-SARS-CoV-2 contacts were estimated for each studied antibody for the snapshots of the last 1 ns of MD simulations. It resulted in 10 values of energies for each antibody; their average value provided a single energy score. The ranking of studied antibodies was implemented based on the energy scores.

Estimation of energy scores with test dataset

The 3D models of anti-SARS-CoV-2 antibodies in complex with RBD S protein were obtained from the PDB. The energy scores for all 10 anti-SARS-CoV-2 antibodies were estimated according to Eq. (1). As antibody-protein structures were available we didn’t run MD simulations for these complexes.

For 10 anti-Flu antibodies we downloaded antibody’s 3D models from PDB and followed SARS-AB pipeline as described above: ran MD simulations and estimated energy scores.

Estimation of energy scores with validation datasets

For independent validation of SARS-AB we followed the full SARS-AB pipeline: we built models of all studied antibodies using SWISS-MODEL software, ran MD simulations and ranked antibodies by their energy scores.

Antibody expression and purification

The Regdanvimab and its analogs DNA sequences variable regions were synthesized in GenScript. The heavy and light chain variable regions were inserted into pfuse2ss-chig-hg1 and pfuse2ss-clig-hl2 (InvivoGen) vector accordingly. Antibodies were produced using ExpiCHO cells (Thermo Fisher). 8 days after transfection, antibodies were purified from medium using Protein A/G resin (Thermo Fisher). After dialysis in PBS, antibodies were used in neutralization assay. Antibodies’ concentrations were measured by BCA (Bicinchoninic Acid) Protein Assay (Thermo Fisher).

Pseudovirus generation

SARS-CoV-2 pseudovirus were generated following published paper [89]. HIV-1 pseudovirus coated with SARS-cov2 Spike protein were produced in 293T cells by co-transfecting pcDNA3.1-Spike with HIV-1 NL4-3 ΔEnv ΔVpr Luciferase Reporter Vector. Mutant pseudovirus were produced by mutagenesis on wild-type pcDNA3.1-Spike vectors. The mutation are as follows: UK Variant B.1.1.7 (501Y.V1) 6970del + N501Y; South Africa Variant B.1.351 (501Y.V2) K417N + E484K + N501Y; Brazil Variant P.1 (501Y.V3) K417T + E484K + N501Y.

Pseudovirus neutralization assay

The pseudovirus neutralization assays were performed using ACE2-Expressing Huh-7. Huh-7 cells (100 μL, 5 × 103 in DMEM) were added to 96 well-plate for overnight incubation. Various concentrations of mAbs (4-fold serial dilution with starting point at 30 μg/mL, 50 μL aliquots, triplicates) were mixed with the same volume of SARS-CoV-2 pseudovirus in a 96 well-plate. The mixture was incubated for 1 h at 37 °C, supplied with 5% CO2. No virus control wells were supplied with 100 μL DMEM (1% (v/v) antibiotics, 25 nM HEPES, 10% (v/v) FBS). Virus only control wells were supplied with 50 μL DMEM and 50 μL pseudovirus. After 1 h, medium was removed from Huh-7 cells, and then 100 μL pseudovirus and antibody mixture were added into Huh-7 cells containing plates. The 96-well plates were incubated for 48 h at 37 °C supplied with 5% CO2. After the incubation, supernatants were removed, and 100 μL Nano-Glo® Luciferase Assay Reagent (Promega) (1:1 diluted in PBS) was added to each well and incubated for 5 mins. After the incubation, the Luminescence was measured using CLARIOstar Plus Microplate Reader (BMG labtech). The relative luciferase unit (RLU) was calculated by normalizing Luminescence signal to virus only control group. IC50 were determined by a four-parameter nonlinear regression using GraphPad Prism 9.0 (GraphPad Software Inc.).

Results

Binding prediction of SARS-CoV-2 with monoclonal antibodies

We developed a computational approach based on MD simulations to estimate the binding energy of protein-antibody complex and predict interactions of S protein SARS-CoV-2 with monoclonal antibodies. The spike (S) protein is the major surface antigen of SARS-CoV-2 [22], which mediates viral entry into host cells by binding to the host receptor ACE2 through the RBD [4]. As of July 30, 2021, more than 3,000 SARS-CoV-2 monoclonal antibodies have been reported in 47 different studies [20], [23], [24], [25], [26], [27], [28], [29], [30], [31], [32], [33], [34], [35], [36], [37], [38], [39], [40], [41], [42], [43], [44], [45], [46], [47], [48], [49], [50], [51], [52], [53], [54], [55], [56], [57], [58], [59], [60], [61], [62], [63], [64], [65], [66], [67], [68], [69]. 753 antibodies were shown to bind RBD [20], [23], [24], [25], [26], [27], [28], [29], [30], [31], [32], [33], [34], [35], [36], [37], [38], [39], [40], [41], [42], [43], [44], [45], [46], [47], [48], [49], [50], [51], [52], [53], [54], [55], [56], [57], [58], [59], [60], [61], [62], [63], [64], [65], [66], [67], [68], [69], 427 of which had a measurable neutralization activity (Table S1, data taken from CoV-AbDab [70]) [20], [23], [24], [25], [26], [27], [28], [29], [30], [31], [32], [33], [34], [35], [36], [37], [38], [39], [40], [41], [42], [43], [44], [45], [46], [47], [48], [49], [50], [51], [52], [53], [54], [55], [56], [57], [58], [59], [60], [61], [62], [63], [64], [65], [66], [67], [68], [69]. We analyzed the V gene usage in these antibodies (Fig. S1) and found that IGHV3-53, IGHV3-66 and IGHV1-2 are the most dominant heavy chain genes, which is in line with previous studies [22], [71]. Currently, the structures of 43 human antibodies in complex with SARS-CoV-2 are available in PDB. Yuan et al. have shown that antibodies can be clustered based on the interactions with 3 different RBD binding sites: receptor binding site (RBS), CR3022 cryptic site and S309 proteoglycan site [71]. The antibodies that bind the RBS can be further sub-divided into sub-groups according to the angle of approach to RBD and their relative disposition on the surface of RBS [71]. The crystallographic studies of antibodies encoded by IGHV3-53 and IGHV3-66 revealed that they bind the same epitopes of RBS, while representatives of IGHV1-2 family bind in the close proximity. These antibodies represent 36% of all discovered by now antibodies with SARS-CoV-2 neutralization activity and will be used as template structures for our computational analysis. A structural and crystallographic analysis of CB6, a representative IGHV3-53 antibody (pdb id 7C01), revealed that CB6 recognizes an epitope that overlaps with ACE2 binding site in the RBD-SARS-CoV-2, and thereby interface with virus-receptor interactions by both steric hindrance and direct competition for interface residues [23]. In this work, we used 7C01 as a template structure to model antibody–S protein interactions.

The SARS-AB method is as follows (Fig. 1): First, we have built antibody’s 3D model using SWISS-MODEL software [72]. It is a fully automated protein structure homology-modeling server, which takes sequence of antibody as input and outputs its 3D model. The 3D model represents a structure of a molecule with x,y,z coordinates of the atom positions. Next, we performed a superposition of 3D model of antibody with the selected template structure using the Chimera software [73]. The purpose of superposition is to place antibody on the same binding site of RBD S protein as in the template structure. All files for the MD simulations were prepared as described in Methods. Simulation studies of 10 nano-seconds (ns) for all the selected antibody-RBD-SARS-CoV-2 complexes were performed. The energy of antibody–RBD contacts was estimated as implemented in AutoDock Tools [74], [75] and in the LigEnergy module [14]. The AutoDock free energy scoring function is based on the AMBER force field and was parameterized using a large number of protein-inhibitor complexes [74].

Fig. 1

The SARS-AB pipeline. The pipeline starts with building 3D model of an antibody using SWISS-MODEL software. Next, we do superposition of built antibody with a template structure (pdb 7C01) using Chimera software. Then, we prepare files for MD simulation, which includes minimization of antibody-RBD SARS-CoV-2 structure, solvation and ionization, heating and equilibration of structure in VMD and NAMD. Next, we run 10 ns of MD simulation in NAMD. The last step is to estimate energy of antibody-RBD contacts using AutoDock Tools.

Performance evaluation of SARS-AB on a dataset with crystal structures

To evaluate the ability of SARS-AB to distinguish true RBD binders from the non-binders, we generated a dataset containing 20 antibodies, with half known to neutralize SARS-CoV-2 and the other half effective against human influenza virus. Due to the highly specific interactions of these antibodies with their targets, the anti-flu antibodies are not expected to be cross-reactive with SARS-CoV-2, and they were considered as negative control. The PDB accession numbers of all the 20 structures are available in Table S2. We ran SARS-AB to estimate the binding energy of each of the 20 antibody-RBD complexes. The purpose of this step is to test if estimated energy, as a single predictor, is sufficient to distinguish true RBD binders from the non-binders. As expected, anti-SARS-CoV-2 antibodies showed higher affinity to S protein than anti-Flu antibodies. The 10 anti-SARS-CoV-2 antibodies had a negative energy between −30.59 and −22.25 kcal/mol while energy of anti-flu antibodies with S protein were unanimously higher than −15.5 kcal/mol (Table S2, Fig. 2A).

Fig. 2

SARS-AB ability to distinguish SARS-CoV-2 binders from non-binders. (A) Energy score distribution across 20 antibody-RBD-SARS-CoV-2 complexes displayed as violin plots. Binders correspond to 10 anti-SARS-CoV-2 antibodies, non-binders are 10 anti-Flu antibodies. (B) ROC curve measuring the performance of SARS-AB in distinguishing SARS-CoV-2 binders from non-binders for a dataset of 20 antibodies.

As expected, the antibody with the highest affinity to SARS-CoV-2 is CB6, which is the template structure used in SARS-AB. It was shown that this antibody specifically recognize RBD and block the binding of the SARS-CoV-2 RBD to ACE2. The antibody’s heavy chain dominates the interaction with RBD by all three complementarity-determining region (CDR) loops forming hydrophobic interactions and polar contacts, while light chain forms limited contacts with RBD (Fig. S2). The structural and energetic analysis implemented for anti-Flu antibodies showed that these antibodies fail to form strong interactions between the their CDR loops and S protein. These results indicated that the energy estimation in SARS-AB faithfully reflect the structural interactions. Therefore, we used the estimated energy to predict the binding status, with strong SARS-CoV-2 binders having interaction energy lower then −22 kcal/mol and non-binders having energy higher than −15.5 kcal/mol. Using binding energy as a single predictor, SARS-AB reached an AUC of 100% (Fig. 2B).

Validation of SARS-AB on an independent dataset without crystal structures

To evaluate if SARS-AB can predict RBD binding antibodies without knowing the structure, we created a validation dataset containing equal amount of binders and non-binders of SARS-CoV-2 from recent literature. Cao et al. reported S protein neutralizing antibodies identified by high-throughput single-cell RNA and VDJ sequencing of antigen-enriched B cells from 60 convalescent patients [28]. From this work we selected all the neutralizing antibodies with full-length sequences available. We also randomly selected antibodies with diverse reactivity from PDB database as negative control, but left out the structural information. In total, our dataset contained 15 binders and 15 non-binders of SARS-CoV-2 (Table S3). We implemented SARS-AB to estimate the interaction energy for all 30 antibodies with RBD, blind to the binder class labels. The distribution of estimated binding energy score showed near perfect separation between binders and non-binders (Fig. 3A). The cluster of SARS-CoV-2 binders has median value of −23.23 kcal/mol, whereas the non-binders has significantly higher median energy of −11.24 kcal/mol (p = 2.58⋅10−8). To date, in this dataset, BD-501 failed the downstream neuralization assay, as it had high KD and IC-50 values [28]. SARS-AB correctly assigned a higher binding energy and predicted it to be a non-binder.

Fig 3

Predictive performance of SARS-AB. (A) Energy score distribution of 30 antibodies displayed as violin plots. (B) ROC curve measuring the performance of SARS-AB on an independent dataset of 30 antibodies.

Further, SARS-AB predicted interaction energies are consistent with those estimated from the experimental data. All predicted strong binders have low values of interaction energy (<−22 kcal/mol) that correspond to their high binding affinity to the RBD (KD < 2 nM). All 15 non-binders, except BD-495, have much higher interaction energy (>−15.5 kcal/mol). Only BD-495 has lower value of interaction energy, −18.44 kcal/mol, that may be a result from its low RBD binding affinity, but still measurable neutralizing activity, to RBD SARS-CoV-2 (KD > 50 nM, IC50 = 18.1 μg/mL). These data strongly suggested that the energy estimation by SARS-AB truthfully reflect the experimentally measured binding affinity. Indeed, when using the binding energy as a single predictor, SARS-AB can separate anti-SARS-CoV-2 antibodies from non-binders with 99.6% AUC (Fig. 3B) for this independent test cohort.

Interestingly, several true-binders from this dataset differ by only a few amino acids [28], which allowed us to contrast the RBD binders and non-binders with high resolution. Specifically, we selected two antibodies – a true binder, BD-505, and non-binder, BD-495 (Fig. 4A,B, Table S3), and investigated their predicted docking structures. These two antibodies are encoded by the same VH, JH and VL genes with a difference in JL combination. The sequences of CDR3H loops are different by 4 amino acids and CDR3L sequences are almost identical (Fig. 4C). Our modeling suggests that the heavy chain of antibody BD-505 forms overall stronger H-bonds, VDW and hydrophobic interactions with RBD in contrast to BD-495. Close investigation over the contacts formed between CDR3H of antibodies with RBD revealed different interactions of the two antibodies. R99 of BD-505 formed H-bond with Q493 RBD,V100 and V101 create hydrophobic interactions with L455 (Fig. 4A). Unlike BD-505, BD-495 was unable to make H-bonds between atoms of CDR3H and RBD, and forms only hydrophobic interactions between I100 and L455, F456. The light chain of BD-495 forms weaker VDW interactions with RBD compared to BD-505 antibody. The BD-495 loose H-bond between D28 of antibody and Y505 RBD (Fig. 4B). This analysis identified critical residues on the RBD region for antibody recognition.

Fig. 4

The structures of BD-505-RBD-SARS-CoV-2 and BD-495-RBD-SARS-CoV-2. (A) H-bonds and hydrophobic interactions between CDR3H BD-505 and RBD-SARS-CoV-2 are shown. RBD is coloured blue, light chain of BD-505 is coloured green and heavy chain is yellow. (B) H-bonds and hydrophobic interactions between CDR3H BD-495 and RBD-SARS-CoV-2 are shown. (C) CDR3 sequence comparison and VDJ combination for BD-505 and BD-495 antibodies. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Binding prediction of wild type SARS-CoV-2, 501Y.V2 and P.1 variants with monoclonal antibody Regdanvimab (CT-P59) and its analogs

Some of the recently emerged variants of SARS-CoV-2, including 501Y.V2, P.1 and B.1.1.529, are more contagious and potentially more fatal [15]. Concerns have been raised that the new variants could escape antibody protection in immunized individuals, or the existing antibody drugs. Regdanvimab is a potent monoclonal antibody, which received conditional marketing authorization from the Korean Ministry of Food and Drug Safety [76], and is under discussion with U.S. FDA for emergency use authorization. It was recently announced [77] that antibody’s neutralizing activity was weak against 501Y.V2 and P.1 variants, and completely lost against Omicron [21]. Therefore, we next applied SARS-AB to investigate how it interacts with the mutant Spike protein of the 501Y.V2, P.1 and B.1.1.529 variants.

The crystal structure of Regdanvimab in complex with SARS-CoV-2 has been obtained from literature (PDB: 7CM4) [78]. The 501Y.V2/P.1 variants contain 3 amino acid changes from the original SARS-CoV-2 (N501Y, K417N/T and E484K). The latest severe acute respiratory syndrome coronavirus 2 variant Omicron (B.1.1.529) has more than 30 mutations in its spike (S) protein, 15 of which occur in the receptor binding domain (RBD). In order to evaluate how 501Y.V2, P.1 and B1.1.529 mutations affect the Regdanvimab binding, we ran SARS-AB to predict the interaction energy of Regdanvimab binding to either 501Y.V2, P.1, B1.1.529 or to the wild type (WT) SARS-CoV-2 complex from the PDB. Expectedly, the complex with WT SARS-CoV-2 showed stronger interaction energy with Regdanvimab, −22.4 kcal/mol, where for the mutants, the binding energy changed to −17.6 kcal/mol, decreased by 21.4% for 501Y.V2, −20.9 kcal/mol for P.1 and −10.1 kcal/mol for B1.1.529 variant. To understand what structural changes have caused this difference, we compared the contacts formed by antibody with WT SARS-CoV-2, the 501Y.V2 variant, P.1 strain and Omicron variant (Fig. 5A–D). The Regdanvimab forms H-bonds between R109 VH and E484 SARS-CoV-2, and between Y33 VL and E484 of WT SARS-CoV-2. The S32 of VH makes VDW contacts with K417 of WT S RBD (Fig. 5A). The K417N and E484K mutations in 501Y.V2 variant lead to the loss of H-bonds between R109 VH and K484 RBD, which increased the distance between chains to 10.3 Å. The mutations also weakened the interactions between S32 VH and N417 (distance is 5.4 Å, Fig. 5B). The E484K mutation in P.1 strain also abolish the formation of hydrogen bonds between residues K484 and R109, while K484 makes VDW interactions with Y111 (distance is 4.6 Å, Fig. 5C), which helps to stabilize the complex and results in a lower binding energy compared to the 501Y.V2 variant. Based on the above results, we concluded that the loss of two strong hydrogen bonds and weakened VDW interactions induced by E484K and K417N/T can significantly decrease the binding of Regdanvimab with 501Y.V2 and P.1 variants. The loss of two hydrogen bonds between E484A, K417N and Regdanvimab as well as the loss of contacts with Q493R and Q498R of S protein (Fig. 5D) allows Omicron strain to completely escape from neutralization by antibody.

Fig. 5

Comparison of contacts formed by Regdanvimab with SARS-CoV-2 and 501Y.V2 SARS-CoV-2. (A) Structure of Regdanvimab with WT SARS-CoV-2. (B) Structure of Regdanvimab with 501Y.V2 SARS-CoV-2. (C) Structure of Regdanvimab with P.1 SARS-CoV-2. (D) Structure of Regdanvimab with Omicron variant. (E) Structure of Regdanvimab analog (S32E, R109E) with P.1 SARS-CoV-2. (F) Structure of Regdanvimab analog (K55G, D56S, D57G, G102T, R109L, Y111F) with Omicron. RBD is coloured blue, light chain of antibody is coloured green and heavy chain is yellow. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Next we sought to identify minimum modification(s) to the original sequence of Regdanvimab that could restore its interaction with 501Y.V2, P.1 and B.1.1.529. We created a number of Regdanvimab analogs by mutating amino acids on the CDR loops of heavy chain, and investigated their binding affinity using SARS-AB pipeline in order to find antibody with improved affinity to 501Y.V2, P.1 and B.1.1.529 SARS-CoV-2. We found that S32E and R109E mutations in VH of antibody (mut 1) can restore the interaction energy (-21.6 kcal/mol and −21.7 kcal/mol) between antibody and the new variants (501Y.V2 and P.1). This improved energy is a result of additional stabilization of the antibody-S RBD P.1 complex by a formation of H-bonds between E32 VH and T417 RBD (distance is 2.6 Å), and between E109 VH and K484 RBD (distance is 2.7 Å, Fig. 5E). The same binding pose was predicted for a complex of antibody with 501Y.V2 variant. Thus, the proposed modification of Regdanvimab could potentially make it effective against the 501Y.V2 and P.1 variants of SARS-CoV-2 and re-establish antibody protection against the virus. For a number of tested antibodies against Omicron, SARS-AB predicted the Regdanvimab analog (K55G, D56S, D57G, G102T, R109L, Y111F) to be a B1.1.529 binder with interaction energy −16.1 kcal/mol. The antibody forms 2 hydrogen bonds between Y60, S56 of H chain and R493 of S protein, one hydrogen bond between T102 and G485, and additionally stabilized by hydrophobic interactions between F111 and F490 (Fig. 5F). The proposed modifications should restore an antiviral activity of Regdanvimab analog against B.1.1.529 variant.

Next, we aimed to confirm SARS-AB predictions in-vitro. We purified the proteins of the original Regdanvimab and all 5 analogs, and performed neutralization assays against WT SARS-CoV-2, 501Y.V2, P.1, B1.1.7 and D614G pseudo viruses (Fig. 6). The Regdanvimab showed slight decrease of neutralization of P.1 variant (IC50 = 0.149 μg/mL) compared to WT SARS-CoV-2 (IC50 = 0.086 μg/mL). As predicted by SARS-AB, the Regdanvimab showed decreased neutralization against 501Y.V2 variant (IC50 = 0.243 μg/mL) compared to WT and P.1 strains. The Regdanvimab analogs (mut 1 and mut 5, Fig. 6) containing two mutations (32E, 109E and 32Q, 109T) in the VH, showed the strongest neutralization among all studied analogs. These Regdanvimab analogs were predicted to be the strongest binders for 501Y.V2 and P.1 strains by SARS-AB. Their neutralization was increased 3–6-folds, 2.2–2.5 folds and 1.8–2-folds for mut 1and mut 5 against WT, P.1 and 501Y.V2 SARS-CoV-2 variants compared to the original antibody (Fig. 6).

Fig. 6

The pseudovirus neutralization assay. (A) The inhibition curve for Regdanvimab and its analogs in different pseudovirus types. Pseudovirus carry luciferase reporter gene and its titer can be determined with luminescence. The relative luciferase unit (RLU) were measured at 48 h after infection. The neutralization activity of antibodies can be determined by RLU fold change. The x-axis shows the antibody concentrations, and the y-axis shows %RLU. (B) IC50 for Regdanvimab and its analogs are shown in the tables. IC50 values are in ug/mL. Five different analogs have been studied, containing 2 to 4 mutations in heavy chain compared to the original Regdanvimab antibody. Mutants/analogs have the following mutations: mut 1 (32E, 109E), mut 2 (32D, 109E), mut 3 (32E, 56E, 104D, 109E), mut 4 (32Q, 109L), mut 5 (32Q, 109T).

Computational efficiency of SARS-AB

The average computational time to build a 3D model of one antibody with SWISS-MODEL webserver was around 20 min. The MD simulations were run on NVIDIA Tesla P100 GPU card with 56 CPUs on UT Southwestern BioHPC cluster. The 10 ns of MD simulation in NAMD for each antibody-RBD-SARS-CoV-2 complex with 100 ps step took ∼5 h. The analysis of MD trajectories and estimation of energy of antibody-RBD contacts took another 2 h. The abilities of BioHPC cluster allowed to run 4 MD simulations jobs in parallel. As a result, SARS-AB could process 8 antibodies in 10 h, which takes 5 days to analyze 100 antibodies. This performance can be further improved if large cloud-based computing, such as AWS, is enabled.

Discussion

To improve the efficiency of discovering new neutralizing antibodies, we designed SARS-AB, which allows to predict binding of antibodies to RBD-SARS-CoV-2. We tested the method on its ability to distinguish SARS-CoV-2 binders from non-binders using a dataset of antibodies derived from the PDB, and validated SARS-AB on independent set of antibodies with known specificity.

It is known that the S1 subunit of spike protein contains two major structural elements, RBD and N-terminal domain (NTD). While majority of discovered antibodies are directed against RBD, some were shown to have a potent neutralization activity against NTD [20], [43]. The latter does not directly compete with the ACE2 binding site and it remains unclear how it blocks SARS-CoV-2 infection. Currently, SARS-AB approach was designed to predict binding of antibodies only with RBD, thus it is possible that predicted RBD-non-binders may interact with NTD or other undiscovered areas on SARS-CoV-2. That being said, the methodology of SARS-AB is readily transferrable to cover NTD binding prediction, given a proper template structure and sufficient training data.

Another limitation of SARS-AB, as with most other binding prediction software, is the need to use a template structure to indicate the docking position, specifically, the location of the binding site. Consequently, the binding energy estimation of SARS-AB is limited to the same docking position as the template structure. Here we applied the structure with a binding site (RBS-A) that is dominant in all currently known neutralizing antibodies against RBD region (CB6). This approach ensures high specificity, but may leave out true binders taking other docking positions. In principle, we could include other known binding sites to SARS-AB, but it will significantly increase computational burden, which may result in similar time cost as in vitro experiments. Future accelerations on SARS-AB would be required to incorporate more binding sites and improve prediction sensitivity.

As screening large antibody databases with the help of MD simulations requires huge computational power, the current design of SARS-AB represent a good balance of prediction accuracy and computational efficiency. The more accurate binding affinity predictions will require longer MD simulations and the use of time consuming molecular mechanics methods. Currently, many antibody screening studies experimentally profiled 102 to 103 antibodies predicted from large-scale sequencing data. Based on the published results, the fractions of the true neutralizing antibodies is lower than 5%, even with flow sorting [28]. Given its high specificity, SARS-AB can be applied to test out most of the non-binders, leaving a handful of remaining prioritized targets for downstream validation. We believe the future application of SARS-AB to large antibody sequence datasets will significantly accelerate the discovery of ultra-high-affinity neutralizing antibodies against SARS-CoV-2.

It was shown that South African and Brazilian strains may escape from neutralization by first-wave monoclonal antibodies and antibodies being developed for clinical use [19]. The virus engages N501Y, K417N/T and E484K mutations, which helps to escape neutralization by monoclonal antibodies directed at the RBD. It is interesting that virus might sacrifice its hACE2 binding affinity through K417N to survive the attack of neutralizing antibodies [79]. In this work we showed that Regdanvimab undergoes 1.7–2.8-fold decrease in neutralization activity against virus new variants found in South Africa and Brazil compared to WT SARS-CoV-2. By applying SARS-AB, we identified minimum modifications to the original sequence of Regdanvimab that restored its interaction with 501Y.V2 and P.1 strains. The SARS-AB can be applied to most antibody screening studies to identify promising therapeutic candidates against SARS-CoV-2.

Finally, as SARS-CoV-2 is a fast-evolving strain of coronavirus, any developed antibody drugs are at risk of rapid immune evasion caused by mutations at contact residues. Despite the above limitations, we demonstrated that SARS-AB can be used to measure the binding energy changes caused by a few mutations at affordable time complexity, and designed novel antibody analogs with improved neutralization activity against SARS-CoV-2 variants. We believe that our approach will be used as a design platform to evaluate interactions of candidate antibodies with new variants of SARS-CoV-2.

Author statement

BL conceived the project. DB developed SARS-AB pipeline, performed data analysis and wrote the manuscript. YF contributed to the writing of the manuscript. YF, MD, YS, and FD implemented pseudovirus neutralization assays. ZJC and BL co-supervised the project and contributed to the manuscript.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.