The crystal structures of four SARS-CoV-2 RBD targeting nanobodies (30–33) R14 (PDB ID: 7WD1) (30), DL4 (PDB ID: 7F5G) (31), V_H ab6 (PDB ID: 8DLX) (32), and Nanosota9 (PDB ID:9CO9) (33) and SARS-CoV-2 WT (PDB ID: 7WD1) (30), BA.2 (PDB ID: 7ZF8) (36), JN.1 (PDB ID: 8Y5J) (37) RBDs were obtained from the RCSB Protein Data Bank (38). The KP.3/XEC RBD structure was generated by computational mutagenesis using RosettaDesign (39) with JN.1 RBD as template. These nanobodies were subsequently docked with the RBDs using the HDOCK server (40), and favorable binding conformations were selected based on docking score and confidence score (>0.99).

MD simulation was conducted using GROMACS 2022.5 with the CHARMM36 force field (41). The complex solvated in a cubic box using the TIP3P water model, and counter ions were added to neutralize the protein in the aqueous system. Energy minimization was executed through the steepest descent algorithm over 50,000 iterations to alleviate any high-energy contacts. Then, a two-step equilibration process was performed for 100 ps. First, canonical ensemble (NVT, constant number of particles, volume, and temperature) equilibration was conducted, maintaining the temperature at 310.15 K using the v-rescale thermostat, followed by isothermal-isobaric ensemble (NPT, constant number of particles, pressure, and temperature) equilibration with the c-rescale method for isobaric coupling at a reference pressure of 1.0 bar. Throughout the simulation, periodic boundary conditions (PBC) were applied, and long-range electrostatic interactions were calculated using the Particle Mesh Ewald (PME) method with a cutoff radius of 12 Å. Finally, a 200 ns production MD run was performed on the optimized system, with a time step of 2 fs, recording energy and log data every 1000 steps. The MD trajectory was analyzed to calculate the root mean square deviation (RMSD) and the number of hydrogen bonds, providing insights into the stability and interactions of the complexes.

Binding free energy calculations were performed using the Molecular Mechanics/Poisson-Boltzmann (Generalized Born) Surface Area (MM/PB(GB)SA) method. To ensure robust sampling, equilibrium trajectory data were extracted from the production phase of the MD simulations. Specifically, for each system, the last 50 ns of equilibrium trajectory data were used, with snapshots collected at 50 ps intervals. This methodology provided a total of 1,000 frames for binding free energy calculations.

The binding free energy ( Δ G b i n d ) of the Nb-RBD complex can be regarded as the sum of enthalpy term (ΔH) and entropy term (-TΔS):

In this context, the enthalpy change is composed of ensemble average interaction energy E and solvation energy ( Δ G s o l ), in which E can be separated into the electrostatic interaction energy (ΔE_ele ) and van der Waals interaction energy (ΔE_vdW ).

The solvation energy ( Δ G s o l ) is further divided into polar solvation energy ( Δ G p b ) and nonpolar solvation energy ( Δ G S A ):

In this study, a single trajectory approach was adopted due to its lower noise level and its suitability for systems with minimal structural rearrangement during binding. The entropy contribution was generally excluded from the calculations (42), as entropy differences associated with relative binding affinities were expected to be small, with only minor variations arising from mutations. MM/PB(GB)SA calculations were performed using gmx_MMPBSA version 1.6.3 (43), a novel tool developed for endpoint free energy calculations from MD trajectories.

To further elucidate the binding interactions, an energy decomposition analysis was conducted to assess the contribution of each residue to the binding free energy. This analysis highlights key residues that stabilize or destabilize the Nb-RBD complex, identifying potential targets for therapeutic intervention.

The HDOCK generated Nb-RBD complexes were subjected to single-point mutations using DDMut-PPI (44). Single-point mutations were introduced at CDRs and hotspot residues, substituting each residue with the 19 alternative amino acids, excluding itself. The resulting changes in binding affinities were calculated as ΔΔG (ΔΔG=ΔG_WT -ΔG_mutant , kcal/mol), with ΔΔG > 0 indicating increased affinity and ΔΔG < 0 indicating decreased affinity compared to the wild type.

Empirical evidence and previous studies indicate that combinations of individually beneficial mutations are more likely to result in additive or synergistic improvements in binding affinity (45, 46). Therefore, we prioritized the combination of single-point mutations with ΔΔG > 0 for subsequent analysis. Stabilizing mutations (ΔΔG > 0) identified across variant complexes were subsequently combined for multiple mutations analysis using DDMut-PPI’s combinatorial module. Utilizing this approach, an in silico library comprising high-binding and stable single-point and multiple-points mutations in the nanobodies R14, DL4, V_H ab6, and Nanosota9 were generated.

To validate our streamlined computational pipeline, mutations identified through computational affinity maturation were introduced into R14. The resulting mutant structures were docked with the RBD of the KP.3/XEC variant using HDOCK, employing its scoring function to filter for the best-performing mutants. The evaluation focused on predicted binding affinities and structural compatibility with the RBD.

Subsequently, the optimized combinations for R14 were subjected to a 200 ns production MD simulation with WT and three variants (BA.2, JN.1, and KP.3/XEC) complexes. The dynamic stability and interaction profiles of the mutants were analyzed using RMSD and MM/PB(GB)SA method. RMSD was utilized to assess conformational changes in the protein backbone, while MM/PB(GB)SA was employed to calculate binding free energies. These analyses enabled the identification of high-affinity and stable mutants, further refining the selection based on dynamic stability and interaction profiles. This approach ensured a robust assessment of the best-performing mutants.

In this study, we evaluated the antiviral spectrum of four RBD-targeting nanobodies (R14, DL4, V_H ab6, and Nanosota9) against SARS-CoV-2 Omicron variants. Despite the RBD accumulating mutations at 3–5 times the rate of other spike protein domains (47), RBD-targeting antibodies remain clinically dominant due to their direct disruption of ACE2 binding (48). To elucidate their binding characteristics, we conducted a sequence alignment of RBDs from BA.2, JN.1, and KP.3/XEC. The analysis identified three critical mutations (L455S, F456L, and Q493E) within the RBD regions that directly interact with the nanobodies. Epitope mapping ( Figure 2A ) revealed that V_H ab6 and Nanosota9 target relatively conserved epitopes on the Omicron RBD, whereas R14 and DL4 bind to more variable regions. These findings indicate that the binding affinities of V_H ab6 and Nanosota9 exhibit minimal variability across different Omicron subvariants, enabling effective neutralization of the most predominant variants ( Figures 2F-I ).

Furthermore, structural analysis revealed that R14 primarily interacts with KP.3/XEC through main-chain functional groups, suggesting that side-chain variations exert minimal influence on binding ( Figure 2B ). This finding accounts for the partial neutralization efficacy of R14 against KP.3/XEC. In contrast, V_H ab6 and Nanosota9 involve changes in both main-chain and side-chain interactions, leading to a more pronounced influence of side-chain groups on overall antibody affinity ( Figures 2D, E ). These structural insights align with biochemical and virological findings regarding the antiviral efficacy of each nanobody against earlier and recent Omicron subvariants.

To further investigate the interactions between nanobodies and SARS-CoV-2 Omicron variants, we performed MD simulations on sixteen systems, including four nanobodies (R14, DL4, V_H ab6, and Nanosota9) and four RBDs (WT, BA.2, JN.1, and KP.3/XEC). We assigned the complexes formed by R14 with these RBDs as R14-RBD systems (R14-WT, R14-BA.2, R14-JN.1 and R14-KP.3/XEC). The naming convention for the complexes formed by the other nanobodies followed the same pattern. Throughout the simulations, all systems maintained stable binding states, with RMSD fluctuations consistently below 6 Å ( Supplementary Figures S1A-D ). Based on the equilibrium portion of the MD trajectory (150–200 ns), the free energy landscape (FEL) analysis ( Supplementary Figure S2 ) was performed to explore the distribution of conformations and stability of Nb-RBD systems. In general, lower energy states indicate the corresponding conformations possess greater stability. Notably, the complexes formed between JN.1 and KP.3/XEC RBDs and the nanobodies exhibited increased flexibility at the binding interface compared to the complexes formed between the WT and BA.2 RBDs and the nanobodies. This enhanced flexibility suggests more dynamic interactions, potentially affecting binding characteristics and antibody efficacy.

To assess the effect of mutations on binding modes and energy changes, we calculated the binding free energies using the MM/PBSA method. The binding affinities of R14 and DL4 were significantly reduced, with the R14-KP.3/XEC having a binding affinity of -45.41 kcal/mol, representing increases of 36.94 kcal/mol and 39.85 kcal/mol, respectively, compared to their WT systems ( Figures 3B, C ). In contrast, the reduction in binding affinity for V_H ab6 and Nanosota9 was more gradual, with increases of 13.96 kcal/mol and 3.99 kcal/mol, respectively, compared to their WT complexes, consistent with the analysis in section 3.1 regarding the antiviral spectrum of nanobodies ( Figures 3D, E ). These findings were further validated through MM/GBSA cross-verification ( Supplementary Table S1 ), confirming the robustness of our predictions. Although the absolute binding energies varied between methods due to differences in the solvation model, both approaches showed consistent trends in changes to binding affinity across variants.

Overall, the binding affinities of JN.1 and KP.3/XEC variants to nanobodies were markedly lower than those of WT and BA.2 variants. This discrepancy primarily arises from a significant reduction in electrostatic and van der Waals (vdW) energies ( Figure 3A ). The former is attributed to the introduction of positively charged residues by mutations, while the latter results from unstable binding modes between the mutated RBD and nanobodies, leading to fewer atomic contacts. In addition, hydrogen bond network analysis ( Supplementary Figure S3 ) revealed that the JN.1 and KP.3/XEC variants had more low-occupancy hydrogen bonds (occupancy < 70%) and fewer high-occupancy hydrogen bonds (occupancy ≥ 70%), which further impacted the stability of the complex.

Subsequently, we calculated the binding free energy contributions of residues near the binding interface using residue decomposition methods. Residues with an absolute binding free energy difference ≥ 1 kcal/mol compared to the WT system were defined as hotspot residues. In the R14-KP.3/XEC complex, a greater number of hotspot residues were identified than in the R14-JN.1 complex ( Figure 4A , Supplementary Figures S4A, B ), likely due to the unique F456L and Q493E mutations in KP.3/XEC. Notably, the Q493E mutation in the R14-KP.3/XEC complex exhibited a substantial energy change, attributed to the transformation from neutral glutamine to negatively charged glutamate at physiological pH, enhancing electrostatic repulsion between nanobody and RBD. This mutation also disrupted existing hydrogen bonds, further destabilizing the complex ( Supplementary Figure S4 ). Similarly, the charged mutation at residue 493 significantly decreased electrostatic interactions in the DL4, V_H ab6, and Nanosota9 systems ( Figure 3A ). In the DL4-KP.3/XEC complex, residues F489, G495, and R497 exhibited notable energy reductions, clustering in a similar region, suggesting substantial conformational changes. The V_H ab6 and Nanosota9 complexes displayed comparable binding modes, primarily influenced by the Q493E mutation in KP.3/XEC ( Figures 4C, D , Supplementary Figures S4E-H ). However, the charged residues introduced by the mutations did not approach the binding interface closely enough to significantly impact binding affinity compared to R14 and DL4.

The sixteen nanobody-RBD systems were analyzed using DDMut-PPI to predict high-affinity single-point mutations targeting CDRs and hotspot residues ( Supplementary Tables S2 - 3 ). This computational affinity maturation approach generated a virtual library containing 2,508 mutations for R14, 1,824 for DL4, 2,432 for V_H ab6, and 2,508 for Nanosota9. Screening for mutations with predicted ΔΔG > 0 (indicating enhanced affinity and stability, Supplementary Figure S6 ) identified 4 beneficial mutations for R14, 5 for DL4, 15 for V_H ab6, and 9 for Nanosota9 ( Figure 5 , Supplementary Table S4 ).

Moreover, mutations within specific CDRs and hotspot residues not only alter the binding characteristics at the mutation sites but also affect the binding modes of other CDRs and framework regions. Analysis of combinatorial mutations revealed substantial variations in binding energy, with the multi-points mutant library containing 28 high-affinity variants for R14, 40 for DL4, 13,748 for V_H ab6, and 176 for Nanosota9. Notably, V_H ab6 displayed markedly greater mutational plasticity, accommodating significantly more viable multi-points mutation than other nanobodies while maintaining stable binding conformations.

High-affinity substitution mutations identified through preliminary screening using DDMut-PPI, including both single and multiple mutations, were introduced into the nanobodies ( Supplementary Table S5 ). For R14 variants, docking against four SARS-CoV-2 variant RBDs using HDOCK identified three top-ranking mutation combinations exhibiting enhanced interaction ( Supplementary Table S6 ).

To validate the feasibility of this nanobody repurposing pipeline, we conducted molecular dynamics simulations on the complexes of three optimized variants of R14 with four different SARS-CoV-2 variants to assess structural stability and energy changes following mutation. The RMSD was utilized to determine the conformational changes occurring in the protein backbone throughout the simulation, indicating the dynamic stability of the complex. Although the RMSD of the optimized structure was higher than that of the initial structure, reflecting internal structural changes that occurred during the optimization, it remained stable throughout the simulation. The RMSD fluctuations of RBD and R14 mutated nanobody complex remained consistently below 6 Å ( Supplementary Figure S6 ), induced by minimal fluctuations and suggesting stable binding. The increased RMSD post-optimization, while stable, indicates a reorganization of the protein structure that may enhance the interaction with the target, thus aligning with the improved binding energy observed in the analysis.

Binding affinity serves as a key metric for evaluating the strength and specificity of nanobody-antigen interactions. We designate the complexes formed by the original R14 nanobody and its three optimized R14 as R14-ORI, R14-OPT1 (L29W/S52C/A101V), R14-OPT2 (L29W/A101V), and R14-OPT3 (L29W/S52C/A101L), respectively. Comprehensive binding free energy analysis ( Figure 6 ) revealed substantial improvements in binding affinity against all SARS-CoV-2 variants mentioned above in this study, following nanobody optimization. The R14-ORI exhibited progressively weaker binding from WT (-82.35 kcal/mol) to KP.3/XEC (-45.41 kcal/mol), consistent with evolving immune evasion. Optimized Nb-RBD complexes demonstrated markedly enhanced and more uniform affinities. R14-OPT1 showed particularly strong JN.1 binding affinity, with a free energy of -76.88 kcal/mol compared to -47.3 kcal/mol for R14-ORI, while maintaining stable WT interactions at -67.04 kcal/mol. R14-OPT2 achieved exceptional KP.3/XEC binding affinity at -78.09 kcal/mol with consistent performance across variants. R14-OPT3 displayed optimal BA.2 affinity at -79.47 kcal/mol despite a moderate reduction for JN.1 to -55.69 kcal/mol. Notably, the binding free energy for the R14-OPT2-KP.3/XEC system exhibited slight discrepancies compared to the DDMut-PPI and HDOCK predictions, which can be attributed to the stronger stability of this system, resulting in lower overall perturbations in its structure compared to other similar systems.

Furthermore, R14-OPT1 exhibited the most balanced binding characteristics, with affinities ranging from -64.49 kcal/mol (KP.3/XEC) to -76.88 kcal/mol (JN.1), a variation of 16.1% compared to the 44.9% variation observed for R14-ORI (-45.41 to -82.35 kcal/mol). These energy landscapes suggest R14-OPT1 represent the optimal compromise between variant coverage and binding potency, particularly given its superior JN.1 recognition (-76.88 kcal/mol versus -47.3 kcal/mol for R14-ORI). Overall, the binding affinities of the optimized R14 nanobody with these four variants were comparatively favorable. These findings highlight the potential of integrating computational methods for affinity maturation in the development of effective therapeutics against SARS-CoV-2.