Overview

This study aims to design and computationally validate a novel BioPROTAC [1] that efficiently and selectively degrades the key transcription factor HIF-1α under hypoxic conditions, thereby blocking VEGF expression and providing an innovative strategy to alleviate the pathological progression of hemorrhoids. Through a series of computational biology models, we systematically validated the feasibility of this BioPROTAC system at both the structural and kinetic levels.

First, we successfully predicted the complete three-dimensional structure of BioPROTAC using AlphaFold 3.0 [2] and Modeller [3]. Combined with 100 ns Amber [4] molecular dynamics simulations and Ramachandran diagram analysis, we confirmed that the geometric structure of this design is highly rational and stable.

Subsequently, during the construction of the BioPROTAC mediated ubiquitination system complex, molecular docking results revealed the BioPROTAC's own spatial masking effect on the classical ubiquitination sites of HIF-1α (K532, K538, K547) [5]. To address this critical issue, we employed PROPKA 3 [6] for evaluation and identified Lys11 as a novel ubiquitinylation site with high reactivity potential, characterized by a significantly lower pKa (6.96) and high surface exposure (98%).

Further 100 ns unconstrained MD simulation results indicate that the large ubiquitinating complex composed of eight subunits tends toward stability under dynamic conditions. More importantly, the Lys11 site on the catalytic core of HIF-1α and the ubiquitin (Ub) reaction center (Gly76) maintained a favorable distance of approximately 3.3 Å over extended periods, providing strong evidence that the BioPROTAC effectively recruits the substrate to the ubiquitination reaction center.

Ultimately, we conducted a quantitative assessment of the kinetic feasibility of the ubiquitination reaction. Analysis of MD trajectories revealed the existence of a transient, short-range attack conformation with an interatomic distance of approximately 3 Å. This thermodynamically non-favorable conformation is nevertheless a critical prerequisite for nucleophilic attack. Based on this, we employed the Gaussian 16 QST3 [7] method to search for the transition state at the B3LYP/6-31G(d) level in conjunction with the SMD solvent model. Its validity was confirmed through IRC calculations. The final calculated reaction energy barrier was determined to be 8.20 kcal/mol. This value falls within the typical range for enzyme-catalyzed reaction energy barriers, providing strong theoretical support. It demonstrates that the designed BioPROTAC is not only geometrically feasible but also possesses the potential for efficient degradation of HIF-1α in terms of reaction kinetics.

Model 1: Design of BioPROTAC

In preliminary studies and literature reviews, we found that the development of hemorrhoids is closely associated with local inflammatory responses and abnormal vascular proliferation. Among these, vascular endothelial growth factor (VEGF) is regulated by hypoxia-inducible factor HIF-1α [8]. As a key transcription factor in cellular responses to hypoxia, HIF-1α activates multiple genes related to inflammation and angiogenesis, thereby promoting the maintenance of the inflammatory microenvironment. Therefore, selective degradation of HIF-1α holds promise for blocking inflammatory pathways at their source and inhibiting angiogenesis. Based on this, we designed and constructed a BioPROTAC aimed at efficiently degrading HIF-1α under hypoxic conditions, thereby alleviating local inflammation and abnormal vascular proliferation in hemorrhoids.

Building upon the aforementioned research foundation, we designed and constructed a BioPROTAC with targeted degradation capabilities. This compound is engineered to selectively degrade HIF-1α under hypoxic conditions, thereby alleviating local inflammation and abnormal vascular proliferation associated with hemorrhoids.

This BioPROTAC primarily consists of three functional parts: First, the HIF-1α-binding protein specifically recognizes and binds the target protein HIF-1α, which is crucial for achieving selective degradation. Second, the E3 ubiquitin ligase recruitment part binds to the E3 ubiquitin ligase complex, triggering the ubiquitination of HIF-1α and subsequent proteasomal degradation. To ensure sufficient spatial flexibility and coordination between these two structural components, we introduced a moderately sized linker peptide to optimize the overall molecular conformation and minimize structural interference. Furthermore, we incorporated a nuclear localization signal (NLS) [9] into the design to enhance the BioPROTAC's intracellular targeting efficiency and functional activity under hypoxic conditions, thereby improving its overall degradation efficacy.

Model 2: Structure Prediction of BioPROTAC

Since the crystal structure of our currently designed BioPROTAC has not been resolved, we aim to predict its structure through rational bioinformatics methods, thereby providing structural basis and guidance for subsequent functional validation.

First, we obtained the amino acid sequence of VHH212 from the literature [10] and attempted to compare the 3D structures of VHH212 reconstructed using AlphaFold 3.0 [2] and homology modeling based on this amino acid sequence. Following a BLAST search [11] in the NCBI database, we selected templates with an E-value cutoff of 1e-50 and a sequence identity cutoff greater than 70%. These templates were then modeled using Modeller to generate candidate structures. After meticulously comparing prediction quality, we selected the 3D structure predicted by AlphaFold 3.0. We then performed 100 ns molecular dynamics simulations using Amber [4] to eliminate unreasonable structures.

From the RMSD results of molecular dynamics simulations, we observe that VHH212 exhibits overall high stability, with its RMSD consistently maintaining at 8.2 Å. Minor fluctuations in RMSD arise from the flexible CDR regions. We extracted stable conformations after 60 ns to serve as the stable structure for VHH212.

After obtaining a reasonable VHH212 structure, we retrieved the X-ray crystal structure of pVHL from the RCSB PDB database [12]. Using both structures as references, we integrated the NLS nuclear penetration tag and linker. We obtained the 3D structure of the BioPROTAC and analyzed it to generate its Ramachandran plot.

From the Ramachandran plot, it can be observed that nearly all amino acid residues in the BioPROTAC structure fall within the permissible regions, with the majority residing in the favorable zone. Only one residue lies outside these regions, indicating the structure's validity. Nevertheless, given that molecular dynamics simulations will be conducted later, we still consider this a highly reliable structural prediction result.

Model 3: Molecular Docking Simulation

Under physiological conditions in the human body, the ubiquitinylation reaction,as a precisely regulated protein modification process,relies primarily on the correct assembly of the ubiquitin ligase complex (E3 ligase complex) for its orderly execution. The E3 complex not only plays a central role in substrate recognition but also directly determines whether ubiquitin molecules can approach the substrate's ubiquitinylation site with the appropriate spatial conformation, thereby enabling efficient ubiquitin transfer. Based on this mechanism, we attempted to rationally model the BioPROTAC-mediated ubiquitinylation system using molecular docking methods. This approach aims to reveal key spatial constraints at the structural level and provide guidance for subsequent functional validation.

First, since the binding site of the substrate ligand VHH212 in our BioPROTAC and the complex structure with HIF-1α remain unresolved, obtaining interaction data between the two is urgently needed. Previous studies have reported that VHH212 specifically binds to the PAS-B domain of HIF-1α [10]. Based on this, we employed HADDOCK [13] for molecular docking to screen for the optimal binding conformation, yielding a complex model of VHH212-HIF-1α-PAS-B domain. The three-dimensional coordinates of the HIF-1α-PAS-B domain were derived from the crystal structure of the HIF-1/ARNT PAS-B complex (PDB ID: 4H6J), while the coordinates of VHH212 were obtained from a reasonable conformation after molecular dynamics simulation.

Subsequently, we analyzed the VHH212 amino acid residues and functional hotspots within 5 Å of the HIF-1α-PAS-B domain using SpotON [14] and PyMol. These sites will serve as a guide for subsequent molecular docking.

Next, we retrieved and obtained the required crystal structures from the RCSB PDB database [12] as modeling references. The BioPROTAC-Elongin B and BioPROTAC-Elongin C complex structures referenced the pVHL/Elongin C/Elongin B crystal structure (PDB ID: 1LQB); For the interactions between Elongin B, Elongin C, and Cullin-2, as well as between Cullin-2 and Rbx1, we referenced the crystal structure of the Cul2-Rbx1-EloBC-VHL complex (PDB ID: 5N4W). The docking structure of UbcH5C with Rbx1 was derived from the highly homologous UbcH5B-Rbx1 complex crystal structure (PDB ID: 6TTU). Using PyMOL, UbcH5C was superimposed onto the coordinates of UbcH5B to obtain a reasonable UbcH5C–Rbx1 complex model. The binding of UbcH5C to ubiquitin was modeled by referencing the crystal structure of the UbcH5C–Ubiquitin complex (PDB ID: 4BVU). It should be noted that in this crystal structure, UbcH5C and ubiquitin are covalently linked via an oxygen ester bond, whereas in the active complex within the human body, they are actually covalently linked via a thioester bond [15]. We speculate this discrepancy may stem from the higher stability of the oxyl ester bond in the crystalline environment, which facilitates crystal preparation and X-ray diffraction measurements. To obtain a model more consistent with physiological conditions, we replaced S85 in UbcH5C with cysteine in PyMOL, thereby generating a sulfur ester bond-linked UbcH5C–Ubiquitin model based on the crystal structure.

Initially, we hypothesized that BioPROTAC could promote modification of HIF-1α at the classical ubiquitinylation sites mediated by VHL (K532, K538, and K547) [5]. However, molecular docking results revealed that BioPROTAC and Elongin B/C block these lysine residues during binding. The longer cyclic Cullin-2 cannot bend sufficiently to guide the ubiquitin molecule to these sites, thereby limiting ubiquitin accessibility. This finding indicates that, in the current conformation, these classical sites are no longer viable as effective ubiquitinylation sites.

Through literature review, we have learned that the application of PROTAC technology often leads to alterations in ubiquitinylation sites [16][17]. This primarily stems from changes in the conformation of the linker domain and increased steric hindrance, rendering the original ubiquitinylation sites inoperative. Therefore, we urgently need to explore new target residues with ubiquitinylation potential.

We utilized GPS-Uber [18] to predict a series of potential ubiquitination sites. According to literature reports, the nature of the ubiquitination reaction can be regarded as a nucleophilic addition process [15]. Under physiological conditions, the ε-amino group of lysine residues (Lys) undergoes deprotonation to convert into the –NH₂ form, thereby acquiring enhanced nucleophilicity. This activated amino group can nucleophilically attack the thioester bond in the ubiquitin-E2 complex, enabling ubiquitin transfer. Concurrently, the deprotonation of lysine is often regulated by its amino acid environment, such as adjacent carboxyl groups or water molecules, thereby influencing its activity and efficiency as a ubiquitinylation site.

Therefore, we employed the PROPKA 3 [6] tool to assess the deprotonation capacity of candidate sites and thereby evaluate their ubiquitinylation potential. We found that the pKa value of candidate site Lys11 was 6.96, significantly lower than the typical value of approximately 10.5 for lysine residues. This indicates that its ε-amino group remains largely deprotonated (–NH₂) at physiological pH (approximately 7.4)[31] at the site of action, thereby favoring the formation of an isopeptide bond with the carboxyl group at the C-terminus of ubiquitin.

The structural environment further supports this conclusion. First, Lys11 exhibits a high surface exposure of 98%, indicating that this residue is spatially accessible and readily recognized by the E2–E3 complex. Second, its strong desolvation effect (–3.46) suggests that this site resides in a relatively low-dielectric-constant environment, effectively lowering the pKa value and stabilizing the deprotonated state. Furthermore, the absence of significant hydrogen bonds prevents potential stabilizing interactions from hindering the reaction. Finally, weak Coulombic interactions with Arg89 and Lys65 (–0.01 and –0.07, respectively) provide additional electrostatic shielding, further facilitating amino deprotonation.

In summary, Lys11 possesses a low pKa, high exposure, and a favorable charge environment, and can therefore be considered a ubiquitinylation site with high reaction potential.

Ultimately, based on the reference structure and the newly identified ubiquitinylation site, we successfully obtained the complete conformation of the ubiquitinylation system using HADDOCK [13] and HDOCK [19]. This result not only validates the spatial feasibility of the selected site but also provides crucial insights into the potential mechanisms underlying substrate recognition and ubiquitin transfer. Consequently, it lays the groundwork for subsequent experimental validation and functional studies.

Model 4: Molecular Dynamics Simulations

Since structures obtained from molecular docking represent only static conformations, they often contain steric conflicts or unreasonable conformational features that may not be stable under physiological conditions. Molecular dynamics simulations can eliminate these unnatural tensions under physical force fields, allowing the system to relax into more reasonable states that better approximate real physiological environments. More importantly, molecular systems within cells are inherently dynamic. We need to evaluate through simulation whether the core reaction region can sustain a spatial conformation conducive to reaction occurrence throughout its motion.

Prior to molecular dynamics simulations, we employed the H++ [20] program to determine the protonation states of ionizable residues within the complex, supplemented by manual verification considering local electrostatic interactions and hydrogen bonding environments. For substrate HIF-1α's K 11, we adopted its deprotonated state under physiological conditions in the simulation. To handle glycine and cysteine linked by thioester bonds, we obtained electrostatic potentials based on HF/6-31G* level quantum chemical calculations and fitted charges using AMBER's built-in RESP method [21]. Subsequently, the complex system was solvated in a truncated octahedral box of TIP3P water molecules [22] under the Amber19SB force field [23], with a solvent layer thickness of 20 Å.

After system construction, two-step energy minimization was performed to eliminate unbound conformations: first, 10,000 steps of minimization under solute-constrained conditions, followed by 20,000 steps of unconstrained minimization. All molecular dynamics simulations were performed using the Amber software suite under NPT ensemble conditions with periodic boundary conditions. The simulation workflow comprised: (1) NVT temperature ramping (1 ns, from 100 K to 310 K); (2) NPT density equilibration (2 ns, under identical constraints); (3) Stepwise NPT equilibration with gradually released constraints (four steps, 2 ns each, with backbone constraint weights decreasing sequentially from 10.0 to 0.5 kcal/mol·Å²); (4) 100 ns unconstrained NPT production simulation.

Electrostatic interactions were calculated using the particle-mesh Ewald [24] method, with hydrogen bond lengths constrained by the SHAKE algorithm [25]. The simulation employed an integration time step of 2 fs, a cutoff radius of 8.0 Å for non-bonding interactions, and non-bonding pairs updated every 25 steps. The system temperature was controlled to 300 K via Langevin dynamics (collision frequency 5.0 ps⁻¹), while pressure was maintained at 1 atm (1 atm = 101.3 kPa) using a Monte Carlo barometer. The pressure coupling time constant was set to 1.0 ps.

Molecular dynamics simulations indicate that the constructed BioPROTAC-Cullin2-Rbx1-E2-Ub-HIF-1α complex gradually relaxes and stabilizes during the simulation process. The RMSD curves for the overall system and subunits (Figure 10) reveal significant conformational adjustments occurring within the first approximately 20 ns, followed by relative stability throughout the remaining 80 ns of simulation. Due to the large size of the system and its high computational demands, only a 100 ns production phase simulation was run. Although RMSD results show some fluctuations, our complex consists of eight components, making it complex and large with numerous flexible regions. Minor RMSD fluctuations are unavoidable. Moreover, the complex RMSD generally remained around 10 Å, leading us to conclude that the complex approached a stable state during the simulation.

More critically, the catalytic core region maintained favorable spatial geometry throughout the simulation. Lys11 on HIF-1α and Gly76 on Ub remained within a distance range of approximately 4 Å for an extended period (Figure 11), indicating that the designed BioPROTAC effectively recruits the substrate to the ubiquitinylation reaction center and maintains an orientation conducive to ubiquitin transfer. The overall results validate the dynamic stability of this complex under near-physiological conditions, supporting the feasibility of the designed system to promote efficient ubiquitination reactions in cellular environments. This lays a solid foundation for subsequent reaction energy barrier analysis.

Model 5: Energy Barrier Assessment

Whether ubiquitination occurs depends not only on spatial conformation but also on the energy required for the reaction. To verify its kinetic feasibility, we plan to calculate the reaction energy barrier for lysine attacking the Ub-E2 thioester bond to determine whether the ubiquitination reaction is feasible under physiological conditions.

First, we need to identify representative conformations from the molecular dynamics trajectory to serve as initial coordinates for subsequent QM region calculations [15]. The QM region must encompass the most critical reactive residues. Therefore, we analyzed the state of the ε-amino group of lysine and the thioester bond throughout the entire molecular dynamics simulation.

We analyzed the dynamic behavior and conformational transition characteristics of the distance between the ϵ-amino group and the carbon atom in the system. We discovered that this system exhibits a bistable conformational space, maintaining rapid dynamic equilibrium between two states. The distance trajectory in Figure 12(a) demonstrates that throughout the entire 100-nanosecond simulation period, this distance continuously undergoes conformational jumps between two distinct macroscopic states. These two primary states are defined by interatomic distances of approximately 4 Å and 10 Å, respectively. Figure 12(b) illustrates the distance variation during a representative sampling episode, where the distance fluctuates around 4 Å—the primary pre-attack conformation. Some of these fluctuations may represent high-energy transients attempting but failing to cross the energy barrier to the 10 Å state. The distance probability density distribution in Figure 12(c), visualized via Kernel Density Estimation (KDE), unequivocally confirms the system's bistable conformational preference. The first peak centered around 4 Å represents the short-range conformation, while the second peak clusters around 10 Å. The deep valley between these peaks at approximately 6 Å signifies the free energy barrier between the two stable conformations. We emphasize here that transition state theory and studies of enzyme catalytic mechanisms consistently indicate that the attacking conformation in many enzymatic reactions is not the system's globally most stable state, but rather an instantaneous minority state formed through spatial confinement by the enzyme or molecular dynamics fluctuations [26]. Similarly, Plechanovová et al. demonstrated in their structural study of RNF4 [27] that lysine's nucleophilic attack on the Ub–E2 thioester bond depends on precise attack angles and distances, with this pre-reaction conformation often requiring dynamic fluctuations to achieve. Consistent with this, the approximately 4 Å short-range state revealed here, though thermodynamically non-favorable, possesses geometric features that precisely satisfy the critical attack conditions required for the ubiquitinylation reaction. This result indicates that the pre-reaction state of ubiquitinylation inherently belongs to a thermodynamically non-favorable conformation, yet its transient existence provides the necessary spatial prerequisites and kinetic feasibility for the reaction to proceed.

Next, we extracted the most favorable attack conformation from the aforementioned 4 Å states and constructed a reasonable QM region to prepare for energy barrier mapping. We first isolated the core residues involved in the reaction: K11 and the thioester bond residue of E2-Ub. The bond was broken exclusively at the carbon-carbon link, with hydrogen atoms added at the ends to ensure the system's charge balance remained intact. Subsequently, we constructed the product conformation according to the reaction mechanism: the thioester bond cleaved to form a thiol group, and a new isopeptide bond formed. Next, the initial conformations of reactants and products were optimized using Gaussian 16 software [7]. Calculations employed the B3LYP density functional with the 6-31G(d) basis set [28], incorporating D3BJ dispersion correction [30]. Solvent effects were introduced using the SMD water solvent model [29]. Following geometric optimization, frequency analysis was performed to confirm the system resides at a genuine potential energy minimum. The transition state structure was searched using the QST3 method in Gaussian 16, employing the same parameters as in the optimization section. Reagents, products, and initial transition state structures were input simultaneously during the calculation. Gaussian performed an interpolation search and optimized the structure to the saddle point. Frequency analysis revealed a unique imaginary frequency (−115.8 cm⁻¹). The corresponding vibrational mode corresponds to the nucleophilic attack of the ε-amino group of lysine onto the Ub–E2 thioester bond, further confirming that the obtained structure represents the transition state for the target reaction.

Table 1. Vibrational Frequencies of the Transition State

Finally, to validate the validity of the transition state, we performed an intrinsic reaction coordinate (IRC) calculation. The calculation employed the B3LYP/6-31G(d) method, incorporating D3BJ dispersion corrections and the SMD aqueous solvent model. At the transition state point, the Hessian was computed to determine the initial direction of the reaction pathway.

Next, we analyzed the results obtained after IRC and drew the following conclusions: Figure 13(a) depicts the variation of the total energy of the system along the IRC path, showing an inverted U-shaped energy profile. The peak of the curve occurs at IRC = 0.00, indicating that this stationary point represents the maximum energy along the path connecting reactants and products, confirming its status as a transition state. Energy decreases monotonically from the transition state toward both sides, converging toward more stable configurations. Complementing this, Figure 13(b) shows the variation of the root-mean-square gradient norm along the IRC path. The gradient norm reaches its minimum at IRC = 0.00, fully consistent with the theoretical requirements for a transition state, where net forces in all directions approach zero except for imaginary-frequency vibrations along the reaction coordinate. As the IRC coordinate moves away from the transition state in both directions, the slope of the potential energy surface increases rapidly, causing the gradient norm to rise sharply in both directions. In summary, the results of the total energy and gradient norm variations along the IRC path mutually validate each other, sufficiently demonstrating the validity of the transition state explored by QST3. This finding can guide subsequent calculations of the energy barrier.

We then performed precise single-point energy calculations for these three states using Gaussian and plotted the energy barrier diagram based on the obtained data.

We ultimately determined the energy barrier for this reaction to be 8.20 kcal/mol. This value falls within a reasonable range, indicating that our designed BioPROTAC not only achieves effective proximity in geometric conformation but also demonstrates reactivity feasibility. This provides theoretical support for its potential biological applications.

Initially, we planned to employ the CP2K–GROMACS coupled QM/MM method to calculate energy barriers with higher precision and further validate the feasibility of the reaction. This method possesses robust capabilities, simultaneously accounting for reaction characteristics in the quantum region and the influence of the surrounding protein environment. However, constrained by computational resources and time limitations, we were unable to complete this portion of the study. In subsequent work, we will continue to pursue this approach to obtain results closer to the actual system.