key: cord-0714067-82x242oc
authors: Yang, Jingyi; Zhang, Zhao; Yang, Fengyuan; Zhang, Haiwei; Wu, Haibo; Zhu, Feng; Xue, Weiwei
title: Computational design and modeling of nanobodies toward SARS‐CoV‐2 receptor binding domain
date: 2021-05-13
journal: Chem Biol Drug Des
DOI: 10.1111/cbdd.13847
sha: 26322b25f3d7ec67ea8d0733685fc91a467583a7
doc_id: 714067
cord_uid: 82x242oc

The ongoing pandemic of coronavirus disease 2019 (COVID‐19) caused by severe acute respiratory syndrome coronavirus 2 (SARS‐CoV‐2) has become a global health concern and pose a serious threat to humanity. There is an urgent need for developing therapeutic drugs and (or) biologics to prevent the spread of the virus. The life cycle of SARS‐CoV‐2 shows that the virus enters host cells by first binding to angiotensin‐converting enzyme 2 (ACE2) through its spike protein receptor‐binding domain (RBD). Therefore, blocking the binding between of ACE2 and SARS‐CoV‐2 RBD can inhibit the virus infection in the host cells. In this study, by grafting the complementarity‐determining regions (CDRs) of developed SARS‐CoV, MERS‐CoVs specific neutralizing antibodies (nAbs) include monoclonal antibodies (mAbs) as well as SARS‐CoV‐2 mAbs onto a known stable nanobody (Nb) scaffold, and a total of 16 Nbs sequences were designed. Five Nbs, namely CS01, CS02, CS03, CS10, and CS16, were selected based on the free energy landscape of protein docking verified by the recently reported Nb‐RBD cocrystal structures. CS01, CS02, and CS03 occupied the ACE2 binding site of RBD, while CS10 and CS16 were proposed to inhibit the interaction between RBD and ACE2 through an allosteric mechanism. Based on the structures of the five Nbs in complex with RBD, seven brand‐new Nbs with enhanced binding affinities (CS02_RD01, CS03_RD01, CS03_RD02, CS03_RD03, CS03_RD04, CS16_RD01, and CS16_RD02) were generated by redesign of residues on the interface of the five Nbs contact with SARS‐CoV‐2 RBD. In addition, the identified “hot spots” on the interface of each complex provide useful information to understand the binding mechanism of designed Nbs to SARS‐CoV‐2 RBD. In sum, the predicted stabilities and high binding affinities of the 11 (re)designed Nbs indicating the potential of the developed computational framework in this work to design effective agents to block the infection of SARS‐CoV‐2.

The ongoing outbreak of novel coronavirus disease 2019 caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection has led to 136,291,755 confirmed cases and 2,941,128 deaths in 223 countries, areas or territories as of April 14th, 2021 (https://www.who.int/emerg encie s/disea ses/novel -coron aviru s-2019). The pandemic of the disease has become a global health concern and pose a serious threat to humanity . To meet this challenge, great efforts have been paid for the development of therapeutic approaches against SARS-CoV-2 . Due to the rapid spread of COVID-19, there is an urgent need to develop highly potent and broad-spectrum and cost-effective anti-SARS-CoV-2 drugs and vaccines.

Like the other two CoVs, SARS-CoV and Middle East respiratory syndrome CoV (MERS-CoV), SARS-CoV-2 belongs to the beta-CoV genera in the family Coronaviridae (Zhou, Yang, et al., 2020) . Life cycle of CoV shows that the spike protein plays an essential role in viral attachment, fusion, entry, and transmission . During infection, CoV enters host cells by first binding to their respective cellular receptors angiotensin-converting enzyme 2 (ACE2) through the virus receptor-binding domain (RBD) of the spike protein . Therefore, development of neutralizing antibodies (nAbs) blocking the interaction between RBD and ACE2 plays crucial roles in inhibiting the infection of pathogenic CoVs in the host cells .

Nanobodies (Nbs) are single-domain antibodies (sdAbs) derived from camelids and sharks showing a large sequence identity with the human VH gene family III (Muyldermans, 2013) . The small size (~15 kDa), thermostability, high binding specificity, and low immunogenicity of Nbs making them suitable for many biotechnology and medicine applications (Muyldermans, 2013) . To fight viruses and prevent their spread, Nbs can interfere at different levels of the viral replication cycle (Steeland et al., 2016) . Currently, nanobody maturation technology was deployed to develop several Nbs targeting SARS-CoV-2 spike protein . Crystal structure showed that these molecules block the interaction between RBD and ACE2 and their neutralizing activity against SARS-CoV-2 suggested that Nbs may serve as useful therapeutics during CoVs outbreaks (Wrapp et al., 2020) .

The present work focused on design of Nbs specifically binding to SARS-CoV-2 RBD by using an integrated computational approach. Based on the developed SARS-CoV, premier released SARS-CoV-2 and MERS-CoVs specific nAbs include monoclonal antibodies (mAbs; Chen et al., 2017; Hwang et al., 2006; Li et al., 2015; Pak et al., 2009; Prabakaran et al., 2006; Walls et al., 2019; Wang et al., 2015 Wang et al., , 2018 Ying et al., 2015; Yu et al., 2015; Yuan et al., 2020; Zhang et al., 2018; Zhou et al., 2019) , their functional antigen-binding fragment (Fab), and the single-chain variable region fragment (scFv), and a series of novel Nbs were first designed by grafting the complementarity-determining regions (CDRs) of nAbs targeting the two CoVs heavy chain onto a known nanobody framework (Kang et al., 2019) . Then, the structures of the designed Nbs were predicted by comparative modeling, and the possibility of Nbs-RBD complexes formation was assessed through protein-protein docking followed by molecular dynamics (MD) simulations. Considering the structural changes of SARS-CoV-2 RBD during infection, multiple conformations of RBDs binding to ACE2 were extracted from the crystal structure or microsecond level MD simulation of a chimeric RBD in complex with ACE2 . Computational redesign of protein-protein interaction affinity and specificity was further applied to generate brand-new Nbs. Finally, the results suggest that eleven of the (re)designed Nbs may bind the SARS-CoV-2 RBD and potentially neutralize the virus.

There are 16 nAbs binding to SARS-CoV, MERS-CoV, or SARS-CoV-2 RBD which the atomic coordinates deposited in RCSB PDB database (Berman et al., 2000) at the time of analysis are collected in Table S1 . To design Nbs targeting SARS-CoV-2 RBD, only the RBD of three CoVs in complex with the "single-domain antibodies (sdAb)" (referred to the heavy chain of the corresponding nAbs) were saved for all of the 16 structures. It is known that each sdAb (Nb) has a highly conserved framework with three CDRs of variable sequence composing the paratope (Wilton et al., 2018) . Therefore, the 16 prepared "sdAb" sequences originated from the heavy chain of SARS-CoV, MERS-CoV, or SARS-CoV-2 nAbs were further divided into seven parts, including four fragment regions (FRs) and three CDRs parts ( Figure S1a ). Then, computational sequence design of 16 Nbs (CS01 to CS16) was performed in two steps. First, choosing the four FRs of a camelid Nbs ( Figure S1b ) as a scaffold (Hamers-Casterman et al., 1993) and grafting the three CDRs of each prepared "sdAb" onto the scaffold (Wagner et al., 2018;  Table S2 ). Second, to make sure the correct folding of designed Nbs, several amino acids of the camelid Nbs, especially the residues located on the edge of the four FRs, were replaced by the sequence of prepared "sdAb".

The structures of designed Nbs were predicted by comparative modeling method. The BLAST algorithm (Schaffer et al., 2001) in NCBI was used to search templates. For each Nb, three different structures were selected as templates (Table S3) for modeling using default parameter in Modeller (v.9.24; Sali & Blundell, 1993) and 10 models were predicted. The model for designed Nb was selected by picking the structure with the best DOPE assessment score considering the Lennard-Jones potential and GBSA implicit solvent interaction (Shen & Sali, 2006 ).

The docking of designed Nbs to RBDs was performed using the protein docking module of Rosetta (v.2019.31; Alford et al., 2017) . All of the structures calculated by Rosetta were scored base on the Rosetta Energy Function2015 (REF2015; Alford et al., 2017; Chaudhury et al., 2011) . The coordinates of Nb and RBD structures were first formed through the script of clean_pdb.py and refined by running the Rosetta relax protocol (Nivon et al., 2013) . To guide protein-protein docking, the PyMOL software (The PyMOL Molecular Graphics System, Version 2.0 Schrödinger, LLC) was applied to construct the structures of Nbs contacted with RBDs based on the 16 prepared crystal structures of "sdAb" in complex with SARS-CoV, SARS-CoV-2, or MERS-CoV. The details of the constructed Nbs-RBD complexes are listed in Figure S2 . The RBDs in different conformational states were extracted from the crystal structure or microsecond level MD simulation of a chimeric RBD in complex with ACE2 . To ensure low-energy starting side-chain conformations for docking, prepacking of the constructed Nbs-RBD complexes was conducted by Rosetta prepack protocol . For each complex, the prepacked low-energy conformation was used as a starting point for several rounds protein-protein docking to generate 1,000 or 10,000 decoys by running Rosetta docking protocol (Gray et al., 2003) with the Monte Carlo (MC) refinement method (Gray et al., 2003) . For each docking trajectory, the RMSD was calculated from the heavy atoms of the interface residues using each pose of the top five scorers as a reference structure (Chaudhury et al., 2011) using Rosetta InterfaceAnalyzerMover (Fleishman et al., 2011) . The docking funnel was then identified through plotting I_score against I_rmsd. The top scoring structure with the lowest I_rmsd was selected as the successful docking conformation of the Nb-RBD complex. Finally, five out of the 16 designed Nbs were found to be able to bind on RBD. To verify the reliability of selecting designed Nbs-RBD complex's near-native conformation from free energy landscape of protein docking, three recently reported Nbs (SR4 (Li, Cai, et al., 2020) , H11-D4 , VHH-72 (Wrapp et al., 2020) ) occupying different binding sites of RBD ( Figure S3a ) were used to redocking study with the same setup as mentioned above.

The starting coordinate used to run MD simulation was nearnative conformation of Nbs-RBD extracted from Rosetta protein-protein docking trajectory. Before the MD simulation, AMBER ff14SB (Maier et al., 2015) force field parameters were assigned to two proteins, and missing atoms were added by applying LEaP (Case et al., 2005) to each complex. Then, an appropriate number of counterions was added to maintain the electro-neutrality of the studied system, and each complex was immersed into a rectangular periodic box of TIP3P (Hornak et al., 2006) water molecular with an edge of 10.0 Å. All molecular dynamics simulations were carried out by running a GPU-accelerated simulation (Case et al., 2005; Gotz et al., 2012) in the PMEMD module of AMBER14 software. First, two steps of energy minimization were performed, the whole system was minimized by a harmonic restraint with a 10.0 kcal mol −1 Å −2 force constant in the NVT simulation followed by second minimization with no restraint. Next, the system was gradually heated to 100.0K in 2,500 steps and then heated to 310.0 K in 5,000 steps also with the force constant 10.0 kcal mol −1 Å −2 to the complex. Subsequently, the equilibrate simulation at constant pressure (Case et al., 2014) with isotropic scaling simulation at 310.0 K in 500 ps. Finally, a 200 ns production run carried out without any restraint on all studied systems at a temperature of 310.0 K controlled by the Langevin thermostat with a collision frequency (Case et al., 2014) of 1 ps −1 and a pressure of 1 atm under the control of Monte Carlo barostat (Case et al., 2014; barostat = 2) . During the whole simulation, bond lengths concluding hydrogens were constrained using SHAKE (Case et al., 2014) , a 2 fs time step was set and electrostatics and Lennard-Jones cutoff radius of 10.0 Å (Case et al., 2014) was considered. The trajectory analysis was done with the CPPTRAJ (Roe & Cheatham, 2013) module of AMBER14 (Case et al., 2014) and PyMOL (The PyMOL Molecular Graphics System, Version 2.0 Schrödinger, LLC) software packages were used to visualize the structure extracted from the MD trajectory.

The binding free energy (∆G tol ) for each complex was calculated by the molecular mechanics/generalized Born surface area (MM/GBSA) method as below equation:

where ∆E vdW and ∆E ele represent van der Waals and electrostatic energy changes in the gas phase. ∆G pol and ∆G nonpol represent polar and non-polar solvent interaction changes. ∆E vdW and ∆E ele were calculated using AMBER ff14SB (Maier al., 2015) in the gas phase. The polar part ∆G pol is calculated by generalized Born (GB) model. Dielectric constants of 8.0 and 80.0 are used for solute and solvent, respectively. The non-polar solvation contribution is computed using the solvent accession surface area (SASA) with LCPO method in AMBER molsurf module, according to:

where γ and β represent the surface tension and constant were set to 0.005 kcal/(mol·Å 2 ) and 0 . SASA was used to calculate the non-polar component with a probe radius of 1.4 Å (Yang et al., 2011) . To identify important residues critical for Nb binding to RBD, energy decomposition analysis was performed by the MM/GBSA method at a per-residue basis with the SASA obtained from icosahedron ICOSA method (Weiser et al., 1999) .

The protein-protein interface redesign was carried out by RosettaDesign (Mandell et al., 2009; Wang et al., 2007) in two steps. First, Rosetta fast relax (Nivon et al., 2013) was performed on the PDB files of Nb bound RBD complex taken from the MD simulations. The chain ID of the complex was deleted before running the relax protocol. Then, a resfile which specifies the residues need to be redesigned was used to run the fixbb protocol (Huang et al., 2011) in the Rosetta. For each Nb-RBD complex, sequence alignment shows that the interface redesign may generate the same Nb sequences, and Nbs with different sequences among the top 5 redesigned structures were selected for further MD simulations and binding free energies calculation analysis using the same setup as mentioned above.

A total of 16 Nbs were designed by grafting the three CDRs as well as the FRs edge of "sdAbs" (Figure S1a) onto the camelid Nb scaffold ( Figure S1b ). The disulfide bonds of the camelid Nb scaffold were reserved during Nbs design. The Clustal Omega (v.1.2.4; Madeira et al., 2019) was used to create the multiple sequence alignment of the designed Nbs (Figure 1a ). The length and type analysis of the amino acid composing the CDRs indicated the sequence diversity of the designed Nbs. Structural models of the 16 Nbs were predicted by comparative modeling method in Modeller (v.9.24; Sali & Blundell, 1993) with suitable templates. Statistics on the models are summarized in Table S4 . Structural alignment of the Nbs models is shown in Figure 1b . As expected, the four FRs are almost the same while the major difference of the 16 Nbs models is in the three CDRs, especially the CDR3 showing significant structure diversity, which is highly consistent with their sequences (Robert & Gouet, 2014) .

The initial structures of Nbs in complex with RBD were constructed using the 16 prepared crystal structures of "sdAb" in complex with SARS-CoV, SARS-CoV-2, and MERS-CoV as references ( Figure S2 ). Analysis of the constructed initial structures indicates Nbs contact with the different binding sites of RBD ( Figure S4 ) which could be categorized into three groups by aligning 16 Nbs-RBD complex with ACE2-RBD complex. In group I, CS01, CS02, CS03, and CS04 occupy the ACE2 binding site of RBD ( Figure S4a ). In group II, CS05, CS06, CS07, CS08, CS10, CS12, CS13, and CS15 bind to the near-ACE2 binding site of RBD ( Figure S4b ). In group III, CS09, CS11, CS14, and CS16 bind to non-ACE2 binding site RBD ( Figure S4c ).

Based on the constructed initial structures, Rosetta docking was used to search the near-native conformations of Nbs bound to SARS-CoV-2 RBD. Docking funnel of the trajectory describes the characteristics of the interface score (Rosetta energy units, REU) of each decoy and the interface root-mean-square deviation (I_rmsd) with the initial structure to search for a near-native structure protein-protein binding complex (Lyskov et al., 2013) . Therefore, docking funnel analysis was performed to evaluate the possibility of designed Nbs binding to RBD. As a result, the trajectories of CS01, CS02, CS03, CS10, and CS16 bind to RBD got the ideal docking funnels which contain lowest REU and I_rmsd ≤4 Å (Chaudhury et al., 2011) . For each complex, the structure with the lowest REU was extracted from the docking trajectory for further analysis (Figure 2) . According to the classified information of constructed initial structures, CS01, CS02, and CS03 occupy the ACE2 binding site (Figure 2a-c) , which may directly block the interaction between RBD and ACE2, and the interaction surface of RBD and nanobody is almost in loop regions of (2) ΔG nonpol = SASA + RBD, whereas CS10 and CS16 bind to the near-ACE2 and non-ACE2 binding sites which contain more helix and beta regions of RBD (Figure 2d-e) , respectively, which may inhibit the interaction between RBD and ACE2 through allosteric mechanism. The proposed inhibitory mechanism of design Nbs could be understood that MERS-CoV or SARS-CoV S protein adopted multiple conformation states and nAbs bind to the different epitopes of the antigen (Zavrtanik et al., 2018) .

In addition, based on three recently reported cocrystal structures of Nbs-RBD complexes (PDB ID: 7C8V (Li, Cai, et al., 2020) , 6YZ5 , and 6WAQ (Wrapp et al., 2020) , redocking was performed to verify the reliability of using free energy landscape to find the near-native conformation of protein docking in this study. From Figure S3a ,b, the sequence of three Nbs SR4, HH-72, and VHH-72 is of great difference compared with designed stable Nbs (CS01, CS02, CS03, CS10, and CS16), and the structures of FRs are almost identical while the conformation of CDRs especially CDR3 is different. Interestingly, the three Nbs occupy the ACE2, near-ACE2, and non-ACE2 binding sites on RBD, similar to the five designed Nbs (Figures 2 and S3a) . The results show that all the three Nbs bind to RBD got ideal docking funnels with a lowest REU and I_rmsd ≤4 Å and show similar trend with designed Nbs-RBD complexes ( Figure S3b-d) .

The selected near-native conformation of the three Nbs-RBD complexes protein docking is consistent well with their cocrystal structures, indicating that the selected nearnative conformation of designed Nb-RBD complex from Rosetta docking is reliable and can be used for following analysis.

Starting from the near-native conformations of the five Nbs in complex with RBDs, MD simulations were performed to optimize their binding modes. The RMSD values of Cα atoms of each snapshot relative to the starting coordinates were calculated to monitor the stabilities of Nbs-RBD complexes. As shown in Figure 3a 1.93, and 3.11 Å, respectively, while CS02-RBD complex possesses a large fluctuation (from 1.85 to 8.94 Å) around 90 ns of the simulation. Compared with the initial conformation of CS02-RBD complex, it is found that CS02 associated from the initial site and bound to a new site on RBD after MD simulation (Figure 4b ). For CDRs of the five Nbs, most of the calculated RMSD values range below 2.5 Å (Figure 3e ). Among them, CDR3 region has larger fluctuation in complexes CS03-RBD and CS10-RBD, while CDR1 region has larger fluctuation in complexes CS01-RBD and CS02-RBD, and CDR2 region has larger fluctuation in complex CS16-RBD, indicating that those regions will accommodate to the conformation of RBD during the recognition process.

In addition, the averaged root-mean-square fluctuations (RMSFs) of each residue in Nbs and RBD protein during the last equilibrated 40 ns MD trajectories were calculated. In Figure 3c -d, all of the five complexes share similar trends of dynamics features. As expected, the fluctuation of α-helix or β-sheet regions on RBD is smaller than loop region in Nbs-RBD complex. Interesting, for the loop region of RBD, the values of RMSF on the interface residues interacting with Nbs (CS01-RBD (472-502), CS02-RBD (400-425), CS03-RBD (442-454, 483-502), [436] [437] [438] [439] [440] [441] [499] [500] [501] [502] [503] [504] [505] [506] [507] [508] , and CS16-RBD (368-382)) are almost smaller than 2 Å, suggesting that the loop region on interface of RBD is in stable state (Figure 3c ). The peak RMSF values of RBD in five complexes occur at residues H519 (CS01), N370 (CS02), S371 (CS03), and S477 (CS10 and CS16) located at loop region, all of the residues are far away from the binding site of each complex. For Nbs, the fluctuation of RMSF values on CDR3 of CS01, CS03, and CS16 is more stable than other regions and all the RMSF values are around 1 Å (Figure 3d) , while for CS02 the RMSF with lower values is concentrated in the regions near R38 on FR2 domain (Figure 3d ), and residues focused on CDR2 have a relatively small RMSF value of CS10 (Figure 3d) , indicating that CS01, CS03, and CS16 interact with RBD primarily through CDR3 domain. Apart from CDR3, the recognition of CS10 and RBD also depend on CDR2 in a degree. But for CS02, FR2 is the principal interacting domain which can remind that CS02-RBD are not in a suitable binding state and compared to the other complexes, the relatively higher RMSF values of CS02 bound RBD implying the weak binding of the complex.

We use the MM/GBSA method (Kollman et al., 2000) , which has been widely employed to characterize the thermodynamics properties of different types of protein-protein interactions (PPIs; Du et al., 2020; Tu et al., 2018; Wang et al., 2019) , to estimate the binding affinities between Nbs and RBD. Table 1 lists the total binding free energy (ΔG tol ) as well as their compositing terms including the electrostatic interaction energy (ΔE ele ), van der Waals interaction energy (ΔE vdW ), non-polar solvent energy (ΔG nonpol ), and polar solvent energy (ΔG pol ).

The calculated ΔG tol values of CS01, CS02, CS03, CS10, and CS16 binding to RBD are −67.09, −39.13, −99.57, −63.22, and −64.65 kcal/mol, respectively, which is probably due to the binding site of each Nb bound to the RBD is different as shown in Figure 2a -e. For the binding interface of Nbs-RBD, the number of key residues on the contact surface of Nbs as well as RBD recognized by MM/GBSA method can result in different calculated ΔG tol values. CS02-RBD complex contains four key residues on RBD interface and three key residues on Nb interface while CS03-RBD involves 10 key residues on RBD interface and 10 key residues on Nb interface, most likely the reason why the ΔG tol values between CS02-RBD and CS03-RBD make a huge difference. More information is shown in Table 2 .

For all of the five complexes, the ΔG tol are contributed primary by the ΔE vdW and the ΔE ele . The ΔG nonpol contributed little to ΔG tol , while the ΔG polar seems to be unfavorable to the binding of the Nbs to RBD (Figure 3b ). For CS03-RBD complex, the most favorable ΔE vdW contribution resulting in the highest binding affinities among the studied complexes. Although the CS02-RBD complex found a relative stable binding state (Figure 2b) , the estimated value of ΔG tol further suggests that the binding interface between CS02 and RBD was not well designed. 

To further understand the differentiated binding affinities of the designed Nbs binding to the different epitopes of the antigen, the important residues located at proteinprotein interface were identified by decomposing the ΔG tol at per-residue basis. Herein the residues with a favorable energy contribution more than −2 kcal/mol were defined as "hot spots". As shown in Table 2 , various "hot spots" were discovered for the Nbs-RBDs complexes formation, and for more clearly, we identified "hot spots" on two parts, Table 2a ,b represent the "hot spots" on the RBD interface and Nbs contact face of all studied systems, respectively. The number of "hot spots" on the RBD interface in CS01-RBD, CS02-RBD, CS03-RBD, CS10-RBD, CS16-RBD complexes are 7, 4, 10, 5, and 7, respectively, and the number of "hot spots" on the Nbs interface are 7, 3, 10, 6 and 6, consistent very well with the trends of their predicted binding affinities. Mapping the "hot spot" on the RBDs and Nbs interfaces is shown in Figure S5 . Taking CS01-RBD complex as an example, N487 (−3.86 kcal/mol), F486 (−3.79 kcal/ mol), Y489 (−3.08 kcal/mol), F456 (−2.48 kcal/mol), K417 (−2.31 kcal/mol), L455 (−2.21 kcal/mol) and A475 (−2.04 kcal/mol) form "hot spots" on the RBD interface, and Y106 (−4.75 kcal/mol), F107 (−3.74 kcal/mol), F59 (−3.46 kcal/mol), N57 (−2.61 kcal/mol), R56 (−2.38 kcal/ mol), K65 (−2.36 kcal/mol), and Y60 (−2.01 kcal/mol) compose the "hot spots" of CS01. All of the "hot spots" on RBDs interfaces of CS01, CS02, and CS10 bound complexes belong to loop region ( Figure  S5a-b,d) , and more than half (4/7) "hot spots" on RBD interface of CS16 ( Figure S5e ) bound complex located in loop region. The RBDs interfaces share the same "hot spot" Y489 in CS01 and CS03 bound complexes, demonstrating that Y489 of RBD will play an important role in interacting with neutralizing Nbs. Different from the CS01, CS02, and CS03 bound RBDs, the "hot spots" on the RBD of CS10 and CS16 bound RBDs interfaces represent the allosteric binding sites ( Figure S5d-e) .

The identified "hot spots" of the five designed Nbs are nearly located in the CDRs interfere with RBDs except for CS02 ( Figure S5a-e) . The most favorable residues of the five Nbs are Y106 (CS01, −4.75 kcal/mol, CDR3), S118 (CS02, −3.8 kcal/mol, FR4), E56 (CS03, −4.03 kcal/ mol, CDR2), I101 (CS10, −4.41 kcal/mol, CDR3), and I29 (CS16, −4.35 kcal/mol, CDR1). The "hot spots" occupancies in CDR2 and CDR3 of CS01, CS03, CS10, and CS16 are 4/7, 7/10, 4/6, and 4/6, respectively. For CS02, all key residues are in FRs, which are consistent with the lower binding affinity (−39.13 kcal/mol) of CS02-RBD compared to the other complexes. Therefore, the residues on CDR2 and CDR3 regions especially CDR3 are of great importance for nanobody bound to RBD, enhancing interaction of CDR2 and CDR3 to RBD can do a great help to design Nbs.

In addition, "warm spots" contributing the binding free energies from −0.50 to −2.00 kcal/mol located at Nbs-RBD interfaces were also identified. As shown in Table S5 , residue Q493 of RBD was identified as "warm spots" both in CS01 and CS03 bound complexes. It is interesting that residues L455 and F486 of RBD were identified as "hot spots" in CS01 bound complexes, while they were "warm spots" for CS03 binding. The result further verified the mechanism that CS01 and CS03 occupied the ACE2 binding site on RBD. Figure 5 shows the interfaces of the representative structures extracted from the last equilibrated 40 ns MD trajectories of the five Nbs-RBD complexes. For clearly view, only "interacting residues" were displayed. It should be noted that residues involved in protein-protein interactions are nearly all "hot spots" and "warm spots" (Figure S5 ), proving that per-residue energy decomposition calculation is suitable for recognizing interactions between the Nbs and RBD interface.

Those information located at the Nbs-RBD interfaces not only imply that targeting loop region of RBD with different mechanisms could be an essential way for SARS-CoV-2 neutralizing Nbs design, but also provide useful information for better understanding the results that CS01 and CS03, CS10, and CS16 still bind at the ACE2, near-ACE2, and non-ACE2 binding sites of RBDs (Figure 4a ,c-e), while CS02 found another binding site on RBD (Figure 4b ) after MD simulations.

Through grafting the CDRs of "sdAbs" onto a selected camelid Nb scaffold, five designed Nbs theoretically binding to SARS-CoV-2 RBD have been obtained by Rosetta docking and MD simulation evaluation. Here, RosettaDesign (Mandell et al., 2009; Wang et al., 2007) of the Nb interface residues was performed to enhance the Nbs-RBD interaction specificity while maintaining or improving the binding affinity (Kortemme et al., 2004; Potapov et al., 2008) . Based on the representative conformations of CS01, CS02, CS03, CS10, and CS16 bound complexes, a number of brand-new sequences were generated. Different sequences of top five redesigned structures ( Figure S6 ) were selected for further MD simulations and binding free energies calculations. Here we listed redesigned sequences bind to the RBD with more stability and improved binding affinities after MD simulations, one redesigned sequence of CS02 (CS02_RD01), four redesigned sequences of CS03 (CS03_RD01, CS03_RD02, CS03_RD03, Table S6 . The ΔE vdW is still the main driving force for redesigned Nbs binding to RBDs followed by ΔE ele . Sequence alignments of the redesigned Nbs and their templates are displayed in Figure 6a with the mutations highlighted in white background. Figure 6b shows the Cα RMSD of the seven redesigned Nbs-RBD complexes. CS02_ RD01-RBD and CS16_RDs-RBD become more stable after redesign, and CS03_RDs-RBD still maintain slight RMSD values fluctuation throughout the MD simulations. Figure 6d indicates that the binding affinity of CS02_ RD01-RBD (−79.13 kcal/mol) is almost two times higher than that of CS02-RBD (−39.13 kcal/mol). Compared with the per-residue energy decomposition result of CS02_ RD01 and CS02 after MD simulation (Table 2) , it can be assumed that the binding site of CS02_RD01 on RBD has changed after redesign. This was further supported by several "hot spots" or "warm spots" on RBD interface shared by CS02_RD01 and CS01, such as Y449, Y505, Y489, and Q493. Structural alignments show that the binding site of CS02_RD01 on RBD is close to that of CS01 (Figure 7) , which can directly block the interaction between RBD and ACE2.

After redesign, the total number of "hot spots" on the Nb-RBD interfaces enlarged (Table 2 ) except for CS16_RD01 bound complex. Table 2 also shows that some of the "hot spots" on the RBD (Table 2a) or Nb (Table 2b) interfaces contributing more energies for protein-protein interactions for 12 studied complexes. Compared with CS02-RBD and CS02_ RD01-RBD, the energy contribution of "hot spot" changed from −2.67 kcal/mol (Q415) to −4.61 kcal/mol (F456) on the interface of RBD, the binding site of CS02_RD01 on RBD has changed a lot from CS02 (Figure 7a ) closing to that of CS01. The residue Y449 of RBD is important for the binding of CS03. In all CS03_RDs-RBDs, the energy contribution of Y449 was improved and with CS03_RD04 increased by −1.26 kcal/mol, by aligning CS03-RBD and CS03_RDs-RBD on the basis of Cα backbone of RBD. Figure 8 shows the orientation of benzene ring on Y449 of RBD getting closer to Nb through monitoring the atom between nitrogen atom on residue F114 of Nb and oxygen atom on residue Y449 of RBD. In particular, the change in distance from 8.8 Å in CS03-RBD complex to 2.8 Å in CS03_RD04-RBD F I G U R E 7 Comparison of binding sites between the initial structures and (re)designed Nbs of (a) CS02_Nbs-RBD, (b) CS03_Nbs-RBD, and (c) CS16_Nbs-RBD complexes, respectively. CS01 shown in (a) demonstrates its similar binding site with CS02_RD01. The RBD protein of SARS-CoV-2 is represented in white surface [Colour figure can be viewed at wileyonlinelibrary.com] complex suggests that redesign method can make tight contact and improve the binding affinity of Nb-RBD complex. The "hot spot" with largest energy contribution on RBD was changed from Y378 (−3.59 kcal/mol) to C377 (−4.55 kcal/ mol) and Y378 (−4.33 kcal/mol) in CS16, CS16_RD01, and CS16_RD02, respectively, and this little difference may be due to the Nb conformation readjusted after CS16 redesign (Figure 7c ). Similar to RBD, "hot spots" on the interface of CS02_ RD01 vary greatly. Figure 6e shows that the mutations I54V and M102D lead to the energy contributions from 0 and 0 kcal/mol (CS02) to −3.06 and −2.86 kcal/mol (CS02_ RD01), and these two major changes happened to CS02 reveal the information of the binding epitope of CS02 to RBD has reorganized after redesign and residues at position 54 and 102 deserve more attention. After the redesign, residue Y51 on CS16_RDs emerged and made a difference in the binding of CS16_RDs and RBD than CS16, and the energy of Y51 both is greater than −3.0 kcal/mol indicating the importance of Y51 means to CS16_RDs. Figure 6g demonstrates mutation energy change of CS16 Nbs, and mutations W32F (CS16_RD01, CS16_RD02) and I101Y (CS16_RD02) are beneficial for CS16 binding.

In all four CS03_RDs, the increased energy contributions of mutations R45P, T59Q, and S110Y (Figure 6f ) suggested that the mutant residues may be more suitable to recognize RBD. Apart from that, mutations E46Q, E56V(Q), D64P, and W108T(R) still have great energy values higher than −2 kcal/mol, the residues on the position of 46, 56, 64, and 108 are critical to CS03 Nbs, and those mutations are of benefit to CS03 binding. Mutations K52Y (CS03_RD01, CS03_RD03), S55N (CS03_RD02), S55R (CS03_RD03), S65Q (CS03_RD04), V107M (CS03_RD04), and S109N (CS03_RD02, CS03_RD04) bring improvement in energy contributions. It is should be noted that the energy contributions of the same mutation (S55N, S109N, or S119K) in different redesigned complexes of CS03_RDs show opposite results, such as the energy of mutation S55R gets lower for CS03_RD03 but higher for CS03_RD01 after redesign. Interestingly, in the series of CS03_RDs, the values of the residue's largest energy were changed, but the number of "hot spots" and "warm hot spots" was hardly to be changed. This could be explained by the more conserved binding mode between CS03_Nbs and RBD before and after redesign (Figure 7b) , formed by residue from CDR2 and CDR3 of CS03_Nbs and residues V113-L120, E152-F158, and F165-T168 from RBD. 

In this work, 16 Nbs (CS01-CS16) were designed by grafting CDRs sequence of CoVs' nAbs or mAbs into a known stable nanobody scaffold. And five out of the 16 Nbs (CS01, CS02, CS03, CS10, and CS16) were successfully docked onto SARS-CoV-2 RBD with different binding mechanisms. Among them, CS01, CS02, and CS03 occupy the ACE2 binding site, while CS10 and CS16 bind to the near-ACE2 and non-ACE2 binding sites, respectively. Further MD simulation and binding free energy analysis indicated that, except for CS02, the other four designed Nbs stably bind to the RBD of SARS-CoV-2. Based on the structures of the five Nbs in complex with RBD and seven brand-new Nbs with improved stabilities and binding affinities were generated by redesign of residues on the interface of the five Nbs contact with RBD. Especially for redesigned CS02-RBD complex, the calculated binding free energy increased from −39.13 kcal/mol (CS02-RBD) to −79.13 kcal/mol (CS02_RD01-RBD). As a result, 11 (re)designed Nbs (CS01, CS02_RD01, CS03, CS03_RD01, CS03_RD02, CS03_RD03, CS03_RD04, CS10, CS16, CS16_RD01, CS16_RD02) show stabilities and high binding affinities to SARS-CoV-2 RBD. Moreover, per-residue binding free energy decomposition analysis was performed to identify the "hot spots" between the Nb and RBD interface. For example, residue Y449 of RBD is important for CS03_Nbs-RBD bound complex, and mutations R45P, T59Q, and S110Y of CS03 can do benefit for the binding of CS03_Nbs and RBD. The identified "hot spots" provided useful information for us to understand the binding mechanism of designed Nbs to SARS-CoV-2 RBD. 

The Rosetta all-atom energy function for macromolecular modeling and design

The protein data bank

The Amber biomolecular simulation programs

Benchmarking and analysis of protein docking performance in Rosetta v3.2

Human neutralizing monoclonal antibody inhibition of middle east respiratory syndrome coronavirus replication in the common marmoset

Molecular simulation of oncostatin M and receptor (OSM-OSMR) interaction as a potential therapeutic target for inflammatory bowel disease

RosettaScripts: A scripting language interface to the Rosetta macromolecular modeling suite

Routine microsecond molecular dynamics simulations with AMBER on GPUs. 1. Generalized born

Protein-protein docking with simultaneous optimization of rigid-body displacement and side-chain conformations

Naturally occurring antibodies devoid of light chains

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

Comparison of multiple Amber force fields and development of improved protein backbone parameters

RosettaRemodel: A generalized framework for flexible backbone protein design

Neutralizing nanobodies bind SARS-CoV-2 spike RBD and block interaction with ACE2

Structural basis of neutralization by a human anti-severe acute respiratory syndrome spike protein antibody, 80R

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

COMBINES-CID: An efficient method for de novo engineering of highly specific chemically induced protein dimerization systems

Calculating structures and free energies of complex molecules: Combining molecular mechanics and continuum models

Computational redesign of protein-protein interaction specificity

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

Updated Approaches against SARS-CoV-2

A potent synthetic nanobody targets RBD and protects mice from SARS-CoV-2 infection. BioRxiv, 143438-v2

A humanized neutralizing antibody against MERS-CoV targeting the receptor-binding domain of the spike protein

Serverification of molecular modeling applications: The Rosetta Online Server that Includes Everyone (ROSIE)

The EMBL-EBI search and sequence analysis tools APIs in 2019

ff14SB: Improving the accuracy of protein side chain and backbone parameters from ff99SB

Sub-angstrom accuracy in protein loop reconstruction by robotics-inspired conformational sampling

Nanobodies: Natural single-domain antibodies

A Pareto-optimal refinement method for protein design scaffolds

Structural insights into immune recognition of the severe acute respiratory syndrome coronavirus S protein receptor binding domain

Computational redesign of a protein-protein interface for high affinity and binding specificity using modular architecture and naturally occurring template fragments

Structure of severe acute respiratory syndrome coronavirus receptor-binding domain complexed with neutralizing antibody

Deciphering key features in protein structures with the new ENDscript server

PTRAJ and CPPTRAJ: Software for processing and analysis of molecular dynamics trajectory data

Comparative protein modelling by satisfaction of spatial restraints

Improving the accuracy of PSI-BLAST protein database searches with composition-based statistics and other refinements

Cell entry mechanisms of SARS-CoV-2

Structural basis of receptor recognition by SARS-CoV-2

Statistical potential for assessment and prediction of protein structures

Nanobodies as therapeutics: Big opportunities for small antibodies

Prediction of GluN2B-CT1290-1310/DAPK1 Interaction by protein(-)peptide docking and molecular dynamics simulation

A two-step approach for the design and generation of nanobodies

Unexpected receptor functional mimicry elucidates activation of coronavirus fusion

Protein-protein docking with backbone flexibility

A novel coronavirus outbreak of global health concern

Improved sidechain modeling for protein-protein docking

End-point binding free energy calculation with MM/ PBSA and MM/GBSA: Strategies and applications in drug design

Importance of neutralizing monoclonal antibodies targeting multiple antigenic sites on the middle east respiratory syndrome coronavirus spike glycoprotein to avoid neutralization escape

Evaluation of candidate vaccine approaches for MERS-CoV

Approximate atomic surfaces from linear combinations of pairwise overlaps (LCPO)

sdAb-DB: The single domain antibody database

Structural basis for potent neutralization of betacoronaviruses by single-domain Camelid antibodies

Molecular dynamics simulation and free energy calculation studies of the binding mechanism of allosteric inhibitors with p38alpha MAP kinase

Junctional and allele-specific residues are critical for MERS-CoV neutralization by an exceptionally potent germline-like antibody

Structural basis for the neutralization of MERS-CoV by a human monoclonal antibody MERS-27

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

Structural basis of epitope recognition by heavy-chain camelid antibodies

Structural definition of a unique neutralization epitope on the receptor-binding domain of MERS-CoV spike glycoprotein

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

Structural definition of a neutralization epitope on the N-terminal domain of MERS-CoV spike glycoprotein

A pneumonia outbreak associated with a new coronavirus of probable bat origin

The authors declare no competing financial interests.

W.X. and F.Z. designed the research. J.Y. and Z.Z. performed the research. J.Y., Z.Z., F.Y., H.Z., H.W., and W.X. analyzed the data. J.Y., Z.Z., F.Z., and W.X. wrote the manuscript. All authors reviewed the manuscript.

The data that supports the findings of this study are available in the supplementary material of this article