Structure-Guide Design and Optimization of Potential Druglikeness Inhibitors for TGFβRI with the Pyrrolopyrimidine Scaffold

Among all types of TGFβ signal blockers, small molecule kinase inhibitors (SMKIs) have attracted wide attention due to their economical production, obvious stability, and ease of oral administration. Nevertheless, SMKIs of TGFβRItypically have low druggability so there are none on the market. In this study, structure-based drug design (SBDD) was performed focusing on the pyrrolopyrimidin scaffold of BMS22 to find TGFβRIinhibitors with excellent medical potential. The binding mode, druggability, and target affinity were assessed by molecular docking, ADMET predictions, and molecular dynamics (MD) simulations for the designed TGFβRIinhibitors. Finally, the highly druggable compound W8 was discovered and then synthesized, which inhibited TGFβRIwith an IC50 value of about 10 μM. In addition, the binding free energies (ΔGbind) of W8 (−42.330 ± 3.341 kcal/mol) and BMS22 (−30.560 ± 6.076 kcal/mol) indicate that the high binding affinity is not necessarily accompanied by high inhibitory activity. Last but not least, the per-residue interaction analysis revealed that the contribution energy of ASP351 to binding was the most significant difference between BMS22 and W8, −2.195 kcal/mol and 1.707 kcal/mol, respectively. As a result, increasing the affinity between SMKIs and ASP351 of TGFβRImay effectively improve the inhibitory activity. The insights gained from this study could help with structure-guided optimization in searching for better SMKIs of TGFβRI.


Introduction
Transforming growth factor β (TGFβ) is one of the important members of the TGFβ superfamily [1]. Normally, TGFβ plays a crucial role in early embryonic development and tissue and organ growth, immunosuppression, tissue repair, and maintenance of tissue homeostasis [2][3][4][5][6][7]. However, TGFβs are frequently overexpressed in cancer, fibrosis, inflammation, or other disease states, and this excessive TGFβ production promotes the development of disease [8]. Interestingly, the TGFβ signal plays a dual role in tumorigenesis [9], exhibiting tumor suppressor properties during the early carcinogenesis and oncogenic properties during the late stage of malignancy [10,11]. In the latter, for example, the TGFβ signal stimulates tumor cell proliferation [12] and enhances stem cell properties while suppressing the sensitivity to anticancer drugs [13]. Additionally, the TGFβ signal interferes with the physiological functions of various immune cells, such as T cells and natural killer cells, and attenuates the inhibitory effect of immune cells on tumor cell migration [14][15][16]. Therefore, inhibiting the TGFβ signaling propagation is a key way for anticancer drug development. It is intensively researched in several clinical trials currently underway involving therapies targeting TGFβ signaling while other members of the TGFβ [14][15][16]. Therefore, inhibiting the TGFβ signaling propagation is a key way for anticancer drug development. It is intensively researched in several clinical trials currently underway involving therapies targeting TGFβ signaling while other members of the TGFβ superfamily are under-represented in current trials [17]. Hence, here we focus on the TGFβ signaling pathway.
There are several potentially active SMKIs of TGFβRI, such as Galunisertib [29,30], TEW-7197 [30,31], LY364947 [32,33], LY2109761 [34,35], and SB431542 [36,37] are undergoing clinical trials for treating cancers. Unfortunately, none of the inhibitors of TGFβRIhave been brought to market. It is necessary to further develop inhibitors for this target with potential druggability and high activity. Based on the findings of our previous studies [38][39][40], the pharmacophore model construction for the BMS22-TGFβRI complex (PDB ID: 6B8Y) [41], structure-activity relationship analysis, and scoring function evaluation, it is determined that: pyrrolotriazine TGFβRI inhibitors wedge with gatekeeper SER280 residue to weaken the stereo rejection in the pocket, and the inhibitors display a Y-shaped skeleton to coincide with the Y-type binding pocket ( Figure 2a); the PMF and Ludi 1 are the top two functions for identifying the best ligand pose. In addition, the structural similarity map (Figure 2b) between BMS22 and the five clinical TGFβRI inhibitors reveals they share the Y-shaped molecular shape and some partially identical molecular fragments, such as a pyridine ring [42], and suggests that BMS22 brings together the active fragments of the five molecules in clinical studies, having clinical application potential and commercial worth. Further detailed BMS22 involves the majority of the marketed drugs' fragments, such as pyrrolopyrimidine, pyridine and rifluoromethyl [42]. Based on the structure of this compound, it is promising to explore TGFβRI inhibitors with high druggability. Hence, BMS22 was chosen as an initial compound based on structure-based drug design to generate and optimize TGFβRI inhibitors. (PDB ID: 6B8Y) [41], structure-activity relationship analysis, and scoring function evaluation, it is determined that: pyrrolotriazine TGFβRⅠ inhibitors wedge with gatekeeper SER280 residue to weaken the stereo rejection in the pocket, and the inhibitors display a Y-shaped skeleton to coincide with the Y-type binding pocket ( Figure 2A); the PMF and Ludi 1 are the top two functions for identifying the best ligand pose. In addition, the structural similarity map ( Figure 2B) between BMS22 and the five clinical TGFβRⅠ inhibitors reveals they share the Y-shaped molecular shape and some partially identical molecular fragments, such as a pyridine ring [42], and suggests that BMS22 brings together the active fragments of the five molecules in clinical studies, having clinical application potential and commercial worth. Further detailed BMS22 involves the majority of the marketed drugs' fragments, such as pyrrolopyrimidine, pyridine and rifluoromethyl [42]. Based on the structure of this compound, it is promising to explore TGFβRⅠ inhibitors with high druggability. Hence, BMS22 was chosen as an initial compound based on structure-based drug design to generate and optimize TGFβRⅠ inhibitors.
(a) (b) Figure 2. A: A schematic view of the binding characteristics of BMS22 in the co-crystal structural domain depicting three colored regions: the hinge region (red), the selective hydrophobic pocket (yellow), and the solvent channel area (blue). As a typeⅠkinase inhibitor, BMS22 accommodates the DGF-in state of TGFβRⅠ. B: The similarity map of clinical trials testing TGFβRⅠinhibitors and BMS22 was predicted by Rdkit (http://rd.sourcceforge.net) based on MorganFingerprint, and the latter was the reference drug in this study (the greener the atom, the greater the similarity contribution). The source of the data literature is indicated in parentheses in Figure 2B.
In this study, based on the analyses of the active sites on TGFβRⅠ, we envisioned structural modifications of the pyrrolopyrimidine scaffold in BMS22 ( Figure 3B), which resulted in three series of modified molecules: H-series (hinge region fragments modified), S-series (selective region fragments modified), and W-series (solvent channel fragments modified). Then, for screening out the optimal and potential druggable TGFβRⅠinhibitors, the previous compounds were subjected to molecular docking, ADMET prediction, and dynamics simulation. Meanwhile, the best molecule was synthesized and its inhibitory activity was evaluated. The work chart is shown in Figure 3. Our work could help to make tangible progress in the direction of the discovery of better drug-like SMKIs of TGFβRⅠ.

Figure 2. (a)
A schematic view of the binding characteristics of BMS22 in the co-crystal structural domain depicting three colored regions: the hinge region (red), the selective hydrophobic pocket (yellow), and the solvent channel area (blue). As a typeIkinase inhibitor, BMS22 accommodates the DGF-in state of TGFβRI. (b) The similarity map of clinical trials testing TGFβRIinhibitors and BMS22 was predicted by Rdkit (http://rd.sourcceforge.net) based on MorganFingerprint, and the latter was the reference drug in this study (the greener the atom, the greater the similarity contribution). The source of the data literature is indicated in parentheses in Figure 2b.
In this study, based on the analyses of the active sites on TGFβRI, we envisioned structural modifications of the pyrrolopyrimidine scaffold in BMS22 ( Figure 3B), which resulted in three series of modified molecules: H-series (hinge region fragments modified), S-series (selective region fragments modified), and W-series (solvent channel fragments modified). Then, for screening out the optimal and potential druggable TGFβRIinhibitors, the previous compounds were subjected to molecular docking, ADMET prediction, and dynamics simulation. Meanwhile, the best molecule was synthesized and its inhibitory activity was evaluated. The work chart is shown in Figure 3. Our work could help to make tangible progress in the direction of the discovery of better drug-like SMKIs of TGFβRI.

The Binding Mode of the Co-Crystal of TGFβRⅠ
Detailed interactions between protein and inhibitor as well as key residues in the protein-inhibitor binding interface could provide substantial information for the discovery and design of TGFβRⅠinhibitors. LigPlot + is a program to show 2D representations of protein-ligand complexes from standard Protein Data Bank file input [43]. According to LigPlot + interaction analysis (Table S1), the drug-target binding pattern of 26 distinct TGFβRⅠinhibitors displays the common features: (1) the hinge region HIS283 acts as a hydrogen bond donor to interact with nitrogen or oxygen in heterocyclic segments of ligands, which may be a key residue for inhibition; (2) the gatekeeper residue SER280 and hydrophobic residues such as ILE, LEU, and ALA establish a significant number of hydrophobic contacts with the ligand, contributing to increased ligand-acceptor affinity; (3) ASP351 forms a hydrogen bond with the ligand in the solvent channel and regulates the kinase stability as one of the characteristic residues of the DLG-motif (in the TGFβRⅠstructure, PHE has been mutated to LEU in the DFG motif). Additionally, LYS232 is also a highfrequency amino acid residue that forms a hydrogen bond. In the active pocket of 6B8Y (Figure 2A), the binding mechanism of BMS22 is consistent with the aforementioned general characteristics of TGFβRⅠinhibitors. On the one hand, BMS22 forms hydrogen bonds with three essential residues, HIS283, LYS232, and ASP351; on the other hand, BMS22 also generates hydrophobic contacts with ILE211, VAL219, ALA230, LEU278, ALA350, LEU340 in the selective active pocket. Furthermore, there is almost no substantial interaction between BMS22 and the P-loop in the solvent channel. In conclusion, improving the

The Binding Mode of the Co-Crystal of TGFβRI
Detailed interactions between protein and inhibitor as well as key residues in the protein-inhibitor binding interface could provide substantial information for the discovery and design of TGFβRIinhibitors. LigPlot + is a program to show 2D representations of protein-ligand complexes from standard Protein Data Bank file input [43]. According to LigPlot + interaction analysis (Table S1), the drug-target binding pattern of 26 distinct TGFβRIinhibitors displays the common features: (1) the hinge region HIS283 acts as a hydrogen bond donor to interact with nitrogen or oxygen in heterocyclic segments of ligands, which may be a key residue for inhibition; (2) the gatekeeper residue SER280 and hydrophobic residues such as ILE, LEU, and ALA establish a significant number of hydrophobic contacts with the ligand, contributing to increased ligand-acceptor affinity; (3) ASP351 forms a hydrogen bond with the ligand in the solvent channel and regulates the kinase stability as one of the characteristic residues of the DLG-motif (in the TGFβRIstructure, PHE has been mutated to LEU in the DFG motif). Additionally, LYS232 is also a high-frequency amino acid residue that forms a hydrogen bond. In the active pocket of 6B8Y (Figure 2a), the binding mechanism of BMS22 is consistent with the aforementioned general characteristics of TGFβRIinhibitors. On the one hand, BMS22 forms hydrogen bonds with three essential residues, HIS283, LYS232, and ASP351; on the other hand, BMS22 also generates hydrophobic contacts with ILE211, VAL219, ALA230, LEU278, ALA350, LEU340 in the selective active pocket. Furthermore, there is almost no substantial interaction between BMS22 and the P-loop in the solvent channel. In conclusion, improving the interactions with HIS283, LYS232, ASP351, SER280, and the P-loop region is critical for the design of TGFβRIinhibitors.

Results of Molecular Modification
Structure-based drug design is crucial in the development of several FDA-approved small-molecule kinase inhibitors, such as the Bcr-abl kinase inhibitor dasatinib [44]. As has been reported, most kinase inhibitors can form hydrogen bonds with the hinge region residues via nitrogen-containing heterocycles, which is a conserved mode of hydrogen bonding interaction for the SMKIs [45]. In response, the FDA-approved six-membered ring system as well as fused five-and six-membered ring systems, such as aminopyridines, aminopyrimidines, pyrrolopyridines, pyr-rolopyrimidines, purines, and imidazo triazine [23], were used to replace the hinge binder, pyridyl, in BMS22. Additionally, the volume of the gatekeeper residue determines ligand access to the active hydrophobic pocket in the kinase domain, which is closely related to kinase selectivity [46]. Our earlier research demonstrated that pyrrolotriazine inhibitors form a Y-shape binding with SER280, the gatekeeper residue in the ATP cleft of TGFβRI [39]. As a result, this Y-shape binding may improve the selectivity of TGFβRI inhibitors. Even though the pyrrole ring in BMS22 extends to the solvent channel, there is a weak connection between the pyrrole ring and the P-loop ring, therefore opening the pyrrole ring to branch chains for improving the utilization of the P-loop region. Figure 3B depicts the three modifications made to BMS22 based on the preceding analysis: i. substituting ring A with some hinge binders; ii. introducing the polar groups such as -SF 5 , -NH 2 , and -SONH 2 to ring B to promote interaction with the SER280 as well as maintaining a Y-shape binding; iii. opening ring C and inserting flexible side chains such as -CH 2 CH 3 OH, -CH 2 CH 2 CH 2 CH 2 OH, -CH 2 C(CH 3 ) 2 OH, or other five-membered tiny ring fragments, etc.
In accordance with the substantial understanding of BMS22's binding mode and the confirmation of key residues in the active pocket, BMS22 was modified to produce 288 new molecules, including 124 H-series molecules, 19 S-series, and 145 W-series molecules ( Figure  S1), none of which were found to be entity molecules in the Scifinder database. All were continued to the next phase of molecular docking research.

Molecular Docking Analysis
Molecular docking is used to analyze the conformation and orientation of small molecules into the active site of a bio-macromolecular target [47], which was conducted by the CDOCKER module [48] in this study. And all modified molecules were docked into the active site on TGFβRI (PDB ID:6B8Y) to obtain the docking score, binding affinity and estimated K i value. Before this step, the docking pose of BMS22 and co-crystal conformation of 6B8Y were overlapped to validate the docking accuracy of the CDOCKER protocol. On the other hand, the ability of scoring functions in molecular docking to distinguish active and inactive molecules was evaluated in accordance with calculating the receiver operating characteristic curve (ROC). Then, the results show the root-mean-square deviation (RMSD) value between the docking pose and eutectic pose of BMS22 in the active pocket of 6B8Y was 0.2797 ( Figure 4), which means the CDOCKER protocol can be applied effectively to obtain a reliable docking pose of TGFβRI inhibitors. Whereas, the AUC of the CDOCKER ENERGY is 0.539, which is an incredible indicator to distinguish the active pose of these docking ligands. Herein, the credible scoring function must be filtered out for this target. As shown in Table S2, the ideal scoring function for the TGFβRIdocking system to accurately obtain the most active conformation of the ligand is PMF with an AUC of 0.921 (>0.9 indicates great accuracy).  The results of molecular docking revealed that 30 modified compounds had excellent PMF values that were higher than 112.00 comparable to the reference molecule BMS22 (PMF = 113.32), including 11 molecules from the H-series, 11 molecules from the S-series, and 8 molecules from the W-series (Table S3-1and Table S3-2). Except for H5 and S4, which formed only two hydrogen bonds with the active pocket, the other molecules formed three or four hydrogen bonds with key residues LYS232, HIS283, ASP351, SER280, or TYR249, as confirmed by manual review and had good hydrophobic interactions with non-polar amino acids such as ILE211, VAL219, LEU260, LEU278, and ALA350. This is consistent with most TGFβRⅠinhibitors in the crystal structure having the same mode of action. Table S3-1 summarized the PMF scores, binding affinity energies, estimated Ki values, and the key residues at the active site. The findings demonstrated that the compounds H1-H10 showed higher binding affinity and Ki values than BMS22. Meanwhile, H1 and H2, the top two PMF scoring molecules, may have potential TGFβRⅠinhibitory activity with substantially lower Ki values than the reference molecule (165.85 nM) of 4.87 nM and 8.50 nM, respectively. However, only S1 (−9.81 kcal/mol, 64.45 nM), S3 (−9.87 kcal/mol, 58.24 nM), S5 (−9.75 kcal/mol, 71.32 nM), and S6 (−9.68 kcal/mol, 80.26 nM) of the S-series showed a better binding affinity and Ki than BMS22. Furthermore, the binding modes revealed that the introduction of the -SONH2 group may effectively enhance the hydrogen bonding of S1, S3, and S6 to SER280, which would potentially improve its selectivity for TGFβRⅠ. For the W-series molecules, W1, W6, and W8 displayed strong columns of potential TGFβRⅠ inhibition with predicted Ki values of 2.55 nM, 1.1 nM, and 2.16 nM, respectively. Subsequently, based on the above analysis, the molecules H1-H10, S1, S3, S5, S6, W1, W6, and W8 were subjected to ADMET property prediction to assess the drug-like properties.

ADMET Prediction
In silico prediction of pharmacokinetic properties is an important component of The results of molecular docking revealed that 30 modified compounds had excellent PMF values that were higher than 112.00 comparable to the reference molecule BMS22 (PMF = 113.32), including 11 molecules from the H-series, 11 molecules from the S-series, and 8 molecules from the W-series (Table S3-1 and Table S3-2). Except for H5 and S4, which formed only two hydrogen bonds with the active pocket, the other molecules formed three or four hydrogen bonds with key residues LYS232, HIS283, ASP351, SER280, or TYR249, as confirmed by manual review and had good hydrophobic interactions with non-polar amino acids such as ILE211, VAL219, LEU260, LEU278, and ALA350. This is consistent with most TGFβRIinhibitors in the crystal structure having the same mode of action. Table S3-1 summarized the PMF scores, binding affinity energies, estimated K i values, and the key residues at the active site. The findings demonstrated that the compounds H1-H10 showed higher binding affinity and K i values than BMS22. Meanwhile, H1 and H2, the top two PMF scoring molecules, may have potential TGFβRIinhibitory activity with substantially lower K i values than the reference molecule (165.85 nM) of 4.87 nM and 8.50 nM, respectively. However, only S1 (−9.81 kcal/mol, 64.45 nM), S3 (−9.87 kcal/mol, 58.24 nM), S5 (−9.75 kcal/mol, 71.32 nM), and S6 (−9.68 kcal/mol, 80.26 nM) of the S-series showed a better binding affinity and Ki than BMS22. Furthermore, the binding modes revealed that the introduction of the -SONH 2 group may effectively enhance the hydrogen bonding of S1, S3, and S6 to SER280, which would potentially improve its selectivity for TGFβRI. For the W-series molecules, W1, W6, and W8 displayed strong columns of potential TGFβRI inhibition with predicted K i values of 2.55 nM, 1.1 nM, and 2.16 nM, respectively. Subsequently, based on the above analysis, the molecules H1-H10, S1, S3, S5, S6, W1, W6, and W8 were subjected to ADMET property prediction to assess the drug-like properties.

ADMET Prediction
In silico prediction of pharmacokinetic properties is an important component of pharmaceutical research and development and has been central to guiding hit-to-lead and lead-optimization efforts in the early phases of drug discovery [49]. In this research, the pharmacokinetic characteristics of modified compounds with high scoring, rational binding pose, and strongly predicted K i were calculated by the SwissADME [50] and ADMETlab 2.0 web servers [51]. Table 1 shows the results of ADMET prediction for 17 selected molecules. All modified compounds were accepted by the Lipinski Rules and only H1, W1, W6, and W8 showed a negative value to P-gp substrates, while BMS22 did not, which indicated that they might degrade the antitumor resistance. What's more, the introduction of cyclic fragments such as 1H-pyrrolo[3,2-c] pyridine (H3), pyrimidine-2 amine (H8), and other cyclic fragments could improve the predicted GI absorption, whereas the introduction of -SONH 2 did not significantly increase GI absorption. Additionally, most modified molecules did not have the potential BBB permeability for minimal CNS side effects and low carcinogenicity, while only the S3 and W8 had low potential cardiotoxicity. Overall, ADMET predictions revealed that W8 had outstanding pharmacokinetic properties and high druggability as the most potential TGFβRIinhibitor, with a synthetic accessibility score of 3.10 (from 1(easy) to 10 (difficult)). Other research showed that the side chain, 1-amino-2-methyl-2-propanol, of W8, contained two common medicinal fragments of isopropyl alcohol and propan-1-amine, and 3-fluoropyridin-4-yl and 6-(trifluoromethyl)pyridin-2-yl also appeared frequently in drug molecules, so it was inferred that W8 had good potential medicinal properties [42]. Based on the molecular docking analysis, we retained H1, S1, W1, and the best drug-like molecule W8 to further estimate the target affinities by MD simulations. In addition, the (6-trifluoromethyl) pyridine in W8 was substituted with the pyridine-3-sulfinamide of molecule S1 to fuse the compound S1W8. With four hydrogen bonding contacts made with HIS283, ASP351, ASN337, and ASN338, S1W8 outperformed BMS22 in terms of a PMF score of 124.65, even though it had a binding affinity of −8.74 kcal/mol and K i of 392.23 nM. The ADMET prediction showed that its LogP was lower than that of BMS22 and its other pharmacokinetic properties were equivalent to those of BMS22. In MD simulations, it was tested alongside H1, S1, W1, and W8. For the MD simulations, the input conformations of the aforementioned small molecules were listed in Figure 5.

MD Trajectories Analysis
In molecular biology and drug discovery, molecular dynamics (MD) simulations d matically capture the behavior of proteins and other biomolecules in full atomic detail a at a very fine temporal resolution [52]. In this study, the best docking conformation H1, S1, W1, W8, S1W8, and BMS22 were submitted to GROMACS 2020.3 software [53 further assess the binding stability and figure out the key interactions between poten SMKIs and TGFβRⅠduring the dynamical progress. The root mean square deviat (RMSD), the root mean square fluctuation (RMSF), and hydrogen bond occupancy w also calculated on the basis of 100 ns trajectory for the best docking conformation of e system.
Firstly, the convergences of MD trajectories were evaluated by RMSDs of Cα ato of the protein and the heavy atoms of ligands. As shown in Figure 6, the complex syste were stable as evidenced by the RMSD values of the TGFβRⅠ and the five modified co pounds, which fluctuated less below 0.30 nm. Furthermore, from the monitoring of RMSD value of the heavy atoms of the ligands, we can justify roughly that the ligand bind to the protein in a stable pose. Meanwhile, it has been seen that H1, S1, W1, W8, a S1W8 reached relative stabilities faster than the reference BMS22 system, with RMS convergent at about 0.10 ~ 0.20 nm. The trajectories were performed for conformatio analysis and energy calculation.

MD Trajectories Analysis
In molecular biology and drug discovery, molecular dynamics (MD) simulations dramatically capture the behavior of proteins and other biomolecules in full atomic detail and at a very fine temporal resolution [52]. In this study, the best docking conformations of H1, S1, W1, W8, S1W8, and BMS22 were submitted to GROMACS 2020.3 software [53] to further assess the binding stability and figure out the key interactions between potential SMKIs and TGFβRIduring the dynamical progress. The root mean square deviation (RMSD), the root mean square fluctuation (RMSF), and hydrogen bond occupancy were also calculated on the basis of 100 ns trajectory for the best docking conformation of each system.
Firstly, the convergences of MD trajectories were evaluated by RMSDs of C α atoms of the protein and the heavy atoms of ligands. As shown in Figure 6, the complex systems were stable as evidenced by the RMSD values of the TGFβRI and the five modified compounds, which fluctuated less below 0.30 nm. Furthermore, from the monitoring of the RMSD value of the heavy atoms of the ligands, we can justify roughly that the ligand can bind to the protein in a stable pose. Meanwhile, it has been seen that H1, S1, W1, W8, and S1W8 reached relative stabilities faster than the reference BMS22 system, with RMSDs convergent at about 0.10~0.20 nm. The trajectories were performed for conformational analysis and energy calculation. The differences in the flexibility of the residues were analyzed by main-chain root mean square fluctuation (RMSF) versus the residue number of Cα based on the last 100 ns trajectory for each system, the result of which is given in Figure 7. For all six systems, it can be seen that the RMSF values of most residues have similar trends, and the RMSF values of the residues in active sites, such as P-loop, hinge region, and DLG-in motif, are all significantly more rigid than other regions. Moreover, the inhibition of all modified molecules in the P-loop region is stronger than that of BMS22, which not only enhances the stability of the P-loop, but also intensifies its contact with the active residues of the Ploop.  The differences in the flexibility of the residues were analyzed by main-chain root mean square fluctuation (RMSF) versus the residue number of C α based on the last 100 ns trajectory for each system, the result of which is given in Figure 7. For all six systems, it can be seen that the RMSF values of most residues have similar trends, and the RMSF values of the residues in active sites, such as P-loop, hinge region, and DLG-in motif, are all significantly more rigid than other regions. Moreover, the inhibition of all modified molecules in the P-loop region is stronger than that of BMS22, which not only enhances the stability of the P-loop, but also intensifies its contact with the active residues of the P-loop. The differences in the flexibility of the residues were analyzed by main-chain root mean square fluctuation (RMSF) versus the residue number of Cα based on the last 100 ns trajectory for each system, the result of which is given in Figure 7. For all six systems, it can be seen that the RMSF values of most residues have similar trends, and the RMSF values of the residues in active sites, such as P-loop, hinge region, and DLG-in motif, are all significantly more rigid than other regions. Moreover, the inhibition of all modified molecules in the P-loop region is stronger than that of BMS22, which not only enhances the stability of the P-loop, but also intensifies its contact with the active residues of the Ploop. The hydrogen bond is a critical nonbonding interaction in biological molecules' spatial conformation and functions, thus the hydrogen bonds between the ligands and 6B8Y were monitored over the whole 100 ns MD simulation. Figure S3 shows the number of hydrogen bonds in each of the six TGFβRⅠ complex systems during the simulation. It has been seen that the number of hydrogen bonds formed by molecules H1 and S1 with the active residues remained stable, while that in the TGFβRⅠ-S1W8 system fluctuated greatly. The hydrogen bond is a critical nonbonding interaction in biological molecules' spatial conformation and functions, thus the hydrogen bonds between the ligands and 6B8Y were monitored over the whole 100 ns MD simulation. Figure S3 shows the number of hydrogen bonds in each of the six TGFβRI complex systems during the simulation. It has been seen that the number of hydrogen bonds formed by molecules H1 and S1 with the active residues remained stable, while that in the TGFβRI-S1W8 system fluctuated greatly. Although S1W8 can generate up to six hydrogen bonds, the hydrogen interactions were less stable. The numbers of hydrogen bonds between the ligand and protein in W1 and W8 complexes were essentially stable in two or three. Stronger hydrogen bonds are indicated by a higher hydrogen bond occupancy. As shown in Table S4 and Figure 8, the occupancy rate of hydrogen bond between H1 and HIS283 was 59.55%, and the other compounds established at least one strong hydrogen bond in the 6B8Y hinge region, especially with the residue HIS283, which had a greater than 80% occupancy. This indicates that all molecules can bind effectively to hinge region residues, retaining the general properties of a potent kinase inhibitor; additional hydrogen bonding contacts may contribute to affinity due to the low solvent accessibility of the internal hydrophobic pocket in the 6B8Y active pocket. The -SONH 2 group of S1 and S1W8 established the H-bond with the -OH of SER280, and the hydrogen occupancy rates were 53.47% and 10.07% respectively, both of which may increase the selectivity of TGFβRI. In this regard, enhancing the interactions with the gatekeeper residue SER280 in this area might improve the activity. W1 and W8 can even form more than two hydrogen bonds with ASP351. The top hydrogen bond occupancy between W1 and ASP351 was 94.67%, which was higher than that of BMS22 (58.87%), indicating that opening the pyrrole ring C and introducing the flexible side chains could enhance the hydrogen bonding interaction with ASP351. The largest hydrogen bond occupancy of the drug-like compound W8 with ASP351 was close to that of the reference molecule, which was 49.02%. indicated by a higher hydrogen bond occupancy. As shown in Table S4 and Figure 8, the occupancy rate of hydrogen bond between H1 and HIS283 was 59.55%, and the other compounds established at least one strong hydrogen bond in the 6B8Y hinge region, especially with the residue HIS283, which had a greater than 80% occupancy. This indicates that all molecules can bind effectively to hinge region residues, retaining the general properties of a potent kinase inhibitor; additional hydrogen bonding contacts may contribute to affinity due to the low solvent accessibility of the internal hydrophobic pocket in the 6B8Y active pocket. The -SONH2 group of S1 and S1W8 established the H-bond with the -OH of SER280, and the hydrogen occupancy rates were 53.47% and 10.07% respectively, both of which may increase the selectivity of TGFβRⅠ. In this regard, enhancing the interactions with the gatekeeper residue SER280 in this area might improve the activity. W1 and W8 can even form more than two hydrogen bonds with ASP351. The top hydrogen bond occupancy between W1 and ASP351 was 94.67%, which was higher than that of BMS22 (58.87%), indicating that opening the pyrrole ring C and introducing the flexible side chains could enhance the hydrogen bonding interaction with ASP351. The largest hydrogen bond occupancy of the drug-like compound W8 with ASP351 was close to that of the reference molecule, which was 49.02%. Figure 8. The 3D diagram of hydrogen bonds between TGFβRⅠ and each potential inhibitor: a, b, c, d, e, f are BMS22, H1, S1, W1, W8 and S1W8, respectively.

Binding Energy Analysis and Energy Decomposition
Although molecular docking is the gold standard in silico approach for identifying novel hit compounds that are active against the desired target receptor, there is still a need to improve the ability to distinguish true active ligands from inactive compounds [54]. In this context, using binding free energy evaluation approaches is a profitable strategy for revalidating docked ligand-protein complexes based on more reliable estimations of ligand-protein binding affinities than those obtained with simple scoring functions. To explore the differences in target affinity between each compound from the energy, the MM/PBSA calculation method was applied to calculate the binding free energies based on the last 20 ns equilibrium trajectories. The detailed information on binding free energies is listed in Table 2. Total binding free energies ΔGbind) were lower for W1, H1, W8, and Figure 8. The 3D diagram of hydrogen bonds between TGFβRI and each potential inhibitor: (a-f) are BMS22, H1, S1, W1, W8 and S1W8, respectively.

Binding Energy Analysis and Energy Decomposition
Although molecular docking is the gold standard in silico approach for identifying novel hit compounds that are active against the desired target receptor, there is still a need to improve the ability to distinguish true active ligands from inactive compounds [54]. In this context, using binding free energy evaluation approaches is a profitable strategy for revalidating docked ligand-protein complexes based on more reliable estimations of ligand-protein binding affinities than those obtained with simple scoring functions. To explore the differences in target affinity between each compound from the energy, the MM/PBSA calculation method was applied to calculate the binding free energies based on the last 20 ns equilibrium trajectories. The detailed information on binding free energies is listed in Table 2. Total binding free energies ∆G bind ) were lower for W1, H1, W8, and S1W8 than for BMS22 (−30.560 ± 6.076 kcal/mol), barring the S1 (−28.466 ± 4.026 kcal/mol). In descending order, values of W1, H1, W8, and S1W8 were −46.992 ± 3.902 kcal/mol, −44.797 ± 4.962 kcal/mol, −42.330 ± 3.341 kcal/mol, −36.399 ± 5.052 kcal/mol, respectively. Besides, the energy component analysis revealed that the van der Waals energies and electrostatic interaction energy play a dominant role in the contribution to the total binding affinity of TGFβRI potential agents. Moreover, ∆G vdw demonstrated a fairly high correlation toward ∆G bind (Pearson's r=0.614); although the value of ∆G SA itself was small, the influence of ∆G SA on the binding free energy could not be canceled (Pearson's r = 0.681). Therefore, it has been indicated that increasing the van der Waals contacts with TGFβRI and the non-polar solvent interaction can substantially boost the binding affinity to the receptor to a certain extent. Furthermore, it is critical to verify the interactions of modified molecules toward TGFβRI and to confirm the key residues in the binding interface. The results of binding free energy residue decomposition ( Figure 9) showed that W1 and W8 could enhance the interaction with the P-loop region compared with the reference BMS22, and LYS232 in the active site contributed the most to the total binding energy. Obviously, H1 greatly increases the binding free energy with LYS232. Despite their weak interactions with LYS232, S1 and S1W8 both demonstrated good contributions to the total binding energy with the gatekeeper residue SER280 that were absent from the original complex. This implied that the addition of -SONH 2 would improve the binding to the SER280, which might be advantageous for increasing the selectivity for TGFβRI. For the hinge-binding region residues, there was no significant difference in the contribution of binding energy in the hinge region. However, at DLG-motif, H1 improved the binding interaction of ASP351 significantly, and W1 bound to ASP351, similar to the original ligand BMS22, while the drug-like compound W8 displayed poor binding to ASP351. In the APE-motif region, there were no outstanding contributions or significant differences to the total binding for each molecule. Finally, because of its strong drug-like qualities, the compound W8 was further synthesized and kinase activity estimated.

Kinase Active Validation
ADP-Glo is a brand-new homogenous, bioluminescent assay for monitoring biological activities that produce ADP; as a result, it is ideally suited for determining enzyme activity when employing a variety of substrates [55]. It is a versatile test that can be utilized with ATPases, lipid kinases, sugar kinases, protein kinases, and many other types of kinases. In this study, the ADP-Glo technique was used to test the inhibition of the compound W8 on TGFβRⅠand its IC50 was evaluated to be 10.77 μM, indicating that it has a certain inhibitory impact on TGFβRⅠkinase and can be further examined as a lead chemical for TGFβRⅠinhibition.

Molecular Docking
The multiple poses generated by the professional search algorithm are ranked according to the scoring functions [58]. Therefore, the top docking pose for each ligand is scored by the optimal function. To confirm the accurate score function, the 50 active molecules were considered as the positive dataset and 300 inactive molecules of TGFβRⅠ were identified as a negative dataset, which were randomly selected from the DUD-e database [59] (http://dude.docking.org) using the Find Diverse Molecules protocol of DS2020. All above compounds were docked into the ATP binding cleft of 6B8Y. Next, the twelve scor-

Kinase Active Validation
ADP-Glo is a brand-new homogenous, bioluminescent assay for monitoring biological activities that produce ADP; as a result, it is ideally suited for determining enzyme activity when employing a variety of substrates [55]. It is a versatile test that can be utilized with ATPases, lipid kinases, sugar kinases, protein kinases, and many other types of kinases. In this study, the ADP-Glo technique was used to test the inhibition of the compound W8 on TGFβRIand its IC 50 was evaluated to be 10.77 µM, indicating that it has a certain inhibitory impact on TGFβRIkinase and can be further examined as a lead chemical for TGFβRIinhibition.

Molecular Docking
The multiple poses generated by the professional search algorithm are ranked according to the scoring functions [58]. Therefore, the top docking pose for each ligand is scored by the optimal function. To confirm the accurate score function, the 50 active molecules were considered as the positive dataset and 300 inactive molecules of TGFβRI were identified as a negative dataset, which were randomly selected from the DUD-e database [59] (http://dude.docking.org) using the Find Diverse Molecules protocol of DS2020. All above compounds were docked into the ATP binding cleft of 6B8Y. Next, the twelve scoring functions, including CDOCKER ENERGY, CDOCKER INTERACTION ENERGY, PLP1, PLP2, PMF, PMF04, Jain, LigScore1, LigScore2, Ludi Energy Estimate 1 (Ludi 1), Ludi Energy Estimate 2 (Ludi 2), Ludi Energy Estimate 3 (Ludi 3), were used to evaluate the ligand binding pose in the active pocket. Then, the score functions were estimated by the area under the curve (AUC) of the receiver operating characteristic (ROC) curve using the Calculation ROC Curve tool of DS 2020 program. The AUC > 0.9 means excellent accuracy, 0.7-0.9 indicates moderate precision, while an AUC of 0.5 represents a random event [60].
The CDOCKER, a semi-flex docking method based on the full CHARMm forcefield in DS2020, was used to identify the optimal binding conformations of ligands and the common interactions between ligands and TGFβRI. Before docking, the binding site was edited according to the original ligand, BMS22, in 6B8Y by the "From Current Selection" protocol. Next, the "Top Hits" was set to "10", "Pose Cluster Radius" was set to"0.5" and the other docking parameters were default. Finally, the optimal pose of each molecule with an excellent docking score was retained for the next analysis.
The estimated K i is a useful value as it is more related to usually measured experimental parameters, as compared to the affinity. The estimated K i and binding affinity of the designed inhibitors that were hit by CDOCKER were calculated using the Auto dock 4 protocol of AMDock 1.4 [61]. Initially, the "simple docking" option was selected to predict the binding mode of a single protein-ligand complex. Then, on the geometric center of the current ligand, a flexible box with a grid center (x = 5.3186, y = 8.8746, z = 5.0928) and optimum dimensions (x = 53, y = 53, y = 53) was defined using the "Center to hetero" option. The other parameters were set to their default values. Following docking, the conformation of each molecule that was consistent with the CDOCKER results was retained.

ADMET Prediction
The molecules with the high scoring, rational binding pose, and strong predicted K i were evaluated by the SwissADME online server [50] and ADMETlab 2.0 [51], the pathway to access absorption, distribution, metabolism, excretion, and toxicity properties, such as water solubility (Log P), gastrointestinal absorption (GI), blood-brain barrier permeability (BBB), carcinogenicity, and hERG inhibition. Meanwhile, the rule-based filters from Lipinski's rules were used to predict the drug-likeness, and the synthetic accessibility of molecules was measured on a scale of 1 to 10 [62]. On this basis, we preserved the compounds with the optimal drug-likeness for molecular dynamics study.

Molecular Dynamics Simulation
GROMACS 2020.3 software was utilized to evaluate the stability and the binding circumstance of a protein-ligand complex [53,[63][64][65]. The initial structures and topologies for protein 6B8Y and its ligands were considered according to the results of molecular docking. Firstly, the topologies and parameters of each molecule were derived from SwissParam, a fast force field generation tool (https://www.swissparam.ch/) [66]. The protein 6B8Y and water were then parameterized by CHARMm36 force field [67] and TIP3P [68], respectively. The complex system was performed in an NVT ensemble at 300 K for 100 ps and subsequently in an NPT constant simulation at 1 bar for 100 ps. After both equilibration processes, all-atom systems were conducted for 100 ns MD simulation. The particle mesh Ewald method was used to evaluate the electrostatic interactions with a cut-off radius of 1.2 nm, and van der Waals interactions were switched between 1.0 and 1.2 nm. The H-bonds were constrained with LINCS, in the order of 4. The relative stable trajectories were used for further analysis, such as the root mean square distance (RMSD) analysis and per-residue root mean square (RMSF) analysis, and so on.

Binding Free Energy Calculation and Decomposition
Molecular mechanics/Poisson-Boltzmann surface area (MM/PBSA) is one of the precise approaches to evaluate binding free energies from the MD trajectories [69]. This computational theory is based on the following equations: Where ∆G bind is the free energy of binding, ∆E MM is the difference between the internal energies of molecules in a vacuum, and ∆G sol is the difference between the free energies of solvation. ∆S represents the entropic change, while T represents the thermodynamic temperature. In addition, due to high computational costs and poor prediction accuracy, the conformational entropy change brought on by ligand binding (−T∆S ) is canceled out. ∆E MM can be subdivided into several terms, which is consistent with van der Waals energy (∆E vdw ), electrostatic interaction energy (∆E ele ) and (∆E ini ). However, the ∆E ini , such as a bond, angle, and torsion energy, is ignored usually. Moreover, ∆G sol is the total of two energies: the polar energy (∆G polar also called ∆G PB ) determined by the continuum solvent Poisson-Boltzmann model (PB) and the nonpolar energy (∆G apolar also called ∆G SA ) determined by the change of the solvent accessible surface area (SASA). Therefore, ∆G bind can be calculated using the following equation: In this study, the free energy of binding was evaluated by exacting the last 20 ns trajectory from the whole MD process, using the g_mmpbsa tool [70] according to the above formulas. Furthermore, the per residual average binding free energy contribution was also quantified by the ∆E vdw , ∆E ele , ∆G polar and ∆G apolar which were used to differentiate crucial residues in the protein-ligand interaction.

General Procedure of Synthetic Reaction General Information
All evaporations were carried out in a vacuum with a rotary evaporator. Analytical samples were dried under a vacuum (1-5 mmHg) at room temperature. Thin layer chromatography (TLC) was performed on silica gel plates, and spots were visualized by UV light (214 and 254 nm). Purification by column and flash chromatography was carried out using silica gel (200-300 mesh). Solvent systems were reported as mixtures by volume. All NMR spectra were recorded on a Bruker 400 (400 MHz) spectrometer. 1 H chemical shifts were reported as δ values in ppm with the deuterated solvent as the internal standard. Data are reported as follows: chemical shift, multiplicity (s = singlet, d = doublet, t = triplet, q = quartet, br = broad, m = multiplet), coupling constant (Hz), and integration. LCMS spectra were obtained on an Agilent 1200 series 6110 or 6120 mass spectrometer with electrospray ionization and excepted as otherwise indicated, the general LCMS condition was as follows: Waters X Bridge C18 column (50 mm × 4.6 mm ×3.5 um), Flow Rate: 2.0 mL/min, the column temperature: 40 • C.
The compound W8 was synthesized based on a previously described synthesis method [71] as depicted in Scheme 1. The NMR spectra and LCMS spectra of synthesized compounds were shown in Figure S4-S6 and Report 1-3 in the supporting file.
Where ΔGbind is the free energy of binding, ΔEMM is the difference between the internal energies of molecules in a vacuum, and ΔGsol is the difference between the free energies of solvation. ΔS represents the entropic change, while T represents the thermodynamic temperature. In addition, due to high computational costs and poor prediction accuracy, the conformational entropy change brought on by ligand binding (−TΔS ) is canceled out. ΔEMM can be subdivided into several terms, which is consistent with van der Waals energy (ΔEvdw), electrostatic interaction energy (ΔEele) and (ΔEini). However, the ΔEini , such as a bond, angle, and torsion energy, is ignored usually. Moreover, ΔGsol is the total of two energies: the polar energy (ΔGpolar also called ΔGPB) determined by the continuum solvent Poisson-Boltzmann model (PB) and the nonpolar energy (ΔGapolar also called ΔGSA) determined by the change of the solvent accessible surface area (SASA). Therefore, ΔGbind can be calculated using the following equation: In this study, the free energy of binding was evaluated by exacting the last 20 ns trajectory from the whole MD process, using the g_mmpbsa tool [70] according to the above formulas. Furthermore, the per residual average binding free energy contribution was also quantified by the ΔEvdw, ΔEele, ΔGpolar and ΔGapolar which were used to differentiate crucial residues in the protein-ligand interaction. All evaporations were carried out in a vacuum with a rotary evaporator. Analytical samples were dried under a vacuum (1-5 mmHg) at room temperature. Thin layer chromatography (TLC) was performed on silica gel plates, and spots were visualized by UV light (214 and 254 nm). Purification by column and flash chromatography was carried out using silica gel (200-300 mesh). Solvent systems were reported as mixtures by volume. All NMR spectra were recorded on a Bruker 400 (400 MHz) spectrometer. 1 H chemical shifts were reported as δ values in ppm with the deuterated solvent as the internal standard. Data are reported as follows: chemical shift, multiplicity (s = singlet, d = doublet, t = triplet, q = quartet, br = broad, m = multiplet), coupling constant (Hz), and integration. LCMS spectra were obtained on an Agilent 1200 series 6110 or 6120 mass spectrometer with electrospray ionization and excepted as otherwise indicated, the general LCMS condition was as follows: Waters X Bridge C18 column (50 mm × 4.6 mm ×3.5 um), Flow Rate: 2.0 mL/min, the column temperature: 40 °C.   The Synthesis of 6-(6-(Trifluoromethyl)Pyridin-2-yl)-1,3,5-Triazine-2,4(1H,3H)-Dione (2 in Scheme 1)

Kinase Active Validation
The ADP-Glo kinase assay was adopted to evaluate the inhibition effect of W8 toward TGFβRI in a 1× kinase buffer (50 mM Tris pH 7.5, 10 mM MgCl 2 , 200 nM NaCl, 0.1% BSA, 1mM DTT) at room temperature. The compound was dissolved in 100% DMSO to prepare the working solution at a concentration of 1.0 mM. The ADP-Glo kinase assays were conducted in wells of a 384-well white plate in the 20 µL final solution. Firstly, the 50 nL compound solution and 2.5 µL 2x kinase solution (TGFβRI) were added in the compound wells while the positive wells only added 2.5 µL 2x kinase solution (TGFβRI). Then, the negative control wells added the 2.5 µL 1x kinase buffer. Next the 2.5 µL 2x substrate (ATP) was added to each well and the reactions were incubated for 120 min. Furthermore, the 5 µL of ADP-Glo R1 agent was added for 180 min reaction time. Subsequently, the 10 µL of ADP-Glo R2 agent was added to incubate for 60 min. All experiments were carried out following the manufacturer's instructions (Promega). The data were collected in Envision and the IC 50 values were calculated in Graph Pad Prism 5.

Conclusions
The highly active SMKIs of TGFβRIhave been involved in various cancer treatments, including non-small cell lung carcinoma, hepatocellular carcinoma, glioma, pancreatic cancer, and scirrhous gastric cancer targeted therapies, and so on. However, no TGFβRI-listed medicine has been effectively created today. In this study, the structural domain of TGFβRI kinase was divided into the hinge region, hydrophobic selective region, and solvent channel region to guide the modification conducted on BMS22 based on the pyrrolo-pyrimidine backbone. The critical residues of TGFβRI were also identified. Then, the preferential binding conformations, drug-like properties, and the target affinity of the modified compounds were evaluated by computational simulations such as molecular docking, ADMET predictions, and molecular dynamics simulations. As a result, four potential TGFβRI inhibitors W1, H1, W8, and S1W8 having better binding free energies than that of BMS22 were obtained. Among these molecules, the ADMET prediction of W8 exhibited high water solubility, good GI absorption, and sufficient drug safety (non-carcinogenic, non-teratogenic, and non-cardiotoxic). Further experiments proved that the target's half maximal inhibitory concentration value against TGFβRI was about 10 µM, which was somewhat less inhibitory than that of the reference molecule BMS22 but has an exceptional drug likeness and may be applied clinically. It defies the ∆G bind prediction and illustrates that high binding energy does not necessarily result in a substantial inhibitory action. Meanwhile, the binding free energy on the per-residue level showed that the contribution of ASP351 in the DLG-motif of W8 changed from −2.194 kcal/mol to 1.706 kcal/mol compared with the BMS22, which may prove that ASP351 has a vital influence in the active site of TGFβRI.
Therefore, we can conclude that the high binding affinity may not necessarily lead to high inhibitory activity, and that interaction with certain amino acids, such as ASP351, is the crucial factor of activity. The subsequent modification of the TGFβRI inhibitor should focus on improving the interaction with ASP351, which is what we will study in the next stage. All the above may give structure-guided optimization in searching for the better SMKIs of TGFβRI.