Introduction

This study establishes an integrated framework for the rational design and safety optimization of camelid VHH antibodies targeting VEGF165. The framework couples machine learning–driven prediction of binding free energy changes (ΔΔG) using a low-rank bilinear model trained on the AbBind dataset with structural quality evaluation (Clashscore) via a random forest model. By integrating these parameters into a weighted, rank-based selection module, we systematically identify VHH variants exhibiting improved affinity and conformational stability.

To further ensure therapeutic precision and minimize systemic toxicity, we incorporated a Probody-based conditional activation strategy. A cyclic masking peptide was de novo designed using deep generative models (RFdiffusion and ProteinMPNN) to sterically block the antibody’s active site. This peptide is tethered through a rigid linker containing an MMP3-cleavable sequence, enabling selective activation in pathological microenvironments characterized by elevated MMP3 expression.

Comprehensive molecular simulations-including Steered Molecular Dynamics (SMD) and Umbrella Sampling (US)-were conducted to delineate the dissociation pathway and quantify the peptide’s binding free energy. The results demonstrated stable masking during systemic circulation and efficient site-specific activation upon cleavage, confirming the designed system’s structural integrity and functional controllability.

ΔΔG Prediction: Assessing Affinity Change

Binding free-energy change (ΔΔG) is a central indicator of how mutations affect antibody-antigen affinity. Negative ΔΔG suggests affinity enhancement, while positive ΔΔG indicates reduced binding. Thus, ΔΔG is a critical energetic reference in antibody engineering.

However, datasets for VHH or tandem VH antibodies are limited, making direct training infeasible. To address this limitation, we considered the evolutionary and structural homology between camelid VHH frameworks and conventional VH domains. Prior comparative analyses have shown that the framework regions (FRs) of VHHs are highly similar to those of human and murine VH domains, with only a few compensatory mutations that stabilize VHHs in the absence of a light chain [4-6]. This homology provides a strong rationale for knowledge transfer: VH-dominant energetic contributions are expected to be conserved across VH-VL and VH-only architectures.

Based on this, we proposed the heavy–light separability assumption: the main contribution of VH to binding free energy can be modeled independently of VL, and the absence of VL does not invalidate the predictive framework. In practice, the light-chain embedding was standardized and replaced with a zero vector during inference, eliminating light-chain terms and reducing the bilinear energy function to a heavy-chain–only predictor. This treatment not only reflects the biological reality of VHH and tandem VH antibodies but also ensures consistency between training and inference, enabling us to extend models trained on AbBind (VH–VL) to VH-only systems.

We trained our predictive model on the AbBind dataset [31], which provides antibody–antigen complexes with single-point mutations and experimentally measured changes in binding free energy (ΔΔG).

Feature construction

Each antibody sequence was divided into VH and VL domains, and sequence embeddings were extracted separately using ESM-1b. For each mutation , both wild-type and mutant embeddings were obtained. We then constructed difference vectors to explicitly encode mutational effects:

$$\Delta H = H_{\text{mut}} - H_{\text{wt}}, \quad \Delta L = L_{\text{mut}} - L_{\text{wt}}$$

where ΔH and ΔL denote the differences in heavy-chain and light-chain embeddings, respectively.

Low-rank bilinear energy function

We modeled mutational effects with a low-rank bilinear energy function:

$$E(H,L)= w_H^\top H + w_L^\top L + \sum_{r=1}^{R} \langle U_r, H \rangle \langle V_r, L \rangle$$

Here, w_H and w_L are linear weights capturing chain-specific main effects, while U_r and V_r are rank-R projection vectors modeling VH-VL interactions. This low-rank approximation reduces parameter complexity while preserving nonlinear expressiveness.

Prediction mapping

The scalar energy score was mapped to the predicted ΔΔG value through a regression layer:

$$\widehat{\Delta\Delta G} = f(E(\Delta H,\Delta L))$$

Heavy–light separability assumption

During inference on VHH or tandem VH antibodies, the light-chain embedding was standardized and replaced with a zero vector, effectively removing light-chain terms:

$$E(H,0)= w_H^\top H$$

This reduces the predictor to a heavy-chain–only model, consistent with the biology of VHH antibodies.

Calibration

Predicted values were further calibrated using isotonic regression:

$$\widehat{\Delta\Delta G}_{\text{calib}} = g\!\left(\widehat{\Delta\Delta G}\right)$$

where g(⋅) is a monotonic transformation that preserves rank ordering while aligning predictions to the experimental ΔΔG scale.

Conversion to fold-change affinity

To aid interpretation, we converted calibrated ΔΔG predictions into fold-change affinity improvements, using the thermodynamic relation:

$$\text{Fold Improvement} = \exp\!\left(-\frac{\widehat{\Delta\Delta G}_{\text{calib}}}{RT}\right)$$

where R is the gas constant and T is the temperature (set to 298 K).

We evaluated the performance of our ΔΔG predictor on the AbBind dataset. The rank-96 bilinear interaction energy model achieved an out-of-fold (OOF) Spearman correlation of approximately 0.5, with a root mean squared error (RMSE) of about 2.9 and a absolute error (MAE) of about 1.7. While this level of accuracy does not reach experimental precision, it is sufficient for high-throughput mutation screening and prioritization.

To further evaluate performance, we analyzed OOF results across mutation groups. Figure 1 shows the distribution of per-group Spearman correlations, indicating that most groups fall into the moderate-to-good range with positive median values. This suggests that the model can stably capture ranking signals within groups, though certain groups are limited by sample size or higher measurement noise.

Figure 1. OOF performance across groups (Spearman correlation).

We also visualized predicted vs experimental ΔΔG after filtering out extreme experimental conditions (ΔΔG outside [-6, 6]) for clarity. As shown in Figure 2, predictions cluster around the diagonal, reflecting that the model successfully captures primary energetic trends while a few outliers remain. This filtered view highlights the model’s capability to provide robust prioritization despite limited VHH-specific training data.

Figure 2. Predicted vs Experimental ΔΔG, with extreme values filtered (ΔΔG ∈ [-6, 6]).

Clashscore: Random-forest structural quality prediction

In protein design, optimizing only binding free energy (ΔΔG) is insufficient, as some mutations may increase affinity but introduce severe steric clashes, reducing structural plausibility and expressibility. ClashScore is a widely used structural validation metric that quantifies atomic clashes in molecular models. Lower ClashScore values indicate higher structural quality and stability.

Therefore, in our pipeline, joint evaluation of ΔΔG and ClashScore is essential: ΔΔG reflects binding affinity changes, while ClashScore constrains structural plausibility. Together, they allow us to prioritize mutations that balance both affinity and structural quality.

We trained a Random Forest regressor to predict the MolProbity ClashScore, which measures the degree of steric clashes in protein structures. The feature set was derived from ESM embeddings extracted from antibody–antigen complexes, representing high-dimensional sequence–structure information.

To improve robustness and reduce redundancy, PCA (100 dimensions) was applied before fitting the Random Forest model. The regressor was trained using 5-fold cross-validation, with hyperparameters set to 700 trees and unrestricted maximum depth to balance bias and variance. Duplicate entries in the dataset (same PDB, mutation, chain) were aggregated by averaging labels, while keeping the first embedding record.

Performance metrics including R², MAE, RMSE, Pearson and Spearman correlations were calculated on the held-out test set. The Random Forest model achieved an R² of 0.50 with a Pearson correlation of 0.73 and a Spearman correlation of 0.63, indicating moderate-to-strong predictive ability in estimating ClashScore (MAE = 3.99; RMSE = 7.19). Feature importance analysis further revealed key contributors to prediction, and multiple diagnostic plots were generated.

Figure 3 (True vs Predicted) shows that most points fall near the diagonal, indicating alignment between predicted and experimental ClashScores, though variance increases at higher values.

Figure 3. True Clashscore vs Predicted Clashscore.

Figure 4 (Residual Distribution) reveals that most errors cluster around zero, suggesting low systematic bias, with a few outliers indicating over- or under-prediction in specific cases.

Figure 4. Residual Distribution.

Figure 5 (Bland–Altman Plot) confirms that the prediction errors are symmetrically distributed around the mean difference (≈0.9), with 95% of points falling within ±14 units, highlighting consistency while also capturing uncertainty.

Figure 5. Bland–Altman Plot.

These results confirm that the model achieves screening-level reliability: while not precise enough for fine-grained structural refinement, it robustly distinguishes between mutations that lead to structurally favorable vs. unfavorable conformations. When combined with ΔΔG predictions, ClashScore serves as an orthogonal filter, enabling Pareto optimization of affinity and structural stability in antibody design.

VHH Prediction

For VHH antibodies, both binding affinity and structural stability are essential determinants of functionality. Binding affinity changes can be quantified through ΔΔG, where negative values indicate stronger antigen binding. Structural stability can be assessed using MolProbity Clashscore, which evaluates steric clashes in the protein structure. Optimizing affinity without maintaining stability risks producing unstable or misfolded antibodies. Therefore, a dual-objective evaluation combining ΔΔG and Clashscore is necessary to ensure rational prioritization of candidate mutations.

We refined the scoring system to ensure that higher scores directly reflect better overall performance. For each mutation, ΔΔG and Clashscore were first converted into ranks and then normalized to a 0–1 scale. To emphasize the stronger biological relevance of ΔΔG as an indicator of binding affinity change, we applied a higher weight to ΔΔG during score calculation. The final score was defined as:

$$ \text{Score} = \alpha \cdot (1 - r_{\text{ddG}}) + \beta \cdot (1 - r_{\text{Clash}}) $$

where \(r_{\text{ddG}}\) and \(r_{\text{Clash}}\) are the normalized ranks of ΔΔG and Clashscore, respectively, and \(\alpha > \beta\) assigns more importance to ΔΔG (e.g., \(\alpha = 2, \beta = 1\)). By construction, higher scores correspond to mutations that achieve greater affinity improvements while maintaining structural stability. This rank-based weighted scoring avoids problems of unit inconsistency and provides a more interpretable prioritization of variants.

We evaluated single-point VHH mutations under a dual-objective view of affinity improvement and structural quality. The x-axis represents affinity improvement (WT Kd / Mut Kd), and the y-axis is the MolProbity Clashscore. Panel A in Figure 6 shows the full mutation space: most variants cluster near the wild type (WT), while a subset falls into the “better quadrant”, where affinity improves (right-shift) and Clashscore decreases (down-shift). Panel B zooms into this quadrant and highlights the Top-K candidates.

To produce a single ranking consistent with both objectives, we used a weighted rank-based score in which higher values indicate better overall performance. Specifically, ΔΔG and Clashscore predictions were converted to ranks, normalized to [0,1], and then combined as:

$$ \text{Score} = \alpha \cdot (1 - r_{\Delta\Delta G}) + \beta \cdot (1 - r_{\text{Clash}}), $$

with

$$ \alpha > \beta \, (\alpha = 2, \beta = 1) $$

This design prioritizes affinity (ΔΔG) while still penalizing excessive steric clashes (Clashscore), yielding an interpretable, unit-free ordering. The Top-K table below reports each candidate’s predicted ΔΔG, Clashscore, affinity fold change, their individual ranks, and the final weighted score. Variants at the top of the list therefore reflect strong affinity gains together with acceptable or improved structural quality, making them suitable for downstream validation.

Note: “Weighted Score” is computed from normalized ranks (higher is better). The qualitative labels (high / moderate) are placeholders for the numeric values produced by the scoring script.

Masking Peptide Design

Why do we need to design ?

During the directed evolution of anti-VEGF antibodies, we observed that VEGF receptors are widely distributed throughout the body. Direct use of high-affinity anti-VEGF antibodies could pose significant safety risks. Fortunately, the recent emergence of Probody technology [9] offers a viable solution to this challenge. This approach introduces masking peptide sequences on antibodies that can be cleaved by specific proteases, rendering the antibodies “inactive” in normal tissues. Only in pathological regions, where protease activity increases, are they cleaved and regain their function. Significantly, matrix metalloproteinase 3 (MMP3) exhibits relatively high expression in hemorrhoidal lesions [10], providing the ideal condition for Probody's localized, specific activation. Therefore, we plan to design masking peptides for anti-VEGF antibodies based on the Probody strategy, enabling activation of the anti-VEGF antibody only at the hemorrhoidal lesion site through cleavage. This approach aims to ensure therapeutic efficacy while significantly reducing the risk of systemic side effects.

How will it be implemented ?

We plan to design masking peptides targeting the binding epitopes of anti-VEGF antibodies and VEGF to achieve controllable modulation of antibody activity. Considering their small molecular weight, structural stability, and ease of expression, we selected cyclic peptides as the design target for masking peptides. Since no natural binding proteins meeting the requirements have been identified to date, we decided to adopt a de novo design strategy to develop novel cyclic peptide molecules. Leveraging the recent advances in RFdiffusion [11] for cyclic peptide design and its robust structure generation capabilities, we aim to establish a rational protein design pipeline using this method. This approach will enable efficient cyclic peptide design and optimization targeting VEGF epitopes.

Since the crystal structure of the anti-VEGF antibody-VEGF-A complex has not yet been determined, we employed HADDOCK [16] for blind docking analysis to identify potential binding epitopes for subsequent targeted design of cyclic peptides mimicking epitope structures. Among the docking results, we selected the conformation cluster with the highest number of clusters and the most favorable HADDOCK composite score for focused analysis.

From the energy component perspective, both van der Waals energy (−33.1 ± 7.2) and electrostatic energy (−107.8 ± 8.5) exhibit significant attractive forces, indicating that the interaction between anti-VEGF antibodies and VEGF-A is primarily electrostatic in nature, accompanied by favorable hydrophobic contacts. The solvation energy (−0.4 ± 2.0) is close to zero, indicating good polarity compensation at the interface and minimal influence of solvent effects on binding stability. The buried surface area (BSA = 1251.6 ± 224.5 Ų) is large, consistent with the binding characteristics of a typical antibody–antigen complex.

Based on this, we performed 100 ns molecular dynamics simulations on the selected complexes. The results indicate (Figure 7) that the complex maintained overall structural stability throughout the simulation, with the RMSD consistently remaining around 2.5 Å. Trajectory analysis revealed that the two components remained tightly bound throughout the entire process. Minor fluctuations in the RMSD primarily stemmed from conformational adjustments in the loop region rather than dissociation at the interface. Comprehensive analysis suggests that the epitope obtained from this docking exhibits stable conformation and high binding rationality, serving as a structural foundation for subsequent epitope-targeting peptide design.

Figure 7. RMSD diagram of the VEGF-A and anti-VEGF antibody complex.

After obtaining a reasonable complex structure, we further visualized key interacting residues at the antibody-antigen interface using PyMOL [32] to identify potential hotspot regions. This interface residue information provides crucial reference for subsequent RFdiffusion-based directed design, enabling the setting of hotspot residues in model inputs to guide generated structures toward achieving more precise binding at this epitope region.

Figure 8. VEGF-A with anti-VEGF antibody interaction.

To design cyclic peptides capable of specifically masking epitopes at the binding interface between VEGF and anti-VEGF antibodies, we employed the RFpeptides extension protocol of RFdiffusion [13]. This method generates large cyclic peptide structures bound to target proteins at atomic resolution.

During the design process, we first used the resolved VEGF-antiVEGF complex as the input structure and identified key hotspot residues at its binding interface as design constraints. Subsequently, we invoked RFpeptides' macrocyclic binder design workflow, specifying macrocyclic parameters during execution to ensure the generated peptide backbone forms a stable ring structure. Concurrently, we incorporated key amino acid residues from the VEGF epitope as constraint inputs to guide the generated macropeptide toward interacting with this region.

In this process, we generated a total of 50 candidate cyclic peptide structures and further screened them based on masking efficacy and peptide size, prioritizing those that blocked as many epitopes as possible while maintaining a smaller size. The cyclic peptide structure with the optimal masking effect was ultimately selected, laying the foundation for subsequent sequence design.

Figure 9. Schematic of Epitope Masking by Cyclopeptides.
A: Front view of the cyclic peptide scaffold structure; B: Top view of the cyclic peptide scaffold structure.Gray represents the anti-VEGF antibody, pink represents VEGF-A, and blue represents the cyclic peptide scaffold structure.

After obtaining the backbone structure of the cyclic peptide, we further optimized its amino acid sequence using ProteinMPNN [14]. To achieve our objectives, the cyclic peptide must not bind too tightly to anti-VEGF, ensuring it remains attached even after the linker is cleaved by matrix metalloproteinase-3. Simultaneously, the binding must not be too weak, as this would compromise its masking efficacy. We therefore adjusted several parameters of ProteinMPNN. Specifically, we employed a higher-noise model to generate more dispersed, less constrained sequences, thereby avoiding overly strong binding. Simultaneously, we increased the sampling temperature to enhance sequence diversity and reduce the model's fixation on strong interface residues. We generated 50 sequence variants for the candidate cyclic peptide backbone. Using Modeller [15], we mapped these sequences onto cyclic peptide structures and selected the best-mapped structure for batch docking. We employed HDOCK [23] for batch docking of the cyclic peptides with anti-VEGF, yielding 500 conformations from which we screened for reasonably docked structures. Subsequently, Rosetta InterfaceAnalyzer [17] analyzed these structures to identify sequences exhibiting non-compact binding with feasible separation potential.

Figure 10. Atomic-Filling Interface and Binding Free Energy Distribution Diagram.

The reasonable range for packstact is 0.6–0.8, indicating the degree of binding. Peptides within the green band may potentially detach.

After obtaining the reasonable sequence and structure, we performed a preliminary analysis of the protein-protein complex interface using PRODIGY [18]. The results indicate that the binding free energy of the complex is −9.3 kcal/mol, corresponding to a dissociation constant Kd ≈ 2.6×10⁻⁷ M (approximately 260 nM), suggesting moderate affinity between the two components. Interfacial interactions were primarily dominated by hydrophobic (apolar–apolar, 16%) and polar–apolar (10%) contacts, while charged–charged (5%) and polar–polar (1%) interactions were less prevalent. Furthermore, the nonpolar surface area at the interface (NIS apolar, 41.0%) was significantly higher than the charged surface area (NIS charged, 23.1%), suggesting that this binding relies on a hydrophobic core rather than a strong electrostatic network. Overall, this interaction pattern is consistent with a moderately affinity-bound yet reversible complex, featuring a relatively open and flexible interface that may be susceptible to dissociation under physiological conditions.

Why we need linker ?

Due to the limited binding affinity of the masking cyclic peptide, we aim to design a rigid linker incorporating an MMP3-specific cleavage site. This design will maintain cyclic peptide sealing throughout the body while enabling localized, specific activation only at the hemorrhoidal lesion site upon MMP3 cleavage.

How will it be implemented ?

To identify suitable MMP3 cleavage sites, we first used the CleaveNet [19] model to predict candidate sites and selected the top 50 highest-scoring sequences. Subsequently, by integrating MMP3 substrate preference information from the MEROPS database [20], we selected the optimal sequence from these 50 as the cleavage site.

Next, we performed structural prediction of this peptide sequence using PEP-FOLD3 [21] and conducted molecular docking with the MMP3 active site using AutoDock Vina [22]. Given that MMP3's catalytic mechanism primarily involves Zn ions promoting hydrolysis reactions, with the reaction core located within the Zn ion coordination pocket, the docking process focused primarily on this pocket. Ultimately, we obtained a best-fit affinity of -8.195 kcal/mol, indicating that this enzyme cleavage site possesses good binding potential.

Figure 11. Enzyme Cleavage Sites and MMP3 Interaction Map.

From the results figure, we can derive the following analysis: Matrix metalloproteinase 3 (MMP3) is a zinc-dependent enzyme that recognizes specific substrate peptides through hydrogen bonds, hydrophobic interactions, and electrostatic forces, thereby promoting their degradation. In Figure 11(A), the substrate peptide (Chain A) extends linearly, forming key interactions with residues in the enzyme's active site (Chain B): Asn668 forms a hydrogen bond with the carbonyl group of Ala3, His705 forms an amide bond with Leu4, Glu702 forms a salt bridge with Val6, and Val698 and and Leu664 with Leu7. Additionally, Pro656 and Tyr723 constrain the peptide into a β-turn conformation, exposing the cleavage site Ala5-Leu6 to the zinc-coordinated catalytic center. Figure 11(B) shows the substrate peptide (magenta) embedded in the S2–S3' active site of the full-length MMP3 structure (orange), with N-terminal Tyr10 and Pro9 sandwiched by Tyr722 and Pro721, and C-terminal Gly8 and Leu7 located in the S1 pocket. The substrate-induced conformational change reduces solvent exposure and promotes the catalytic process. Overall, the selected cleavage sites demonstrate excellent cleavage potential.

Finally, we incorporated the identified MMP3 cleavage sites into the linker peptide. Flexible (GGS) sequences were employed at both ends of the linker to minimize interference with the flanking structures, while a rigid (PA) sequence was used in the central region. The cleavage sites were inserted within this rigid segment, thereby maintaining overall conformation stability while providing a controllable activation mechanism for the cyclic peptide.

Figure 12. Schematic of Anti-VEGF Antibody with Masking Peptide.

Dissociation Capability Analysis

In Model 1, we successfully designed a masking peptide and screened for a complex with moderate affinity, suggesting its potential dissociation capability.

To accurately reconstruct the free energy profile of the masking peptide dissociation process, we designed a rigorous computational workflow covering three core steps: pre-simulation preparation, Steered Molecular Dynamics (SMD) simulation, and Umbrella Sampling with WHAM analysis.

The goal of pre-simulation preparation is to generate a structurally sound and well-equilibrated protein complex, and to construct a coordinate system and simulation box suitable for steered simulation. This involves system equilibration and affine transformation matrix calculation.

System Equilibration

A five-step steepest descent energy minimization was employed to ensure that the system's energy is minimized and equilibrated with respect to temperature, pressure, and density, thereby obtaining stable initial coordinates. This equilibration was performed using AMBER [30] with the ff14SB [26] protein force field and TIP3P [27] water model, which are widely used and validated for protein simulations.

Calculation of the Affine Transformation Matrix

To simplify the setup of restraints and the pulling simulation in GROMACS [28], we calculate a 4x4 affine transformation matrix M to perform a coordinate transformation on the system. This aligns the antibody core with the origin and places the peptide chain along the Z-axis for pulling.

I. Define the core vector

The pulling force direction vector is defined by calculating the vector from the antibody core c_ab to the peptide centroid c_pep.

$$v_{pull} = c_{pep} - c_{ab}$$

The goal is to align its coordinates with the Z-axis unit vector.

$$z = \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix}$$

II. Solving for Rotation and Translation

For any point p, after the coordinate transformation, p' is:

$$p' = R(p - c_{ab})$$

Where the rotation matrix R satisfies:

$$R \frac{v_{pull}}{||v_{pull}||} = z$$

III. Constructing the Affine Matrix

The translation vector t is defined as:

$$t = -R c_{ab}$$

The affine matrix M is:

$$M = \begin{pmatrix} R & t \\ 0^T & 1 \end{pmatrix}$$

This matrix can be directly applied to the homogeneous coordinates of all particles in the system to perform translation and rotation.

Simulation Box Construction

To avoid periodic boundary artifacts, the simulation box must be carefully set up. The overall size and center coordinates of the system are measured in VMD, then the box size is calculated accordingly.

I. X and Y-axis dimensions

The dimensions in X and Y directions are determined by the molecular size plus buffer:

$$L_{box,X} = L_{mol,X} + 2 \times D_{buffer,XY}$$

$$L_{box,Y} = L_{mol,Y} + 2 \times D_{buffer,XY}$$

Where:

• \(L_{mol,X}\) and \(L_{mol,Y}\) are the molecular lengths along X and Y axes.
• \(D_{buffer,XY} = 12.0Å\) is the buffer distance set on X and Y axes.

II. Z-axis dimension

The Z-axis is the pulling direction; its length must accommodate the steered simulation. The final Z-length is chosen as:

$$L_{box,Z} = \max(L_{mol,Z} + \Delta z_{pull} + 2 \times D_{safe}, 2 \times \Delta z_{pull})$$

Where:

• \(L_{mol,Z}\) is the molecule's initial Z-length.
• \(\Delta z_{pull} = 30.0 Å\) is the total pulling distance.
• \(D_{safe} = 12.0 Å\) is the safe buffer distance.

Steered molecular dynamics (SMD) simulations were performed in GROMACS to generate the dissociation trajectory of the masking peptide, providing initial conformations for umbrella sampling. Key parameters include a pulling rate of 1 nm/ns and a harmonic force constant of 500 kJ·mol^-1·nm^-2. The SMD simulations employed the AMBER99SB-ILDN [29] force field and TIP3P [27] water model, following established tutorials for theoretical and practical guidance [24].

The system was solvated and neutralized with Na⁺ and Cl^- ions (0.1 M), followed by energy minimization and pre-equilibration. A 500 ps SMD pulling simulation was conducted to generate the dissociation pathway, ensuring near-equilibrium pulling.

• pull_coord1_geometry = distance: Reaction coordinate is the distance between the centers of mass of the two defined groups.
• pull_coord1_dim = N N Y: Pulling along Z-axis only.
• pull_coord1_k = 500 kJ·mol^-1·nm^-2: Harmonic potential constant.

The trajectory shows successful separation of the masking peptide and antibody along the intended pulling direction.

Figure 13. Animated trajectory of the steered molecular dynamics simulation.

To fully sample the dissociation pathway, we performed the following steps:

1. Window Selection: 63 snapshots along the reaction coordinate (1.90-5.00 nm) were selected at 0.05 nm intervals.
2. Umbrella Sampling: Each window: 100 ps pre-equilibration + 1 ns production MD with harmonic bias to restrain COM distance.
3. WHAM Analysis: Processed all windows using WHAM to reconstruct the PMF curve, removing bias from umbrella potentials.

Figure 14. Potential of Mean Force (PMF) curve for masking peptide dissociation.

4. Binding Free Energy Calculation: The global minimum at ξ = 1.90 nm (PMF = -2.04 kcal·mol⁻¹) corresponds to the stable bound conformation. Plateau region at ξ > 3.9 nm has average PMF = 7.76 kcal·mol⁻¹.

   $$\Delta G_{bind} = G_{bound} - G_{dissociated} = -2.04 - 7.76 = -9.80 \text{ kcal/mol}$$

5. Dissociation Constant Calculation: The relationship between binding free energy and the dissociation constant K_d is as follows:

   $$\Delta G_{bind} = RT\ln(K_d)$$

   According to this formula, at T=310 K, the calculated K_d is 1.22 × 10^-7 M = 122 nM.

6. Verification of Sampling Sufficiency: Histogram overlap confirms effective sampling, validating PMF reliability.

Figure 15. Umbrella sampling histograms showing overlap between adjacent windows.

Umbrella sampling combined with WHAM allowed quantitative determination of binding thermodynamics. The calculated binding free energy is -9.80 kcal·mol⁻¹ with K_d = 122 nM, indicating a stable yet reversible interaction. The nanomolar-level affinity confirms that the designed complex fulfills the goals of conditional dissociation: stability in circulation with potential activation at the target site. Histogram overlap demonstrates sampling reliability, providing a solid quantitative foundation for further structure-based optimization of Probody molecules.

To facilitate reproducibility and provide practical guidance, we also developed a detailed PDF tutorial presenting a comprehensive, optimized computational workflow and automated scripts. This resource is designed to assist iGEM teams in efficiently calculating binding free energies of protein-ligand or peptide complexes using Steered Molecular Dynamics (SMD) and Umbrella Sampling (US) methods. The tutorial addresses common challenges in molecular dynamics simulations-including system orientation, precise simulation box construction, and automation of multi-step equilibration-lowering the technical barrier for teams applying computational approaches in biological research.

Stability Analysis

In a normal assembly, the linker of the Probody is expected to stably hold the masking peptide and the antibody together, thereby preventing non-specific dissociation of the masking peptide and maintaining the antibody in a "dormant" state prior to activation. Verifying the binding stability between the masking peptide and the linker is therefore critical for the rational design of the Probody. To this end, we employed Hydrogen Mass Repartitioning (HMR) technology to accelerate molecular dynamics simulations and assess system stability under physiological conditions.

System Preparation with HMR

For long-timescale simulations, a complete system consisting of the antibody, masking peptide, and linker was prepared. Hydrogen mass repartitioning (HMR) was applied during topology generation. This method transfers a portion of the mass from heavy atoms to bonded hydrogen atoms, lowering the highest vibrational frequencies and enabling a time step of 4 fs. This allows an approximately 1.6-fold acceleration in sampling while conserving energy, enabling efficient long-timescale molecular dynamics [25].

Long-Timescale MD Simulation

Following energy minimization and equilibration, a 500 ns production MD simulation was conducted. To quantitatively evaluate the stability of the complex, we performed two main trajectory analyses:

I. Overall Conformational Stability Analysis (RMSD)

As shown in the RMSD plot (Figure 16), the backbone RMSD of the system rises rapidly during the first 50 ns and subsequently reaches a plateau of 11–12 Å. Although this seemingly high RMSD might suggest structural instability, visual inspection of the trajectory revealed that most fluctuations arise from the highly flexible signal peptide loop region on the antibody. After excluding this flexible segment, the core region of the complex-including the antibody, masking peptide, and linker regions that interact with them-remains highly stable throughout the 500 ns simulation. This demonstrates that the core structure of the designed Probody exhibits excellent conformational stability.

Figure 16. RMSD of the complex over time.

II. Masking Peptide-Antibody Distance and Linker Function Analysis

The distance between the center of mass of the antibody and the masking peptide (Figure 17) provides further evidence of system stability. During the initial 0–100 ns, the distance decreases from ~23 Å and stabilizes within 19–21 Å. This indicates that the masking peptide, constrained by the linker, maintains close contact with the antibody surface, forming a relatively stable conformation in the core region. Throughout the 500 ns simulation, no significant dissociation events were observed, confirming that the linker effectively prevents premature detachment of the masking peptide.

Figure 17. Distance between antibody and masking peptide over time.

The 500 ns MD simulation demonstrates that the Probody maintains a stable overall structure under the constraint of the linker. While the flexible signal peptide loop contributes to an elevated overall RMSD, the core region of the complex remains highly stable. Importantly, the linker effectively secures the masking peptide on the antibody surface, preventing dissociation during the simulation. These results verify the rationality and stability of the Probody design at the atomic level, confirming that the linker preserves the intended "dormant" or "masked" state under physiological conditions.

Although the 500 ns simulation provides preliminary validation, this timescale may not fully capture rare conformational changes or potential dissociation events. Longer simulations, on the microsecond (µs) scale, generally provide higher statistical and physical significance. Future work will involve extended simulations and the application of enhanced sampling methods, such as Replica Exchange MD or Metadynamics, to further verify the stability and dynamic behavior of the linker-constrained Probody system.