Homology Modeling-Based in Silico Affinity Maturation Improves the Affinity of a Nanobody

Affinity maturation and rational design have a raised importance in the application of nanobody (VHH), and its unique structure guaranteed these processes quickly done in vitro. An anti-CD47 nanobody, Nb02, was screened via a synthetic phage display library with 278 nM of KD value. In this study, a new strategy based on homology modeling and Rational Mutation Hotspots Design Protocol (RMHDP) was presented for building a fast and efficient platform for nanobody affinity maturation. A three-dimensional analytical structural model of Nb02 was constructed and then docked with the antigen, the CD47 extracellular domain (CD47ext). Mutants with high binding affinity are predicted by the scoring of nanobody-antigen complexes based on molecular dynamics trajectories and simulation. Ultimately, an improved mutant with an 87.4-fold affinity (3.2 nM) and 7.36 °C higher thermal stability was obtained. These findings might contribute to computational affinity maturation of nanobodies via homology modeling using the recent advancements in computational power. The add-in of aromatic residues which formed aromatic-aromatic interaction plays a pivotal role in affinity and thermostability improvement. In a word, the methods used in this study might provide a reference for rapid and efficient in vitro affinity maturation of nanobodies.


Introduction
CD47 molecule, a ubiquitous cell-surface receptor that promotes immune evasion by interacting with signal-regulatory protein alpha (SIRPα) [1], is a member of immunoglobulin (Ig) superfamily and is considered as a promising cancer biomarker [2]. Targeted blocking CD47-SIRPα interaction for engaging macrophages to attack cancer cells represents a potentially promising immunotherapeutic strategy. Several antibody agents have been applied in clinical trials and display a positive effect [3,4]. However, there is only one approved nanobody medicine [5] and none for targeting CD47 ext . We had screened a panel of single anti-CD47 nanobodies from a semi-synthetic library via Phage display technology (submitted work), and the optimal one (Nb02) was chosen as the parental VHH in this study.
Apart from monoclonal antibodies (mAbs), small antibody fragments of different formats have become increasingly popular binders in research and diagnostic applications due to their exceptionally high stability and selectivity [6][7][8]. Some strategies about affinity optimization were undertaken previously to improve the binding affinity and stability of biotherapeutics [9][10][11][12]. Nevertheless, given the time-consuming process of conventional approach and improved computational capabilities, structure-based computational methods involved the in silico selection of promising candidates with a high likelihood of improving binding affinity has aroused considerable interest [13][14][15][16][17].
was chosen as the optimal one, which was model number 19. The SAVES evaluation is 100% pass, and Ramachandran plots were further used to assess and evaluate Nb02 (Figure 1b and Figure S2). Figure 1b shows the Ramachandran plot of the main residues except for G and P. Main residues located in most favored regions and additional allowed regions by 92.8% and 7.2%, respectively. In Figure S2, all 15 G residue and 3 P residues located in allowed regions. These indicate that the structure of Nb02 is rational and high-quality. Ramachandran plots were further used to assess and evaluate Nb02 (Figure 1b and Figure S2). Figure  1b shows the Ramachandran plot of the main residues except for G and P. Main residues located in most favored regions and additional allowed regions by 92.8% and 7.2%, respectively. In Figure S2, all 15 G residue and 3 P residues located in allowed regions. These indicate that the structure of Nb02 is rational and high-quality. On the Nb02 structure, we ran 50 ns MD (Molecular Dynamics) simulation to eliminate unreasonable clashes and minimize energy, and analyze the whole production run and choose the lowest internal potential energy for subsequent analysis [34]. Compared to the modeled structure, significant conformational rearrangements were observed in their CDRs (Complementarity determining regions), especially with the CDR3 (as in Figure 2a). Figure S3 shows that the modeled VHH is stable over time. Along the simulated time, the backbone root means square deviation (RMSD, Figure S3a) remains smaller than 0.2 nm. Moreover, for the root mean square fluctuation (RMSF, Figure S3b), as expected, it was greater than 0.2 nm only in the predicted CDR loops, which was free to explore different conformations. As shown in Figure 2, there are also some rearrangements in backbone regions; however, the RMSD and RMSF value shows it has a much smaller impact. Moreover, nanobodies usually bind to its antigen via CDRs, especially CDR3; thus, it may be neglected in the analysis. On the Nb02 structure, we ran 50 ns MD (Molecular Dynamics) simulation to eliminate unreasonable clashes and minimize energy, and analyze the whole production run and choose the lowest internal potential energy for subsequent analysis [34]. Compared to the modeled structure, significant conformational rearrangements were observed in their CDRs (Complementarity determining regions), especially with the CDR3 (as in Figure 2a). Figure S3 shows that the modeled VHH is stable over time. Along the simulated time, the backbone root means square deviation (RMSD, Figure S3a) remains smaller than 0.2 nm. Moreover, for the root mean square fluctuation (RMSF, Figure S3b), as expected, it was greater than 0.2 nm only in the predicted CDR loops, which was free to explore different conformations. As shown in Figure 2, there are also some rearrangements in backbone regions; however, the RMSD and RMSF value shows it has a much smaller impact. Moreover, nanobodies usually bind to its antigen via CDRs, especially CDR3; thus, it may be neglected in the analysis.
As experiments show that there is a competition between the VHH and the anti-CD47 mAb B6H12.2, we take the corresponding CD47 ext interaction sites as the CD47 ext active residues in this research, that is to say, the residues 29, 31, 34, 35, 37, 39, 46, 97 and 99-104 [35]. The active residues of VHH were determined as the CDRs. We then dock the Nb02 final configuration obtained after MD simulation to the selected CD47 ext binding sites. First, we performed a short-time pre-docking to check parameters and determine conformation rationality. After analysis of pre-docking results, we changed the docking time parameter and started production docking process. The cluster analysis and estimate results using docking score and Van der Waals energy are shown in Figure S4. We first selected three clusters with the lowest docking score and Van der Waals energy: number 1, 2 and 6. Then, the cluster with the lowest RMSD value, number 1, was chosen for subsequent analysis, as shown in Figure 2b. Comparison between different complex cluster structures indicates different binding modes of Nb02-CD47 ext (as shown in Figure S4d). unreasonable clashes and minimize energy, and analyze the whole production run and choose the lowest internal potential energy for subsequent analysis [34]. Compared to the modeled structure, significant conformational rearrangements were observed in their CDRs (Complementarity determining regions), especially with the CDR3 (as in Figure 2a). Figure S3 shows that the modeled VHH is stable over time. Along the simulated time, the backbone root means square deviation (RMSD, Figure S3a) remains smaller than 0.2 nm. Moreover, for the root mean square fluctuation (RMSF, Figure S3b), as expected, it was greater than 0.2 nm only in the predicted CDR loops, which was free to explore different conformations. As shown in Figure 2, there are also some rearrangements in backbone regions; however, the RMSD and RMSF value shows it has a much smaller impact. Moreover, nanobodies usually bind to its antigen via CDRs, especially CDR3; thus, it may be neglected in the analysis. Side view of Nb02-CD47 ext complex. The residues associated with the B6H12.2 binding sites (yellow) on the CD47 ext (green) and those of the CDRs (red) on the VHH (blue) is highlighted. The gray region presents the potential binding domain.

Identifying Key Residue Positions for Affinity Maturation
For affinity maturation, the identification of critical residues of the interaction between Nb02 and CD47 ext was crucial [10,36,37]. In the ADAPT computational approach, the first vital step is an exhaustive single-point virtual mutagenesis along the entire CDR sequence [30]. This process involved in whole CDR residues and increased the later calculation difficulty, so it is necessary to narrow down the initial interaction hotspots by a credible method, which is described in the Materials and Methods section.
As described in the Materials and Methods section, approximately 50~60% hotspots of CDRs have constant DNA sequence, which usually also encodes the antibody's antigen-binding sites [38,39]. Moreover, the VHH has a more extended CDR3 than other antibody fragments; thus, mutagenesis aims at AGY/RGYW (R = A/G, Y = C/T and W = A/T) mutational hotspot in the CDR3 are more likely to increase the affinity. As shown in Table 1, eventually, there are six binding sites identified. Additionally, to avoid missing some residues and improving accuracy and reliability, the interface residues of Nb02 within 5 Å of CD47 ext were counted using PyMol, the results are also shown in Table 1. In a word, the active residues of the Nb02 are identified: residues S55, V58, G107, T108, S109, F110, and are used as initial mutation residues for the subsequent ADAPT analysis.

ADAPT Cycle-multiple Mutants
In the present study, mutants were narrowed down using in silico screening and computational modeling, thus reducing in vitro analytical effort. There is a disulfide bond between C24 and C99 of Nb02, which is a hallmark of nanobody and helps to stabilize the extended CDR3 loop configuration. Therefore, each residue is mutated in turn to 17 other possible natural residues (C and P were excluded) and totaling 102 single-point mutants. In Figure S5, we show the proximity of the selected sites to the CD47 ext .
After ADAPT round, 102 single-point mutations covered six positions that might alter antigen-binding affinity were computationally evaluated (as shown in Table S1). As described above, the top 50 single-point mutants with improved antigen-binding affinity were predicted, while utilizing a three force-filed-based scoring platform. The results are shown in Table 2. The sub-selection process is based on residue types, site diversity, and visual examination. This process narrowed down the selection to 13 mutants from the 50 top-scoring mutants for further evaluation. The CD47 ext -binding affinities of these 13 single-point mutants were then calculated and compared with the parental VHH. The K D cal and improvements in binding affinities relative to the parental VHH are presented in Figure 3. Almost all the mutants demonstrate improved binding affinities and mutations have a positive effect, excepted V58D showing no improvement. The mutant S109V is showing the highest 21.54-fold increase, with another two mutants (T108H and S109E) an over 10-fold increase. The location of these three mutations relative to the CD47 ext epitope is shown in Figure S5b. Additionally, mutant S109G with the substitution of residues G led to small improvement than with V and E, and T108K shows a minor increase than T108H. Mutation G107T and replacement of S55 shows only marginal gains compared to parental VHH. From the above, in this first round of ADAPT screening, there was no clear trend as to which kind of residues was more useful to the binding affinity improvement. For instance, considering the best two single-point mutants, the substitution of S109E and S109V belongs to the acidic residue and aliphatic residue, respectively, whereas the replacement of S109G had only small improvement. Six leading single-point mutations from round 1, G107I, G107W, T108H, S109E, S109V and F110A, which possessing greater than five-fold affinity improvement in K D cal , were selected to round 2 affinity maturation. Despite a slight affinity improvement, the S55Q was also included in round 2 since its structural proximity to CD47 ext ( Figure S5) and the charge-neutral mutation property. A combination of these mutant sites could lead to a total of 20 double mutants. These 20 mutants considered in this round were calculated computationally about the change in binding affinity within the ADAPT protocols, and then compared predicted affinity improvements with the component single-point mutations and parental VHH (marked as M0). There were three affinity-decreased mutants among the 20 double mutants and combined with analysis of other mutants, G107W and S109E were discarded. Consider the mutation site diversity, and after manual visual inspection of the virtual mutants, the mutations improved more than five-fold were selected to the next step. The next analysis leads to five triple mutants and two quadruple mutants (some mutants with lower improvement did not show in the figure). The seven selected multiple mutants were marked M1-M7, as shown in Figure 3. Figure S6 shows the experimentally determined and calculated K M0 D /K MX D ratios. The experimentally determined ratios different from the calculated ones in numerical value, and this might result from the uncertainty of calculation criteria. However, the order of this value is the same, and it did not affect the selection of the best mutant. Six leading single-point mutations from round 1, G107I, G107W, T108H, S109E, S109V and F110A, which possessing greater than five-fold affinity improvement in KD cal , were selected to round 2 affinity maturation. Despite a slight affinity improvement, the S55Q was also included in round 2 since its structural proximity to CD47 ext ( Figure S5) and the charge-neutral mutation property. A combination of these mutant sites could lead to a total of 20 double mutants. These 20 mutants considered in this round were calculated computationally about the change in binding affinity within the ADAPT protocols, and then compared predicted affinity improvements with the component single-point mutations and parental VHH (marked as M0). There were three affinity-decreased mutants among the 20 double mutants and combined with analysis of other mutants, G107W and S109E were discarded. Consider the mutation site diversity, and after manual visual inspection of the virtual mutants, the mutations improved more than five-fold were selected to the next step. The next analysis leads to five triple mutants and two quadruple mutants (some mutants with lower improvement did not show in the figure). The seven selected multiple mutants were marked M1-M7, as shown in Figure 3. Figure S6 shows the experimentally determined and calculated ⁄ ratios. The experimentally determined ratios different from the calculated ones in numerical value, and this might result from the uncertainty of calculation criteria. However, the order of this value is the same, and it did not affect the selection of the best mutant.

Purification of VHHs
A HisTrap™ column was used for purifying the VHHs. The results of SDS-PAGE indicated that all of the purified VHHs showed a single band of ~15 kDa ( Figure 4). Moreover, we measured the protein yield based on the purified VHHs. The yield of M0 was 1.10 mg/L, while the yield of M1-M7 were 3.20 mg/L, 2.78 mg/L, 1.92 mg/L, 3.84 mg/L, 1.46 mg/L, 1.05 mg/L and 0.60 mg/L, respectively.

Purification of VHHs
A HisTrap™ column was used for purifying the VHHs. The results of SDS-PAGE indicated that all of the purified VHHs showed a single band of~15 kDa ( Figure 4). Moreover, we measured the protein yield based on the purified VHHs. The yield of M0 was 1.10 mg/L, while the yield of M1-M7 were 3.20 mg/L, 2.78 mg/L, 1.92 mg/L, 3.84 mg/L, 1.46 mg/L, 1.05 mg/L and 0.60 mg/L, respectively.

SPR
The KD values were determined by SPR measurements, and the results are presented in Figure  5. By BIAcore (Biomolecular Interaction Analysis), the binding affinity of M0 was low (KD M0 = 278 ± 9.3 nM). To get affinity matured, we selected seven mutants by homologous modeling-based ADAPT platform. Compared with M0, the best mutation shows a 162.58-fold calculated improvement of M0 cal

SPR
The K D values were determined by SPR measurements, and the results are presented in Figure 5. By BIAcore (Biomolecular Interaction Analysis), the binding affinity of M0 was low (K D M0 = 278 ± 9.3 nM). To get affinity matured, we selected seven mutants by homologous modeling-based ADAPT platform. Compared with M0, the best mutation shows a 162.58-fold calculated improvement of binding affinity (K D M0 /K D cal ). Logically, the SPR was used to verify the K D values and improvement folds. As shown in Figure 5, the seven mutants bound CD47 ext with K D of 74.5 ± 5.5 nM, 11.4 ± 2.1 nM, 69.5 ± 6.1 nM, 57.0 ± 7.3 nM, 25.3 ± 3.7 nM, 58.2 ± 4.5 nM and 3.2 ± 1.0 nM, respectively, compared to M0. All mutants displayed a significant increase in binding affinity relative to M0, which suggests that our virtual selection protocol enables robust identification of affinity-improved mutants.

SPR
The KD values were determined by SPR measurements, and the results are presented in Figure  5. By BIAcore (Biomolecular Interaction Analysis), the binding affinity of M0 was low (KD M0 = 278 ± 9.3 nM). To get affinity matured, we selected seven mutants by homologous modeling-based ADAPT platform. Compared with M0, the best mutation shows a 162.58-fold calculated improvement of binding affinity (KD M0 /KD cal ). Logically, the SPR was used to verify the KD values and improvement folds. As shown in Figure 5, the seven mutants bound CD47 ext with KD of 74.5 ± 5.5 nM, 11.4 ± 2.1 nM, 69.5 ± 6.1 nM, 57.0 ± 7.3 nM, 25.3 ± 3.7 nM, 58.2 ± 4.5 nM and 3.2 ± 1.0 nM, respectively, compared to M0. All mutants displayed a significant increase in binding affinity relative to M0, which suggests that our virtual selection protocol enables robust identification of affinity-improved mutants. With a KD of 3.18 nM (Chi2 = 3.64), mutant M7 exhibited the greatest improvement in binding affinity of around 87.4-fold.   Figure S7 shows a direct comparison of affinity versus stability (T M -K D ) for M0 and seven mutants. mutants and parental VHH was carried out with qPCR. The results are shown in Figure 6. It shows that the melting curves of M0 and the TM value is 43.38 °C, which is averaged from five independent experiments. The other seven mutants with TM of 57.47 °C, 60.59 °C, 60.01 °C, 43.08 °C, 51.16 °C, 48.91 °C and 50.74 °C, respectively. With a TM of 60.59 °C, mutant M2 shows the highest enhancement in thermal stability of 17.21 °C, while the M7 enhanced 7.36 °C. Figure S7 shows a direct comparison of affinity versus stability (TM-KD) for M0 and seven mutants.

Indirect-ELISA by Lead Variant
The best mutant M7 obtained at the end of the virtual screening and validated after SPR and qPCR was used to test binding activity and thermal stability with M0 through ELISA, as shown in Figure 7. Targeting CD47 ext in vitro binding affinity data as a function of VHH concentration are shown in Figure 7a. The best affinity-matured mutant M7 exhibited a stronger binding signal than M0 at all levels. In this assay, both VHHs show high bind ability to CD47 ext with the decreasing of the sample concentration. It is interesting to note that M7 can reach CD47 ext inhibition levels around 120.2% (at 360 μg/mL) compared to M0, whereas the value is 157.5% at 5.625 μg/mL. The BSA (Bovine Serum Albumin) was also tested as a negative control. As shown in Figure 7b, the thermostability of VHHs with different temperature treated 10 min was measured by indirect-ELISA. The data show that it still has a great binding activity to CD47 ext after sample treatment, and all the OD450 values are more than 1.0 (except OD450 of M0 is 0.936 at 70 °C). In particular, the OD450 value of M7 at 70 °C is around 45.0% compared with at 4 °C, while the value is just 32.7% for M0. These results indicate that the thermal stability of M7 is stronger than M0, and it also echoes the results of the qPCR assessment.

Indirect-ELISA by Lead Variant
The best mutant M7 obtained at the end of the virtual screening and validated after SPR and qPCR was used to test binding activity and thermal stability with M0 through ELISA, as shown in Figure 7. Targeting CD47 ext in vitro binding affinity data as a function of VHH concentration are shown in Figure 7a. The best affinity-matured mutant M7 exhibited a stronger binding signal than M0 at all levels. In this assay, both VHHs show high bind ability to CD47 ext with the decreasing of the sample concentration. It is interesting to note that M7 can reach CD47 ext inhibition levels around 120.2% (at 360 µg/mL) compared to M0, whereas the value is 157.5% at 5.625 µg/mL. The BSA (Bovine Serum Albumin) was also tested as a negative control. As shown in Figure 7b, the thermostability of VHHs with different temperature treated 10 min was measured by indirect-ELISA. The data show that it still has a great binding activity to CD47 ext after sample treatment, and all the OD 450 values are more than 1.0 (except OD 450 of M0 is 0.936 at 70 • C). In particular, the OD 450 value of M7 at 70 • C is around 45.0% compared with at 4 • C, while the value is just 32.7% for M0. These results indicate that the thermal stability of M7 is stronger than M0, and it also echoes the results of the qPCR assessment.

Discussion
A month ago, Caplacizumb was the first in class biopharmaceutical launched in both Europe and the U.S. for the treatment of aTTP (acquired Thrombotic Thrombocytopenic Purpura) [40]. Impressive progress in the application and economic outlook of nanobody has revived the

Discussion
A month ago, Caplacizumb was the first in class biopharmaceutical launched in both Europe and the U.S. for the treatment of aTTP (acquired Thrombotic Thrombocytopenic Purpura) [40]. Impressive progress in the application and economic outlook of nanobody has revived the researcher's interest in developing nanobody-based biological drugs [41]. The nanobodies screened by ribosome display or phage display always possess a weaker binding affinity than mAbs, which limits its potential application in drug discovery [42][43][44]. In vivo, experimental methods such as display-based screening and enrichment based on kinds of libraries dominated the affinity maturation field [45]. With these methods, the affinity improvement could reach around 10~1000-fold [46][47][48]. However, its weakness such as being time-consuming and its randomness are also apparent and limit full application. Recent advancements in computational power, molecular dynamics (MD) trajectories and simulation have become a helpful tool to characterize the properties of protein [49]. With the increase of structural data for antibody-antigen complexes and the development and establishment of algorithms for binding affinity prediction, the computer-aided structure-based rational approach has been an attractive choice for antibody affinity maturation. Currently, the dominating limitation of ADAPT is the need for a pre-existing crystal structure of the target antibody-antigen complex, even though the ADAPT platform could bring an improvement around 10-100-fold range based on several limited applications. Ideally, it could be replaced by a high-quality model based on homologous structure relative to the parental sequence.
The primary consideration of this study is to develop an efficient, engineered, and universal protocol for the in vitro affinity maturation of nanobodies to lay the foundation for drug development.
The parental VHH of this study was selected via three rounds of panning of a synthetic nanobody library, which display a weak binding affinity targeting CD47 ext (K D = 278 nM). Initial concerns emerged due to the uncertainty of created strictly homologous structure compared to the actual fabric. Here, we build its 3D structure by MODELLER based on five homologous crystal nanobody structures. Ramachandran plots (Figure 1b and Figure S2) and MD simulation results (Figure 2 and Figure S3) show that the structures of Nb02 were rational and high-quality. The Nb02 was then utilized in ADAPT-based virtual affinity maturation, with the replacement of crystal structure, and the experimental validation of mutants shows a significant improvement in binding affinity and thermostability.
One limitation of the ADAPT protocol is the inadequate conformational sampling of the binding between the mutants and the target protein, essential for an accurate estimation of the binding affinity. This affinity maturation protocol approximates that the changes in the binding conformation after a single mutation are negligible. To observe these conformational changes, an MD simulation of each new complex was performed. The structure aligns results of M1~M7 show that their binding does not suffer great rearrangements.
Rational structure-based affinity maturation contributes to an understanding of the structural basis for the improvement of binding affinity, which is one of its advantages. Figure 8 displays the structural basis and interactions predicted for quadruple mutant M7 (G107W, T108H, S109V, F110A) with CD47 ext [50]. The substitution of residue T108 by H introduce a salt bridge with D51 on CD47 ext (Figure 8a). Furthermore, two sets of ionic interaction formed between T108H and D51, E35 and R103 of VHH. On the G107 site, mutation to W generate a hydrophobic interaction with Y37 of antigen with a distance 2.74 Å (Cα-Cα); relatively, the distance of previous Y37 hydrophobic interaction with I105 of VHH is 5.51 Å (Figure 8b). Additionally, the mutation of G107W constitutes hydrogen bonds: the main chain-side chain H-bound between K39 and G107W, the side chain-side chain H-bound between G107W and T99 (Figure 8c). The new intramolecular H-bound interactions introduced by W at position 107 are also likely contributed to the observed improvement of binding affinity. In Figure 8e, it shows there still formed a cation-Pi interaction between G107W and K39 (4.63 Å) and K41 (5.13 Å). This study reveals several key factors which may affect the efficiency and robustness during nanobody in silico affinity maturation. First, these results illuminated that multiple mutations are necessary to achieve a considerable improvement in binding affinity for the anti-CD47 VHH. Previous works demonstrate that multi-single point affinity-matured mutants result in a moderate increase in binding affinity of multiple mutations [10,25,26,30], while our finding is generally consistent with it. It is possible and relatively straightforward to identify and combine several affinity-enhanced single mutations; nevertheless, the additive effects of multiple mutations on nanobody affinity are complex and incompletely understood thus far. Additionally, it sometimes hurts binding affinity [55,56]. When ADAPT platform was used in Fab affinity maturation, the first step is an exhaustive single-point virtual mutagenesis along the entire CDR sequence. It is timeconsuming and abnormal to generate all possible combinations of single-point mutations, which led Another critical factor is the add-in of aromatic residue, which is beneficial to protein structure stabilization [51]. The aromatic residue W interacted with cation provided by electro-positive residue K and formed a stronger bond than a salt bridge, which could also be likely responsible for affinity improvement. Additionally, the mutation of G107 to W introduces a new aromatic-aromatic interaction with Y37 (4.65 Å) of CD47 ext (Figure 8e). The structure model shows that the F110 formed an aromatic-aromatic interaction with W62 neared CDR2 (Figure 8f). From above, we can see that even on the same site, there could be different mechanisms of improving binding affinity. Hence, the result figures out the importance of aromatic residue in the enhancement of affinity and stability. Inoue et al. carried out an in silico virtual analysis of a CDR3-grafted VHH and obtained a mutant that had a~20-fold improvement over the original protein as measured by SPR [10]. Their work and some others are all substitutions by charged residues, such as R and D. In contrast, the designed mutants in this study have informative ways in affinity improvement as they mutated to several kinds of residues, such as aromatic G107W and aliphatic F110A, and did not solely rely on electrostatics (alkaline T108H).
Taken together, our finding demonstrates that high-quality homology structure can afford ADAPT-guided binding affinity improvements. These data also show that ADAPT-guided affinity maturation of VHH for improved antigen-binding ability can also translate into an enhancement of thermostability. The combination of improved affinity and thermostability is commonly tricky during the process of nanobody affinity maturation. Prior research has shown that increased affinity results from mutations of CDRs sites always accompanied by structure destabilizing [26,30]. This destabilization may attribute to the change of strain of framework results from the residue replacement and corresponding structure modification. Another likely reason might be the change of chemical property of antigen-binding sites selected for affinity maturation [37]. Encouragingly, our six selected multiple mutants displayed addition stability (except M4) and the highest affinity increased variant M7 with 7.36 • C enhancement of T M value was identified by qPCR (T M of 50.74 • C relatives to 43.38 • C of M0). In detail, the M4 and M7, which include F110A mutation, have a weaker stability enhancement than others, the T M of M4 even smaller than M0. The most probable cause is that the six antigen-binding sites for affinity maturation are all located in the CDRs region, and its substitutions scarcely affect the whole stability. Another reason, as above, is the aromatic residue. The side chains of aromatic residue, which possess directional Pi-systems, always stabilize protein interiors and interface. Thus, substitution by aromatic residue could enhance the structural stability of protein [52,53], and in this study is F and W. Nevertheless, antibodies aggregate at high temperatures and thus may influence the accurate determination of their T M value [54]. Admittedly, the stability improvement observed by thermal denaturation might occur due to delayed onset of aggregation.
This study reveals several key factors which may affect the efficiency and robustness during nanobody in silico affinity maturation. First, these results illuminated that multiple mutations are necessary to achieve a considerable improvement in binding affinity for the anti-CD47 VHH. Previous works demonstrate that multi-single point affinity-matured mutants result in a moderate increase in binding affinity of multiple mutations [10,25,26,30], while our finding is generally consistent with it. It is possible and relatively straightforward to identify and combine several affinity-enhanced single mutations; nevertheless, the additive effects of multiple mutations on nanobody affinity are complex and incompletely understood thus far. Additionally, it sometimes hurts binding affinity [55,56]. When ADAPT platform was used in Fab affinity maturation, the first step is an exhaustive single-point virtual mutagenesis along the entire CDR sequence. It is time-consuming and abnormal to generate all possible combinations of single-point mutations, which led to a large amount of affinity evaluation and multiple-rounds of screening. Moreover, nanobody has a more extended CDR1 and CDR3 than the traditional antibody. These two factors determined that it is challenging to utilize ADAPT screening platform with whole CDR. As described in Materials and Methods section, the binding sites were identified by a Rational Mutation Hotspots Design Protocol (RMHDP). Each has different predicted binding sites, and we selected the overlapped parts as a credible and rational result. This method immensely narrowed down the initial mutagenesis amounts (17 mutations per CDR site across six sites in this study). These sites were utilized in virtual affinity maturation and eventually led to positive and encouraging results, as shown above.
Technically, virtual screening filtered the best mutant among all virtual candidates in the development to affinity-improved mutant and hence moved ahead of the affinity validation study. In this study, we undertook the preparation of a nanobody structure by homologous modeling. Compared to crystallography-determined structure data, the use of homologous modeling considerably lowering the hurdle to conduct affinity maturation in vitro to some degree. However, comparable data is still scarce, and further studies are still necessary. Furthermore, additional experiments will be based on synthetic biology concepts to rationally design bispecific nanobodies to enhance the affinity and targeting the ability of nanobody-based biotherapeutics. Moreover, the effect of the drug increases the phagocytic capacity of macrophages to cancer cells will be tested, while the antitumor effect on an animal model of carcinoma in situ.

Homology Modeling
The sequence of Nb02 (GenBank accession number MK780744) and CD47 ext are shown in Figure S8. The sequence of the Nb02 used for experimental validation in this study has a C-terminal His6 tag that was not included for in silico studies. We used the MODELLER v9.19 to build the structural models of Nb02. The experimentally determined sequences were searched in the protein databank with BLAST (http://blast.ncbi.nlm.nih.gov/Blast.cgi), E-value cutoff 10.0 and sequence identity cutoff 90%. The resulting nanobody sequences were aligned with CLC Sequence Viewer 7.8.1 (CLC bio, Finlandsgade, Aarhus N, Denmark). The highest resolution and the lowest B-factor of the results, namely 5LMW [57], 5LZ0 [57], 5HDO [58], 2X1O [59] and 5LMJ [57], were chosen as templates. First, 50 models were built and the optimum one was selected based on the use of three different criteria (molpdf, DOPE and GA341). After that, it was evaluated and optimized using the webserver SAVES (http://services.mbi.ucla.edu/SAVES/), and then assessed using the Ramachandran plot. Eventually, an MD simulation was performed to eliminate unreasonable structure clash and minimize the energy.

MD Simulation
All MD simulations and their analysis were run as implemented in the GROMACS package (http://www.gromacs.org) [60] with the AMBER99SB-ILDN42 force field [61] and SPC water. The Nb02-CD47 ext complexes were first minimized without solvent, and then placed the molecules in a cubic box with a water layer of 1.2 nm and the system was neutralized by adding two Cl − as counterions. Subsequently, a second minimization of the system was performed before starting the MD simulation. After minimization, the equilibrium was performed for 100 ps in the NVT ensemble, followed by 20 ns NPT production run. The temperature was controlled at 300 K by a modified Berendsen thermostat [62], with a time constant of 0.5 ps. The pressure was controlled at 1 bar by isotropic Parrinello-Rahman method [63] with coupling constants of 1.0 ps.
Moreover, a Periodic boundary condition was used for all MD simulations. The iteration time step was set to 2 fs with the Verlet integrator and LINCS (Linear Constraint Solver) constraints [64]. The MD simulations were then performed for 50 ns.

Molecular Docking
Through computational analysis combined with analysis of solvent-exposed hydrophobic regions, the active residues of the Nb02 are identified: residues 55, 58, 107-110. Additionally, passive residues were automatically defined. Molecular docking using AutoDock Vina (http://vina.scripps.edu/) [65] and its standard protocols were conducted to predict the Nb02-CD47 ext complex structure. The atomic coordinates of the CD47 ext were taken from the structure of the SIRPα-CD47 complex crystallized at pH 8.5 (PDB ID: 4CMM) [66]. Active residues were defined for each VHH as mentioned above, while CD47 ext active residues were those known to be in contact with the commercially available anti-CD47 mAb B6H12.2: residues 29, 31, 34, 35, 37, 39, 46, 97 and 99-104 [35]. While other parameters were set to the default values, the size of the grid box was 40 Å × 40 Å × 40 Å, and the Lamarckian genetic algorithm was used to search the complex pose. Among the docking results obtained, we chose the docking cluster with the lowest docking score and binding free energy for the top conformations.

Rational Mutation Hotspots Design Protocol
Two different criteria are used to identify and pre-select the interface residues which might affect the affinity. First, the CDR residues within a contact distance with the antigen might affect the binding affinity. Thus, we counted the surface residue residues of parental VHH within a 5 Å distance of CD47 ext using PyMol (http://www.pymol.org/). Second, mutations of variable regions in the affinity maturation frequently belong to some constant CDR sequences, namely "Hot Spots", which residues generally were AGY or RGYW (R = A or G, Y = C or T,W = A or T) [38,39,67,68]. These residues are also often located at antigen-binding sites for antibodies, and mutations aim at these sites is more likely than others to improve the binding affinity. The overlapped parts are selected as a credible and rational result of rational mutation hotspots design. The mutations were done with Swiss-PDB Viewer 4.1 (http://spdbv.vital-it.ch/) [69], and all obtained models underwent an MD simulation.

In Silico Affinity Maturation
It has been reported that a few residues' changes could significantly improve the binding affinity of proteins in a previous publication [10]. In the present study, ADAPT was used to optimize the binding affinity of Nb02. A 3D structure model of the antibody-antigen complex is necessary for the ADAPT calculation. In this study, the initial sequence of the Nb02 was screened from a semi-synthetic library via Phage display technology (submitted work). The protein structure of Nb02 is predicted by homology modeling using the protocol described above and were used as a starting point for virtual affinity maturation.
Some structural preprocessing had to be done before the computational screening. Initially, water molecule and counterions of the VHH modeling structure were removed. Secondly, hydrogen atoms were added and adopted to protonation states at neutral pH. Thirdly, structural refinement was carried out by energy-minimization using the AMBER force-field, with a distance-dependent dielectric (4r ij ) and an infinite distance cutoff for non-bonded interactions. The non-hydrogen atoms were restrained at their initial positions with a harmonic force constant of 10 kcal/(mol·A 2 ).
In the first round of affinity optimizations, six positions (S55, V58, G107, T108, S109, F110) within the CDRs of the parental sequence were mutated to 17 other possible natural residues (Cys and Pro were excluded). 102 single-point scanning mutagenesis was carried out in this step. Three protocols, SIE-SCWRL, FoldX and Rosetta, were used to building and repacking the structures and evaluating the energies of these single-point mutations. Each component energy function would give a Z-score about the binding affinity of these mutants, and a consensus Z-score was further applied for scoring these mutants to alleviate the influence of the median and median absolute deviation. For the ratio of improved binding affinity, K M0 D /K cal D , we employed HADDOCK webserver and Schrodinger to calculate the K D values. Further technical and implementation details of this consensus approach and its component methods can be found in Sulea et al. [26]. Consider the mutation site diversity, and after manual visual inspection of the virtual mutants, approximately 15 mutants from the 50 top-scoring mutants are selected for further optimizations. In the second round of optimization, combine the favorable single-point mutations with Z-score using the same computational protocol above, a series of double-and triple-point VHH mutants were generated, and so do the quadruple-point mutants' generation.

Protein Expression and Purification
The DNA sequences of VHH mutants were synthesized commercially by GENEWIZ (Beijing, China), subcloned into the pET-32a (+) expression plasmid, with N-terminal RBS/TATA box and C-terminal His6 tag, and were subsequently expressed in TransB (DE3) E. coli. For expression, 2 mL of the overnight culture was inoculated into LB media (200 mL) supplemented with 100 µg/mL ampicillin in 500 mL flasks with air top seals and grown until an OD 600~0 .6-0.8. Cultures were then induced with 200 µL IPTG (500 mM) and grown for 16 h at 18 • C with shaking at 180 rpm. After that, the interest proteins were extracted by sonication and centrifugation, and the supernatants get filtered through 0.45µM filters. The purification of VHHs was carried out by employing Ni 2+ -affinity chromatography (IMAC) using HisTrap™ HP (GE Healthcare, Tianjin, China). The binding buffer containing PBS pH 7.4, 500 mM NaCl, 80 mM imidazole, and the elution buffer containing PBS pH 7.4, 500 mM NaCl, 500 mM imidazole. The recombinant monomeric protein was visualized by SDS-PAGE.

SPR
The IMAC-purified VHHs were subjected to desalting and soluble in PBS (pH 7.4) before the SPR binding experiments. Binding interactions were conducted with an SPR-based Biacore 3000 instrument (GE Healthcare) at 25 • C. Protein concentrations were quantified by 280 nm absorbance with a Nanodrop2000 spectrometer (Thermo Scientific, Shanghai, China). Typically, approximately 350 resonance units (RUs) of enzymatically biotinylated CD47 ext (ab174029, Abcam, Shanghai, China) were directionally immobilized on a CM5 sensor chip (GE Healthcare). An unrelated biotinylated protein was immobilized with an RU value matching that of the reference surface to control for nonspecific binding. Measurements were made with serial VHH concentrations in Biacore running buffer (10 mM HEPES, 150 mM NaCl, 3 mM EDTA, 0.05%(v/v) P20, pH 7.4; GE Healthcare) at a flow rate of 40 µL/min. Running buffer without VHH was then passed over the chip to allow spontaneous dissociation at the same flow rate. The CD47 surface was regenerated by three 60 s injections of 10 mM Glycine-HCl buffer (pH 2.1) after each sample injection. All data were analyzed with the Biacore 3000 evaluation software (GE Healthcare) with a 1:1 Langmuir binding model. All VHHs were run in duplicates.

Thermal Stability Measurements
To investigate the thermal stability changes of VHHs, real-time fluorescence quantitative PCR (qPCR) was used to measuring the melting curves of the parental and mutant VHHs, and the melting temperatures (T m ) were acquired by taking a negative first derivative transformation [70]. qPCR was carried out in a Roche LightCycler ® 480 II real-time PCR instrument (country Roche, Beijing, China). A total volume of 20 µL in each tube of 96-well plates was used. 1.0 µL HEPES (1 M) and 0.04 µL TCEP (0.25 M) was added first. Samples were diluted in PBS at a final concentration, after mixing, of 2.0 µM. 17.46 µL 1× PBS was then added. SYPRO ® Orange (Roche, Beijing, China) was diluted 25-fold from the 5000× concentrated stock to the working dye solution in PBS, and 0.5 µL was added to the mixture just before the experiment. Thermal denaturation temperature was ranged from 25 • C to 95 • C at a rate of 0.01 • C/s. Fluorescence intensity, with excitation at 465 nm and emission at 580 nm, was collected at 1 • C intervals and analyzed with the correlative software Exor4 (Roche Applied Science) and LightCycler Thermal Shift Analysis. Each sample was measured in quintuplicate and the average arithmetic value was calculated.

Indirect ELISA
The specificity of M0 and selected indirect ELISA assessed the best mutant with indirectly coated CD47 ext /VHH complex in 96-well plates, where the bound VHH was detected using an HRP-conjugated anti-His tag mAb. The CD47 ext and BSA were diluted to 30 µg/mL and coated overnight at 4 • C. After PBST wash, add 200 µL blocking liquid and incubated with 2.5 h at 37 • C. For binding affinity test, add 100 µL VHH and incubated with 2.5 h at 37 • C. As for thermal sensitivity experiment, add 100 µL VHH with diluted concentration 30 µg/mL and incubated with 2 h at 37 • C. After anti-His tag mAb added, the plates were then stained with TMB (Solarbio ® , Beijing, China), and the OD 450 values were immediately measured in Infinite M Nano + (Tecan, Tianjin, China). For the thermal stability assessment, the VHHs were treated with a different temperature for 10 min rather than diluted to serial concentrations. BSA was used as negative control. Origin 2018 was used for the representation of all ELISA data.

Conclusions
In summary, the combination of in silico screening of VHH mutants and experimental measurements may be used for the estimation of the relative binding affinity maturation in similar protocols for the design and optimization of VHHs as binders of protein targets. The fast and efficient approach designed in this study, which uses homologous modeling structure in ADAPT platform to screen multiple mutations, is useful for increasing the binding affinity of a phage display VHH while enhancing stability. These results will need to be assessed for other VHHs to evaluate their commonality. With the development of computational power and algorithms' optimization, there is still room for further improvement in the accuracy and usability of virtual mutation screening. The structure features of nanobody guaranteed in silico affinity maturation done in a fast and convenient way. Moreover, the RMHDP utilized in optimizing binding sites identification was efficient and rational, which turns out to be valid to affinity improvement. It reduced unnatural residues and additional calculation, to overcome in silico affinity maturation period to 3~4 weeks from more than two months. These two points form the foundation of our computational approach. The platform described in this study for in silico nanobody affinity maturation based on homology modeling, the RMHDP and ADAPT will be further optimized. The beneficial properties of nanobody provide rapid development of nanobody-based high-affinity binders and offer new ideas for the treatment of "orphan" diseases and explosive infectious diseases.

Conflicts of Interest:
The authors declare no conflict of interest.