Structural Basis for TGF-β Mimetic Peptide-Induced Signaling Activation Through Molecular Dynamics Simulations

TL;DR

This study uses simulations to explain how two TGF-β mimetic peptides interact differently with receptors, showing one is more effective in activating a key signaling pathway.

Contribution

The study reveals the structural dynamics of TGF-β mimetic peptides and identifies a key interaction for signaling activation.

Findings

TB2 stabilizes a complex with TβRII and promotes TβRI engagement through a hydrogen bond.

The TB2–TβRII–TβRI trimer has higher binding affinity and stability compared to TB1.

TB2 activates the TGF-β/Smad pathway, as shown by Smad3 expression and phosphorylation.

Abstract

Transforming growth factor-β (TGF-β) mimetic peptides offer significant therapeutic potential due to their superior pharmacological properties over the native cytokine. Our previous work identified two such peptides, TB1 and TB2, which bind to the type II TGF-β receptor (TβRII) yet elicit distinct cellular responses. To uncover the mechanistic basis for the functional divergence, we employed integrated molecular dynamics (MD) simulations with the AlphaFold3-predicted structures. Our analytical results indicated that TB2 stabilizes a dynamic complex with TβRII and is predicted to facilitate type I receptor (TβRI) engagement possibly involving a critical hydrogen bond between TB2-Gly11 and TβRI-Phe60. The resulting trimeric assembly (TB2–TβRII–TβRI) exhibits a higher relative binding affinity (−67.76 ± 7.70 kcal/mol) and structural stability. In contrast, the TB1–TβRII complex fails to productively engage TβRI. These computational results were experimentally validated. Western blot analysis confirmed that TB2, but not TB1, activates the canonical TGF-β/Smad pathway by enhancing the expression and phosphorylation of Smad3. This study will elucidate the dynamic structural basis for the activity of TGF-β mimetic peptides and suggest TB2 as a promising lead candidate for the rational design of tissue-regenerative therapeutics.

Figures

Full text

1. Introduction

Transforming growth factor-beta (TGF-β) is a critical cytokine that regulates diverse cellular processes, including proliferation, differentiation, apoptosis, migration, and immune responses [1,2,3,4]. Its activity is mediated through binding to transmembrane receptors, primarily the high-affinity TGF-β type II receptor (TβRII) and the low-affinity TGF-β type I receptor (TβRI) [5,6,7]. Due to its biological activity, TGF-β holds considerable therapeutic potential in chronic wound healing, osteoarthritis treatment, and tissue repair [8,9,10,11,12,13,14]. During the regeneration of tissues, TGF-β promotes repair by stimulating collagen synthesis and extracellular matrix (ECM) deposition [13,15].

TGF-β signaling is mediated by both Smad-dependent and Smad-independent pathways [1,16]. In the canonical Smad-dependent pathway, TGF-β ligand binding to TβRII leads to the recruitment of TβRI, which subsequently phosphorylates Smad2/3 at their C-terminal domains [17,18]. Phosphorylated Smad2/3 form a complex with Smad4 and translocate into the nucleus to regulate the expression of target genes [19,20] (Scheme 1). In contrast, Smad-independent pathways trigger signaling cascades such as MAPK (ERK/JNK/p38) and PI3K-Akt, which modulate cytoskeletal dynamics, apoptosis, and immune responses [2,21,22,23]. The interplay between these pathways enables TGF-β signaling to respond dynamically to changes in the cellular microenvironment [2].

Current drug development targeting the TGF-β pathway exhibits a distinct trend toward inhibitors—such as neutralizing antibodies, ligand traps, and kinase inhibitors—which aim to block pathway hyperactivation in cancer and fibrosis [2,24]. In contrast, therapeutic agonists capable of promoting tissue repair are scarce. While recombinant proteins like Avotermin (rhTGF-β3) have shown promise in clinical trials for scar improvement [25], their widespread application faces challenges inherent to protein therapeutics, such as high production costs, complex production processes, immunogenicity and poor tissue penetration [26]. Peptide-based therapeutics have attracted increasing interest due to their stability, ease of synthesis, and efficient tissue delivery [27,28,29,30]. To develop peptides that mimic the biological functions of TGF-β and activate the TGF-β/Smad3 pathway, our laboratory employed phage display technology [31,32] targeting the extracellular domain of TβRII, which led to the identification of two peptides: TB1 (LDSQAPMGRINH) and TB2 (KLHHHLHVPRGP). In vitro experiments demonstrated that the two peptides bind to TβRII but exhibit distinct biological activities (these data are derived from a separate, ongoing study): TB1 promotes proliferation of 3T3 cells, whereas TB2 not only stimulates 3T3 cell proliferation but also enhances the expression of type I and type III collagen. However, the molecular mechanisms and binding profiles for these different functions remain unclear. Elucidating these mechanisms will benefit the development of tissue repair by the peptide-based therapeutics.

In this study, we combined AlphaFold3 [33,34,35] predictions with molecular dynamics (MD) simulations to investigate the structural basis of TB1 and TB2 binding to TβRII and the structural basis for their capacity to recruit TβRI. Although AlphaFold3 significantly improves the accuracy of protein–ligand interaction modeling [35], it primarily yields static structural snapshots. Biological molecules exert their functions in vivo through dynamic behaviors rather than static structures [36]. To capture the conformational fluctuations essential for receptor activation, we performed extensive MD simulations [37,38,39] on both dimeric (TB1–TβRII, TB2–TβRII) and trimeric (TB1–TβRII–TβRI, TB2–TβRII–TβRI) complexes. We further employed MM-PBSA [40,41,42] calculations to evaluate binding free energies and residue-specific contributions. Finally, Western blot analysis was conducted to experimentally validate the activation of the TGF-β/Smad pathway by the peptides. Our integrated computational and experimental approach reveals that TB2 stabilizes a dynamic complex with TβRII and actively recruits TβRI, a process facilitated by a critical hydrogen bond between TB2-Gly11 and TβRI-Phe60. In contrast, TB1 is predicted to show limited interaction with TβRI, which may partly account for its weaker downstream signaling activity. These findings not only elucidate the mechanistic basis for the functional divergence between TB1 and TB2 but also highlight TB2 as a promising lead for the rational design of TGF-β-based regenerative therapeutics.

2. Results

2.1. AlphaFold3-Driven Modeling

AlphaFold3 was employed to construct the dimeric complexes TB1-TβRII and TB2-TβRII. As shown in Figure 1, in the two generated structures (Figure 1a,c), most regions of the expected position error (EPE) matrix plots (Figure 1b,d) displayed a relatively uniform dark tone, indicating that the expected positional errors of these residues were small and that the model’s structural predictions for these regions were relatively accurate. For the TB1 model, the pTM score was 0.88 and the ipTM score was 0.60, while for the TB2 model, the pTM score was 0.86 and the ipTM score was 0.46. From the results, the pTM scores are all above 0.8, indicating that the overall structures are relatively reasonable [34,43]. However, for the ipTM score, the TB2 model exhibited a lower ipTM of 0.46, suggesting potential imperfections at the interaction interface. A low ipTM does not inherently invalidate the model; for instance, the TB2 peptide interacts with TβRII primarily through partial N-terminal residues, resulting in fewer atomic contacts and a smaller interface area compared to TB1.

Structural alignment of the predicted complexes with the X-ray crystallographic structure of the TGF-β receptor complex (PDB: 2PJY) revealed partial overlap of the epitopes where TB1 and TB2 interact with TβRII. Notably, the predicted TB2 binding site on TβRII appears to partially resemble that of TGF-β3. In the TB2-TβRII complex, residues of TβRII located within 5 Å of TB2 included Leu27–Asp32 and Ser49–Glu55 (residue numbering based on PDB 2PJY, with TβRII sequence spanning Ala19–Phe126). This region substantially overlapped with the binding domain of TGF-β3 on TβRII in the 2PJY structure. In contrast, in the TB1-TβRII system, residues within 5 Å of TB1 included Leu27–Asp32, Ser49–Glu55, Val77–Asp80, and Pro84–Asp87. The predicted similarity of TB1 and TB2 binding sites to those of TGF-β3 may account for their TGF-β3-like biological activities in promoting cell proliferation (Figure 2a–c).

To further investigate whether the peptide–TβRII complexes could facilitate the recruitment of TβRI, homology modeling was performed using the crystal structures of the TβRI–TβRII complex (PDB: 2L5S and 1PLO) and the TGF-β receptor complex (PDB: 2PJY). Based on these templates, two hypothetical trimeric models were constructed: TB1-TβRII-TβRI and TB2-TβRII-TβRI (Figure 2d–f).

2.2. Stability of Molecular Dynamics Simulations

For each of the four systems (TB1-TβRII, TB2-TβRII, TB1-TβRII-TβRI, and TB2-TβRII-TβRI), three independent replicates of unrestrained production molecular dynamics simulations were performed. After analysis, the results showed good consistency, and a representative trajectory from each system was selected for subsequent analysis. The stability of the four MD trajectories was analyzed in terms of root mean square deviation (RMSD) [44,45,46], root mean square fluctuation (RMSF) [47], radius of gyration (Rg) [48] and the number of hydrogen bonds.

2.2.1. RMSD Analysis

Root mean square deviation (RMSD) quantifies the average displacement between corresponding atoms of two structures, typically using backbone atoms to evaluate conformational alignment and assess the stability of an MD system over time. In the present study, RMSD analysis was applied to assess the equilibrium states and structural stability of the peptide–receptor complexes across the 200 ns trajectories. The results demonstrated that the simulation time was sufficient to capture the dynamic properties of all four systems, each of which eventually reached global equilibrium.

As shown in Figure 3a, TB1-TβRII and TB2-TβRII equilibrated after approximately 70 ns. The average RMSD values were 0.03376 nm for TB1-TβRII and 0.03654 nm for TB2-TβRII. In contrast, TB1-TβRII-TβRI reached global equilibrium only after ~100 ns and exhibited larger fluctuations, with an average RMSD of 0.05011 nm. These results suggest that although the trimeric TB1-TβRII-TβRI system ultimately achieved equilibrium, its conformation continued to undergo substantial dynamic transitions. By comparison, TB2-TβRII-TβRI equilibrated much earlier, at ~40 ns, with minimal fluctuations. The average RMSD was 0.0359 nm, indicating a more stable structural configuration during the simulation.

2.2.2. Rg Analysis

The radius of gyration (Rg) reflects the overall compactness of a macromolecule and is defined as the root mean square distance of atoms from their common center of mass. A smaller Rg value indicates greater structural compactness and stability. To minimize interference from side-chain fluctuations, Rg values were calculated using only the protein backbone atoms.

As shown in Figure 3b, the Rg values of the TB1–TβRII and TB2–TβRII complexes reached equilibrium at approximately 40 ns and subsequently remained stable around 1.4 nm, with fluctuations confined within ~0.01 nm. This stability is consistent with the RMSD trends observed for these two systems. In contrast, the Rg trajectories of the trimeric complexes TB1–TβRII–TβRI and TB2–TβRII–TβRI exhibited distinct behaviors. For TB1–TβRII–TβRI, equilibrium was achieved at around 70 ns, but the Rg values displayed relatively large fluctuations, consistent with the dynamic instability revealed by RMSD analysis.

In the TB2–TβRII–TβRI system, however, the Rg values showed a sharp decline during the initial 0–40 ns, decreasing from 2.25 nm to 2.05 nm. Beyond 40 ns, the system entered a stable equilibrium phase that persisted until 200 ns, with minimal fluctuations. This rapid reduction in Rg indicates that the trimeric complex became more compact within the first 40 ns. These simulation results suggest that the specific binding pattern of TB2 to TβRII enhances its cooperative interaction with TβRI, thereby reinforcing the overall compactness and stability of the system.

2.2.3. RMSF Analysis

Root mean square fluctuation (RMSF) is a key metric for assessing structural flexibility, as it quantifies the average deviation of each atom from its mean position over the course of the simulation. To evaluate residue-specific flexibility and positional variations, RMSF values of the Cα atoms were calculated for all four complexes. For consistency, the analysis was performed using the relatively stable equilibrium phase (100–200 ns) of each trajectory.

As shown in Figure 3c–f, the RMSF profiles of TβRII across the four complexes exhibited a high degree of consistency. The primary flexible region was located at the N-terminus (Ala19–Gln26), which constitutes the interaction interface with TβRI. The flexibility observed in this region may reflect lower predicted stability at the modeled interface of TβRI. Additional flexible segments of TβRII were identified at Thr37–Ser43, Glu55–Gln58, Asn68–Ile72, and Lys104–Gly107. Notably, the RMSF values of TβRII in the TB2–TβRII–TβRI system (orange curve in Figure 3c) were lower than in the other three systems, suggesting that the combined presence of TB2 and TβRI exerts a stabilizing effect.

For TβRI, the RMSF analysis of the two trimeric complexes revealed four major flexible regions: Cys14–Phe22, Ser33–Ser45, Pro55–Thr74, and the C-terminal residues Lys84–Leu87 (Figure 3d).

Figure 3e,f illustrate the RMSF values of TB1 and TB2 in both dimeric and trimeric complexes. For both peptides, the N- and C-terminal residues displayed the highest flexibility, a feature commonly observed for short peptides. The RMSF values of TB1 were higher in the trimeric complex compared to the dimeric system, suggesting that incorporation into the trimer did not restrict TB1’s conformational mobility. This observation may be attributed to the absence of direct interactions between TB1 and TβRI (Figure 3e).

Remarkably, it was observed that in the TB2–TβRII complex, TB2 exhibited reduced flexibility at its N-terminus but pronounced fluctuations at the C-terminus, with the RMSF of the terminal residue reaching 0.2943 nm (Figure 3f). This result implies that the binding of TB2 to TβRII is primarily mediated by its N-terminal region, whereas the C-terminus does not engage in stable interactions, which may also be the reason for the lower ipTM score of the AlphaFold3 TB2–TβRII model. Most notably, in the TB2–TβRII–TβRI system, the RMSF values of the TB2 C-terminus were markedly reduced compared to those in the dimeric complex, with the terminal residue showing an RMSF of only 0.1117 nm. This reduction indicates that the C-terminal region of TB2, due to the presence of TβRI, has its flexibility constrained, and it is likely that the C-terminal region of TB2 may form stabilizing interactions with TβRI.

These findings suggest that although both TB1 and TB2 act as binding peptides for TβRII, their dynamic behaviors differ significantly. The contrasting flexibilities and interaction patterns of TB1 and TB2 may help explain the observed differences, although experimental confirmation is required.

2.2.4. Number of Hydrogen Bonds

The number of hydrogen bonds is a critical parameter for evaluating the strength and stability of intermolecular interactions, with a higher number of hydrogen bonds generally indicating stronger binding affinity. We first analyzed the hydrogen bonds formed between the peptides and TβRII, where the peptides were defined as ligands and TβRII as the receptor. The average number of hydrogen bonds between TB1 and TβRII was 4, whereas TB2 formed an average of 6 hydrogen bonds with TβRII (Figure 3g). These results suggest that TB2 may establish a more stable interaction with TβRII compared with TB1.

In the trimeric systems TB1–TβRII–TβRI and TB2–TβRII–TβRI, hydrogen bond analysis was performed with TβRI considered as the ligand and the TB–TβRII complexes as receptors. The average number of hydrogen bonds formed between TB1–TβRII and TβRI was 4.7, while the TB2–TβRII–TβRI system exhibited a higher value of 6.6 (Figure 3h). The greater number of hydrogen bonds observed for TB2 indicates stronger interactions and enhanced stability of the TB2–TβRII–TβRI complex compared with the TB1 counterpart.

2.3. Inter-Subunit Distances

The above analyses preliminarily suggested that in the TB2–TβRII–TβRI trimeric complex, TβRI interacts with the TB2–TβRII complex, whereas in the TB1–TβRII–TβRI trimer, such interactions appear absent. To further validate these observations, we analyzed the time-dependent distance changes between different subunits in all four systems.

As shown in Figure 4a, the distances between TBs and TβRII in the two dimeric complexes remained stable throughout the simulations. Similarly, in the trimeric systems (Figure 4b), the distances between TBs and TβRII were also stable. Analysis of the distances between TβRII and TβRI in the two trimeric models suggested relatively stable trends as well (Figure 4c). The average inter-receptor distance in the TB1 system was 0.175 nm, while in the TB2 system it was slightly shorter at 0.170 nm.

However, pronounced differences were observed in the distances between the peptides and TβRI (Figure 4d). In the TB2 trimer, the TB2–TβRI distance rapidly decreased from ~0.35 nm to ~0.19 nm within the first 20 ns and subsequently remained stable, with a fluctuation range of only 0.012 nm. In contrast, the TB1–TβRI distance exhibited large variations over time, without reaching a stable state, with an average distance of 0.790 nm and an average fluctuation of 0.1990 nm.

Taken together, these results indicate that in the TB2–TβRII–TβRI trimer, TB2 establishes a clear and stable interaction with TβRI, rapidly approaching and maintaining close contact. In contrast, TB1 shows little to no association with TβRI, consistent with the absence of stable interactions in the TB1–TβRII–TβRI system.

2.4. Free Energy Landscape (FEL)

Free energy landscape (FEL) analysis based on principal component analysis (PCA) is a powerful approach for characterizing the energetic stability and conformational states of protein–ligand complexes [49]. PCA was performed for the four systems—TB1–TβRII, TB2–TβRII, TB1–TβRII–TβRI, and TB2–TβRII–TβRI—using the first two principal components (PC1 and PC2) to calculate the Gibbs free energy of each conformation and construct three-dimensional free energy landscapes (3D FELs). PCA captures the dominant motions of proteins during simulations, while the FELs, projected onto the PC1–PC2 space, map Gibbs free energy variations through surface height and color gradients, reflecting transitions from high-energy to low-energy states.

As shown in Figure 5, deeper blue regions correspond to low-energy conformations, representing the most stable states. The TB1–TβRII complex exhibited multiple energy basins (Figure 5a), suggesting that this system explored several conformational substates during the simulation, indicative of conformational diversity. In contrast, TB1–TβRII–TβRI (Figure 5b), TB2–TβRII (Figure 5c), and TB2–TβRII–TβRI (Figure 5d) each displayed well-defined energy minima, consistent with strong thermodynamic stability. Collectively, the FEL analysis strongly supports our findings that TB2, unlike TB1, forms stable conformations upon binding to TβRII. Furthermore, the TB2–TβRII–TβRI complex maintains a stable conformation throughout the simulations, supporting the computational hypothesis of TB2-mediated receptor engagement.

2.5. TB2 Peptide Movement

To further elucidate the direction of movement of TB2 within the TB2–TβRII–TβRI system and clarify whether it facilitates the recruitment of TβRI, we extracted the most energetically favorable representative conformation of the TB2–TβRII–TβRI system from the energy basin identified in the FEL analysis. Structural superposition of the initial conformation with this representative conformation was then performed to analyze the direction of TB2 peptide movement. As shown in Figure 6, during the dynamic transition of TB2–TβRII–TβRI, the TB2 peptide underwent a significant positional shift. This movement was particularly pronounced in the C-terminal residues, with Pro9, Arg10, Gly11, and Pro12 all exhibiting displacements exceeding 10 Å. Most notably, the Cα atom of Pro12 displayed a displacement of 20.761 Å (Table 1). Upon binding to TβRII, the N-terminus of TB2 exhibited a certain degree of movement but remained associated with TβRII. In contrast, the C-terminus underwent a substantial shift toward TβRI, consistent with the results obtained from RMSF, Rg, and distance analyses. In addition to TB2, conformational adjustments were also observed in both TβRI and TβRII.

Taken together, these findings further substantiate our earlier conclusions: TB2 may behave similarly to TGF-β3 in terms of predicted receptor engagement. Specifically, after the N-terminus of TB2 binds to TβRII, its C-terminus progressively and rapidly approaches TβRI, ultimately anchoring to TβRI and facilitating the formation of a stable TB2–TβRII–TβRI trimeric complex. This structural recruitment of TβRI by TB2 is likely not only a structural phenomenon but may also be reflected in the interaction network and energy landscape of the system.

2.6. Dynamic Cross-Correlation Matrix (DCCM)

The Dynamic Cross-Correlation Matrix (DCCM) provides insights into the coordinated motions of residues within a protein complex [50]. To further elucidate residue-specific correlated movements and internal dynamics of the TB2–TβRII–TβRI complex, we calculated the covariance matrix based on the trajectories of Cα atoms and generated the corresponding DCCM.

As shown in Figure 7, DCCM values range from +1 to −1, where +1 (deep red) represents fully correlated motions and –1 (deep blue) represents fully anti-correlated motions. The intensity of the color indicates the degree of correlation, with darker colors reflecting stronger associations. Initially, DCCM analysis was performed for the entire TB2–TβRII–TβRI system, comprising 199 residues (Figure 7a).

To focus on the interfacial interactions, we defined a specific subgroup consisting of residues within 5 Å of inter-subunit atomic contacts (Cα atoms only): TβRII (residues 21–27, 51–53, 118–119), TβRI (residues 29–33, 54–67), and TB2 (residues 1–12). The DCCM of this subgroup is shown in Figure 7b. Notably, TβRI contributed more residues (19 residues) to the interfacial interactions than TβRII (12 residues). Residue 58 of TβRI displayed dynamic correlations with multiple residues of TβRII during the MD simulation. TB2 residues 1–3 exhibited dynamic correlations with TβRI residues 30–33 as well as several residues of TβRII; residues 4–7 correlated with TβRII residues 51–53; residues 6–9 correlated with TβRII residues 24–27; and residues 10–12 correlated with TβRI residues 58–63.

These results suggest that nearly all residues of TB2 are actively involved in interactions with both TβRII and TβRI, suggesting a potential role in bridging the two receptors and stabilizing the trimeric complex.

2.7. Key Hydrogen Bond Interactions

Hydrogen bond formation plays a critical role in stabilizing protein structures and mediating protein recognition processes. In the TB2–TβRII–TβRI complex, we observed variations in hydrogen bond formation among the three subunits (Table 2). Using an occupancy threshold of ≥40% during MD simulations, a total of 12 hydrogen bonds were identified between the subunits.

Notably, a highly stable hydrogen bond was formed between the hydrogen atom of Gly11 in TB2 and the oxygen atom of Phe60 in TβRI, exhibiting an occupancy of 77% and a distance of 0.293 ± 0.014 nm. This interaction is crucial for stabilizing the TB2–TβRI interface.

TB2 also established multiple hydrogen bonds with TβRII, among which several were particularly significant: (1) the nitrogen atom of His5 in TB2 formed a stable hydrogen bond with the hydrogen atom of Thr51 in TβRII (occupancy = 92.9%, distance = 0.307 ± 0.014 nm); (2) the hydrogen atom of Leu6 in TB2 formed a hydrogen bond with the oxygen atom of Thr51 in TβRII (occupancy = 88.4%, distance = 0.302 ± 0.017 nm); and (3) the oxygen atom of Glu119 in TβRII formed three hydrogen bonds with the hydrogen atoms of Arg10 in TB2, with occupancies of 49.1%, 45.5%, and 44.3%, respectively.

Similarly, multiple stable hydrogen bonds were observed between TβRI and TβRII, with four exhibiting occupancies greater than 90%: (1) the hydrogen atoms of Arg58 in TβRI formed two hydrogen bonds with the oxygen atom of Asp118 in TβRII, showing occupancies of 98.6% (distance = 0.277 ± 0.010 nm) and 97.0% (distance = 0.282 ± 0.011 nm); (2) the hydrogen atom of Arg58 in TβRI formed a hydrogen bond with the oxygen atom of Pro25 in TβRII (occupancy = 95.9%, distance = 0.285 ± 0.011 nm); and (3) the hydrogen atom of Cys76 in TβRI formed a hydrogen bond with the oxygen atom of Val22 in TβRII (occupancy = 91.2%, distance = 0.298 ± 0.015 nm).

These hydrogen bond interactions may collectively contribute to the stabilization of the TB2–TβRII–TβRI trimeric complex and likely underlie the high binding affinity and functional cooperation observed in this system.

2.8. Salt Bridge Interactions

In protein complexes, salt bridges are formed between oppositely charged amino acid residues when the distance between their side-chain centroids falls within a defined threshold and the geometric angle satisfies specific criteria. As an important stabilizing force for protein tertiary structures and molecular recognition, salt bridges are crucial for understanding biomolecular interactions. A salt bridge is considered stable when the distance between the charged side-chain centroids remains consistently within 3.5–5.0 Å. The dynamic behavior and stability of these salt bridges directly influence the overall stability of the complex.

Using the salt bridge analysis tool in VMD [51], we examined the TB2–TβRII–TβRI complex and identified several salt bridge interactions during the interaction between TβRI and the TB2–TβRII complex (Figure 8). Key salt bridges included TβRI_Asp57–TB2_Lys1, TβRII_Asp118–TβRI_Arg58, TβRII_Asp118–TβRI_Lys67, TβRII_Asp122–TβRI_Lys67, and TβRII_Glu119–TβRI_Lys67. We further analyzed the temporal variation of the Cα–Cα distances for these salt bridges to evaluate their stability over the simulation.

The salt bridge between TβRII_Asp118 and TβRI_Arg58 (yellow) maintained an extremely close distance throughout the simulation, with over 95% of the time below 4 Å, consistent with the characteristics of a strong salt bridge. Additionally, the TβRII_Glu119–TβRI_Lys67 interaction (blue) also remained relatively stable, primarily fluctuating within 2–8 Å, although transient dissociation events were observed. The dynamic variation in these salt bridge distances indicates the presence of persistent ionic interactions between TβRI and TβRII, which may significantly enhance inter-subunit binding through charge complementarity. The other three salt bridge pairs mostly remained above 4 Å during the MD simulations and did not exhibit pronounced salt bridge characteristics.

2.9. Relative Binding Free Energy Calculation

To investigate the relative binding affinities of TB peptides to TβRII and their potential to recruit TβRI, we performed MM-PBSA [52,53] calculations on the four systems. As shown in Table 3, TB2–TβRII exhibited a relative binding free energy of −35.91 ± 9.29 kcal/mol, while TB1–TβRII showed −20.35 ± 7.89 kcal/mol, indicating that both peptides can bind TβRII, with TB2 displaying higher affinity, consistent with previous experimental observations. For the trimeric complexes, the binding free energies between TβRI and the TB–TβRII complexes were −67.76 ± 7.70 kcal/mol for TB2–TβRII–TβRI and −38.12 ± 4.07 kcal/mol for TB1–TβRII–TβRI. Decomposition of the binding energy revealed that in TB2–TβRII–TβRI, TβRI interacts with TβRII (ΔG_binding = −47.19 ± 5.94 kcal/mol) and TB2 (ΔG_binding = −19.58 ± 4.00 kcal/mol), whereas in TB1–TβRII–TβRI, TβRI exhibits negligible binding to TB1 (ΔG_binding = −0.18 ± 0.11 kcal/mol) and primarily binds TβRII (ΔG_binding = −37.94 ± 4.06 kcal/mol).

Further analysis of energy components showed that van der Waals interactions (ΔG_vdw) and electrostatic interactions (ΔG_ele) are the major contributors to the stability and affinity of TB2–TβRII–TβRI, with ΔG_vdw = −76.12 ± 5.09 kcal/mol and ΔG_ele = −205.01 ± 34.61 kcal/mol, compared to ΔG_vdw = −37.05 ± 5.79 kcal/mol and ΔG_ele = −84.78 ± 21.33 kcal/mol in TB1–TβRII–TβRI. Although the polar solvation energy partially opposes binding, strong van der Waals and electrostatic interactions compensate for this unfavorable contribution, whereas the effect of nonpolar energy is minor.

Taken together, the computational results indicate that formation of the TB2–TβRII complex not only stabilizes the peptide–receptor interaction but also actively recruits and stably binds TβRI, with all three components contributing positively to complex formation. In contrast, TB1 binding to TβRII does not promote TβRI recruitment, highlighting the differential biological activities of the two peptides. To reduce computational costs, the MM-PBSA calculations in this study omitted entropic contributions. This represents a notable limitation, particularly for flexible peptide systems, as neglecting entropy may overestimate receptor–ligand binding affinities. Consequently, the resulting relative binding free energies are suitable only for comparative analysis of ligand binding within similar systems and should not be interpreted as absolute binding free energies.

2.10. Per-Residue Energy Decomposition

To elucidate the contributions of individual amino acid residues to the relative binding free energy in the TB2–TβRII–TβRI complex, we performed per-residue energy decomposition at the interaction interface between TβRI and the TB2–TβRII complex, as well as for all residues in TB2. The analyzed residues include TβRI: 29–33, 54–67; TβRII: 21–27, 51–53, 118–119; TB2: 1–12 (Figure 9).

In TβRI, residues Cys30, Phe31, Pro55, Arg58, Pro59, Phe60–Pro64, Ser66, and Lys67 contributed significantly to binding, with per-residue energies below −0.3 kcal/mol. In TβRII, residues Ala21, Val22, Phe24, Pro25, Leu27, Ile53, Asp118, and Glu119 showed notable contributions (energy < −0.3 kcal/mol). In TB2, the primary energy-contributing residues were located at the C-terminal region of the peptide, especially Leu6, His7, Val8, Gly11, and Pro12.

Among these, TβRII Asp118 (−5.833 kcal/mol), TβRI Arg58 (−7.045 kcal/mol) and Phe60 (−4.756 kcal/mol), as well as TB2 Gly11 (−1.69 kcal/mol) and Pro12 (−2.779 kcal/mol), play crucial roles in facilitating the recruitment of TβRI by the TB2–TβRII complex. Conversely, some residues at the interface, such as TβRI Asp57 and TB2 Pro9 and Arg10, exhibited positive energy contributions (>0 kcal/mol), suggesting a minor inhibitory effect on the binding of TB2–TβRII to TβRI, although this effect is limited.

By integrating hydrogen bond, salt bridge, and per-residue energy decomposition analyses, we identified key residues that facilitate the recruitment of TβRI by the TB2–TβRII complex (Table 4). In TβRI, Arg58 forms three stable hydrogen bonds with TβRII Asp118 and Pro25 (occupancies 98.6%, 97%, and 95.9%) and a salt bridge with Asp118, contributing strongly via electrostatic (−46.333 kcal/mol) and van der Waals (−4.67 kcal/mol) interactions. Phe60 forms a hydrogen bond with TB2 Gly11 (occupancy 77.8%) and contributes −5.804 kcal/mol via van der Waals interactions, while Ser66 hydrogen bonds with TβRII Asp118 (occupancy 46.7%) with an electrostatic contribution of −8.091 kcal/mol. In TβRII, Asp118 and Glu119 form multiple hydrogen bonds with TB2 and TβRI, contributing significantly through electrostatic interactions (−35.103 and −27.26 kcal/mol, respectively), and Val22 and Pro25 also form stabilizing hydrogen bonds with TβRI. On the TB2 peptide, Gly11 forms a critical hydrogen bond with Phe60, contributing −1.69 kcal/mol overall, and Pro12 contributes mainly via polar solvation and van der Waals interactions. Collectively, these residues constitute key interaction hotspots that stabilize the TB2–TβRII–TβRI complex and enable efficient recruitment of TβRI.

2.11. Experimental Validation of Peptide Activation of the Smad Pathway

Western blot (WB) analysis was performed to experimentally validate the peptide effects on the TGF-β/Smad signaling pathway and corroborate molecular simulation predictions. It is important to note that this method confirms the functional divergence between TB1 and TB2, rather than directly assessing receptor binding or recruitment events. TB1, TB2, and their corresponding random-sequence controls (RTB1: RPMIANQGDSHL; RTB2: VLGHKHPHRHLP) were applied to 3T3 cells for 48 h, with untreated cells serving as the control. As shown in Figure 10a–c, TB1 treatment significantly reduced Smad3 protein expression compared with both the control and RTB1 groups, while the level of phosphorylated Smad3 (p-Smad3) decreased, but the change did not reach statistical significance. These results indicate that TB1 suppresses Smad3 expression, reflecting an inhibitory effect on the TGF-β/Smad signaling pathway. This observation aligns with our computational analyses, which suggest that although TB1 binds the extracellular domain of TβRII, its binding site may hinder TβRI recruitment by TβRII.

In contrast, TB2 treatment markedly increased both Smad3 and p-Smad3 protein levels relative to the control and RTB2-treated groups (Figure 10d–f), confirming that TB2 effectively activates the TGF-β/Smad signaling pathway. The random-sequence control peptide RTB2 did not produce similar effects, emphasizing the sequence-specific activity of TB2. These experimental findings are in line with our molecular dynamics (MD), principal component analysis (PCA), dynamic cross-correlation matrix (DCCM), and MM-PBSA analyses, which suggest that TB2 not only forms a stable complex with TβRII but also facilitates TβRI recruitment, thereby enhancing Smad3 expression and triggering downstream signaling. This mechanistic insight likely underlies previous observations that TB2 promotes the secretion of type I and type III collagen.

3. Discussion

In this study, we investigated the biological functional differences between two TGF-β-mimetic peptides targeting TβRII, TB1 and TB2, through a multidimensional approach encompassing structural analysis, dynamic behavior, and experimental validation. The main findings were as follows: (1) TB2 forms a complex with TβRII via its N-terminal region, and its binding domains (Leu27–Asp32 and Ser49–Glu55) highly overlap with the TGF-β binding sites on TβRII. (2) MD simulations predict that the TB2–TβRII complex further recruits TβRI to form a stable heterotrimer, with a calculated relative binding free energy of −67.76 ± 7.70 kcal/mol. (3) Key interaction sites include the C-terminal residues of TB2 (Arg10, Gly11 and Pro12) and the N-terminal flexible region of TβRII (Ala19–Leu26). Notably, a stable hydrogen bond is formed between the H atom of TB2 Gly11 and the O atom of TβRI Phe60 (occupancy: 77.8%), and between TβRI Arg58 and TβRII Asp118 (occupancy: 97%). Additionally, a stable salt bridge is observed between TβRII Glu119 and TβRI Lys67. (4) Western blot experiments confirmed that TB2 activates the Smad3-dependent TGF-β/Smad signaling pathway.

Taken together, our computational modeling suggests that TB2 is predicted to mimic TGF-β by binding to TβRII and recruiting TβRI, thereby stabilizing the trimeric receptor complex, activating the Smad-dependent TGF-β/Smad signaling pathway, and promoting collagen secretion. In contrast, TB1 binds TβRII but appears unable to recruit TβRI, resulting in ineffective pathway activation. These findings not only validate the computational predictions but also suggest TB2 as a promising peptide template meriting further experimental evaluation for the rational design of therapeutics targeting the TGF-β/Smad signaling pathway and tissue repair. Of course, all results and conclusions remain subject to future experimental validation, such as through receptor-binding assays (e.g., SPR or ITC) and structural characterization techniques, to confirm the predicted molecular interactions.

In summary, by integrating AlphaFold 3 modeling, molecular dynamics (MD) simulations, relative binding free energy calculations, and experimental validation, this study proposes a possible molecular mechanism by which the peptide TB2 mimics TGF-β function and activates the TGF-β/Smad signaling pathway. The approaches and findings presented herein provide a valuable reference for the development of novel therapeutics for tissue repair.

4. Materials and Methods

4.1. Modeling of TB1–TβRII and TB2–TβRII Using AlphaFold3

AlphaFold3 [33,34,35], developed by DeepMind (Google DeepMind, 2024, London, UK), represents the latest generation of AI-based models for biomolecular structure prediction, leveraging deep learning algorithms to achieve high-accuracy structural forecasts. Initial structure predictions were performed using the AlphaFold3 online server (https://alphafoldserver.com/) with the default system parameters, which have been extensively optimized based on experimental data to balance computational efficiency and predictive accuracy in most scenarios.

For template searching, the default parameters were applied to ensure broad and effective homology template selection. The number of iterative predictions was set to 20, a parameter widely validated in previous studies to provide reliable structural predictions while maintaining computational efficiency. In addition, internal random seeds were used to enable reproducibility and account for stochastic variation in model predictions. The Pairformer module of AlphaFold3 was employed to process sequence alignments and generate the initial three-dimensional structures of the complexes. Following multi-step optimization of atomic coordinates and iterative refinement, final complex models were obtained and visualized using Discovery Studio v4.5 (DS45) [54]. A comprehensive evaluation strategy combining energy-based scoring and structural similarity comparison was applied.

The quality of the predicted structure was evaluated using the interface predicted Template Modeling score (ipTM) and the predicted Template Modeling score (pTM) to assess the accuracy of the interaction interfaces between different subunits (chains), as well as the Expected Position Error (EPE) of overall three-dimensional folding structure of the entire complex.

4.2. Construction of TB1–TβRII–TβRI and TB2–TβRII–TβRI Hypothetical Models

Based on the hypothesis proposed in this study, the TGF-β mimetic peptides are expected to simulate the biological function of TGF-β by binding to TβRII and further enhancing the affinity between TβRII and TβRI. Accordingly, we constructed the hypothetical models of TB1–TβRII–TβRI and TB2–TβRII–TβRI by incorporating the TβRI structural region into the previously obtained TB1–TβRII and TB2–TβRII complexes.

Model construction was performed using the Homology Modeling [55,56,57] module in Discovery Studio v4.5 (DS45) [54]. The crystal structures of the TGF-β–receptor complex (PDB ID: 2PJY, containing a TβRI–TβRII–TGF-β3 trimer) [16] and the individual receptors TβRI (PDB ID: 2L5S) [58] and TβRII (PDB ID: 1PLO) [59] were employed as templates. The Align Structures function in DS45 was used to spatially superimpose the AF3-generated TB1–TβRII and TB2–TβRII complexes onto the templates. This alignment is based on multiple sequence alignment (MSA), with a gap opening penalty set to 3.0.

Subsequently, the DS_Modeler module [56] was used to build the initial structures of the TB1–TβRII–TβRI and TB2–TβRII–TβRI models. Energy minimization was then performed in the GBSW implicit solvent model using the native DS protein force field via the DS_CHARMM [60] module, with a convergence criterion of 0.4184 kJ·mol^−1^·nm^−1^. These minimized structures were used as starting points for the subsequent molecular dynamics simulations.

4.3. Molecular Dynamics (MD) Simulations

Molecular dynamics (MD) simulations provide a powerful approach for modeling protein conformational changes under physiological conditions, thereby offering critical insights into the dynamic interactions between ligands and receptors [37,39]. To investigate the dynamic properties of TB1–TβRII, TB2–TβRII, TB1–TβRII–TβRI, and TB2–TβRII–TβRI, we performed 200 ns long-term molecular dynamics (MD) simulations using the GROMACS (Version 2024.2) package [61,62] with the RSFF2C forcefield [63,64] to examine atomic-level interactions and dynamic energetic behavior. Each initial model was solvated in a cubic box containing TIP3P [65] water molecules, and physiological ion conditions were established by adding Na^+^ and Cl^−^ ions to a concentration of 0.15 M, ensuring overall charge neutrality. Then, the protein molecules were subjected to 1000 steps of steepest descent (SD) energy minimization to remove steric clashes. Subsequently, a 2 ns restrained MD simulation of water molecules and ions was performed (1 ns each in the NVT and NPT ensembles) [66,67,68]. Finally, three independent replicates of unrestrained production MD simulations were conducted for 200 ns under the NPT ensemble at 300 K and 1 atmospheric pressure (atm) with periodic boundary conditions (PBC) [69]. Electrostatic interactions were calculated using the Particle-Mesh Ewald (PME) method, all the hydrogen bond lengths were constrained with the LINCS algorithm [70], the integration time step was set to 2 fs, and configurations were recorded every 200 ps.

The resulting MD trajectories were analyzed for Root Mean Square Deviation (RMSD) [44,45,46], Root Mean Square Fluctuation (RMSF) [47], Radius of Gyration (Rg) [48], and the number of hydrogen bonds to assess system stability. Additionally, inter-subunit distance analysis and principal component analysis (PCA)-based free energy landscape (FEL) calculations [49] were performed to further characterize the dynamic behavior of TB1 and TB2 during the MD simulations.

4.4. Dynamical Cross-Correlation Matrix (DCCM) Analysis

The Dynamical Cross-Correlation Matrix (DCCM) analysis [50] is an essential approach for investigating coordinated motions in functional biomacromolecules, such as proteins. DCCM captures the correlated motion between individual residues within a protein, allowing identification of residues that exhibit similar dynamic behavior during molecular motions. In receptor–ligand systems, conformational changes in the receptor and ligand are often interdependent, exhibiting coordinated motion. DCCM quantifies these correlations, providing insights into how ligand binding affects protein dynamics and highlighting key residues at the interaction interface that contribute to binding, thereby elucidating the mechanistic basis of protein–ligand interactions. For a system comprising N atoms, the DCCM is calculated as follows:

$$[eqn]$$

Among them, the correlation coefficient C_ij_ ranges from −1 to 1, and Δr_i_ = r_i_(t) − <r_i_>, where r_i_(t) is the position vector of the atom at time t, and <r_i_> is the average of the atom’s position vector over the entire trajectory time range; the angle brackets < $[eqn]$ > denote the average operation over the entire trajectory time. We used the covar tool in the GROMACS software (Version 2024.2) package to calculate the covariance of pairwise amino acid residues. In the representative conformations of the Free Energy Landscape (FEL), we selected the alpha carbon atoms of amino acid residues at the interaction interfaces between TB2 and TβRII, between TβRII and TβRI, and between the TB2 polypeptide and TβRI, to obtain the correlation of motion patterns between pairwise residues.

We utilized the GROMACS covar tool to compute the covariance matrix for paired amino acid residues. For the representative conformation extracted from the Free Energy Landscape (FEL) analysis, we selected Cα atoms at the following interaction interfaces: between TB2 and TβRII, between TβRII and TβRI, and between the TB2 peptide and TβRI. The motions of these paired residues were analyzed to quantify correlations in their dynamic behavior.

4.5. Relative Binding Free Energy Calculation by MM-PBSA

Relative binding free energy was calculated using the MM-PBSA (Molecular Mechanics Poisson–Boltzmann Surface Area) method, which combines molecular mechanics (MM) energies, polar solvation energy (PB), and nonpolar solvation energy (SA). PB and SA describe solvent effects, while MM energy accounts for van der Waals and electrostatic interactions between the receptor and ligand. A negative energy indicates spontaneous binding, with lower values corresponding to stronger binding affinity.

Fifty representative snapshots were extracted from the equilibrium phase of each MD trajectory for relative binding free energy calculation and per-residue energy decomposition using the gmx_MMPBSA (Version 1.63) software package [52,53]. In the MM-PBSA [40,41,42,71] approach, the enthalpic contribution of the system is evaluated using molecular mechanics (MM). The polar component of solvation free energy is calculated by solving the Poisson-Boltzmann (PB) equation, whereas the nonpolar contribution is estimated based on the solvent-accessible surface area (SA). The overall relative binding free energy can be expressed as follows:

$$[eqn]$$

where ΔG_bind_ represents the relative binding free energy, ΔE_MM_ denotes the molecular mechanical energy difference in vacuum, comprising electrostatic interactions (ΔE_ele_) and van der Waals interactions (ΔE_vdw_); ΔG_sol_ corresponds to the solvation free energy difference, which includes the polar solvation contribution (ΔG_PB_), calculated by solving the finite-difference Poisson-Boltzmann equation, and the nonpolar solvation contribution (ΔG_SA_), estimated from the solvent-accessible surface area (SASA). T is the absolute temperature, ΔS is the entropy change, and to save computational costs and improve computational efficiency, the entropic term (TΔS) was omitted in MM-PBSA calculations.

4.6. Peptide Synthesis and Western Blot (WB) Analysis

The peptides TB1 and TB2 and their corresponding random-sequence controls, RTB1 and RTB2 (sequences synthesized by Sangon Biotech, Shanghai, China; purity > 95%), were dissolved in sterile PBS to prepare 1 mM stock solutions. To avoid artificial charge interactions, the peptides were synthesized with N-terminal acetylation and C-terminal amidation. Mouse embryonic fibroblasts (NIH/3T3 cells) were used as an in vitro model system. Cells were seeded in 6-well plates at a density of 5 × 10^5^ cells per well and cultured for 24 h prior to treatment. The experimental design included the following groups: 1. Control group: Cells maintained in complete culture medium without peptide treatment. 2. Treatment groups: Cells treated with TB1 (0.2 mg/mL), TB2 (0.2 mg/mL), RTB1 (0.2 mg/mL), or RTB2 (0.2 mg/mL) for 48 h. After treatment, cells were harvested by trypsinization and subjected to Western blotting to evaluate protein expression. GAPDH (glyceraldehyde-3 phosphate dehydrogenase) was used as the internal loading control. Primary antibodies against Smad3, phospho-Smad3 (p-Smad3), and GAPDH were applied according to the manufacturer’s instructions. Protein bands were quantified using ImageJ (Version v1.54r), and the relative protein expression levels were normalized to GAPDH. Statistical analysis was performed using one-way ANOVA in GraphPad Prism 9.0 (GraphPad Software Inc., Boston, MA, USA), with p < 0.05 considered statistically significant.