key: cord-0744766-pgr1i90r
authors: Martí, Didac; Alsina, Marc; Alemán, Carlos; Bertran, Oscar; Turon, Pau; Torras, Juan
title: Unravelling the molecular interactions between the SARS-CoV-2 RBD spike protein and various specific monoclonal antibodies
date: 2021-10-25
journal: Biochimie
DOI: 10.1016/j.biochi.2021.10.013
sha: 7719300d005c607fc6b33a019a0a55e2b8f65e6e
doc_id: 744766
cord_uid: pgr1i90r

Vaccination against SARS-CoV-2 just started in most of the countries. However, the development of specific vaccines against SARS-CoV-2 is not the only approach to control the virus and monoclonal antibodies (mAbs) start to merit special attention as a therapeutic option to treat COVID-19 disease. Here, the main conformations and interactions between the receptor-binding domain (RBD) of spike glycoprotein of SARS-CoV-2 (S protein) with two mAbs (CR3022 and S309) and the ACE2 cell receptor are studied as the main representatives of three different epitopes on the RBD of S protein. The combined approach of 1 μs accelerated molecular dynamics (aMD) and ab-initio hybrid molecular dynamics is used to identify the most predominant interactions under physiological conditions. Results allow to determine the main receptor-binding mapping, hydrogen bonding network and salt bridges in the most populated antigen-antibody interface conformations. The deep knowledge on the protein-protein interactions involving mAbs and ACE2 receptor with the spike glycoprotein of SARS-CoV-2 increases background knowledge to speed up the development of new vaccines and therapeutic drugs.

Late 2019, an outbreak of pneumonia (coronavirus disease 2019, was detected in Wuhan (China) caused by an unknown etiology that is believed to be related to a pathogen coming from a wildlife reservoir that was able to migrate to human host [1] . This pathogen was identified as a new coronavirus named SARS-CoV-2 (Severe Acute Respiratory Syndrome coronavirus 2) [2] . The new virus has spread very quickly, recently through more infective mutations, stressing health systems around the world to the limit and causing a tremendous impact in terms of deaths and economic losses.

Vaccination against SARS-CoV-2 have started in most of the countries with 40.1% of the world population which has received at least one dose after 5.41 billion doses administered globally (September 3rd, 2021) [3] , showing a markedly uneven trend among the vaccinated population. However, the development of specific vaccines against SARS-CoV-2 is not the only approach to control the virus and other strategies are still needed considering the potential for mutation of the coronavirus. Monoclonal antibodies (mAbs) merit special attention as a therapeutic option against COVID-19 [4] . Extensive work is being done to explore and find the best mAbs to minimize the short term impact of the pandemic and complement the prophylactic vaccination plan [4, 5] . Therefore, the understanding of how SARS-CoV-2 interact with mAbs is still of outmost importance.

The densely glycosylated spike protein of SARS-CoV-2 (hereafter, S protein) has attracted the attention of researchers as it plays a crucial role in virus entry to a cell and mediates its transmission and infectivity [6] . Recently, the crystallographic structure of S glycoprotein has been resolved [7, 8] . Spike glycoprotein is made up of three protomers, each one of them with two subunits (i.e., S1 and S2). S1 subunit is the most interesting for this work as contains an N-terminal subunit responsible for virus-receptor binding. It is mainly structured in two domains, the N-terminal domain (NTD) and the receptor-binding domain (RBD). The latter is the responsible to interact with the angiotensin-converting enzyme 2 (ACE2), the cell receptor [9e11] . Interestingly, S1-RBD/ACE2 interaction triggers conformational changes in the S2 subunit that facilitate the virus fusion with the cell membrane and it allows the penetration of the virus inside the cell [12, 13] . Consequently, the research and selection of mAbs that are able to block or distort such interactions merit special attention as potential treatments against the infection.

The initial research about feasible mAbs was focused on those that were already known to be able to bind SARS-CoV and MERS-CoV, considering the relatively high similarity of such virus with SARS-CoV-2 in terms of RBD. As a result of the sampling, the CR3022 and S309 mAbs were identified as potential therapeutic candidates [14] . However, others (i.e., m336, m396, CR3014) failed in targeting SARS-CoV-2 spike implying that dissimilarities among coronavirus impair the cross-reactivity of some neutralizing antibodies [15] . Moreover, CR3022 and S309 together with ACE2 have been identified as representatives of three well differentiated interaction zones with RBD [16] .

CR3022 epitope does not overlap with the ACE2 binding site when binding to SARS-CoV-2 RBD as does not show any competition with ACE2 for such site [15, 17] . Controversial results have been reported on the CR3022 mAb in vivo effectiveness. More specifically, Yuan et al. [18] , among others [19e22] , reported that CR3022 failed to neutralize SARS-CoV-2. Structural modeling demonstrated that the binding epitope could only be targeted by CR3022 when at least two RBDs on the trimeric S protein were slightly rotated and in the "open/up" conformational state. The cryptic epitope targeted by CR3022 is deemed important to achieve a high probability against resistance mutations in this region as it stabilizes the "closed" state of the prefusion trimer. However, although the in vitro assays showed that CR3022 does not neutralize SARS-CoV-2, the possibility that this epitope could confer in vivo protection was accepted. More recently Huo et al. [23] reported that CR3022 has neutralizing activity against SARS-CoV-2 with IC50 of~0.114 mg/mL in vitro and its role as a neutralizing agent in combination with other mAbs is still under research.

On the other hand, S309 mAb was identified from memory B cell screening using peripheral blood mononuclear cells of a SARS-CoV infected patient in 2003 but collected in 2013 [14, 24] . S309 targets the S glycoprotein of SARS-CoV and potently neutralizes SARS-CoV-2 by recognizing an epitope which is containing a fucosylated glycan at position N343 [14] . This glycan seems to play an important role on the interaction RBD/S309. It is worth noting that S309 mAb shows ability to bind in both "up" and "down" RBD conformations and it does not compete with ACE2 binding. Therefore, S309 binds a highly conserved epitope shared by the SARS-CoV and SARS-CoV-2, reducing the likelihood of mutational escape. Currently, an engineered antibody (VIR-7831) based on S309 is under assessment in clinical trials (i.e. NCT04545060) [25, 26] .

Until October 2020, more than 800 SARSCoV-2-targeting mAbs, obtained from COVID-19 patients, have been disclosed [16] . These mAbs were grouped considering their interaction with three different RBD binding sites: i) mAbs targeting the receptor-binding motif (RBM) of ACE2, ii) mAbs targeting CR3022 cryptic site, which is the most frequent epitope targeted by cross-neutralizing antibodies (i.e. COVA1-16, H014, EY6A and ADI-56046), and iii) mAbs targeting S309 binding site. It is worth to noting that mAbs that are able to bind the CR3022 and S309 sites, usually cross-reacted with the RBM of SARS-CoV, as compared to the ones that bind to the RBM of ACE2, which show high specificity for SARS-CoV-2 [16] .

In this work the main interactions between the RBD of spike glycoprotein of SARS-CoV-2 and the main systems that recognize the three main RBD epitopes described by Yuan et al. [16] (i.e., the mAbs CR3022 and S309) and the human ACE2 receptor, have been studied (see Fig. 1 ). Several theoretical works on the RBD interface can be found in the bibliography, mainly addressing S/ACE2 interactions, although the theoretical works with the CR3022 and S309 mAbs are scarce [27e31] . Most of these works present a classic approach to the interactions either from classical molecular dynamics (cMD) point of view or using docking methodologies. This work aims to carry out an exhaustive sampling on the different conformations of the protein-protein interface under physiological conditions using accelerated molecular dynamics (aMD), thus analyzing the most favored poses under these conditions by means of its two-dimensional potential of mean force (2D-PMF) considering the binding energy and buried surface interface, to determine the predominant interactions throughout the conformational sampling. Conformational sampling can be achieved using cMD and aMD as it is shown in the bibliography [32, 33] . However aMD simulations present a greater efficiency when compared to MD simulation as it was shown by Pierce et al. where 500 ns of aMD can be equivalent to about 1 million ns at cMD time scale [34] . Finally, the most energetically favorable set of poses is subjected to relaxation by means of the quantum mechanics/molecular mechanics molecular dynamics (QM/MM MD) approach, which is taking into account the environment effect to obtain a better description of the physical-chemical interactions. The proposed methodology allows to determine the main interactions in antibody-antigen complexes in a more realistic way in comparison with the conventional docking and/or molecular dynamics techniques [35] .

The objective of developing a cross-coronavirus vaccine is of interest for the future especially as it remains active and ready to mutate and other coronavirus strains remain a potential pandemic thread. For that reasons the deep knowledge on the protein-protein interactions involving mAbs and ACE2 receptor with the spike glycoprotein of SARS-CoV-2 increases the odds to be better prepared for the next novel coronavirus outbreak, giving relevant background to speed the development of new vaccines and therapeutic drugs.

Three initial structures were built considering the complex between the receptor-binding-domain (RBD) of SARS-CoV-2 spike (S) glycoprotein and some potential anchorage proteins. Specifically, two mAbs (CR3022 and S309) and the Angiotensin Converting Enzyme 2 (ACE2 receptor) were considered (hereafter S/CR3022, S/ S309 and S/ACE2, respectively)(see Fig. 1 ). All initial protein structures were obtained from the Protein Data Bank: (a) The S/CR3022 complex was build using a crystal structure with a resolution of 1.95 Å. (PDB ID code 6YLA) [23] . (b) The S/S309 complex was taken from the crystal structure with a resolution of 3.10 Å (PDB ID code 6WPT) [14] . (c) The S/ACE2 complex was built from a crystal structure with a resolution of 2.45 Å (PDB ID code 6M0J) [10] . S glycoprotein is extensively covered with N-linked glycans which is thought to be related for a proper folding [36] and controlling the accessibility to neutralizing antibodies and proteases [37, 38] . The RBD's crystallographic structure of the SARS-CoV-2 S glycoprotein presents some glycans that were also included and modeled in the final complex systems.

Missing protein residues were incorporated using the Modeller algorithm [39] as implemented in the UCSF Chimera 1.14 program [40] and the Z-DOPE (Discrete Optimized Protein Energy), which is an statistical potential based for the choice of best model (i.e. that with the lowest Z-DOPE) for each conformational state, among the generated ones. All PDB files were edited and processed before the simulations adjusting the protonation states to pH 7 and adding the hydrogen atoms to generate Amber topology files and coordinates files. All glycan molecules present at the crystallographic structure were preserved in the final simulation structure. The protein system was placed in a orthorhombic box of water molecules represented by the TIP3P model [41] , containing about 20 Å of water buffer surrounding the protein. The solvated protein was subsequently neutralized and filled with a concentration of 0.15 M of NaCl to represent a more typical biological environment. Accordingly, the models used to represent the three protein interfaces, i.e., CR3022, S309 and ACE2, with the RBD of SARS-CoV-2 S glycoprotein presented 286389, 203384 and 291368 explicit atoms, respectively.

All classical simulations were performed using the AMBER 18 simulation package [42] . AmberTools Leap program was used to set up input files for MD simulations with Amber. Thus, initial system parameter files were generated using the Amber ff14SB force field for the protein atoms [43] in combination with the Glycam06 force field for the glycan molecules at the S protein [44] and TIP3P parameters for water [41] whilst the parameter for the solvated free ions were obtained from the Li, Merz and co-workers parameterization [45] . Before starting the system equilibration, an initial relaxation on the protein regions that were filled with the UCSF Chimera program, was performed. The protein coordinates were restrained to the crystal coordinates, except for the new added residues, which were relaxed by means of the Limited-memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) quasi-Newton algorithm [46] as implemented in the AMBER package [42] . Next, the system was minimized with all the protein atoms restrained to the crystal coordinates to remove close contacts, and the restrained system was gradually heated up to 298 K using an NVT ensemble along 60 ps? Covalent bonds involving hydrogen atoms were constrained using the SHAKE algorithm [47] , and long-range electrostatic interactions were treated with particle-mesh Ewald using a realspace cutoff of 10 Å [48] . The protein restraints were released following the next 1 ns of simulations using a NPT ensemble at 1 atm and 298 K up to constant density. An integration time step of 2 fs was used in all MD simulations. Finally, the simulation was extended for 250 ns using a NVT ensemble at 298 K with coordinates recorded every 40 ps?

The accelerated molecular dynamics (aMD) technique [49] allows an enhanced-sampling on the protein conformational space by artificially reducing energy barriers separating different states of a system. Therefore, a faster and more extensive exploration of new configurations is easily achieved compared to cMD simulations.

Specifically, aMD is based on applying a modification to the potential energy of the system, such that at the moment when the potential energy of the system is below a defined threshold (E) a new summand modifying the potential energy will be added [49] .

where DV(r) is the boost potential, and is responsible for modifying the potential profile being calculated as follows: 

where a is a parameter that determines the extent to which we want to modify the potential profile. Values of this parameter close to or equal to 1 will modify the potential profile in such a way that it becomes completely flat, disabling the molecular dynamics; while extremely large values will change the profile very little, to the point of having no influence on the potential energy. The parameters E and a were calculated from the averaged dihedral energy obtained at the end of 100 ns of cMD simulations, following the methodology described by Miao et al. [50] The final parameter pairs (E; a) in kcal/mol used were (10746.60; 520.0), (6857.23; 340.0) and (13528.60; 638.8) for the S/CR3022, S/S309 and S/ACE2 systems, respectively. aMD simulations were performed along 50 ns using an NVT ensemble and starting from 20 snapshots. Initial structures were obtained from the last 150 ns of cMD at each system equally spaced and with a fairly large separation interval to avoid correlation.

Accordingly, a total of 1 ms of aMD trajectories were produced for each spike-antibody complex (50 ns  20 snapshots). All aMD productions were carried out with the NAMD molecular simulation package [51] . In all cases, the simulations were run using the same force field, parameters, and protocols than previous cMD using Amber in order to be comparable among them. The Langevin dynamics method [52] with a relaxation time for the coupling of 5 ps À1 for the temperature, was used. Similarly to the previous cMD simulation using a NVT ensemble, the rest of the simulation parameters were taken accordingly.

The potential of mean force (PMF) were obtained from reweighting aMD simulations to recover the original free energy profiles of studied protein complexes with the RBD of SARS-CoV-2 spike glycoprotein. Energetic reweighting was conducted by using a cumulant expansion to the second order, which is able to recover the most accurate free energy profiles within statistical errors [53] . After studying several variables in 1D-and 2D-PMF profiles to obtain the best representative of protein-protein interface interactions, two different variables were considered to build the two-dimensional PMF profiles:

i. The binding energy between the protein (antibody or ACE2) and the RBD of S protein (BE).

This value was calculated using the MMPBSA approach [54] as implemented within AMBER package [42] .

ii. The interface buried surface between the protein (antibody or ACE2) and the RBD of S protein (BS). The BS is the area that is masked by the solvation effect between the two proteins. It is calculated as follows:

where SASA refers to the accessible solvent surface area (i.e. the area of the molecule that is exposed to the solvent. SASA values are determined using the Connolly algorithm [55] as implemented in the cpptraj software of AMBER package [42] ). Finally, 2D-PMF profiles (i.e. PMF ¼ f(BE, BS)) in all the studied protein complexes were obtained.

Thirty snapshots were obtained from the global minimum and adjacent areas (within a PMF radius of 0.1 kcal/mol) of twodimensional PMF plots for each subsystem. These conformations will serve as starting points for ab initio QM/MM MD simulations to relax the interface among the RBD of protein S and either the studied antibodies or the ACE2 entry receptor. Indeed, the system is divided in two regions, one involving the main residues that are closely linked within the interface, which are treated at quantum level (QM region), and the rest of the system that is treated classically (MM region) under the QM/MM MD approach. This methodology has been previously used in biomolecular systems that a further description on the chemical interactions are needed, such as enzymatic processes and large metalloproteins, among others [56e58]. More specifically, the QM region was derived from a contact analysis observing those atoms pairs lying at both sides of the interface with a distance equal to or less than 3 Å and with a probability of occurrence equal to or greater than 90%. The QM regions of three studied systems contained a total amount of 288, 199, and 279 quantum atoms located at the interface between the RBD of SARS-CoV-2 S glycoprotein and the proteins CR3022, S309 and ACE2 entry receptor, respectively. The interface with CR3022 was modeled using two different QM regions following the maz-QM/MM MD approach [59] to reduce the simulation time, which contained 231 and 57 atoms, respectively, whereas the other two models only involved one QM region to model the interface.

QM/MM MD simulations were conducted by means of the AMBER-PUPIL-Gaussian tandem [60] , where the PUPIL program [61, 62] interfaces the classical engine (AMBER 18 [42] ) with the quantum engine (Gaussian 16 [63] ) calculating and managing all the QM/MM coupling within the QM/MM MD approach. All starting conformations were initially relaxed along 0.25 ps of QM/MM MD trajectory. An NVT ensemble at 298 K with a temperature regulation through Langevin dynamics (collision frequency of 10 ps À1 ) and a time step of 0.5 fs, were used. However, the rest of simulation parameters and force fields used for the QM/MM MD simulations were the same as those used for cMD protocol. Finally, a QM/MM MD production trajectory of 1.0 ps was obtained. The coordinates in the last 0.75 ps at any trajectory were saved for subsequent analysis. In total, 30 trajectories of 1.25 ps (2500 steps) were produced for each one of the starting structures of each system, thus obtaining a total of 37.5 ps of QM/MM MD production for each antibody and the ACE2 receptor. The MM region was parameterized using ffSB14 Amber force field [43] whilst the QM region was modeled by means of Density Functional Theory (DFT) using the M06e2X functional [64] in combination with a minimum basis set (STO-3G) for all quantum atoms except the sulphur atom, which was holding an additional polarization orbital (STO-3G*). It should be noticed that medium-range (i.e. 5 Å) non-covalent interactions, such as conventional (e.g. NeH / O and HeO/H) and non-conventional (e.g. CeH/O) hydrogen bonds, are better described with M06e2X functional than with usual DFT functionals [65] . Long-range electrostatics on the QM regions were considered by means of Ewald summations with a real-space cutoff of 25 Å around each QM region for the electrostatic interactions within the QM/MM coupling methodology, whilst the reciprocal-space cutoff D. Martí, M. Alsina, C. Alem an et al.

Biochimie xxx (xxxx) xxx value was set to 6 (maximum integer translation of the reciprocal lattice).

3.1. RMSD/RMSF analysis of protein-protein complex Figure S1 shows the Root Mean Square Distance (RMSD) along 250 ns of cMD of three protein complexes with the RBD region of spike glycoprotein of SARS-CoV-2. The variable domain on the Fab regions of the CR3022 (S/CR3022) and S309 (S/S309) antibodies, and the human ACE2 entry receptor protein (S/ACE2), were considered. The three studied systems were equilibrated after 100 ns of production. In fact, the average RMSD values in the last 150 ns of production take values of 3.3 ± 0.5 Å, 4.2 ± 0.3 Å and 2.9 ± 0.3 Å for the S/CR3022, S/S309 and S/ACE2 systems, respectively. In all cases the standard deviation is below 0.5 Å along all the fraction of production trajectory observed for statistical studies. Thus, the stability of the three examined protein-protein complex conformations shows that S/S309 system is a little less stable than the others two, which present a similar instability when compared with the crystallographic coordinates. Protein flexibility has attracted the attention of the scientific community due to the importance of the intrinsically disordered proteins, mostly associated with essential biological functions such as regulation and signaling [66] . However, flexibility is also important for their influence on the protein complexes formation as it facilitates the formation of asymmetric interfaces among different protein domains [67] . Root mean square fluctuation (RMSF) measures the amplitude of atom motions during cMD simulation and thus, the system flexibility. Figure S2 illustrates the three RMSF computed considering all atoms fluctuations of each residue on the three protein complexes studied along 250 ns of cMD, starting from crystallographic coordinates. A larger fluctuation is observed on the RBD of S protein compared with those observed in the Fab domain of antibodies and ACE2 protein. Besides of RBD protein ends, the sections with the highest residue fluctuation are similar on both complexes with mAbs (residues 475e500 of RBD), whereas the S/ ACE2 system presents a different region of greater fluctuation (residues 370e385 of the RBD).

Moreover, in order to study the flexibility of those residues involved in the interface between the two proteins, an analysis of the contacts between pairs of amino acid residues of each side of the interface was performed. Thus, the amino acids which presents C a atoms with a proximity less than or equal to 10 Å with another amino acid located on the other side of the interface and a population higher than 90% along the last 100 ns of the cMD trajectory, were obtained. Figure S2 shows the areas belonging to the detected protein-protein interface highlighted with a red ellipse. As expected, a lower fluctuation of the interface atoms is obtained compared to the rest of the protein. However, in order to better study the fluctuations at the interface, Fig. 2 shows the RMSF only for those amino acids with greater proximity at the protein-protein interface. Most of the residues present fluctuations lower than 2 Å. However, S/309 system registers a higher mobility at the interface showing 6 residues with a fluctuation higher than 2 Å (N334, R346, K444, V445, G446 and G447) compared with S/CR3022 and S/ACE2 systems that present only 3 (N370, R408, G519) and 2 (S19, C325) residues with fluctuations >2 Å, respectively. Again, S/309 complex shows a higher flexibility not just in the whole protein structure but also on the protein-protein interface.

Conformational scanning on the protein complexes with the RBD of SARS-CoV-2 spike glycoprotein were conducted by means of accelerated molecular dynamics (aMD) approach [49] . This technique allows a quick and extended biological conformation sampling that remains inaccessible to standard cMD due to the necessary time scale required to address such a large systems using cMD approach. Indeed, millisecond-timescale events in both globular and membrane proteins can be captured in just a hundreds-of-nanosecond aMD simulation [34, 68] . Fig. 3 shows the two-dimensional potential of mean force (2D-PMF) of S/CR3022, S/ S309 and S/ACE2 protein-protein complexes, obtained from energetic reweighting along 250 ns of aMD simulations to recover original free energy profiles. Binding energy (BE) and buried surface area (BS) variables were chosen to represent the two dimensional free energy profiles, (see methods section).

It is known that aMD leads to a large conformational rearrangement on the protein-protein interfaces, which might help to explore a large number of potential poses. The large interface rearrangement due to aMD simulation on current systems under simulation can be observed on the large dispersion of BE and BS The comparison of the 2D-PMF profiles shows the S/ACE2 system with the lower conformational variability, which is holding a reduced 2D-PMF surface dispersity. The S/S309 system presents a greater conformational dispersion both for the global system and for the protein-protein interface, which coincides with the greater flexibility observed above. However, despite the greater flexibility, the absolute free energy minimum on the 2D-PMF S/S309 profile is the one closest to the PMF coordinates of the crystallographic structure (green circle, Fig. 3b ). Taking into account that both conformations, i.e., the crystallographic structure and the absolute minimum structure of the 2D-PMF profile can be represented as a single point in the 3D surface of Fig. 3 , the Euclidean distance between projected 3D-points over the plane (BE,BS) of Fig. 3 is obtained. This comparison leads to Euclidian distances of 10.0, 11.8 and 26.3 on the (BE, BS) plane for the S/S309, S/ACE2 and S/CR3022 systems, respectively. Table 1 lists the binding energy and the buried surface among the three protein-protein interfaces studied but considering different stage of system relaxation under physiological conditions and starting from crystal coordinates. Interestingly, the inclusion of solvation under physiological conditions increases the energetic interactions between RBD domain and CR3022 antibody, decreases in ACE2 and remains very similar in the S/S309 system with respect to those observed in the crystalline system (Table 1) . Indeed, energy affinity of protein-protein complexes obtained after a conformational sampling with aMD are, in decreasing order, S/CR3022 > S/ S309 > S/ACE2. However, it is experimentally observed that the order of affinity would correspond to S/S309 (K D ¼ 0.001 nM [14] ) > S/ACE2 (K D ¼ 1. [23] ). This is because in the S/CR3022 system, despite having the most favorable BE, the RBD protein glycosylation on the region of the epitope recognized by the CR3022 antibody is almost completely shielded in the open conformation and not accessible at all in the closed conformation [38, 70] . Since the RBD must have an "open/up" state in order for the antibody-antigen complex to form (see Figure S3 ) [23] , the protein glycosylation might reduce the likelihood of this happening, in agreement with the lower observed experimental affinity. Moreover, it is observed that as the binding energy decreases, the buried surface of the complexes increases, thus favoring a greater energetic interaction between the amino acids buried in the interface. Fig. 4 shows the contact map of the conformations around the global minimum of the 2D-PMF profiles obtained above (Fig. 3) . All distances between alpha carbons of the amino acids found on both sides of the interface that have a population of more than 90% have been sampled. The contact map was prepared using cpptraj program as implemented in AMBER [42] and results were processed by an in-house program to color map the information as pleased. As can be seen, the S/CR3022 system has a greater number of shorter contacts that extend along all six complementarity-determining regions (CDRs), i.e., from H1 up to L3, similarly to the interactions observed in the crystal structure determined by Yuan et al. [18] A detailed description on the amino acid sequence of different CDRs can be found at Table S1 of supplementary information. The positions of the residues belonging to the CDRs are labeled according to the Kabat numbering scheme and obtained from SabDab structural antibody database [71, 72] . S/S309 complex shows a lower number of contacts than S/CR3022, which are mainly localized on the CDR-H3 region and in less extension on the light chain (CDR-L1, L2, and L3). Moreover, S/ACE2 complex is showing the lowest number of contacts in this work, mainly localized on the 50 first residues as well as those close to the residue 350. This reduction on the contact number (S/CR309 > S/S309 > S/ACE2) was expected considering the order on the binding energy previously obtained in this work.

Close inspection to the residues that present a greater number of contacts with distances less or equal to 4 Å among the most populated contacts on the conformations around the global minimum of 2D-PMF profile ( Figure S4 ) was done. Thus, the mAb CR3022 paratope shows the highest number of contacts on residues W33 (CDR-H1), Y52 (CDR-H2), S100 (CDR-H3), I34 (CDR-L1) and in Table 1 Binding energy (BE, Kcal/mol) and buried surface area of interface (BS, Å 2 ) among different protein-protein complexes derived from crystallographic coordinates (cryst), along last 150 ns of cMD simulation (cMD), and averaged among all the conformations close to the absolute aMD minimum (FBE ¼ f(BE,BS) < 0.1 kcal/mol). Standard error is also shown. less extension on the Y38 (CDR-L2) residue, leaving the L3 region out of the most populated interactions with the RBD of protein S. Moreover, epitope of CR3022 is clearly identified attending Figure S4a . Indeed, most of the most stable conformations in physiological conditions show the most populated interactions on three regions, i.e., 369e392, 428e430, and 515e518 residues. Those results are in agreement with the experimental crystallographic structure of S/CR3022 system [18, 23, 69] . However, observing a highest sampling of conformational structures derived from aMD in physiological conditions, a mapping of the most conserved residue interactions along the protein-protein interface can be shown. Specifically, the residues from the 377 up to 381 constitute the epitope sub-region that presents the highest number of close contacts. Also, residues K378, T430, and L517 are those more populated within each one of the previously localized epitope regions, respectively.

The paratope of mAb S309 shows a greater number of contacts in the residues along the CDR-H3 region, such as Y100 to a greater extent, and residues W105, F106, E108, L110 and I111 to a lesser extent ( Figure S4b) . Also, some contacts are detected involving residues from CDR-L1 and L2 regions but with a lower weight within the most populated interactions with the RBD of protein S. Interestingly, the residue Y50 between L1 and L2 regions is accumulating an important number of contacts with the RBD protein.

On the other hand, the epitope of mAb S309 presents one main region of interactions along residues 334e345, which also involves N343 glycosylation site. Indeed, E340 residue and the a-L-fucopyranose component of the oligosaccharide at N343 position are the residues that present the highest number of contacts from the epitope region, being the residue K356 the one with the highest number of conserved contacts on those latter epitope regions along conformational space studied. All observed regions are similar to those reported from crystallographic data [14] , albeit aMD data allows to consider those interactions that seems to be conserved along the conformational study.

Crystallography observations of Lan et al. [10] on the interface between the human ACE2 receptor and the RBD region of the S protein show about 17 residues of RBD protein interacting with 20 residues from the ACE2 (cut-off distance of 4 Å), where most of these interactions appear in the N-terminal helix of ACE2. The present study considers an important set of solvated conformations around the absolute minimum of 2D-PMF profile, where two important contact region are shown, i.e., on the N-terminal helix of ACE2 (residues 24e42) and the second one starting from residue 353 up to 357 ( Figure S4c) . Interestingly, the two residues with the highest number of contacts in each of these two regions are Y41 and K353 residues, which are very close (d Ca-Ca~7 .5 Å). Indeed, the mentioned residues are interacting with the same region of the RBD (residues 498e505) which supports a greater number of interactions with ACE2 than the N-terminal.

Once the most stable conformations under physiological conditions of the protein-protein complexes studied were obtained, a subset of 30 different conformations were selected (see methodology section) and were allowed to relax under the QM/MM-MD approach. Thus, the residues along the interface with the largest number of prevalent contacts between the conformations around the minimum of the 2D-PMF profile were selected (cut-off distance of 3 Å and population !90%) to be considered as the quantum region of the simulation (QM region). The rest of the system was modeled under a classical molecular mechanics approach (MM region) to take into account the influence of the environment in a more realistic way. This subsequent treatment of the more stable conformations will allow us to see not only those interactions of the interface that might persist under the application of thermal energy, but also we will obtain more realistic interaction distances between residues on both sides of the interface. In fact, M06e2X density functional used in this simulation is taking into account not only electrostatic interactions but also van der Waals interactions and short-range dispersion forces that are not considered under previous cMD simulation [65] .

The different interactions at the crystal structure have been described elsewhere [10, 14, 23] . However, a comparison was made among the different approaches carried out in this work for a deep insight on the main interactions of hydrogen bond network and salt bridges that will persist under physiological conditions at room temperature. Tables S2, S3, and S4 (supplementary information) summarize the major interactions of hydrogen bonds (HB) and salt bridges (SB) between the mAbs CR3022 and S309, and the human ACE2 receptor with the S protein of SARS-CoV-2 (RBD domain), respectively. The comparison includes five sets of different conformations: First, the interactions collected from the literature for crystalline structures (I xRay ); Second, the interactions obtained along the cMD trajectory (I cMD ); Third, the main interactions of the conformational set around the absolute minimum of the 2D-PMF profile (I raMD ); Fourth, the interactions at the absolute minimum snapshot of the 2D-PMF profile (I maMD ); Fifth, the main interactions obtained from the QM/MM MD trajectories of a selected conformational set around the minimum 2D-PMF profile is relaxed (I QMMM ) . In all cases, only those interactions that are conserved throughout the conformational set with a population of 60% or higher are considered.

Comparing the set of crystallographic interactions (I xRay ) with those of the absolute minimum of the 2D-PMF profiles (I maMD ), the latter usually present a larger number of interactions due to the exhaustive conformational scanning that takes place around the starting crystallographic structure. However, the coincidence of interactions is low, and despite of new interactions, most of them retain the major residues involved but with different heavy atoms. This conformational study does not incorporate the complete glycan chains of the S protein, but it allows us to see the predominant interactions at the protein-protein interface studied. We would like to point out that when comparing the set of I raMD interactions around the absolute minimum, the cardinality of this set is much lower, hinting only at those predominant interactions for each system and remaining throughout most conformations. The case of the S/S309 system deserves special attention, where the number of persistent interactions in I raMD is the lowest of all (one single SB), which agrees with the above observed higher flexibility of the RBD in the S/S309 complex. Flexible RBD amino acids are those that are close to or directly involved on the I maMD interactions set (i.e. N334, R346, K444, V445, G446 and G447). Thus, the introduction of thermal energy into the simulation will cause the number of perdurable interactions to decrease in all simulated systems due to system flexibility and physiological conditions (i.e. a lower cardinality of I cMD and I QMMM compared with the I xRay and I maMD , as observed in Tables S2-S4).

The interactions between CR3022 mAb and RBD of S protein have been described with an important hydrophobic component [18, 23] . However, HB and SB type interactions are also quite important (Fig. 5) and would be interesting to study to what extend are affected by the flexibility of the system when introducing temperature and solvation effects. The epitope residues involved in the most important interactions described by Yuan et al. [18] (F377, K378, C379, G381, P384, K386, and T430) are also found in the set of residues holding the most populated interactions among the most stable conformations of the S/CR3022 complex (I raMD , Figure S4a ), whilst among paratope residues only Y52 (within the variable domain of heavy chain, V H ) was detected within the most important set of interactions of the protein-protein complex. In fact, considering only HB and SB interactions, F377(N):Y52(OH) HB interaction is the only one that coincides in both experimental work of Yuan et al. [18] and Huo et al. [23] by crystallography. This is in agreement with the main interactions found in this work, where this interaction is one of the strongest HB (d QMMM ¼ 2.60 Å) and lasting throughout all the conformations studied at the QM/MM MD level (Table S2) 

Similarly to the S/CR3022 interface, the S/S309 system also presents an increased set of interactions in the absolute energy minimum of the 2D-PMF profile (I maMD ) when compared to those reported at the experimental level (I xRay ) but with a high coincidence on the S309 paratope residues involved. Indeed, I maMD presents a significant number of HB and SB interface interactions. However, the interactions found in the I maMD set are very labile due to the great flexibility of both the glycan chain and some of the amino acids within the S309 paratope. In fact, in the conformational study around the absolute energy minimum, no interaction was detected that occurred in more than 60% of the conformational set (Table S3) . However, by relaxing I raMD under physiological conditions at room temperature and at the ab initio level using the QM/MM MD approach, 4 HBs and 3 SBs were recovered, which persist in most of the conformations obtained. It should be noted the strongest interactions that are maintained throughout the I QMMM conformational set are three HBs involving three residues located on the CDR-H3 and interacting with two residues in the Nterminal RBD domain (Fig. 5b) (Table S3 ).

An important point to be considered in this interface is the role that glycan residues play in the interaction with the variable region (Fv) of S309 mAb. Residue N343 is highly fucosylated (98% of glycans bear a fucose residue) [38] and the glycan attached interacts mostly with CDR-H3 (residue D115) and CDR-L2 (residues Y50, S53 and R55, see Fig. 5b and Table S3 ) because it is deposited on the pocket formed between both CDRs, as observed in experimental crystallographic structures [14] . Interestingly, a significant HB interaction with the CDR-H1 (S53(O):MAN1324(O3) d QMMM ¼ 2.74 Å) was also detected, thus showing a glycan chain sandwiched between H1 and L2, and laying along H3 and L2 (Fig. 5 and Table S3 ).

The most significant contributions to the HB and SB network were observed between the N-terminal helix of ACE2 (Q24, D30, K35, E37, D38 and Y41) and the RBD binding loop (residues 449 to 505) of the SARS-CoV-2 (Fig. 5c ). In fact, the residue range 487e505 is the one that concentrate most of the RBD contacts ( Figure S4c ). However, when introducing temperature and relaxing the system with cMD only the SB interaction D30(OD2):K417(NZ) remains bound, maintaining a more sporadic behavior the rest of HB along cMD trajectories (Table S4) . A more detailed study on the most probable S/ACE2 conformations using aMD approach, it was observed that at the energy minimum of 2D-PMF profile (I maMD ) an important fraction of the interactions previously observed at the experimental level [10] were recovered. However respectively. The latter HB is not observed at the experimental level [10] but remains throughout all simulation sets under physiological conditions. In addition, R393(NH2):Y505(OH) salt bridge interaction observed both at the crystallographic level and on the minimum of 2D-PMF profile (I maMD ) seems not to be so important at physiological conditions and is lost after relaxing at QM/MM MD level (I QMMM ) at room temperature.

The study of the BE decomposition allows to better understand the main contributions to the stability of each protein-protein complex, that is, the role played by each residue and the importance of the different SB, HB and hydrophobic interactions. Thus, the BE decompositions between the absolute minimum of the 2D-PMF profiles (ImaMD) with the set of minimum energy conformations that is around them (IraMD) are compared. Fig. 6 shows the main RBD residues that contribute to the interaction with mAbs and the human ACE2 receptor, as well as the surface area of energy contribution to the BE for each complex.

On the other hand, Tables S5, S6, and S7 show the main contributions to the binding free energy on a per-residue basis for SARS-CoV-2 RBD with CR3022, S309 and ACE2, respectively, using MMPBSA methodology [54] as implemented in AMBER [42] . Similarly, Tables S8, S9, and S10 show the main interaction energies among residue pairs located at both sides of the interface for each system, and derived through BE decomposition on a pairwise perresidue basis. The comparison between the residues that contribute the most to the absolute minimum I maMD and those of their surrounding region can give us another point of view of those that have a greater weight in the most probable conformational evolution in a physiological environment, at room temperature. In fact, I34(L1), Y52(H2), I102(H3), I30(heavy chain) residues of the CR3022 paratope, ordered from highest to lowest energy contribution and with values higher than 4 kcal/mol, are the most contribute to the binding energy of the S/CR3022 complex. Interestingly, the residue E57(H2) loses much prominence compared to that shown in the conformation of the absolute minimum of energy (I maMD ). Regarding the epitope of the S/CR3022 system, the K378 residue is the one that presents a greater energy contribution by far, mainly due to four important interactions, one of them being a strong ionic interaction (SB) with the CR3022 (Table S8) . However, residues L517, K386, Y380, F377 and C379 are also important with contributions >4 kcal/mol. Also, within the set of hydrophobic interactions, the Y306 and T430 residues gain considerable prominence in the I maMD conformational set compared with low profile exhibited in I raMD .

The BE decomposition of S/S309 system shows an important contribution from the H3 paratope region mostly, and to a lesser extent from L1 and L2 CDR regions. This agrees with what was previously observed when studying the contacts between the interface residues, as well as the main and most persistent HB and SB interactions along conformational study. Comparing the absolute minimum with the conformational set around I raMD , a greater reinforcement of the interactions with H3, L1 and L2 CDR regions is observed, while the energy contribution of residues located on the H2 and L3 regions decreases. More specifically, residues W105; F106 and I111 of the CDR-H3 are those that present a greater energetic contribution (|DBE|>4 kcal/mol) while the best contributions of L1 and L2 present values of |DBE|<2 kcal/mol to the BE. On the other hand, the residues of the RBD epitope that contribute the most are T345, A344, P337, N343, R346, T333 and V445 (|DBE| >2 kcal/mol). Interestingly, most of these residues are located around the glycosylated residue N343. In fact, glycans represent an important weight contribution to BE, showing significant interactions with the S309 mAb, and thus, favoring the interface stabilization. Glycans are mainly interacting with the CDR-H3 (D115 and Y100), CDR-H1 (Y32), and with CDR-L2 (R55, S54, S53 and Y50), all of them with a pairwise per-residue basis contributions of |DBE| >2 kcal/mol, (see Table S9 ). It is worth mentioning that the E340 residue, which presents a high number of pairwise per-residue basis interactions, its averaged energy contribution throughout the entire I raMD set is not very high (|DBE| ¼ À1.02 kcal/mol). This is because in the different averaged conformations, the histogram of energy contribution (not shown) is fairly balanced between the attractive and repulsive contributions to the BE. In comparison, residue T345 (|DBE| ¼ À3.85 kcal/mol) presents an energy balance with clearly attractive contributions, although when decomposing in pairwise per-residue it does not present a large number of energetically favorable interactions.

From an energetic point of view, the main residues of the S/ACE2 complex are pretty similar with what has been seen above on the study by contacts and by its HB/SB network at the protein-protein complex interface. Thus, we want to highlight D355, Y41, K353, Q24, T27, S19 and D38 residues with energetic values of |DBE | > 2 kcal/mol. All of them, holding attractive energy contributions but S19 residue, which contributes to BE with repulsive energies. D355 shows the highest contribution, and as seen above, presents the strongest HB within the I raMD conformational set, which is enhanced by relaxing the conformational set using a higher calculation level (I QMMM ). Something similar happens with K353 residue. On the other hand, besides of those residues belonging to the b-sheet/turn region of human ACE2 receptor (residues 353e357), the largest set of residues that contribute energetically to the BE is observed in the N-terminal helix of ACE2. Among them, we would like to highlight D30 residue, which despite not being within the most energetically contributing residues in the I raMD set, it holds the most energetic predominant interaction D30(ACE2):K417(RBD) which corresponds to a SB (Tables S4 and  S10) . Similarly to what was observed with residues K355 and K353, the SB observed by this residue is also reinforced when using a higher calculation level (I QMMM ). Regarding the RBD region of the spike glycoprotein of SARS-CoV-2, its region binding motif (RBM) was reported to be in the residue region 338e509 [10] . However, observing Table S7 most of the energy contributions are located in the region of residues 486e505, in agreement to what was previously observed in this work studying contacts and HB network. In fact, 10 (67%) of the 15 residues that contribute the most energy to BE belong to this region, meanwhile 14 (93%) of these residues are located in the 449e505 range (Table S7 ). This can be explained by the observed well steric-coupling between 449 and 505 amino acids region of the RBD with the N-terminal helix of human ACE2 receptor, so favoring this region to have a high number of polar and hydrophobic interactions with ACE2.

After the outbreak of COVID-19 numerous studies appeared concerning either new or known monoclonal antibodies (mAbs) with a potential affinity to bind to the region binding domain (RBD) of the spike protein of SARS-CoV-2. However, the molecular mechanism of recognition and the structural origin of binding affinities between the amino acids belonging to the epitope and paratope of RBD-mAb complexes are under-explored. In this work, we have studied the complex interface between the RBD region of SARS-CoV-2 with two mAbs, the CR3022 and S309 that showed different binding affinity with the S glycoprotein. The exploration of RBD-mAb interface was compared with the interface of RBD complexed with the angiotensin-converting enzyme 2 (ACE2), which is known as the entry point of SARS-CoV-2 to the human host cells, using the same methodology.

A combined in silico strategy consisting of accelerated molecular dynamics (aMD) and Quantum Mechanics-Molecular Mechanics Molecular dynamics (QM/MM-MD) simulation was applied to elucidate the structural origin of binding affinities of RBD-mAb complexes of SARS-CoV-2. A broad conformational space of protein-protein interface through all-atoms accelerated MD simulations was explored. This methodology can be associated to a kind of flexible docking, which combined with two-dimensional PMF profile involving interface descriptors such as the binding energy and buried surface, allowed us to wisely achieve a selected set of poses surrounding the absolute minimum of 2D-PMF profile but taking into account the flexibility of the protein-protein interface. A later structural relaxation through high level ab-initio QM/MM-MD approach leaded to a much improved binding interaction geometries, thus elucidating the most prominent interaction residue-pairs under physiological conditions at room temperature. The obtained affinity of protein-protein complexes are, in decreasing order, S/CR3022 > S/S309 > S/ACE2. The S/CR3022 binding affinity was overestimated since no protein glycosylation on the region of the epitope recognized by the CR3022 antibody was considered. However an interactions mapping on the S/CR3022 interface was conducted and the most persistent and energetic interactions of the paratope and epitope were determined. S/S309 complex interface is the most flexible due mainly to the epitope flexibility and the N433 glycosylation involved in the proteinprotein interface, which is playing a very important role in stabilizing it, mainly due to existence of pairwise interactions involving CDR-H3 and L2 regions of S309 mAb. S/ACE2 complex was compared with previous experimental and theoretical results showing good agreement thus validating the current used methodology. Besides a strong salt bridge (SB) involving K417, most of the region binding motif (RBM) are located in the 449e505 range of residues of SARS-CoV-2 RBD when is attached to human ACE2 receptor.

Overall, the three studied systems of SARS-CoV-2 RBD bound with CR3022, S309 and ACE2 proteins, and the key binding residues and interactions identified in this work provide new insights into understanding mechanisms of SARS-CoV-2 infection of host cells and on the development of novel immunoassays based in antibodymediated immobilization, which could facilitate to develop new treatments to facing the ongoing COVID-19 pandemic and potential future coronavirus outbreaks. 

The authors declare no conflict of interest either financial, personal or other.

Group of the International Committee on Taxonomy of, the species Severe acute respiratory syndrome-related coronavirus: classifying 2019-nCoV and naming it SARS-CoV-2

Coronavirus Pandemic (COVID-19), OurWorldInData.org

An update to monoclonal antibody as therapeutic option against COVID-19

Rapid isolation and profiling of a diverse panel of human monoclonal antibodies targeting the SARS-CoV-2 spike protein

Neutralizing antibodies against SARS-CoV-2 and other human coronaviruses

Cryo-EM structure of the 2019-nCoV spike in the prefusion conformation

Structure, function, and antigenicity of the SARS-CoV-2 spike glycoprotein

Structure of the SARS-CoV-2 spike receptor-binding domain bound to the ACE2 receptor

SARS-CoV-2 cell entry depends on ACE2 and TMPRSS2 and is blocked by a clinically proven protease inhibitor

The spike protein of SARS-CoV d a target for vaccine and therapeutic development

MERS-CoV spike protein: a key target for antivirals

Crossneutralization of SARS-CoV-2 by a human monoclonal SARS-CoV antibody

Potent binding of 2019 novel coronavirus spike protein by a SARS coronavirus-specific human monoclonal antibody

Recognition of the SARS-CoV-2 receptor binding domain by neutralizing antibodies

A human monoclonal antibody blocking SARS-CoV-2 infection

A highly conserved cryptic epitope in the receptor binding domains of SARS-CoV-2 and SARS-CoV

A human neutralizing antibody targets the receptor-binding site of SARS-CoV-2

Structural basis for the neutralization of SARS-CoV-2 by an antibody from a convalescent patient

Neutralization of SARS-CoV-2 by destruction of the prefusion spike

Receptor-binding domain-specific human neutralizing monoclonal antibodies against SARS-CoV and SARS-CoV-2

Anti-SARS-CoV-2 neutralizing monoclonal antibodies: clinical pipeline

Enhanced receptor binding of SARS-CoV-2 through networks of hydrogen-bonding and hydrophobic interactions

Is the rigidity of SARS-CoV-2 spike receptor-binding motif the hallmark for its enhanced infectivity? Insights from all-atom simulations

Why does the novel coronavirus spike protein interact so strongly with the human ACE2? A thermodynamic answer

Silico analysis of interaction between full-length SARS-CoV2 S protein with human Ace2 receptor: modelling, docking, MD simulation

On the interactions of the receptor-binding domain of SARS-CoV-1 and SARS-CoV-2 spike proteins with monoclonal antibodies and the receptor ACE2

Improved modeling of peptideprotein binding through global docking and accelerated molecular dynamics simulations

Using Accelerated Molecular Dynamics Simulation to elucidate the effects of the T198F mutation on the molecular flexibility of the West Nile virus envelope protein

Routine access to millisecond time scale events with accelerated molecular dynamics

The latest automated docking technologies for novel drug discovery

The viral spike protein is not involved in the polarized sorting of coronaviruses in epithelial cells

Unexpected receptor functional mimicry elucidates activation of coronavirus fusion

Site-specific glycan analysis of the SARS-CoV-2 spike

Comparative protein structure modeling using MODELLER

UCSF Chimerada visualization system for exploratory research and analysis

Comparison of simple potential functions for simulating liquid water

Simmerling, ff14SB: improving the accuracy of protein side chain and backbone parameters from ff99SB

GLYCAM06: a generalizable biomolecular force field, Carbohydrates

Systematic parameterization of monovalent ions employing the nonbonded model

On the limited memory BFGS method for large scale optimization

Settle: an analytical version of the SHAKE and RATTLE algorithm for rigid water models

Towards an accurate representation of electrostatics in classical force fields: efficient implementation of multipolar interactions in biomolecular simulations

Accelerated molecular dynamics: a promising and efficient simulation method for biomolecules

Accelerated molecular dynamics simulations of protein folding

Scalable molecular dynamics with NAMD

Langevin stabilization of molecular dynamics

Improved reweighting of accelerated molecular dynamics simulations for free energy calculation

MMPBSA.py: an efficient program for end-state free energy calculations

Analytical molecular surface calculation

QM/MM molecular dynamics studies of metal binding proteins

Mechanistic insights into the reaction of chlorination of tryptophan catalyzed by tryptophan 7-halogenase

Massive quantum regions for simulations on bionanomaterials: synthetic ferritin nanocages

Multiple active zones in hybrid QM/MM molecular dynamics simulations for large biomolecular systems

A versatile AMBER-Gaussian QM/MM interface through PUPIL

Software integration in multi-scale simulations: the PUPIL system

Chapter one -PUPIL: a software integration system for multi-scale QM/MM-MD simulations and its application to biomolecular systems

The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals

Assessment of the performance of the M05À2X and M06À2X exchange-correlation functionals for noncovalent interactions in biomolecules

Predicting intrinsic disorder in proteins: an overview

Protein flexibility facilitates quaternary structure assembly and evolution

Mapping of allosteric druggable sites in activation-associated conformers of the M2 muscarinic receptor

A Cryptic Site of Vulnerability on the Receptor Binding Domain of the SARS-CoV-2 Spike Glycoprotein

Beyond shielding: the roles of glycans in the SARS-CoV-2 spike protein

SAbDab: the structural antibody database

Thera-SAbDab: the therapeutic structural antibody database

Supplementary data related to this article can be found at https://doi.org/10.1016/j.biochi.2021.10.013.