Complementary Roles of Physics-Based Approaches in Predicting VHH Thermal Stability

Complementary Roles of Physics-Based Approaches in Predicting VHH Thermal Stability | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Article Complementary Roles of Physics-Based Approaches in Predicting VHH Thermal Stability Yusuke Tomimoto, Rika Yamazaki, Hiroki Shirai This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-7555873/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract VHHs are frequently employed in protein therapeutics; however, enhancing their thermal stability (Tm) remains a significant challenge. We previously developed a two‑step in silico strategy to enhance Tm by identifying mutation sites (first round) and selecting favorable substitutions (second round) using dStability, a Gibbs free energy-based score. While dStability performed reasonably in the second round ( r = -0.64), it failed to predict the effect of G97F mutation, ranking it as most stabilizing despite its strong destabilizing impact on experimental Tm. In this study, we re-evaluated this dataset using high-temperature molecular dynamics (MD) simulations, employing Q-values, originally proposed by Bekker et al., as a quantitative metric to assess the extent of structural degradation. This method successfully identified G97F as least stable, demonstrating utility for detecting destabilizing mutations for the first round, though its second‑round performance was weaker ( r = 0.44). We also developed a temperature‑increase MD protocol, ramping simulation temperature from 300 K to 1000 K over 4 ns to completely degrade protein. Despite being 25-fold faster than fixed-temperature simulations, this approach retained comparable predictive performance. Overall, combining these two physics-based approaches, dStability and Q‑value analysis with optimized MD protocols, enables efficient identification of stabilizing mutations. Biological sciences/Biophysics Biological sciences/Computational biology and bioinformatics Biological sciences/Drug discovery Physical sciences/Mathematics and computing Figures Figure 1 Figure 2 Figure 3 Introduction Single-domain antibodies (sdAbs), also known as VHHs (Variable domain of Heavy chain of Heavy-chain antibody), are attracting significant attention as therapeutic and diagnostic agents owing to their small size, high solubility, and capacity to refold after heat denaturation [ 1 , 2 ]. Nonetheless, thermal stability - usually quantified by the melting temperature (Tm) - remains a critical determinant of their developability, especially under stress conditions during manufacturing, storage or transportation [ 3 ]. Various approaches have been explored to enhance the thermal stability of sdAbs, including disulfide bond engineering [ 4 – 9 ], CDR grafting [ 10 – 12 ], consensus sequence optimization [ 13 – 15 ], and mutagenesis [ 16 – 18 ]. While some of these methods have shown success in raising Tm, they often come at the cost of impaired antigen affinity [ 12 ] or decreased refolding efficiency [ 8 ]. In this context, computational approaches offer a promising and cost-effective complement to experimental screening. In our previous study [ 19 ], we proposed a two-step in silico strategy to improve the thermal stability of a model anti-lysozyme VHH, D3-L11. This approach combined energy-based ranking with targeted experimental validation. In the first round, we computed dStability scores (MOE [ 20 ], a molecular modeling software package, used to estimate ΔΔG) for all possible single mutations and selected a limited number of candidate sites predicted to yield the greatest Tm improvements. Experimental evaluations were then performed using differential scanning fluorometry (DSF), which identified A79 as the most promising mutation site. Surface plasmon resonance (SPR) and enzyme-linked immunosorbent assay (ELISA) were used to evaluate binding affinity. In the second round, we exhaustively substituted the most promising site, A79, with all other amino acids to identify optimal stabilizing mutations. This strategy successfully identified mutants such as A79I, which improved Tm by more than 5°C while retaining antigen-binding affinity. dStability performed reasonably well in second round (Pearson’s correlation coefficient, r = -0.64), but in first round, its predictive performance was limited. For example, the G97F mutation - which was ranked as the most stabilizing candidate based on dStability - resulted in a significant decrease in experimental Tm, representing a critical misprediction. This failure highlighted a key limitation of relying solely on static energy-based estimates for mutation site selection. A distinct physics-based approach was proposed by Bekker et al. [ 21 ], who assessed the thermal stability of sdAbs using high-temperature molecular dynamics (MD) simulations and the Q value, a structural metric quantifying the retention of native atomic contacts over time. Their study demonstrated a strong correlation ( r = 0.79) between Q values and experimental Tm values across five VHHs and two VH domains. Although their analysis focused exclusively on wild type structures, they proposed a set of stabilizing single- or double-point mutations based on per-residue Q value analysis. While these predictions were not experimentally tested in their own study, a subsequent investigation by another group later confirmed that the designed mutations indeed improved thermal stability [ 22 ] - supporting the utility of the Q value framework for rational stabilization design. To address the limitations of dStability observed in our previous work - particularly its poor performance in mutation site selection - we investigated the applicability of the Q value-based approach to mutant ranking. We found that the Q-value-based approach is better suited for the first round, whereas dStability performs better in the second round. In experimental methods, Tm measurements involve completely melting the antibody, whereas in this isothermal high-temperature MD, it was confirmed that the system remained only in a partially unfolded state. Therefore, in order to fully ‘melt’ the antibody in silico as in experiments, we additionally performed simulations under a condition where the temperature was ramped from 300 K to 1000 K over 4 ns. This approach was found to improve computational efficiency without significantly compromising predictive accuracy. Results and Discussion In this study, we aim to evaluate the performance of Q-value measurements based on molecular dynamics (MD) simulations [ 21 ] by directly using the experimental data reported in our previous study [ 19 ], and to compare the results with those obtained from dStability analysis. Therefore, we begin with a brief review of our previous work. In the first round of our prior approach [ 19 ], dStability scores were computed for all single amino acid substitutions, and the five top-ranked sites were selected: A24L, G26I, A79F, G97F, and G118F. Among these, G118F could not be tested due to cloning failure. While A24L, G26I, and A79F showed moderate increases in experimental Tm (≥ 1.5°C), G97F - despite being ranked as the most stabilizing mutation by dStability - exhibited a substantial reduction in thermal stability, highlighting a critical misprediction. In the second round, A79 was selected for residue scanning based on its strong performance (A79F: +7.25°C). Sixteen additional mutations were introduced at this position, several of which (e.g., A79I, A79W, A79Y, A79C) also led to Tm increases greater than 5°C. The dStability scores for these A79 mutants correlated well with experimental Tm values ( r = 0.64), confirming the utility of the method for residue prioritization once a mutation site is specified. Taken together, these results demonstrated that while dStability is useful for fine‑tuning substitutions at a given site (second round), its site‑selection performance in the initial round was suboptimal - particularly for outliers such as G97F. To explore a complementary structural metric for predicting thermal stability, we revisited the Q‑value method originally proposed by Bekker et al. [ 21 ] for our dataset. In their study, high‑temperature MD simulations at 400 K were used to calculate Q values, defined as the fraction of native atomic contacts retained during simulation. Evaluation of Residue Contact Schemes (hydrophilic-hydrophilic vs hydrophilic-all) in Isothermal MD Simulations (400 K, 100 ns) Bekker et al. classified amino acids into three physicochemical groups - hydrophobic, hydrophilic, and small - and computed Q values across all ten possible residue‑contact pair types (hydrophilic–hydrophilic, hydrophilic-hydrophobic, hydrophilic-small, hydrophilic-all, hydrophobic-hydrophobic, hydrophobic-small, hydrophobic-all, small-small, small-all, and all-all). They showed that the hydrophilic-all Q value exhibited the highest correlation with experimental Tm, followed by all-all contacts. In this study, we first examined which of the two top contact schemes – hydrophilic-all versus all-all provides better predictive power. To this end, we employed a small set of representative mutants: A79F, the most stabilizing mutation in the first round; G97F, the most destabilizing; and the wild type. For each mutant, we performed triplicate MD simulations (n = 3) at 400 K for 100 ns, and computed Q values over the last 30 ns of the trajectories using both hydrophilic–all and all–all contact definitions. Contrary to the trend reported by Bekker et al., our results indicated that the all-all Q values showed a stronger correlation with the experimentally determined Tm in our dataset (data not shown). These findings suggest that the optimal contact scheme for Q value calculation may vary depending on the structural context and mutation set. We then extended this comparison to the second-round mutants. Specifically, we selected five VHH mutants: A79I and A79W, which exhibited the largest increases in Tm, A79R and A79Q, which showed the most pronounced decreases; and the wild type. For each mutant, we performed MD simulations under the same conditions (400 K, 100 ns, n = 3) and computed Q values over the final 30 ns using both hydrophilic-all and all-all contact definitions. Consistent with the first-round results, the all-all Q values again showed a stronger correlation with experimental Tm than the hydrophilic-all Q values in this second-round mutant set (data not shown). This reinforced the applicability of the all-all contact definition for analyzing our VHH system and provided the basis for applying the same contact definition in subsequent temperature-increase MD simulations. Isothermal MD Simulations for first-round Mutants Based on our preliminary evaluation of residue contact schemes, we adopted the all-all contact definition for subsequent Q value calculations. We then performed a more comprehensive analysis of the full first-round mutant set, which comprised four single-point mutants (A24L, G26I, A79F, and G97F) along with the wild type VHH. For each of these five proteins, we conducted molecular dynamics simulations at 400 K for 100 ns in ten independent replicates (n = 10). Q values were calculated using the all-all contact definition over the final 30 ns of each trajectory. Figures 1 a-e show the Q value trajectories over time, with Fig. 1 a showing the wild type alone, and Figs. 1 b-e overlaying each mutant with the wild type. All proteins exhibited a gradual decrease in Q values, indicative of progressive unfolding. Between 70 and 100 ns, Q values generally stabilized within the range of 0.7 to 0.9, suggesting partially unfolded conformations. Figure 2 a presents the distributions of Q values as box-and-whisker plots, which reveal substantial variation among the different mutants. Notably, mutants such as A79F and G97F showed broader distributions than the others. These observations suggest that mutations involving larger amino acid residues may introduce greater variability in the unfolding behavior during partial denaturation. One-way ANOVA with Tukey’s multiple comparisons indicated that the mean Q value of G97F was significantly lower than those of the other mutants. This aligns with experimental findings showing that G97F was the only first-round mutant to reduce Tm - an outcome that was not predicted by dStability but was clearly captured by Q value analysis. Figure 2 b summarizes the relationship between the average Q value and the experimentally measured Tm, yielding a Pearson correlation coefficient of r = 0.59. In contrast, the correlation between dStability scores and Tm was weaker and in the opposite direction ( r = 0.42). Given that dStability is expected to exhibit a negative correlation with Tm, the observed positive correlation is problematic in itself. While both correlations are based on a limited sample size (n = 5) and should therefore be interpreted with caution, the results suggest that Q values derived from isothermal MD simulations may provide a more reliable indicator of thermal stability than static energy-based metrics for mutation site selection. Isothermal MD Simulations for second-round Mutants We next applied the same protocol to the second-round mutant set, consisting of seventeen A79X mutants along with the wild type, all evaluated under identical simulation conditions (400 K, 100 ns, n = 10). Q values were also computed over the final 30 ns using the all-all contact definition. The distribution of Q values is shown as box-and-whisker plots in Fig. 2 c, illustrating that the degree of variation again differed across proteins. Five mutants - A79F, A79Y, A79V, A79E, and A79R - exhibited particularly large variability. Among these, three mutations (F, Y, and R) involve amino acids with bulky side chains, which may echo the trend observed in first round, where mutations to larger residues were associated with greater variability in partial unfolding. However, this relationship does not hold consistently: for example, tryptophan (W), the largest amino acid, showed only modest variation, while medium-sized residues like valine (V) and glutamate (E) exhibited relatively high variability. Thus, no clear correlation was established between side-chain size and unfolding variability. One-way ANOVA with Tukey’s multiple comparisons revealed that the Q value of the A79E mutant was significantly lower than those of nine other mutants (A79C, A79I, A79W, A79T, A79L, A79S, A79N, A79V, A79H) in the second-round mutants. Similarly, the A79R mutant showed statistically significant differences in Q value compared to six other mutants (79C, A79I, A79W, A79L, A79N, A79H) in the second-round mutants. Both A79E and A79R resulted in decreased experimental Tm values, with A79R being the most destabilizing among all second-round mutants. These findings indicate that Q value analysis can effectively identify markedly destabilizing mutations even in the context of specific site residue scanning. We then assessed the correlation between the average Q value and the experimentally measured Tm across all 17 mutants and wild type. As shown in Fig. 2 d, the Pearson correlation coefficient was r = 0.44, which is notably lower than that observed for the dStability score ( r = − 0.64, inverse correlation). Taken together, these results suggest a complementary relationship between the two metrics: Q value-based MD analysis outperformed dStability in first round (mutation site selection), while dStability showed better predictive power in second round (amino acid selection at a specific site). This indicates that combining both approaches may provide a more balanced and effective strategy for rational antibody stabilization. Temperature-Increase MD Simulations for First-Round Mutants While experimental approaches for detecting Tm, such as DSC, completely disrupt protein structure, isothermal simulations above 400 K correspond only to partially unfolded states. Moreover, their reliance on long isothermal trajectories imposes a significant computational cost. To better capture structural behavior during complete unfolding, we aimed to establish a protocol that both more closely mimics experimental conditions and is more time-efficient while maintaining predictive performance. To this end, we devised a short (4 ns) temperature-ramp MD protocol, in which the simulation temperature is linearly increased from 300 K to 1000 K. This design allows monitoring of native-to-unfolded transitions within a compact simulation window. Q values were calculated across the entire trajectory to assess structural retention throughout the unfolding process. We applied this protocol to the first-round mutant set, comprising four single-point mutants (A24L, G26I, A79F, and G97F) and the wild type. Each mutant was simulated in ten independent runs (n = 10), and Q values were computed using the all-all contact definition. Figures 1 f-j display the Q value trajectories over time, with Fig. 1 f showing the wild type alone, and Figs. 1 g-j overlaying each mutant with the wild type. All proteins underwent complete structural degradation during the simulation as expected. Box-and-whisker plots of the resulting Q values are shown in Fig. 3 a. As in the isothermal case, considerable variability was observed across mutants, particularly for A24L and G97F. However, unlike the isothermal results, no clear relationship was found between side-chain size and Q value variability. This contrast with the isothermal protocol likely reflects the different unfolding regimes captured: under isothermal conditions, partial unfolding may be more sensitive to steric hindrance introduced by bulky side chains, while the temperature-increase protocol induces full denaturation in all proteins, thereby minimizing size-dependent effects. Under these simulation conditions, the G97F mutant still exhibited the lowest average Q value; however, unlike in the isothermal protocol, G97F did not show statistically significant differences from the wild type or other mutants. Figure 3 b presents the correlation between average Q values and experimental Tm, yielding a strong Pearson correlation coefficient of r = 0.85. This result exceeds the correlation observed with isothermal simulations at 400 K ( r = 0.59), despite requiring only 1/25th of the computational time. These findings demonstrate that the temperature-increase protocol offers a highly efficient and accurate approach for mutation site selection. Temperature-Increase MD Simulations for Second-Round Mutants We next applied the same temperature-increase protocol to the second-round mutant set, which included seventeen A79X mutants and the wild type. Each mutant was simulated in ten independent runs (n = 10), and Q values were computed over the full simulation span using the all-all contact definition. Figure 3 c shows the resulting Q value distributions. As observed previously, variability differed among mutants. The most variable mutants were A79N, A79M, and A79Q. However, consistent with first-round results, no clear relationship was observed between side-chain size and variability. Furthermore, unlike the isothermal simulations, no statistically significant differences were observed among the mutants, including the destabilizing A79R mutant. As in first round, Fig. 3 d shows the correlation between average Q values and experimental Tm, yielding a Pearson correlation coefficient of r = 0.40. This is nearly identical to that observed with isothermal MD at 400 K ( r = 0.44), indicating that both protocols perform comparably in the context of residue substitution ranking at a specific site. Notably, the A79R mutant, which exhibited the greatest Tm decrease, also showed the second-lowest Q value in this analysis. Taken together, these findings demonstrate that the temperature-increase MD protocol achieves predictive performance comparable to that of conventional isothermal simulations with significantly reduced computational cost. Importantly, the specific protocol used in this study - linearly increasing the temperature from 300 K to 1000 K over 4 ns and averaging Q values across the full simulation window - represents just one implementation designed to capture the full transition from the native state to complete denaturation. It is likely that further optimization of simulation parameters, including total MD duration, temperature range, and the time window used for Q value calculation, could lead to even more accurate and efficient predictions in future applications. Recently, a curated database of VHH sequences and their Tm values, NbThermo [ 23 ], has been developed. Furthermore, machine learning methods such as NanoMelt [ 24 ] and TEMPRO [ 25 ], which partially utilize NbThermo dataset, have also been introduced. In this study, we demonstrated that the two physics-based approaches, Q-value analysis and dStability, are complementary. Moving forward, we plan to investigate the potential of combining these physics-based methods with AI-based approaches to further enhance predictive performance. Conclusion Our study demonstrates that the Q value-based MD approach is more effective than dStability for mutation site selection (first round), whereas dStability remains superior for evaluating residue substitutions at a specific site (second round) in our dataset. Furthermore, the temperature-increase MD protocol achieved comparable predictive performance to the longer isothermal simulations, despite requiring significantly less computational time. These findings suggest that a combined application of Q value-based MD analysis and dStability scoring - along with carefully optimized MD protocols - could offer a powerful and efficient framework for improving the thermal stability of VHHs. Materials and Methods In the preparations and MD simulations, GROMACS 2018 [ 26 ] was used. VHH mutants previously reported to yield successful Tm measurements were selected [ 19 ]. All VHH mutants were modeled using the crystal structure 6JB9 as a template. The PDB structure of 6JB9 contained no missing residues; however, an additional N-terminal serine (Ser0), resulting from the difference of expression host, was present. This residue was removed, and the N-terminus was capped with an N-methyl (NME) group. The modeling was conducted using the MOE software (Chemical Computing Group) [ 20 ], and the mutations were introduced using the residue mutation function. Side chains at the mutation sites (e.g., A24L, G26I, A79F, G97F in first round; and A79Y, A79W, A79I, A79L, A79V, A79C, A79N, A79Q, A79M, A79H, A79R, A79E, A79K, A79P in second round) were automatically rebuilt by MOE based on its internal rotamer library, and no further energy minimization was applied. For buried mutation sites such as A79, side chain rotamers were manually selected based on the steric compatibility and hydrogen bonding environment when necessary. To prepare the system, hydrogens were attached, and a cubic box was constructed surrounding each protein such that the box was at least 1.0 nm away from the protein and filled with water molecules and 0.153 M of NaCl. Protonation states were assigned assuming a neutral pH (~ 7.0) using the MOE Protonate 3D module. Histidine residues were manually reviewed and assigned as HID (proton on δ-nitrogen), HIE (proton on ε-nitrogen), or HIP (doubly protonated) based on local hydrogen bonding environments. The Amber99-sb-ildn force field [ 27 – 29 ] was used to parameterize the protein, TIP3P [ 30 ] for waters and the monovalent ion parameters [ 31 ] were used for the ions. Our simulation protocol commenced with energy minimization, followed by position restraints on the heavy solute atoms during NVT for the first step and NPT for the second step. 3 or 10 parallel simulations were conducted, each initialized with distinct random seed velocities. To equilibrate velocities, a 100 ps NVT simulation at 298K was performed, which allowed a subsequent 100 ps NPT simulation at the same temperature. For these equilibration simulations, position restraints were used on the heavy solute atoms. For the “isothermal” production run, an unrestrained NVT simulation was performed at 400K, spanning 100 ns. For the “temperature-increase” production run, an unrestrained NVT simulation was performed with temperature increasing from 298K to 1000K, spanning 4 ns. A 1.0-nm cutoff was used for the short-range interactions in combination with the Particle Mesh Ewald method for the long-range electrostatic interactions. A consistent 2 fs timestep was maintained, with constraints managed using LINCS [ 32 ] for the protein and water molecules. Temperature was controlled using the V-rescale thermostat [ 33 ] and pressure was maintained using Parrinello-Rahman barostat [ 34 ] at 1 bar during the NPT simulation. During the production run, structures were saved at 10-ps intervals. Q values were calculated following the definition of Best et al. [ 35 ], using a contact threshold of 0.45 nm between non-hydrogen atoms. Residues within the N-terminal cap (NME) and residues at the flexible C-terminal region (last 3 residues) were excluded from the calculation. Based on the previous study [ 21 ], two types of Q-value calculation methods that showed the highest correlation with Tm were defined: all-all (all residue pairs) and HA (hydrophilic-all), with hydrophilic residues defined as D, E, N, Q, R, K, and H. The average Q-value under the isothermal condition was calculated from the final 30 ns of the entire trajectory, as in the previous study [ 21 ], whereas for the temperature-increase condition, the average Q-value was calculated over the entire 4 ns trajectory. Statistical analysis were performed using Microsoft Excel with Real Statistics Add-In. For Q value comparisons, one-way ANOVA with Tukey’s multiple comparisons was used. Declarations Competing interests statement The authors declare no competing interests. Funding statement This study was carried out with the support of Astellas Pharma Inc. Author Contribution Y.T. and H.S. designed the study. R.Y. conducted the experiments, Y.T., R.Y. and H.S. analyzed the results. Y.T., R.Y. and H.S. wrote and approved the manuscript. Acknowledgement We extend our deepest gratitude to Ms. Nami Tanigawa (WORLD INTEC), whose expertise and dedication in proposing conditions, conducting molecular dynamics simulations, and providing thoughtful data analyses were indispensable to this work. Data Availability The datasets generated and/or analyzed during the current study are available from the corresponding author on reasonable request. References Hamers-Casterman, C. et al. Naturally occurring antibodies devoid of light chains. Nature 363 , 446–448. 10.1038/363446a0 (1993). Könning, D. et al. Camelid and shark single domain antibodies: structural features and therapeutic potential. Curr. Opin. Struct. Biol. 45 , 10–16. https://doi.org/10.1016/j.sbi.2016.10.019 (2017). Goldman, E. R., Liu, J. L., Zabetakis, D. & Anderson, G. P. Enhancing Stability of Camelid and Shark Single Domain Antibodies: An Overview. Front. Immunol. 8 , 865. 10.3389/fimmu.2017.00865 (2017). Hagihara, Y., Mine, S. & Uegaki, K. Stabilization of an immunoglobulin fold domain by an engineered disulfide bond at the buried hydrophobic region. J. Biol. Chem. 282 , 36489–36495. 10.1074/jbc.M707078200 (2007). Saerens, D., Conrath, K., Govaert, J. & Muyldermans, S. Disulfide bond introduction for general stabilization of immunoglobulin heavy-chain variable domains. J. Mol. Biol. 377 , 478–488. 10.1016/j.jmb.2008.01.022 (2008). Arbabi-Ghahroudi, M. et al. Aggregation-resistant VHs selected by in vitro evolution tend to have disulfide-bonded loops and acidic isoelectric points. Protein Eng. Des. Sel. 22 , 59–66. 10.1093/protein/gzn071 (2009). Zabetakis, D., Olson, M. A., Anderson, G. P., Legler, P. M. & Goldman, E. R. Evaluation of disulfide bond position to enhance the thermal stability of a highly stable single domain antibody. PLoS One . 9 , e115405. 10.1371/journal.pone.0115405 (2014). Akazawa-Ogawa, Y., Uegaki, K. & Hagihara, Y. The role of intra-domain disulfide bonds in heat-induced irreversible denaturation of camelid single domain VHH antibodies. J. Biochem. 159 , 111–121. 10.1093/jb/mvv082 (2016). Shinozaki, N., Hashimoto, R., Fukui, K. & Uchiyama, S. Efficient generation of single domain antibodies with high affinities and enhanced thermal stabilities. Sci. Rep. 7 , 5794. 10.1038/s41598-017-06277-x (2017). Saerens, D. et al. Identification of a universal VHH framework to graft non-canonical antigen-binding loops of camel single-domain antibodies. J. Mol. Biol. 352 , 597–607. 10.1016/j.jmb.2005.07.038 (2005). Vincke, C. et al. General strategy to humanize a camelid single-domain antibody and identification of a universal humanized nanobody scaffold. J. Biol. Chem. 284 , 3273–3284. 10.1074/jbc.M806889200 (2009). Zabetakis, D., Anderson, G. P., Bayya, N. & Goldman, E. R. Contributions of the complementarity determining regions to the thermal stability of a single-domain antibody. PLoS One . 8 , e77678. 10.1371/journal.pone.0077678 (2013). Anderson, G. P., Liu, J. H., Zabetakis, D., Liu, J. L. & Goldman, E. R. Thermal stabilization of anti-alpha-cobratoxin single domain antibodies. Toxicon 129 , 68–73. 10.1016/j.toxicon.2017.02.008 (2017). Shriver-Lake, L. C., Zabetakis, D., Goldman, E. R. & Anderson, G. P. Evaluation of anti-botulinum neurotoxin single domain antibodies with additional optimization for improved production and stability. Toxicon 135 , 51–58. 10.1016/j.toxicon.2017.06.002 (2017). Kunz, P. et al. Exploiting sequence and stability information for directing nanobody stability engineering. Biochim. Biophys. Acta Gen. Subj. 1861 , 2196–2205. 10.1016/j.bbagen.2017.06.014 (2017). Harmsen, M. M., van Solt, C. B., van Zijderveld-van Bemmel, A. M., Niewold, T. A. & van Zijderveld, F. G. Selection and optimization of proteolytically stable llama single-domain antibody fragments for oral immunotherapy. Appl. Microbiol. Biotechnol. 72 , 544–551. 10.1007/s00253-005-0300-7 (2006). Turner, K. B., Zabetakis, D., Goldman, E. R. & Anderson, G. P. Enhanced stabilization of a stable single domain antibody for SEB toxin by random mutagenesis and stringent selection. Protein Eng. Des. Sel. 27 , 89–95. 10.1093/protein/gzu001 (2014). Hussack, G. et al. Protease-resistant single-domain antibodies inhibit Campylobacter jejuni motility. Protein Eng. Des. Sel. 27 , 191–198. 10.1093/protein/gzu011 (2014). Tomimoto, Y., Yamazaki, R. & Shirai, H. Increasing the melting temperature of VHH with the in silico free energy score. Sci. Rep. 13 , 4922. 10.1038/s41598-023-32022-8 (2023). Molecular Operating Environment (MOE), 2019.01; Chemical Computing Group ULC, 1010 Sherbooke St. West, Suite #910, Montreal, QC & Canada H3A 2R7, (2021). Bekker, G. J., Ma, B. & Kamiya, N. Thermal stability of single-domain antibodies estimated by molecular dynamics simulations. Protein Sci. 28 , 429–438. 10.1002/pro.3546 (2019). Zabetakis, D., Shriver-Lake, L. C., Olson, M. A., Goldman, E. R. & Anderson, G. P. Experimental evaluation of single-domain antibodies predicted by molecular dynamics simulations to have elevated thermal stability. Protein Sci. 28 , 1909–1912. 10.1002/pro.3692 (2019). Valdes-Tresanco, M. S., Valdes-Tresanco, M. E., Molina-Abad, E. & Moreno, E. NbThermo: a new thermostability database for nanobodies. Database (Oxford) (2023). 10.1093/database/baad021 (2023). Ramon, A. et al. Prediction of protein biophysical traits from limited data: a case study on nanobody thermostability through NanoMelt. MAbs 17, 2442750, (2025). 10.1080/19420862.2024.2442750 Alvarez, J. A. E. & Dean, S. N. TEMPRO: nanobody melting temperature estimation model using protein embeddings. Sci. Rep. 14 , 19074. 10.1038/s41598-024-70101-6 (2024). Abraham, M. J. et al. High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 1–2 . GROMACS https://doi.org/10.1016/j.softx.2015.06.001 (2015). 19–25, doi. Abu-Soud, H. M., Clark, A. C., Francisco, W. A., Baldwin, T. O. & Raushel, F. M. Kinetic destabilization of the hydroperoxy flavin intermediate by site-directed modification of the reactive thiol in bacterial luciferase. J. Biol. Chem. 268 , 7699–7706. https://doi.org/10.1016/S0021-9258(18)53013-1 (1993). Hornak, V. et al. Comparison of multiple Amber force fields and development of improved protein backbone parameters. Proteins Struct. Funct. Bioinform. 65 , 712–725. https://doi.org/10.1002/prot.21123 (2006). Lindorff-Larsen, K. et al. Improved side-chain torsion potentials for the Amber ff99SB protein force field. Proteins Struct. Funct. Bioinform. 78 , 1950–1958. https://doi.org/10.1002/prot.22711 (2010). Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W. & Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79 , 926–935 (1983). Joung, I. S. & Cheatham, T. E. III Determination of alkali and halide monovalent ion parameters for use in explicitly solvated biomolecular simulations. J. Phys. Chem. B . 112 , 9020–9041 (2008). Hess, B. & P-LINCS: A Parallel Linear Constraint Solver for Molecular Simulation. J. Chem. Theory Comput. 4 , 116–122. 10.1021/ct700200b (2008). Bussi, G., Donadio, D. & Parrinello, M. Canonical sampling through velocity rescaling. The J. Chem. physics 126 (2007). Parrinello, M. & Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys. 52 , 7182–7190 (1981). Best, R. B., Hummer, G. & Eaton, W. A. Native contacts determine protein folding mechanisms in atomistic simulations. Proceedings of the National Academy of Sciences 110, 17874–17879 (2013). Additional Declarations No competing interests reported. Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-7555873","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Article","associatedPublications":[],"authors":[{"id":526170236,"identity":"37547567-676b-4ccc-9a6e-d9210565132b","order_by":0,"name":"Yusuke Tomimoto","email":"","orcid":"","institution":"Astellas Pharma Inc","correspondingAuthor":false,"prefix":"","firstName":"Yusuke","middleName":"","lastName":"Tomimoto","suffix":""},{"id":526170237,"identity":"bd3a0e1c-97ea-4936-a7d5-9554c689b899","order_by":1,"name":"Rika Yamazaki","email":"","orcid":"","institution":"Astellas Pharma Inc","correspondingAuthor":false,"prefix":"","firstName":"Rika","middleName":"","lastName":"Yamazaki","suffix":""},{"id":526170238,"identity":"122feaf5-66fc-4ca7-8d18-a6ff51cfce71","order_by":2,"name":"Hiroki Shirai","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA5ElEQVRIiWNgGAWjYBACA2YGBsYGAwiDIaGCgUcCJCxBvJYzxGhhAGmBMhgY2/CrBgNzdt6HD2cUMNiDGB8ezrsnI9nee4DBcgduLZbN7MaGGwwYEncCGRKJ24p5pHnOJTBInsHjsMNsbJIPDBgSgAwGoJYEHjmJHAMGyTbCWuyBDOYfiXOAWuTfEKEF6DDGDUCGRGJDAo+0BA9BLcyGMwwkEkFaLBKOJfBI9uQYHMDrl/PHGB/2/LGxBzKYb/6oSbCXOH7G8LEknhCDArTIOCzZQFALGmD8SLKWUTAKRsEoGMYAAMC+RV0EvBJHAAAAAElFTkSuQmCC","orcid":"","institution":"Astellas Pharma Inc","correspondingAuthor":true,"prefix":"","firstName":"Hiroki","middleName":"","lastName":"Shirai","suffix":""}],"badges":[],"createdAt":"2025-09-07 11:23:11","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-7555873/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-7555873/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":93143395,"identity":"522da7e0-5f5f-4d23-a8d1-6c96cbcae518","added_by":"auto","created_at":"2025-10-09 13:21:12","extension":"docx","order_by":0,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":54841,"visible":true,"origin":"","legend":"","description":"","filename":"Manuscript20250909.docx","url":"https://assets-eu.researchsquare.com/files/rs-7555873/v1/dc2888650a4c7c6d9136ea42.docx"},{"id":93144276,"identity":"08d82ddf-f57c-4592-b620-f641ad8b537c","added_by":"auto","created_at":"2025-10-09 13:29:12","extension":"tif","order_by":1,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":6130076,"visible":true,"origin":"","legend":"","description":"","filename":"Figure1.tif","url":"https://assets-eu.researchsquare.com/files/rs-7555873/v1/a08abc1db6c30d0c74b1fdf2.tif"},{"id":93144274,"identity":"97844f88-d8ae-4b25-aa9c-d73b992bc39b","added_by":"auto","created_at":"2025-10-09 13:29:12","extension":"tif","order_by":2,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":3072202,"visible":true,"origin":"","legend":"","description":"","filename":"Figure2.tif","url":"https://assets-eu.researchsquare.com/files/rs-7555873/v1/9ad94a7d1370e127dec69fda.tif"},{"id":93143400,"identity":"01687141-79f8-4bdb-bd0e-7ca7b1f08711","added_by":"auto","created_at":"2025-10-09 13:21:12","extension":"tif","order_by":3,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":3229462,"visible":true,"origin":"","legend":"","description":"","filename":"Figure3.tif","url":"https://assets-eu.researchsquare.com/files/rs-7555873/v1/701a874dec5f888788e8bab6.tif"},{"id":93143399,"identity":"fb76ee8f-6a9f-429d-b143-0f0d5e3ff42a","added_by":"auto","created_at":"2025-10-09 13:21:12","extension":"json","order_by":4,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":5165,"visible":true,"origin":"","legend":"","description":"","filename":"e0bf0eca7d8848619d1066345f2767a1.json","url":"https://assets-eu.researchsquare.com/files/rs-7555873/v1/1bdf40720037786ca189b150.json"},{"id":93144273,"identity":"aa6903e0-894e-4b3b-91c9-3ac40282e59b","added_by":"auto","created_at":"2025-10-09 13:29:12","extension":"xml","order_by":5,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":88789,"visible":true,"origin":"","legend":"","description":"","filename":"e0bf0eca7d8848619d1066345f2767a11enriched.xml","url":"https://assets-eu.researchsquare.com/files/rs-7555873/v1/4a6729db824d24729180ca63.xml"},{"id":93143411,"identity":"d38ea69b-adeb-41ba-8426-9e3c62cf9d03","added_by":"auto","created_at":"2025-10-09 13:21:12","extension":"tif","order_by":6,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":6130076,"visible":true,"origin":"","legend":"","description":"","filename":"Figure1.tif","url":"https://assets-eu.researchsquare.com/files/rs-7555873/v1/ce2628832be8f11b74f6ba34.tif"},{"id":93143410,"identity":"1837ad0e-0b48-4e6d-ab65-254038a99abb","added_by":"auto","created_at":"2025-10-09 13:21:12","extension":"tif","order_by":7,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":3072202,"visible":true,"origin":"","legend":"","description":"","filename":"Figure2.tif","url":"https://assets-eu.researchsquare.com/files/rs-7555873/v1/685191170aec666a31d582ec.tif"},{"id":93143408,"identity":"6684e627-ed7a-4015-a8f5-a91a63d67a4d","added_by":"auto","created_at":"2025-10-09 13:21:12","extension":"tif","order_by":8,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":3229462,"visible":true,"origin":"","legend":"","description":"","filename":"Figure3.tif","url":"https://assets-eu.researchsquare.com/files/rs-7555873/v1/f2ff35d57cbf463281639817.tif"},{"id":93143404,"identity":"29e74c1d-521b-4b6e-9eec-71d1376866ee","added_by":"auto","created_at":"2025-10-09 13:21:12","extension":"png","order_by":9,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":545928,"visible":true,"origin":"","legend":"","description":"","filename":"OnlineFigure1.png","url":"https://assets-eu.researchsquare.com/files/rs-7555873/v1/63e1b4cc4c5590b7d430636c.png"},{"id":93144571,"identity":"77e3ff4c-3661-4813-8ba9-b22602d1b255","added_by":"auto","created_at":"2025-10-09 13:37:12","extension":"png","order_by":10,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":246948,"visible":true,"origin":"","legend":"","description":"","filename":"OnlineFigure2.png","url":"https://assets-eu.researchsquare.com/files/rs-7555873/v1/70ffd12ef9e314e65a21c053.png"},{"id":93143409,"identity":"4256e119-6baa-4b2e-829a-cf46c073cbb7","added_by":"auto","created_at":"2025-10-09 13:21:12","extension":"png","order_by":11,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":236611,"visible":true,"origin":"","legend":"","description":"","filename":"OnlineFigure3.png","url":"https://assets-eu.researchsquare.com/files/rs-7555873/v1/ba67d2ceed65984939c5551d.png"},{"id":93143407,"identity":"f74fde60-adf2-42f6-95c2-9cd5c74e428c","added_by":"auto","created_at":"2025-10-09 13:21:12","extension":"xml","order_by":12,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":88767,"visible":true,"origin":"","legend":"","description":"","filename":"e0bf0eca7d8848619d1066345f2767a11structuring.xml","url":"https://assets-eu.researchsquare.com/files/rs-7555873/v1/631560278884fdf8d3490338.xml"},{"id":93143402,"identity":"0827694f-b00d-4b31-8b63-eb2474c8ca58","added_by":"auto","created_at":"2025-10-09 13:21:12","extension":"html","order_by":13,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":97914,"visible":true,"origin":"","legend":"","description":"","filename":"earlyproof.html","url":"https://assets-eu.researchsquare.com/files/rs-7555873/v1/59ceaa6d72aa7ce95dce73d0.html"},{"id":93144272,"identity":"71de0ab2-b039-480a-8b98-55da5dde5ec9","added_by":"auto","created_at":"2025-10-09 13:29:12","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":4188822,"visible":true,"origin":"","legend":"\u003cp\u003eQ value trajectories over time of wild type and first-round mutants under the isothermal condition and the temperature-increase condition.\u003c/p\u003e\n\u003cp\u003e(\u003cstrong\u003ea\u003c/strong\u003e) Q value trajectories of wild type under the isothermal condition.\u003c/p\u003e\n\u003cp\u003e(\u003cstrong\u003eb\u003c/strong\u003e) Overlaid Q value trajectories of wild type (grey) and A79F mutant (cyan) under the isothermal condition.\u003c/p\u003e\n\u003cp\u003e(\u003cstrong\u003ec\u003c/strong\u003e) Overlaid Q value trajectories of wild type (grey) and A24L mutant (orange) under the isothermal condition.\u003c/p\u003e\n\u003cp\u003e(\u003cstrong\u003ed\u003c/strong\u003e) Overlaid Q value trajectories of wild type (grey) and G26I mutant (green) under the isothermal condition.\u003c/p\u003e\n\u003cp\u003e(\u003cstrong\u003ee\u003c/strong\u003e) Overlaid Q value trajectories of wild type (grey) and G97F mutant (magenta) under the isothermal condition.\u003c/p\u003e\n\u003cp\u003e(\u003cstrong\u003ef\u003c/strong\u003e) Q value trajectories of wild type under the temperature-increase condition.\u003c/p\u003e\n\u003cp\u003e(\u003cstrong\u003eg\u003c/strong\u003e) Overlaid Q value trajectories of wild type (grey) and A79F mutant (cyan) under the temperature-increase condition.\u003c/p\u003e\n\u003cp\u003e(\u003cstrong\u003eh\u003c/strong\u003e) Overlaid Q value trajectories of wild type (grey) and A24L mutant (orange) under the temperature-increase condition.\u003c/p\u003e\n\u003cp\u003e(\u003cstrong\u003ei\u003c/strong\u003e) Overlaid Q value trajectories of wild type (grey) and G26I mutant (green) under the temperature-increase condition.\u003c/p\u003e\n\u003cp\u003e(\u003cstrong\u003ej\u003c/strong\u003e) Overlaid Q value trajectories of wild type (grey) and G97F mutant (magenta) under the temperature-increase condition.\u003c/p\u003e","description":"","filename":"Figure1.png","url":"https://assets-eu.researchsquare.com/files/rs-7555873/v1/e9558e2badd6dc77ec9003e0.png"},{"id":93143397,"identity":"0638e7e4-8f76-4a16-834e-1ecd2ccfb444","added_by":"auto","created_at":"2025-10-09 13:21:12","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":1140484,"visible":true,"origin":"","legend":"\u003cp\u003eQ value in wild type and each mutant under the isothermal condition.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(a\u003c/strong\u003e) Box-and-whisker plot of wild type and first-round mutants. * shows statistically significant difference (p\u0026lt;0.05).\u003c/p\u003e\n\u003cp\u003e(\u003cstrong\u003eb\u003c/strong\u003e) Correlation between average Q value and melting temperature of wild type and each first-round mutant. Pearson’s correlation coefficient was 0.59.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(c\u003c/strong\u003e) Box-and-whisker plot of wild type and second-round mutants. * and # show statistically significant difference (p\u0026lt;0.05) with A79E mutant and A79R, respectively.\u003c/p\u003e\n\u003cp\u003e(\u003cstrong\u003ed\u003c/strong\u003e) Correlation between average Q value and melting temperature of wild type and each second-round mutant. Pearson’s correlation coefficient was 0.44.\u003c/p\u003e","description":"","filename":"Figure2.png","url":"https://assets-eu.researchsquare.com/files/rs-7555873/v1/3b3d2ab4b97b0ef3a54da4e6.png"},{"id":93143394,"identity":"f8bf3aeb-61d3-4bc3-b71b-cd35991e062f","added_by":"auto","created_at":"2025-10-09 13:21:12","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":1107675,"visible":true,"origin":"","legend":"\u003cp\u003eQ value in wild type and each mutant under the temperature-increase condition.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(a\u003c/strong\u003e) Box-and-whisker plot of wild type and first-round mutants.\u003c/p\u003e\n\u003cp\u003e(\u003cstrong\u003eb\u003c/strong\u003e) Correlation between average Q value and melting temperature of wild type and each first-round mutant. Pearson’s correlation coefficient was 0.85.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(c\u003c/strong\u003e) Box-and-whisker plot of wild type and second-round mutants.\u003c/p\u003e\n\u003cp\u003e(\u003cstrong\u003ed\u003c/strong\u003e) Correlation between average Q value and melting temperature of wild type and each second-round mutant. Pearson’s correlation coefficient was 0.40.\u003c/p\u003e","description":"","filename":"Figure3.png","url":"https://assets-eu.researchsquare.com/files/rs-7555873/v1/2f0672846039b679f3576c68.png"},{"id":96077923,"identity":"155159eb-8024-475b-8a93-765b74a24d42","added_by":"auto","created_at":"2025-11-17 10:55:03","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":8031824,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7555873/v1/7906002b-fea0-4fc1-945e-5e878557ea97.pdf"}],"financialInterests":"No competing interests reported.","formattedTitle":"Complementary Roles of Physics-Based Approaches in Predicting VHH Thermal Stability","fulltext":[{"header":"Introduction","content":"\u003cp\u003eSingle-domain antibodies (sdAbs), also known as VHHs (Variable domain of Heavy chain of Heavy-chain antibody), are attracting significant attention as therapeutic and diagnostic agents owing to their small size, high solubility, and capacity to refold after heat denaturation [\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e, \u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e]. Nonetheless, thermal stability - usually quantified by the melting temperature (Tm) - remains a critical determinant of their developability, especially under stress conditions during manufacturing, storage or transportation [\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e]. Various approaches have been explored to enhance the thermal stability of sdAbs, including disulfide bond engineering [\u003cspan additionalcitationids=\"CR5 CR6 CR7 CR8\" citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e], CDR grafting [\u003cspan additionalcitationids=\"CR11\" citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e], consensus sequence optimization [\u003cspan additionalcitationids=\"CR14\" citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e], and mutagenesis [\u003cspan additionalcitationids=\"CR17\" citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e]. While some of these methods have shown success in raising Tm, they often come at the cost of impaired antigen affinity [\u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e] or decreased refolding efficiency [\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e]. In this context, computational approaches offer a promising and cost-effective complement to experimental screening.\u003c/p\u003e\u003cp\u003eIn our previous study [\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e], we proposed a two-step in silico strategy to improve the thermal stability of a model anti-lysozyme VHH, D3-L11. This approach combined energy-based ranking with targeted experimental validation. In the first round, we computed dStability scores (MOE [\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e], a molecular modeling software package, used to estimate ΔΔG) for all possible single mutations and selected a limited number of candidate sites predicted to yield the greatest Tm improvements. Experimental evaluations were then performed using differential scanning fluorometry (DSF), which identified A79 as the most promising mutation site. Surface plasmon resonance (SPR) and enzyme-linked immunosorbent assay (ELISA) were used to evaluate binding affinity. In the second round, we exhaustively substituted the most promising site, A79, with all other amino acids to identify optimal stabilizing mutations. This strategy successfully identified mutants such as A79I, which improved Tm by more than 5\u0026deg;C while retaining antigen-binding affinity. dStability performed reasonably well in second round (Pearson\u0026rsquo;s correlation coefficient, \u003cem\u003er\u003c/em\u003e = -0.64), but in first round, its predictive performance was limited. For example, the G97F mutation - which was ranked as the most stabilizing candidate based on dStability - resulted in a significant decrease in experimental Tm, representing a critical misprediction. This failure highlighted a key limitation of relying solely on static energy-based estimates for mutation site selection.\u003c/p\u003e\u003cp\u003eA distinct physics-based approach was proposed by Bekker et al. [\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e], who assessed the thermal stability of sdAbs using high-temperature molecular dynamics (MD) simulations and the Q value, a structural metric quantifying the retention of native atomic contacts over time. Their study demonstrated a strong correlation (\u003cem\u003er\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.79) between Q values and experimental Tm values across five VHHs and two VH domains. Although their analysis focused exclusively on wild type structures, they proposed a set of stabilizing single- or double-point mutations based on per-residue Q value analysis. While these predictions were not experimentally tested in their own study, a subsequent investigation by another group later confirmed that the designed mutations indeed improved thermal stability [\u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e] - supporting the utility of the Q value framework for rational stabilization design. To address the limitations of dStability observed in our previous work - particularly its poor performance in mutation site selection - we investigated the applicability of the Q value-based approach to mutant ranking. We found that the Q-value-based approach is better suited for the first round, whereas dStability performs better in the second round.\u003c/p\u003e\u003cp\u003eIn experimental methods, Tm measurements involve completely melting the antibody, whereas in this isothermal high-temperature MD, it was confirmed that the system remained only in a partially unfolded state. Therefore, in order to fully \u0026lsquo;melt\u0026rsquo; the antibody in silico as in experiments, we additionally performed simulations under a condition where the temperature was ramped from 300 K to 1000 K over 4 ns. This approach was found to improve computational efficiency without significantly compromising predictive accuracy.\u003c/p\u003e"},{"header":"Results and Discussion","content":"\u003cp\u003eIn this study, we aim to evaluate the performance of Q-value measurements based on molecular dynamics (MD) simulations [\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e] by directly using the experimental data reported in our previous study [\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e], and to compare the results with those obtained from dStability analysis. Therefore, we begin with a brief review of our previous work.\u003c/p\u003e\u003cp\u003eIn the first round of our prior approach [\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e], dStability scores were computed for all single amino acid substitutions, and the five top-ranked sites were selected: A24L, G26I, A79F, G97F, and G118F. Among these, G118F could not be tested due to cloning failure. While A24L, G26I, and A79F showed moderate increases in experimental Tm (\u0026ge;\u0026thinsp;1.5\u0026deg;C), G97F - despite being ranked as the most stabilizing mutation by dStability - exhibited a substantial reduction in thermal stability, highlighting a critical misprediction.\u003c/p\u003e\u003cp\u003eIn the second round, A79 was selected for residue scanning based on its strong performance (A79F: +7.25\u0026deg;C). Sixteen additional mutations were introduced at this position, several of which (e.g., A79I, A79W, A79Y, A79C) also led to Tm increases greater than 5\u0026deg;C. The dStability scores for these A79 mutants correlated well with experimental Tm values (\u003cem\u003er\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.64), confirming the utility of the method for residue prioritization once a mutation site is specified.\u003c/p\u003e\u003cp\u003eTaken together, these results demonstrated that while dStability is useful for fine‑tuning substitutions at a given site (second round), its site‑selection performance in the initial round was suboptimal - particularly for outliers such as G97F.\u003c/p\u003e\u003cp\u003eTo explore a complementary structural metric for predicting thermal stability, we revisited the Q‑value method originally proposed by Bekker et al. [\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e] for our dataset. In their study, high‑temperature MD simulations at 400 K were used to calculate Q values, defined as the fraction of native atomic contacts retained during simulation.\u003c/p\u003e\u003cp\u003e\u003cb\u003eEvaluation of Residue Contact Schemes (hydrophilic-hydrophilic vs hydrophilic-all) in Isothermal MD Simulations (400 K, 100 ns)\u003c/b\u003e\u003c/p\u003e\u003cp\u003eBekker et al. classified amino acids into three physicochemical groups - hydrophobic, hydrophilic, and small - and computed Q values across all ten possible residue‑contact pair types (hydrophilic\u0026ndash;hydrophilic, hydrophilic-hydrophobic, hydrophilic-small, hydrophilic-all, hydrophobic-hydrophobic, hydrophobic-small, hydrophobic-all, small-small, small-all, and all-all). They showed that the hydrophilic-all Q value exhibited the highest correlation with experimental Tm, followed by all-all contacts.\u003c/p\u003e\u003cp\u003eIn this study, we first examined which of the two top contact schemes \u0026ndash; hydrophilic-all versus all-all provides better predictive power. To this end, we employed a small set of representative mutants: A79F, the most stabilizing mutation in the first round; G97F, the most destabilizing; and the wild type. For each mutant, we performed triplicate MD simulations (n\u0026thinsp;=\u0026thinsp;3) at 400 K for 100 ns, and computed Q values over the last 30 ns of the trajectories using both hydrophilic\u0026ndash;all and all\u0026ndash;all contact definitions. Contrary to the trend reported by Bekker et al., our results indicated that the all-all Q values showed a stronger correlation with the experimentally determined Tm in our dataset (data not shown). These findings suggest that the optimal contact scheme for Q value calculation may vary depending on the structural context and mutation set.\u003c/p\u003e\u003cp\u003eWe then extended this comparison to the second-round mutants. Specifically, we selected five VHH mutants: A79I and A79W, which exhibited the largest increases in Tm, A79R and A79Q, which showed the most pronounced decreases; and the wild type. For each mutant, we performed MD simulations under the same conditions (400 K, 100 ns, n\u0026thinsp;=\u0026thinsp;3) and computed Q values over the final 30 ns using both hydrophilic-all and all-all contact definitions. Consistent with the first-round results, the all-all Q values again showed a stronger correlation with experimental Tm than the hydrophilic-all Q values in this second-round mutant set (data not shown). This reinforced the applicability of the all-all contact definition for analyzing our VHH system and provided the basis for applying the same contact definition in subsequent temperature-increase MD simulations.\u003c/p\u003e\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e\u003ch2\u003eIsothermal MD Simulations for first-round Mutants\u003c/h2\u003e\u003cp\u003eBased on our preliminary evaluation of residue contact schemes, we adopted the all-all contact definition for subsequent Q value calculations. We then performed a more comprehensive analysis of the full first-round mutant set, which comprised four single-point mutants (A24L, G26I, A79F, and G97F) along with the wild type VHH. For each of these five proteins, we conducted molecular dynamics simulations at 400 K for 100 ns in ten independent replicates (n\u0026thinsp;=\u0026thinsp;10). Q values were calculated using the all-all contact definition over the final 30 ns of each trajectory. Figures\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003ea-e show the Q value trajectories over time, with Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003ea showing the wild type alone, and Figs.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eb-e overlaying each mutant with the wild type. All proteins exhibited a gradual decrease in Q values, indicative of progressive unfolding. Between 70 and 100 ns, Q values generally stabilized within the range of 0.7 to 0.9, suggesting partially unfolded conformations.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eFigure \u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ea presents the distributions of Q values as box-and-whisker plots, which reveal substantial variation among the different mutants. Notably, mutants such as A79F and G97F showed broader distributions than the others. These observations suggest that mutations involving larger amino acid residues may introduce greater variability in the unfolding behavior during partial denaturation. One-way ANOVA with Tukey\u0026rsquo;s multiple comparisons indicated that the mean Q value of G97F was significantly lower than those of the other mutants. This aligns with experimental findings showing that G97F was the only first-round mutant to reduce Tm - an outcome that was not predicted by dStability but was clearly captured by Q value analysis.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eFigure \u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eb summarizes the relationship between the average Q value and the experimentally measured Tm, yielding a Pearson correlation coefficient of \u003cem\u003er\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.59. In contrast, the correlation between dStability scores and Tm was weaker and in the opposite direction (\u003cem\u003er\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.42). Given that dStability is expected to exhibit a negative correlation with Tm, the observed positive correlation is problematic in itself. While both correlations are based on a limited sample size (n\u0026thinsp;=\u0026thinsp;5) and should therefore be interpreted with caution, the results suggest that Q values derived from isothermal MD simulations may provide a more reliable indicator of thermal stability than static energy-based metrics for mutation site selection.\u003c/p\u003e\u003c/div\u003e\n\u003ch3\u003eIsothermal MD Simulations for second-round Mutants\u003c/h3\u003e\n\u003cp\u003eWe next applied the same protocol to the second-round mutant set, consisting of seventeen A79X mutants along with the wild type, all evaluated under identical simulation conditions (400 K, 100 ns, n\u0026thinsp;=\u0026thinsp;10). Q values were also computed over the final 30 ns using the all-all contact definition. The distribution of Q values is shown as box-and-whisker plots in Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ec, illustrating that the degree of variation again differed across proteins. Five mutants - A79F, A79Y, A79V, A79E, and A79R - exhibited particularly large variability. Among these, three mutations (F, Y, and R) involve amino acids with bulky side chains, which may echo the trend observed in first round, where mutations to larger residues were associated with greater variability in partial unfolding. However, this relationship does not hold consistently: for example, tryptophan (W), the largest amino acid, showed only modest variation, while medium-sized residues like valine (V) and glutamate (E) exhibited relatively high variability. Thus, no clear correlation was established between side-chain size and unfolding variability.\u003c/p\u003e\u003cp\u003eOne-way ANOVA with Tukey\u0026rsquo;s multiple comparisons revealed that the Q value of the A79E mutant was significantly lower than those of nine other mutants (A79C, A79I, A79W, A79T, A79L, A79S, A79N, A79V, A79H) in the second-round mutants. Similarly, the A79R mutant showed statistically significant differences in Q value compared to six other mutants (79C, A79I, A79W, A79L, A79N, A79H) in the second-round mutants. Both A79E and A79R resulted in decreased experimental Tm values, with A79R being the most destabilizing among all second-round mutants. These findings indicate that Q value analysis can effectively identify markedly destabilizing mutations even in the context of specific site residue scanning. We then assessed the correlation between the average Q value and the experimentally measured Tm across all 17 mutants and wild type. As shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ed, the Pearson correlation coefficient was \u003cem\u003er\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.44, which is notably lower than that observed for the dStability score (\u003cem\u003er\u003c/em\u003e = \u0026minus;\u0026thinsp;0.64, inverse correlation).\u003c/p\u003e\u003cp\u003eTaken together, these results suggest a complementary relationship between the two metrics: Q value-based MD analysis outperformed dStability in first round (mutation site selection), while dStability showed better predictive power in second round (amino acid selection at a specific site). This indicates that combining both approaches may provide a more balanced and effective strategy for rational antibody stabilization.\u003c/p\u003e\n\u003ch3\u003eTemperature-Increase MD Simulations for First-Round Mutants\u003c/h3\u003e\n\u003cp\u003eWhile experimental approaches for detecting Tm, such as DSC, completely disrupt protein structure, isothermal simulations above 400 K correspond only to partially unfolded states. Moreover, their reliance on long isothermal trajectories imposes a significant computational cost. To better capture structural behavior during complete unfolding, we aimed to establish a protocol that both more closely mimics experimental conditions and is more time-efficient while maintaining predictive performance. To this end, we devised a short (4 ns) temperature-ramp MD protocol, in which the simulation temperature is linearly increased from 300 K to 1000 K.\u003c/p\u003e\u003cp\u003eThis design allows monitoring of native-to-unfolded transitions within a compact simulation window. Q values were calculated across the entire trajectory to assess structural retention throughout the unfolding process.\u003c/p\u003e\u003cp\u003eWe applied this protocol to the first-round mutant set, comprising four single-point mutants (A24L, G26I, A79F, and G97F) and the wild type. Each mutant was simulated in ten independent runs (n\u0026thinsp;=\u0026thinsp;10), and Q values were computed using the all-all contact definition. Figures\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003ef-j display the Q value trajectories over time, with Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003ef showing the wild type alone, and Figs.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eg-j overlaying each mutant with the wild type. All proteins underwent complete structural degradation during the simulation as expected.\u003c/p\u003e\u003cp\u003eBox-and-whisker plots of the resulting Q values are shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ea. As in the isothermal case, considerable variability was observed across mutants, particularly for A24L and G97F. However, unlike the isothermal results, no clear relationship was found between side-chain size and Q value variability. This contrast with the isothermal protocol likely reflects the different unfolding regimes captured: under isothermal conditions, partial unfolding may be more sensitive to steric hindrance introduced by bulky side chains, while the temperature-increase protocol induces full denaturation in all proteins, thereby minimizing size-dependent effects. Under these simulation conditions, the G97F mutant still exhibited the lowest average Q value; however, unlike in the isothermal protocol, G97F did not show statistically significant differences from the wild type or other mutants.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eFigure \u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eb presents the correlation between average Q values and experimental Tm, yielding a strong Pearson correlation coefficient of \u003cem\u003er\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.85. This result exceeds the correlation observed with isothermal simulations at 400 K (\u003cem\u003er\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.59), despite requiring only 1/25th of the computational time. These findings demonstrate that the temperature-increase protocol offers a highly efficient and accurate approach for mutation site selection.\u003c/p\u003e\n\u003ch3\u003eTemperature-Increase MD Simulations for Second-Round Mutants\u003c/h3\u003e\n\u003cp\u003eWe next applied the same temperature-increase protocol to the second-round mutant set, which included seventeen A79X mutants and the wild type. Each mutant was simulated in ten independent runs (n\u0026thinsp;=\u0026thinsp;10), and Q values were computed over the full simulation span using the all-all contact definition.\u003c/p\u003e\u003cp\u003eFigure \u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ec shows the resulting Q value distributions. As observed previously, variability differed among mutants. The most variable mutants were A79N, A79M, and A79Q. However, consistent with first-round results, no clear relationship was observed between side-chain size and variability. Furthermore, unlike the isothermal simulations, no statistically significant differences were observed among the mutants, including the destabilizing A79R mutant.\u003c/p\u003e\u003cp\u003eAs in first round, Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ed shows the correlation between average Q values and experimental Tm, yielding a Pearson correlation coefficient of \u003cem\u003er\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.40. This is nearly identical to that observed with isothermal MD at 400 K (\u003cem\u003er\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.44), indicating that both protocols perform comparably in the context of residue substitution ranking at a specific site. Notably, the A79R mutant, which exhibited the greatest Tm decrease, also showed the second-lowest Q value in this analysis. Taken together, these findings demonstrate that the temperature-increase MD protocol achieves predictive performance comparable to that of conventional isothermal simulations with significantly reduced computational cost.\u003c/p\u003e\u003cp\u003eImportantly, the specific protocol used in this study - linearly increasing the temperature from 300 K to 1000 K over 4 ns and averaging Q values across the full simulation window - represents just one implementation designed to capture the full transition from the native state to complete denaturation. It is likely that further optimization of simulation parameters, including total MD duration, temperature range, and the time window used for Q value calculation, could lead to even more accurate and efficient predictions in future applications.\u003c/p\u003e\u003cp\u003eRecently, a curated database of VHH sequences and their Tm values, NbThermo [\u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e], has been developed. Furthermore, machine learning methods such as NanoMelt [\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e] and TEMPRO [\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e], which partially utilize NbThermo dataset, have also been introduced. In this study, we demonstrated that the two physics-based approaches, Q-value analysis and dStability, are complementary. Moving forward, we plan to investigate the potential of combining these physics-based methods with AI-based approaches to further enhance predictive performance.\u003c/p\u003e"},{"header":"Conclusion","content":"\u003cp\u003eOur study demonstrates that the Q value-based MD approach is more effective than dStability for mutation site selection (first round), whereas dStability remains superior for evaluating residue substitutions at a specific site (second round) in our dataset. Furthermore, the temperature-increase MD protocol achieved comparable predictive performance to the longer isothermal simulations, despite requiring significantly less computational time. These findings suggest that a combined application of Q value-based MD analysis and dStability scoring - along with carefully optimized MD protocols - could offer a powerful and efficient framework for improving the thermal stability of VHHs.\u003c/p\u003e"},{"header":"Materials and Methods","content":"\u003cp\u003eIn the preparations and MD simulations, GROMACS 2018 [\u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e26\u003c/span\u003e] was used. VHH mutants previously reported to yield successful Tm measurements were selected [\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e]. All VHH mutants were modeled using the crystal structure 6JB9 as a template. The PDB structure of 6JB9 contained no missing residues; however, an additional N-terminal serine (Ser0), resulting from the difference of expression host, was present. This residue was removed, and the N-terminus was capped with an N-methyl (NME) group. The modeling was conducted using the MOE software (Chemical Computing Group) [\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e], and the mutations were introduced using the residue mutation function. Side chains at the mutation sites (e.g., A24L, G26I, A79F, G97F in first round; and A79Y, A79W, A79I, A79L, A79V, A79C, A79N, A79Q, A79M, A79H, A79R, A79E, A79K, A79P in second round) were automatically rebuilt by MOE based on its internal rotamer library, and no further energy minimization was applied. For buried mutation sites such as A79, side chain rotamers were manually selected based on the steric compatibility and hydrogen bonding environment when necessary. To prepare the system, hydrogens were attached, and a cubic box was constructed surrounding each protein such that the box was at least 1.0 nm away from the protein and filled with water molecules and 0.153 M of NaCl. Protonation states were assigned assuming a neutral pH (~\u0026thinsp;7.0) using the MOE Protonate 3D module. Histidine residues were manually reviewed and assigned as HID (proton on δ-nitrogen), HIE (proton on ε-nitrogen), or HIP (doubly protonated) based on local hydrogen bonding environments. The Amber99-sb-ildn force field [\u003cspan additionalcitationids=\"CR28\" citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e29\u003c/span\u003e] was used to parameterize the protein, TIP3P [\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e] for waters and the monovalent ion parameters [\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e] were used for the ions. Our simulation protocol commenced with energy minimization, followed by position restraints on the heavy solute atoms during NVT for the first step and NPT for the second step. 3 or 10 parallel simulations were conducted, each initialized with distinct random seed velocities. To equilibrate velocities, a 100 ps NVT simulation at 298K was performed, which allowed a subsequent 100 ps NPT simulation at the same temperature. For these equilibration simulations, position restraints were used on the heavy solute atoms. For the \u0026ldquo;isothermal\u0026rdquo; production run, an unrestrained NVT simulation was performed at 400K, spanning 100 ns. For the \u0026ldquo;temperature-increase\u0026rdquo; production run, an unrestrained NVT simulation was performed with temperature increasing from 298K to 1000K, spanning 4 ns. A 1.0-nm cutoff was used for the short-range interactions in combination with the Particle Mesh Ewald method for the long-range electrostatic interactions. A consistent 2 fs timestep was maintained, with constraints managed using LINCS [\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e] for the protein and water molecules. Temperature was controlled using the V-rescale thermostat [\u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e33\u003c/span\u003e] and pressure was maintained using Parrinello-Rahman barostat [\u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e34\u003c/span\u003e] at 1 bar during the NPT simulation. During the production run, structures were saved at 10-ps intervals.\u003c/p\u003e\u003cp\u003eQ values were calculated following the definition of Best et al. [\u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e35\u003c/span\u003e], using a contact threshold of 0.45 nm between non-hydrogen atoms. Residues within the N-terminal cap (NME) and residues at the flexible C-terminal region (last 3 residues) were excluded from the calculation. Based on the previous study [\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e], two types of Q-value calculation methods that showed the highest correlation with Tm were defined: all-all (all residue pairs) and HA (hydrophilic-all), with hydrophilic residues defined as D, E, N, Q, R, K, and H. The average Q-value under the isothermal condition was calculated from the final 30 ns of the entire trajectory, as in the previous study [\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e], whereas for the temperature-increase condition, the average Q-value was calculated over the entire 4 ns trajectory. Statistical analysis were performed using Microsoft Excel with Real Statistics Add-In. For Q value comparisons, one-way ANOVA with Tukey\u0026rsquo;s multiple comparisons was used.\u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003ch2\u003eCompeting interests statement\u003c/h2\u003e\u003cp\u003eThe authors declare no competing interests.\u003c/p\u003e\u003c/p\u003e\u003ch2\u003eFunding statement\u003c/h2\u003e\u003cp\u003eThis study was carried out with the support of Astellas Pharma Inc.\u003c/p\u003e\u003ch2\u003eAuthor Contribution\u003c/h2\u003e\u003cp\u003eY.T. and H.S. designed the study. R.Y. conducted the experiments, Y.T., R.Y. and H.S. analyzed the results. Y.T., R.Y. and H.S. wrote and approved the manuscript.\u003c/p\u003e\u003ch2\u003eAcknowledgement\u003c/h2\u003e\u003cp\u003eWe extend our deepest gratitude to Ms. Nami Tanigawa (WORLD INTEC), whose expertise and dedication in proposing conditions, conducting molecular dynamics simulations, and providing thoughtful data analyses were indispensable to this work.\u003c/p\u003e\u003ch2\u003eData Availability\u003c/h2\u003e\u003cp\u003eThe datasets generated and/or analyzed during the current study are available from the corresponding author on reasonable request.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\u003cli\u003e\u003cspan\u003eHamers-Casterman, C. et al. Naturally occurring antibodies devoid of light chains. \u003cem\u003eNature\u003c/em\u003e \u003cb\u003e363\u003c/b\u003e, 446\u0026ndash;448. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1038/363446a0\u003c/span\u003e\u003cspan address=\"10.1038/363446a0\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (1993).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eK\u0026ouml;nning, D. et al. Camelid and shark single domain antibodies: structural features and therapeutic potential. \u003cem\u003eCurr. Opin. Struct. Biol.\u003c/em\u003e \u003cb\u003e45\u003c/b\u003e, 10\u0026ndash;16. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.1016/j.sbi.2016.10.019\u003c/span\u003e\u003cspan address=\"10.1016/j.sbi.2016.10.019\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2017).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eGoldman, E. R., Liu, J. L., Zabetakis, D. \u0026amp; Anderson, G. P. Enhancing Stability of Camelid and Shark Single Domain Antibodies: An Overview. \u003cem\u003eFront. Immunol.\u003c/em\u003e \u003cb\u003e8\u003c/b\u003e, 865. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.3389/fimmu.2017.00865\u003c/span\u003e\u003cspan address=\"10.3389/fimmu.2017.00865\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2017).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eHagihara, Y., Mine, S. \u0026amp; Uegaki, K. Stabilization of an immunoglobulin fold domain by an engineered disulfide bond at the buried hydrophobic region. \u003cem\u003eJ. Biol. Chem.\u003c/em\u003e \u003cb\u003e282\u003c/b\u003e, 36489\u0026ndash;36495. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1074/jbc.M707078200\u003c/span\u003e\u003cspan address=\"10.1074/jbc.M707078200\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2007).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eSaerens, D., Conrath, K., Govaert, J. \u0026amp; Muyldermans, S. Disulfide bond introduction for general stabilization of immunoglobulin heavy-chain variable domains. \u003cem\u003eJ. Mol. Biol.\u003c/em\u003e \u003cb\u003e377\u003c/b\u003e, 478\u0026ndash;488. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1016/j.jmb.2008.01.022\u003c/span\u003e\u003cspan address=\"10.1016/j.jmb.2008.01.022\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2008).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eArbabi-Ghahroudi, M. et al. Aggregation-resistant VHs selected by in vitro evolution tend to have disulfide-bonded loops and acidic isoelectric points. \u003cem\u003eProtein Eng. Des. Sel.\u003c/em\u003e \u003cb\u003e22\u003c/b\u003e, 59\u0026ndash;66. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1093/protein/gzn071\u003c/span\u003e\u003cspan address=\"10.1093/protein/gzn071\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2009).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eZabetakis, D., Olson, M. A., Anderson, G. P., Legler, P. M. \u0026amp; Goldman, E. R. Evaluation of disulfide bond position to enhance the thermal stability of a highly stable single domain antibody. \u003cem\u003ePLoS One\u003c/em\u003e. \u003cb\u003e9\u003c/b\u003e, e115405. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1371/journal.pone.0115405\u003c/span\u003e\u003cspan address=\"10.1371/journal.pone.0115405\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2014).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eAkazawa-Ogawa, Y., Uegaki, K. \u0026amp; Hagihara, Y. The role of intra-domain disulfide bonds in heat-induced irreversible denaturation of camelid single domain VHH antibodies. \u003cem\u003eJ. Biochem.\u003c/em\u003e \u003cb\u003e159\u003c/b\u003e, 111\u0026ndash;121. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1093/jb/mvv082\u003c/span\u003e\u003cspan address=\"10.1093/jb/mvv082\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2016).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eShinozaki, N., Hashimoto, R., Fukui, K. \u0026amp; Uchiyama, S. Efficient generation of single domain antibodies with high affinities and enhanced thermal stabilities. \u003cem\u003eSci. Rep.\u003c/em\u003e \u003cb\u003e7\u003c/b\u003e, 5794. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1038/s41598-017-06277-x\u003c/span\u003e\u003cspan address=\"10.1038/s41598-017-06277-x\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2017).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eSaerens, D. et al. Identification of a universal VHH framework to graft non-canonical antigen-binding loops of camel single-domain antibodies. \u003cem\u003eJ. Mol. Biol.\u003c/em\u003e \u003cb\u003e352\u003c/b\u003e, 597\u0026ndash;607. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1016/j.jmb.2005.07.038\u003c/span\u003e\u003cspan address=\"10.1016/j.jmb.2005.07.038\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2005).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eVincke, C. et al. General strategy to humanize a camelid single-domain antibody and identification of a universal humanized nanobody scaffold. \u003cem\u003eJ. Biol. Chem.\u003c/em\u003e \u003cb\u003e284\u003c/b\u003e, 3273\u0026ndash;3284. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1074/jbc.M806889200\u003c/span\u003e\u003cspan address=\"10.1074/jbc.M806889200\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2009).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eZabetakis, D., Anderson, G. P., Bayya, N. \u0026amp; Goldman, E. R. Contributions of the complementarity determining regions to the thermal stability of a single-domain antibody. \u003cem\u003ePLoS One\u003c/em\u003e. \u003cb\u003e8\u003c/b\u003e, e77678. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1371/journal.pone.0077678\u003c/span\u003e\u003cspan address=\"10.1371/journal.pone.0077678\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2013).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eAnderson, G. P., Liu, J. H., Zabetakis, D., Liu, J. L. \u0026amp; Goldman, E. R. Thermal stabilization of anti-alpha-cobratoxin single domain antibodies. \u003cem\u003eToxicon\u003c/em\u003e \u003cb\u003e129\u003c/b\u003e, 68\u0026ndash;73. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1016/j.toxicon.2017.02.008\u003c/span\u003e\u003cspan address=\"10.1016/j.toxicon.2017.02.008\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2017).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eShriver-Lake, L. C., Zabetakis, D., Goldman, E. R. \u0026amp; Anderson, G. P. Evaluation of anti-botulinum neurotoxin single domain antibodies with additional optimization for improved production and stability. \u003cem\u003eToxicon\u003c/em\u003e \u003cb\u003e135\u003c/b\u003e, 51\u0026ndash;58. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1016/j.toxicon.2017.06.002\u003c/span\u003e\u003cspan address=\"10.1016/j.toxicon.2017.06.002\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2017).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eKunz, P. et al. Exploiting sequence and stability information for directing nanobody stability engineering. \u003cem\u003eBiochim. Biophys. Acta Gen. Subj.\u003c/em\u003e \u003cb\u003e1861\u003c/b\u003e, 2196\u0026ndash;2205. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1016/j.bbagen.2017.06.014\u003c/span\u003e\u003cspan address=\"10.1016/j.bbagen.2017.06.014\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2017).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eHarmsen, M. M., van Solt, C. B., van Zijderveld-van Bemmel, A. M., Niewold, T. A. \u0026amp; van Zijderveld, F. G. Selection and optimization of proteolytically stable llama single-domain antibody fragments for oral immunotherapy. \u003cem\u003eAppl. Microbiol. Biotechnol.\u003c/em\u003e \u003cb\u003e72\u003c/b\u003e, 544\u0026ndash;551. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1007/s00253-005-0300-7\u003c/span\u003e\u003cspan address=\"10.1007/s00253-005-0300-7\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2006).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eTurner, K. B., Zabetakis, D., Goldman, E. R. \u0026amp; Anderson, G. P. Enhanced stabilization of a stable single domain antibody for SEB toxin by random mutagenesis and stringent selection. \u003cem\u003eProtein Eng. Des. Sel.\u003c/em\u003e \u003cb\u003e27\u003c/b\u003e, 89\u0026ndash;95. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1093/protein/gzu001\u003c/span\u003e\u003cspan address=\"10.1093/protein/gzu001\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2014).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eHussack, G. et al. Protease-resistant single-domain antibodies inhibit Campylobacter jejuni motility. \u003cem\u003eProtein Eng. Des. Sel.\u003c/em\u003e \u003cb\u003e27\u003c/b\u003e, 191\u0026ndash;198. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1093/protein/gzu011\u003c/span\u003e\u003cspan address=\"10.1093/protein/gzu011\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2014).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eTomimoto, Y., Yamazaki, R. \u0026amp; Shirai, H. Increasing the melting temperature of VHH with the in silico free energy score. \u003cem\u003eSci. Rep.\u003c/em\u003e \u003cb\u003e13\u003c/b\u003e, 4922. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1038/s41598-023-32022-8\u003c/span\u003e\u003cspan address=\"10.1038/s41598-023-32022-8\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2023).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eMolecular Operating Environment (MOE), 2019.01; Chemical Computing Group ULC, 1010 Sherbooke St. West, Suite #910, Montreal, QC \u0026amp; Canada H3A 2R7, (2021).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eBekker, G. J., Ma, B. \u0026amp; Kamiya, N. Thermal stability of single-domain antibodies estimated by molecular dynamics simulations. \u003cem\u003eProtein Sci.\u003c/em\u003e \u003cb\u003e28\u003c/b\u003e, 429\u0026ndash;438. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1002/pro.3546\u003c/span\u003e\u003cspan address=\"10.1002/pro.3546\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2019).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eZabetakis, D., Shriver-Lake, L. C., Olson, M. A., Goldman, E. R. \u0026amp; Anderson, G. P. Experimental evaluation of single-domain antibodies predicted by molecular dynamics simulations to have elevated thermal stability. \u003cem\u003eProtein Sci.\u003c/em\u003e \u003cb\u003e28\u003c/b\u003e, 1909\u0026ndash;1912. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1002/pro.3692\u003c/span\u003e\u003cspan address=\"10.1002/pro.3692\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2019).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eValdes-Tresanco, M. S., Valdes-Tresanco, M. E., Molina-Abad, E. \u0026amp; Moreno, E. NbThermo: a new thermostability database for nanobodies. \u003cem\u003eDatabase (Oxford)\u003c/em\u003e (2023). \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1093/database/baad021\u003c/span\u003e\u003cspan address=\"10.1093/database/baad021\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2023).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eRamon, A. et al. Prediction of protein biophysical traits from limited data: a case study on nanobody thermostability through NanoMelt. \u003cem\u003eMAbs\u003c/em\u003e 17, 2442750, (2025). \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1080/19420862.2024.2442750\u003c/span\u003e\u003cspan address=\"10.1080/19420862.2024.2442750\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eAlvarez, J. A. E. \u0026amp; Dean, S. N. TEMPRO: nanobody melting temperature estimation model using protein embeddings. \u003cem\u003eSci. Rep.\u003c/em\u003e \u003cb\u003e14\u003c/b\u003e, 19074. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1038/s41598-024-70101-6\u003c/span\u003e\u003cspan address=\"10.1038/s41598-024-70101-6\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2024).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eAbraham, M. J. et al. High performance molecular simulations through multi-level parallelism from laptops to supercomputers. \u003cem\u003eSoftwareX 1\u0026ndash;2\u003c/em\u003e. \u003cb\u003eGROMACS\u003c/b\u003e \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.1016/j.softx.2015.06.001\u003c/span\u003e\u003cspan address=\"10.1016/j.softx.2015.06.001\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2015). 19\u0026ndash;25, doi.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eAbu-Soud, H. M., Clark, A. C., Francisco, W. A., Baldwin, T. O. \u0026amp; Raushel, F. M. Kinetic destabilization of the hydroperoxy flavin intermediate by site-directed modification of the reactive thiol in bacterial luciferase. \u003cem\u003eJ. Biol. Chem.\u003c/em\u003e \u003cb\u003e268\u003c/b\u003e, 7699\u0026ndash;7706. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.1016/S0021-9258(18)53013-1\u003c/span\u003e\u003cspan address=\"10.1016/S0021-9258(18)53013-1\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (1993).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eHornak, V. et al. Comparison of multiple Amber force fields and development of improved protein backbone parameters. \u003cem\u003eProteins Struct. Funct. Bioinform.\u003c/em\u003e \u003cb\u003e65\u003c/b\u003e, 712\u0026ndash;725. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.1002/prot.21123\u003c/span\u003e\u003cspan address=\"10.1002/prot.21123\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2006).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eLindorff-Larsen, K. et al. Improved side-chain torsion potentials for the Amber ff99SB protein force field. \u003cem\u003eProteins Struct. Funct. Bioinform.\u003c/em\u003e \u003cb\u003e78\u003c/b\u003e, 1950\u0026ndash;1958. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.1002/prot.22711\u003c/span\u003e\u003cspan address=\"10.1002/prot.22711\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2010).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eJorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W. \u0026amp; Klein, M. L. Comparison of simple potential functions for simulating liquid water. \u003cem\u003eJ. Chem. Phys.\u003c/em\u003e \u003cb\u003e79\u003c/b\u003e, 926\u0026ndash;935 (1983).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eJoung, I. S. \u0026amp; Cheatham, T. E. III Determination of alkali and halide monovalent ion parameters for use in explicitly solvated biomolecular simulations. \u003cem\u003eJ. Phys. Chem. B\u003c/em\u003e. \u003cb\u003e112\u003c/b\u003e, 9020\u0026ndash;9041 (2008).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eHess, B. \u0026amp; P-LINCS: A Parallel Linear Constraint Solver for Molecular Simulation. \u003cem\u003eJ. Chem. Theory Comput.\u003c/em\u003e \u003cb\u003e4\u003c/b\u003e, 116\u0026ndash;122. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1021/ct700200b\u003c/span\u003e\u003cspan address=\"10.1021/ct700200b\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2008).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eBussi, G., Donadio, D. \u0026amp; Parrinello, M. Canonical sampling through velocity rescaling. \u003cem\u003eThe J. Chem. physics\u003c/em\u003e 126 (2007).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eParrinello, M. \u0026amp; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. \u003cem\u003eJ. Appl. Phys.\u003c/em\u003e \u003cb\u003e52\u003c/b\u003e, 7182\u0026ndash;7190 (1981).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eBest, R. B., Hummer, G. \u0026amp; Eaton, W. A. Native contacts determine protein folding mechanisms in atomistic simulations. \u003cem\u003eProceedings of the National Academy of Sciences\u003c/em\u003e 110, 17874\u0026ndash;17879 (2013).\u003c/span\u003e\u003c/li\u003e\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":true,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"","lastPublishedDoi":"10.21203/rs.3.rs-7555873/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-7555873/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eVHHs are frequently employed in protein therapeutics; however, enhancing their thermal stability (Tm) remains a significant challenge. We previously developed a two‑step in silico strategy to enhance Tm by identifying mutation sites (first round) and selecting favorable substitutions (second round) using dStability, a Gibbs free energy-based score. While dStability performed reasonably in the second round (\u003cem\u003er\u003c/em\u003e = -0.64), it failed to predict the effect of G97F mutation, ranking it as most stabilizing despite its strong destabilizing impact on experimental Tm. In this study, we re-evaluated this dataset using high-temperature molecular dynamics (MD) simulations, employing Q-values, originally proposed by Bekker et al., as a quantitative metric to assess the extent of structural degradation. This method successfully identified G97F as least stable, demonstrating utility for detecting destabilizing mutations for the first round, though its second‑round performance was weaker (\u003cem\u003er\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.44). We also developed a temperature‑increase MD protocol, ramping simulation temperature from 300 K to 1000 K over 4 ns to completely degrade protein. Despite being 25-fold faster than fixed-temperature simulations, this approach retained comparable predictive performance. Overall, combining these two physics-based approaches, dStability and Q‑value analysis with optimized MD protocols, enables efficient identification of stabilizing mutations.\u003c/p\u003e","manuscriptTitle":"Complementary Roles of Physics-Based Approaches in Predicting VHH Thermal Stability","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2025-10-09 13:21:07","doi":"10.21203/rs.3.rs-7555873/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"eb2db859-b7a4-41db-a7b4-34876f4114d7","owner":[],"postedDate":"October 9th, 2025","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[{"id":55925903,"name":"Biological sciences/Biophysics"},{"id":55925904,"name":"Biological sciences/Computational biology and bioinformatics"},{"id":55925905,"name":"Biological sciences/Drug discovery"},{"id":55925906,"name":"Physical sciences/Mathematics and computing"}],"tags":[],"updatedAt":"2025-11-17T10:54:09+00:00","versionOfRecord":[],"versionCreatedAt":"2025-10-09 13:21:07","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-7555873","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-7555873","identity":"rs-7555873","version":["v1"]},"buildId":"8U1c8b4HqxoKbykW_rLl7","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}