Free energy perturbation (FEP)-guided scaffold hopping.

Scaffold hopping refers to computer-aided screening for active compounds with different structures against the same receptor to enrich privileged scaffolds, which is a topic of high interest in organic and medicinal chemistry. However, most approaches cannot efficiently predict the potency level of candidates after scaffold hopping. Herein, we identified potent PDE5 inhibitors with a novel scaffold via a free energy perturbation (FEP)-guided scaffold-hopping strategy, and FEP shows great advantages to precisely predict the theoretical binding potencies ΔG FEP between ligands and their target, which were more consistent with the experimental binding potencies ΔG EXP (the mean absolute deviations | Δ G FEP - Δ G EXP |  < 2 kcal/mol) than those ΔG MM-PBSA or ΔG MM-GBSA predicted by the MM-PBSA or MM-GBSA method. Lead L12 had an IC50 of 8.7 nmol/L and exhibited a different binding pattern in its crystal structure with PDE5 from the famous starting drug tadalafil. Our work provides the first report via the FEP-guided scaffold hopping strategy for potent inhibitor discovery with a novel scaffold, implying that it will have a variety of future applications in rational molecular design and drug discovery.

Introduction

One of the biggest challenges in drug discovery is to identify high-quality hit and lead compounds. Until now, natural products, high-throughput screening, and combinatorial chemistry have provided a wide range of molecular diversities with multiple privileged scaffolds for drug discovery (Fig. 1)1, 2, 3, 4, 5. However, how to identify privileged scaffolds efficaciously remains a great challenge for organic and medicinal chemists.

Figure 1

Privileged scaffolds in drug discovery.

Scaffold hopping, an effective approach to identify privileged scaffolds, usually refers to a molecule that gains potent bioactivity when its molecular scaffold is replaced with another scaffold, which has a different chemical structure but a similar shape and pharmacophore features, enabling it to interact in the same way with the target as the original molecule. It belongs to the core of drug design methods and is a topic of high interest in medicinal chemistry. Several methods have been developed for scaffold hopping, such as heterocycle replacements, ring opening or closure, computational methods (topological pharmacophore searching,, shape searching, machine learning methods,, chemical similarity searching, and structure-based similarity searching, Fig. 2)11, 12, 13, 14.

Figure 2

Various approaches to scaffold hopping.

Compared with the experimental methods, computational methods can significantly reduce the time and cost of scaffold hopping. However, most of the theoretical approaches do not consider the binding free energies of compounds as a parameter,. Recently, the core hopping FEP method performed relative binding free energy (RBFE) calculations well for limited and minor scaffold hopping. Most scaffold hopping procedures usually involve large topology changes of the entire ligand, and it is necessary to use absolute binding free energy (ABFE) calculations rather than RBFE calculations to predict the ligand binding free energies after scaffold hopping. For both ABFE and RBFE, one critical issue is high computational costs in practical applications. Moreover, the setup of systems for FEP calculations is complicated and requires experience.

To achieve the goal of discovering privileged scaffolds that possess potent affinities, we identified a novel scaffold for phosphodiesterase-5 (PDE5) inhibitors via the FEP-ABFE-guided scaffold hopping strategy (Fig. 3). Based on the similar pharmacophores of the famous PDE5 inhibitor tadalafil and the potent inhibitor LW1607, we performed scaffold hopping to achieve compound 5-(4-chlorobenzyl)-7-methoxy-2,3,4,6-tetrahydro-1H-azepino[5,4,3-cd]indol-1-one (L1) with a novel scaffold. Predicted by FEP calculations and confirmed by bioassay and X-ray crystallography, we performed FEP-guided structural optimizations based on L1. As a result, compound L12 exhibits a potent affinity of 8.3 nmol/L, high selectivity, and favorable pharmacodynamic effects, which provides the first report about the FEP-guided scaffold hopping strategy to discover potent inhibitors against PDE5.

Figure 3

Free energy perturbation (FEP)-guided scaffold hopping and hit-to-lead optimizations of 5-benzyl-2,3,4,6-tetrahydro-1H-azepino[5,4,3-cd]indol-1-ones as potent PDE5 inhibitors.

Results and discussion

FEP-guided scaffold hopping strategy to discover potent PDE5 inhibitors with a novel scaffold

As shown in Fig. 4, both tadalafil and LW1607 share an aromatic ring as an H-bond donor and a hydrophobic aromatic pharmacophore characteristic, which can form an H-bond with the conserved Gln817 and π–π stacking interactions with Phe820 and form π–π stacking interactions with Phe786, respectively. Thus, we designed 2-(4-chlorobenzyl)-9-methoxy-4,5-dihydro-1H-azepino[5,4,3-cd]indol-6(3H)-one (L1) to retain these pharmacophore characteristics. We carried out the FEP-ABFE protocol developed in our previous study to calculate the theoretical binding free energy ΔGFEP between PDE5 and tadalafil or LW1607.

Figure 4

Free-energy perturbation (FEP)-guided scaffold hopping. Alignment of tadalafil (blue sticks), LW1607 (yellow sticks), and L1 (orange sticks) with their common pharmacophore features. The key residues of the binding site are represented by gray sticks. The hydrophobic and hydrogen-donor pharmacophore features are represented by green and magenta grid balls, respectively. The H-bond interactions are represented by red dashes.

As shown in Table 1, the ΔGFEP for the PDE5-tadalafil complex (PDB code: 1XOZ) is very close to the experimental binding free energy ΔGEXP (ΔG EXP ≈ RT × ln IC50) with an absolute deviation (ΔGEXP‒ΔGFEP) of 0.07 kcal/mol, while that of the PDE5‒LW1607 complex is 2.3 kcal/mol higher than the experimental value. Thus, we used the crystal structure of the PDE5‒tadalafil complex as the receptor, and compound L1 was docked into the above receptor by using the Glide molecular docking program,. Then, this structure of the docked PDE5‒L1 complex was used as the initial structure in the subsequent FEP-ABFE calculations. The ΔGFEP value for the PDE5‒L1 complex is −10.98 kcal/mol, demonstrating that L1 exhibits considerable inhibition toward PDE5. Furthermore, organic synthesis of L1 (Supporting Information) followed by bioassay shows that L1's IC50 and ΔGEXP values with PDE5 are 55 nmol/L and −9.89 kcal/mol, respectively, verifying that the compound is a potent PDE5 inhibitor. As shown in the cocrystal structure (Fig. 5), L1 occupies the active pocket with a unique binding pattern. It is worth mentioning that the amide fragment of L1 formed another H-bond with residue Tyr612, which was not observed in the binding pattern of PDE5/tadalafil or PDE5/LW1607.

Table 1

The predicted binding free energies ΔGFEP, ΔGMM-GBSA and ΔGMM-PBSA (kcal/mol) by the FEP, MM-PBSA, and MM-GBSA calculations.

No.	Structure	IC50 (nmol/L)	ΔGEXP (kcal/mol)	ΔGFEP (kcal/mol)	ΔΔGFEP-EXP (kcal/mol)	ΔGMM-PBSA (kcal/mol)	ΔGMM-GBSA (kcal/mol)
Tadalafila	Image 1	1.8 ± 0.1	‒11.92 ± 0.01	‒11.99 ± 0.19	‒0.07 ± 0.02	‒16.16 ± 9.07	‒17.87 ± 8.99
LW1607	Image 2	5.6 ± 0.3	‒11.24 ± 0.03	‒13.54 ± 0.14	‒2.30 ± 0.17	‒15.18 ± 9.32	‒17.31 ± 9.09
L1	Image 3	55 ± 3	‒9.89 ± 0.03	‒10.98 ± 0.19	‒1.09 ± 0.22	‒12.39 ± 9.12	‒17.97 ± 8.84
L2	Image 4	55 ± 1	‒9.89 ± 0.01	‒11.10 ± 0.20	‒1.21 ± 0.21	‒12.47 ± 11.25	‒18.14 ± 10.91
L3	Image 5	346 ± 57	‒8.81 ± 0.10	‒8.42 ± 0.12	0.39 ± 0.22	‒9.80 ± 9.83	‒16.67 ± 9.34
L4	Image 6	150 ± 15	‒9.30 ± 0.06	‒9.05 ± 0.20	0.25 ± 0.26	‒14.17 ± 10.89	‒18.73 ± 10.74
L5	Image 7	30 ± 1	‒10.25 ± 0.02	–	–	14.12 ± 17.80	‒10.31 ± 11.69
L6	Image 8	10 ± 1	‒10.88 ± 0.07	‒9.10 ± 0.27	1.18 ± 0.34	‒8.51 ± 15.29	‒17.21 ± 13.79
L7	Image 9	8.8 ± 1.6	‒10.98 ± 0.11	‒12.28 ± 0.16	‒1.30 ± 0.27	‒9.27 ± 11.78	‒20.03 ± 9.92
L8	Image 10	39 ± 3	‒10.10 ± 0.04	‒8.84 ± 0.18	1.26 ± 0.22	‒10.54 ± 15.84	‒22.43 ± 14.43
L9	Image 11	14 ± 1	‒10.72 ± 0.02	‒12.00 ± 0.24	‒1.28 ± 0.26	‒12.12 ± 10.53	‒18.57 ± 9.12
L10	Image 12	32 ± 1	‒10.21 ± 0.02	–	–	‒9.47 ± 13.33	‒11.17 ± 13.07
L11	Image 13	12 ± 1	‒10.80 ± 0.05	‒11.69 ± 0.21	‒0.89 ± 0.26	‒14.10 ± 11.26	‒16.93 ± 10.89
L12	Image 14	8.3 ± 0.2	‒11.01 ± 0.01	‒12.73 ± 0.18	‒1.72 ± 0.19	‒16.50 ± 10.19	‒18.48 ± 10.00

Under identical assay condition, the reference drug tadalafil with an IC50 of 1.8 nmol/L.

Figure 5

Cocrystal structure of the PDE5‒L1 (PDB ID: 7FAQ) or PDE5‒L12 (PDB ID: 7FAR) complex. (A) Surface model for L1 (yellow sticks) binding; the red dotted lines represent hydrogen bonds. (B) Surface model for L12 (cyan sticks) binding. (C) Structural superposition between the PDE5‒L1 complex (yellow sticks) and the PDE5‒L12 complex (cyan sticks). The dotted lines represent H-bonds.

Subsequently, we carried out two rounds of efficient structural modifications on starting compound L1 (Supporting Information). In each round of modifications, we performed the FEP-ABFE protocol prior to the chemical synthesis and bioassay. As a result, the protocol significantly reduced our synthetic and bioassay efforts.

Chemistry

The targeted compounds were prepared by the synthetic routes reported in Scheme 1, Scheme 2, Scheme 3, Scheme 4, Scheme 5 as follows. Our initial efforts focused on design and syntheses of 5-benzyl-2,3,4,6-tetrahydro-1H-azepino[5,4,3-cd]indol-1-ones as novel PDE5 inhibitors (Scheme 1). The starting material 7-methoxyindole (1) was protected with benzenesulfonyl chloride to afford 2 to deactivation of the pyrrole ring, which allowed selective carboxylation in 4-position and followed by esterification to get 3. Intermediate 3 was treatment with NaOH to afford the acid 4 and esterification in methanol to obtain the ester 5. The ester 5 undergoes selective Vilsmeier formylation and Knovenagael reaction with nitromethane gives the nitrostyrene 7. The double bond is reduced with sodium borohydride by using a THF-methanol mixed solvent while the nitro group is reduced by hydrogenation in the presence of Pd catalyst to afford the spontaneous cyclization product 9. Then, the expected 5-benzyl-2,3,4,6-tetrahydro-1H-azepino[5,4,3-cd]indol-1-ones L1‒3 were synthesized by the alkylated reaction with substituted benzyl bromides to afford compounds 10–12 followed by rearrangement in polyphosphoric acid, respectively.

Scheme 1

Synthesis of compounds L1‒L3. Reagents and conditions: (a) benzenesulfonyl chloride, sodium hydride, DMF, 0 °C to rt, 1 h; (b) i. oxalyl dichloride, aluminum chloride, DCM, 0 °C to rt, 3 h; ii. MeOH, reflux, 3 h; (c) NaOH, EtOH, H2O, reflux, 3 h; (d) MeOH, thionyl chloride, reflux, 4 h; (e) DMF, phosphorus oxychloride, 0 °C to rt, 2 h, 60 °C, 1 h; (f) nitromethane, ammonium acetate, 105 °C, 15 min; (g) sodium borohydride, MeOH, THF, 0 °C, 1 h; (h) Pd/C, H2, MeOH, 40 °C, 24 h; (i) ArCH2Br, sodium hydride, DMF, 0 °C to rt, 1 h; (j) polyphosphoric acid, 90 °C, 3 h.

Scheme 2

Synthesis of compounds L4. Reagents and conditions: (a) vinylmagnesium bromide, THF, −78 °C, 20 min; (b) copper cyanide, 1-methylpyrrolidin-2-one, 200 °C, 4 h; (c) i. 5.0 mol/L NaOH; ii. MeOH, thionyl chloride, reflux, 4 h; (d) DMF, phosphorus oxychloride, 0 °C to rt, 2 h, 60 °C, 1 h; (e) nitromethane, ammonium acetate, 105 °C, 15 min; (f) sodium borohydride, MeOH, THF, 0 °C, 1 h; (g) Pd/C, H2, MeOH, 40 °C, 24 h; (h) 4-chlorobenzyl bromide, sodium hydride, DMF, 0 °C to rt, 1 h; (i) polyphosphoric acid, 90 °C, 3 h.

Scheme 3

Synthesis of compounds L5‒L10. Reagents and conditions: (a) RBr, sodium hydride, DMF, rt to 40 °C, overnight; (b) polyphosphoric acid, 60 or 90 °C, 3 h; (c) NaOH, MeOH, THF, H2O, 40 °C, 3 h; (d) methylamine hydrochloride, HATU, DIPEA, DMF, rt, 2 h.

Scheme 4

Synthesis of compound L11. Reagents and conditions: (a) boron tribromide, DCM, 0–50 °C, 4 h; (b) diethyl (bromodifluoromethyl)phosphonate, KOH, MeCN:H2O = 1:1, 0 °C to rt, 15 min.

Scheme 5

Synthesis of compound L12. Reagents and conditions: (a) i. MeOH, thionyl chloride, reflux, 4 h; ii. (chloromethylene)dibenzene, K2CO3, DMF, 60 °C, 16 h; (b) vinylmagnesium bromide, THF, −40 °C, 20 min; (c) Pd(OH)2/C, H2, EA, rt, 12 h; (d) 2-bromopropane, K2CO3, DMF, 70 °C, 16 h; (e) DMF, phosphorus oxychloride, 0 °C to rt, 2 h, 60 °C, 1 h; (f) nitromethane, ammonium acetate, 105 °C, 30 min; (g) sodium borohydride, MeOH, THF, 0 °C, 1 h; (h) i. Pd/C, H2, MeOH, 40 °C, overnight; ii. NaOH, C2H5OH, reflux, 5 h; (i) 1-(bromomethyl)-4-chlorobenzene, sodium hydride, DMF, 0 °C to rt, 1 h; (j) 1-bromo-3-methoxypropane, sodium hydride, DMF, 0–40 °C, 1 h; (k) polyphosphoric acid, 60 °C, 2 h; (l) diethyl(bromodifluoromethyl)phosphonate, KOH, MeCN:H2O = 1:1, 0 °C to rt, 15 min.

The synthetic approach to obtain 5-benzyl-2,3,4,6-tetrahydro-1H-azepino[5,4,3-cd]indol-1-one L4 was outlined in Scheme 2. The 4-bromo-1-fluoro-2-nitrobenzene (13) undergoes Bartoli reaction with vinylmagnesium bromide and nucleophilic substitution with copper cyanide to afford the 7-fluoro-1H-indole-4-carbonitrile (15),, which was treatment with NaOH to afford the acid and esterification in methanol to obtain the indole 16. Then, the indole 16 undergo the same procedures mentioned above to get product L4. The synthetic approach to obtain 5-benzyl-2,3,4,6-tetrahydro-1H-azepino[5,4,3-cd]indol-1-ones L5‒L10 were outlined in Scheme 3. Herein, the intermediate 10 was treatment with alkyl halides followed by the rearrangement reaction in polyphosphoric acid to afford the compounds L5‒L9 and 28, respectively. Then hydrolysis of the ester 28 to get the carboxylic acid L10.

Cleave the ether linkage of L1 to obtain 29 by treatment with boron tribromide in DCM at reflux temperature (Scheme 4). The phenol 29 was then treated with (bromodifluoromethyl)phosphonate and aqueous KOH in acetonitrile to get the compound L11,. Although the target compound L12 can be synthesized from compound L11, we developed a novel synthesis route and it works more efficiently (Scheme 5). The 4-hydroxy-3-nitrobenzoic acid (30) was esterification in methanol and followed by the alkylation reaction to get the ester 31. The ester 31 was then treated with vinylmagnesium bromide in THF at −40 °C to obtain the indole 32. It is worth noting that the benzhydryl group used as the protected group result in good yield. Finally, the indole 32 undergo the same procedures mentioned above to get the product L12.

Structure‒activity relationships (SARs)

Since the cocrystal structure of PDE5 with bound L1 shows that the inhibitor formed two H-bonds with the conserved Gln817, π–π stacking interaction with Phe820 and the other H-bond with Tyr612, we began the modification at the 5-position or 7-position substitution of the scaffold of L1. We designed L2 and L3 to investigate the structure–activity relationships of 5-position substitution by using 3-chloro-4-methoxybenzyl or 4-bromobenzyl substitution at the 5-position of the 2,3,4,6-tetrahydro-1H-azepino[5,4,3-cd]indol-1-one ring. Furthermore, L4 (fluoro group at the 7-position of the scaffold) and L11 (difluoromethoxy group at the 7-position) were designed to determine the favorable substituents at the 7-position. By using the crystal structure of the PDE5‒L1 complex as the receptor, we docked the above four compounds into the receptor, leading to several initial structures for FEP-ABFE calculations. Meanwhile, we carried out bioassays of these compounds. Both the FEP-ABFE calculations and bioassay results indicate that 4-chlorobenzyl or 4-bromobenzyl substitution at the 5-position and difluoromethoxy substitution at the 7-position on this scaffold are favorable for enhancing the binding potencies.

Subsequently, we designed L7 to study the SAR of 2-position substitution. Initially, based on predictive binding mode, we first calculated the ΔGFEP of L7 with PDE5 is −12.28 kcal/mol. And bioassays showed that its ΔGEXP is −10.98 kcal/mol, which is only 1.30 kcal/mol smaller than its ΔGFEP. Therefore, we designed six additional compounds (L5, L6, L8, L9, L10, and L12) to investigate the influence of the length and charging properties of the 2-position substituents on their activities. Since the ABFE-FEP protocol we used does not contain the correction for the charged ligands, we did not calculate the ΔGFEP values of L5 and L10. In addition, the inhibitory affinities of these six compounds against PDE5 were determined (Table 1). The results indicated that most of the compounds showed improved inhibitory activities compared with L1. Introduction of an appropriate length of oxygen-containing alkane chains, such as 3-methoxypropyl at the 2-position in L7, resulted in the most potent inhibition among the compounds with a methoxy substitution at the 7-position. Thus, L12 with a difluoromethoxy substitution at the 7-position and 3-methoxypropyl at the 2-position was selected as the candidate for subsequent study.

For comparison, we also used the commonly used MM-PBSA, and MM-GBSA methods embedded in our previously developed AutoMD protocol, to calculate the binding free energies ΔGMM-PBSA and ΔGMM-GBSA of the 14 compounds (details of the MM-PBSA and MM-GBSA calculation results can be seen in Supporting Information Table S1). The mean absolute deviations of the FEP-ABFE, MM-PBSA and MM-GBSA calculations are 1.13, 4.21, and 6.84 kcal/mol, respectively. It is worth mentioning that the former is much smaller than the latter two, which suggests that FEP-ABFE shows greater advantages in precisely predicting the receptor-ligand binding affinities between PDE5 and ligands than the commonly used MM-PBSA and MM-GBSA methods. Additionally, the predicted ΔGFEP values of the compounds corresponded closely with the experimental ΔGEXP values (Pearson's r values: 0.72 for the FEP method, Fig. 6). This statistically linear correlation between ΔGEXP and ΔGFEP (r = 0.72), demonstrated that this FEP-guided scaffold hopping method exhibits a remarkable statistical result and will have a larger variety of future applications in drug discovery than the MM-PBSA and MM-GBSA methods.

Figure 6

Linear regression and scatter plots of the FEP (left), MM-PBSA (middle), and MM-GBSA calculation (right) results. MAE refers to mean absolute error.

Cocrystal structure of PDE5 with bound L12

Given the potent inhibition of L12, the cocrystal structure of PDE5 with bound L12 was obtained at 2.4 Å resolution. The liganded PDE5 crystals had the trigonal space group P3121 with unit cell dimensions of a = b = 73.9 and c = 132.2 Å, and the structure was refined to R/Rfree of 0.21/0.25 (Supporting Information Table S2). As shown in Fig. 5, the 2Fo–Fc electron density unambiguously revealed the binding of L12 in PDE5 pocket. The scaffold of L12 formed an H-bond with the conserved Q817 and π–π stacking interaction with the hydrophobic clamp consisting of F820 and F786/V782. The chlorobenzene group fitted well in the Q2 subpocket and interacted with residues M816 and F787, which might be important for its selectivity profile. In addition, the amide fragment of L12 formed two H-bonds, including residues H613 and Y612, which were mediated by a water molecule. Further structural superposition of the crystal structures of PDE5‒L1 and PDE5‒L12 revealed that L12 has a slightly different binding pattern from L1 derived from the hydrophilic substitution at the 2-position, such as excluding the coordinated group and forming stronger H-bond interactions with the residues and waters (Fig. 5C). Therefore, L12 resulted in tighter binding with PDE5 than L1 and exhibited stronger inhibition.

Selectivity of compound L12 across PDE families

The selectivity of compound L12 across PDE families was also measured (Table 2). The inhibition toward PDE2A, PDE3A, PDE7A1, PDE8A1, and PDE9A2 is very weak (IC50 > 10,000 nmol/L). Its inhibitory values against PDE1B, PDE4D2, PDE6C, PDE10A, and PDE11A were 814-, 229-, 6.5-, 176-, and 100-fold higher than that against PDE5A1. As mentioned above, tadalafil has remarkable selectivity versus PDEs except PDE11 (5-fold) which results in the back and muscle pain. And the selectivity profile of sildenafil to the PDEs is similar to that of L12, and its IC50 potency against PDE6A was 5-fold higher than that against PDE5A1. Thus, compound L12 exhibited high selectivity over other PDEs except PDE6.

Table 2

Selectivity profile of compound L12 against PDE families.

PDE isozyme	IC50 (nmol/L)	Selectivity index
PDE5A1 (535–860)	8.3 ± 0.2	–
PDE1B (10–487)	6760 ± 50	814
PDE2A (580–919)	>10,000	>1205
PDE3A (679–1087)	>10,000	>1205
PDE4D2 (86–413)	1900 ± 50	229
PDE6C (1–858)	54.3 ± 3.1	6.5
PDE7A1 (130–482)	>10,000	>1205
PDE8A1 (480–820)	>10,000	>1205
PDE9A2 (181–506)	>10,000	>1205
PDE10A (449–770)	1460 ± 50	176
PDE11A (588–911)	831 ± 13	100

PK and acute toxicity studies

After intravenous injection (i.v.) of 2.5 mg/kg of L12 to male rats, a t1/2 of 1.99 ± 0.39 h and AUC0–24 h of 1125 ± 65 h ng/mL was obtained (Table 3). These results indicate that L12 could be used as a promising lead for further development. The acute toxicity of L12 was evaluated with twenty mice randomly divided into two groups. Single oral dose of 0 or 1.5 g/kg L12 was given on the first day. As a result, lead L12 was well tolerated up to a dose of 1.5 g/kg with no acute toxicity.

Table 3

The PK and safety profile of compound L12.

Content	Value
Pharmacokinetics parameters (i.v., 2.5 mg/kg)
t1/2 (h)	1.99 ± 0.39
AUC0–24 h (h·ng/mL)	1125 ± 65
CL (mL/min/kg)	35.2 ± 2.8
MRT0–24 h (h)	2.34 ± 0.32
Vss (mL/kg)	4918 ± 285
Acute toxicity	>1.5 g/kg

Pharmacodynamics effects of lead L12

To evaluate the pharmacodynamic effects of L12 against pulmonary arterial hypertension (PAH) in vivo, a monocrotaline-induced PAH rat model was adopted. As shown in Fig. 7, the mean pulmonary artery pressure (mPAP) was significantly increased in the model group compared to the control group. At the same time, the same trend was detected in the right ventricle hypertrophy index (RVHI%) and wall thickness percentage (WT%), both of which were significantly increased in the model group compared with those in the control group, which suggested that the model was built successfully after injecting monocrotaline (60 mg/kg) for 3 weeks.

Figure 7

Therapeutic effects of compound L12 and sildenafil citrate on a pulmonary arterial hypertension (PAH) rat model. (A) Mean pulmonary artery pressure (mPAP) of different groups. (B) Right ventricle hypertrophy index (RVHI%) of different groups. (C) Wall thickness percentage (WT%) of different groups. (D) Body weight of different groups. (E) Representative images for each group with hematoxylin and eosin staining. Data are presented as the means ± SEM (n = 6–10 animals per group). ##P < 0.01, ###P < 0.001, ####P < 0.0001, compared with the control group. ∗P < 0.05, ∗∗P < 0.01, and ∗∗∗P < 0.001 compared with the model group.

Rats were treated with lead L12 at a dose of 2.5 mg/kg (ip) and sildenafil citrate at a dose of 10.0 mg/kg (po) daily time for 3 weeks, respectively. Rats treated with compound L12 exhibited a remarkable decrease in mPAP in comparison with the model group, which indicated notable therapeutic effects against PAH. For the RVHI% and wall thickness percentage (WT%), similar trends were also achieved, and both of them decreased significantly compared to the model group. Meanwhile, a similar phenomenon was observed for the positive control sildenafil citrate (mPAP, RVHI%, and WT% decreased significantly), which provided comparable therapeutic effects to compound L12.

Conclusions

In this study, we report the discovery of novel PDE5 inhibitors via the FEP-guided scaffold hopping strategy. Lead L12 has a potent IC50 of 8.3 nmol/L with a totally different scaffold from the starting compounds and exhibits comparable therapeutic effects to sildenafil citrate on rats with PAH. Furthermore, L12 was revealed to show a different binding pattern from tadalafil and LW1607 in their cocrystal structures, which provides structural bioinformatics for the discovery of highly potent PDE5 inhibitors. In summary, our work provides the first study on the FEP-guided scaffold hopping strategy, which was successfully applied to the discovery of new PDE5 inhibitors with a novel scaffold different from tadalafil and LW1607. This strategy shows greater advantages in precisely predicting the receptor-ligand binding affinities ΔGFEP than the commonly used method MM-GBSA, in which the ΔGFEP values are consistent with their experimental binding potencies ΔGEXP (the mean absolute deviations  < 2 kcal/mol). Given these possibilities, the FEP-guided scaffold hopping method may have a variety of future applications in drug discovery and molecular design.

Experimental

Molecular docking

The crystal structure (PDB ID: 1XOZ) of PDE5-tadalafil was used for molecular docking of L1. For L2, L3, L4, L6, L7, L8, L9, and L11, they were docked into the protein of the crystal structure (PDB ID: 7FAQ) of PDE5‒L1. The molecular docking of all compounds is performed in the Glide molecular docking program with default parameters,. In order to verify that glide is suitable for predicting the conformation of small molecules in PDE5, we docked the original ligands in the used co-crystal structure into the corresponding proteins. We found that for the two complexes used (1XOZ and 7FAQ), the RMSD values between the predicted optimal conformations by Glide and the original ligand conformation were both below 0.5 Å (0.45 and 0.41 Å, respectively), which indicates that the Glide is suitable for the PDE5 system.

Free energy perturbation calculation

FEP protocol

To carry out the absolute binding free energy calculation based on free energy perturbation, we use the double doupling method described by the thermodynamic cycle proposed by Gilson and co-workers. For receptor-ligand complex annihilation process, we use 16 λ to sample the probability distribution of the potential energy differences between the adjacent λ leading to free energy difference . While we use 15 λ for the ligand annihilation process. For the restraint addition process, we use the strategy proposed by Karplus. The restrain parameters were determined by a “4 ns”-preliminary molecular dynamics simulation. And we use single-step perturbation to sample the probability distribution of the potential energy differences between the free state and restrained state. And we use the restraint energy distribution (RED) function to fit the probability distribution of the potential energy differences for this process. For the charge annihilation, 6 λ values were used (0.0, 0.2, 0.4, 0.6, 0.8, 1.0). Besides, 10 λ values were used for the van der Waals transformations (0.0, 0.10, 0.20, 0.30, 0.50, 0.70, 0.80, 0.90, 0.95, 1.00). We use Gaussian functions to fit the probability distribution for the charges and van der Waals transformations. The Bennet acceptance ratio (BAR), method was used to calculate the free energy differences between the adjacent window.

Molecular dynamic parameters

For each single λ window, there are four stages of molecular simulations. All ligands are parameterized by the general AMBER force field (GAFF). Use Gaussian 03 program to calculate the restrained electrostatic potential (RESP) charge of the ligands at the HF/6-31G∗ level. The protein parameters are provided by AMBER FF14SB force field. The TIP3P force field model is used for water molecules. If necessary, neutralize the system by adding a counter ion (Na+ or Cl‒). The first stage is minimization by using 5000 cycles of the steepest descent energy minimization. And then Langevin dynamics, was used for temperature coupling to heat the system to 298 K for 100 ps in the NVT ensemble. For third stage, as the position restraints applied, the system is simulated in an NPT ensemble for 500 ps with Parrinello-Rahman pressure coupling. Finally, the system is simulated freely in NPT ensemble for 4 ns. For the last three stages, the H-bonds are constrained by applying the LINCS constraint algorithm. The particle mesh Ewald (PME) algorithm is used to calculate the long-range electrostatic interactions.

MM-GBSA and MM-PBSA calculation

After 8 ns MD simulation with the same parameters above, the MM-GBSA calculation was performed by extracting 100 snapshots of the last 1 ns trajectories. The gas-phase energies were calculated by using the AMBER FF14SB force field as same as that was used in the MD simulations. In the MM-GBSA, Onufriev's GB model was used for GB calculation, and LCPO algorithm was used to calculate the nonpolar desolvation free energy with Eq. (1):

The dielectric constant value we used for the solute was set to 1 and dielectric constant value for the surrounding solvent was set to 80 in GB calculation. In the MM-PBSA calculations, we used the radii optimized by Tan and Luo. And we used molsurf to calculate the nonpolar desolvation free energy based on Eq. (2):

The dielectric constant value we used for the solute was set to 1 and dielectric constant value for the surrounding solvent was set to 80 in both PB and GB calculation. The entropy calculations are performed by normal mode analysis by using MM-PBSA.py.

All chemicals and reagents were bought from several commercial suppliers (Bide, Adamas, Energy, Sigma–Aldrich, and J&K) and tested directly without further purification. Silica gel plates with fluorescence F254 (0.1–0.2 mm) were performed for thin-layer chromatography (TLC) analysis, and chemical HG/T2354-92 silica gel (200–300 mesh) was carried out for column chromatography. Reactions requiring anhydrous conditions were used under argon or a calcium chloride tube. 1H NMR/13C NMR spectra were tested on a Bruker AVANCE III 400 instrument with tetramethylsilane as an internal standard. The following abbreviations are used: s (singlet), br (broad signal), d (doublet), dd (doublet of doublets), dt (doublet of triplets), t (triplet), td (triplet of doublets), q (quartet), and m (multiplet), and coupling constants are reported in Hz. High resolution mass spectra (HRMS) were recorded on a MAT-95 spectrometer. The purity of tested compounds was determined by reverse-phase high performance liquid chromatography (HPLC) analysis confirming to be more than 95%. HPLC instrument: SHIMADZU LC-20AT (detector: SPD-20A UV/Vis detector, UV detection at 254 nm; column: GL science InertSustain C18, 5.0 μm, 4.6 mm × 250 mm; Elution, MeOH in water (70%–90%, v/v); T = 25 °C; and flow rate = 1.0 mL/min).

Synthesis and characterization data of targeted compounds were given in Supporting Information.

Protein expression and purification

The expression and purification of PDE5A were carried out similarly to our previously published protocols. In brief, the catalytic domain coding (535–860) of PDE5A was cloned to vector pET-15b and then the cDNA was transferred to Escherichia coli strain BL21 (CodonPlus, Stratagene) for overexpression. When the cell carrying the plasmid was cultivated in LB medium at 37 °C until OD600 = 0.7, 0.1 mmol/L isopropyl b-d-thiogalactopyranoside (IPTG) was added to induce PDE5A expression for further 40 h growth at 15 °C. PDE5A protein was purified through Ni-NTA column (ϕ = 2.5 cm, 15 mL QIAGEN agarose beads), Q-column (ϕ 2.5 × 8.0 cm, GE Healthcare) and Superdex 200 column (ϕ 2.5 × 45 cm, GE Healthcare). A typical batch cell yielded over 10 mg PDE5A protein from 2L LB medium, with a purity > 95% shown by SDS-PAGE.

The catalytic domains of PDE1B (10–487), PDE2A (580–919), PDE3A (679–1087), PDE4D (86–413), PDE6A (484–817), PDE7A (130–482), PDE8A (480–820), PDE9A (181–506) and PDE10A (449–770) were purified by the similar protocols.

PDE enzymatic assays

Enzymatic activity assays of PDEs were performed similarly to our previously published protocol. The assays were measured by using corresponding 3H-cGMP or 3H-cAMP as substrate in an assay mixture buffer containing 50 mmol/L Tris-HCl (pH 8.0), 10 mmol/L MgCl2 or 4 mmol/L MnCl2, 1 mmol/L DTT. The reaction was carried out at room temperature for 15 min and terminated by adding 0.2 mol/L ZnSO4 and Ba(OH)2. The reaction product was concentrated to the precipitate while the unreacted substrate remained in the supernatant. Radioactivity in the supernatant was measured in 2.5 mL of Ultima Gold liquid scintillation cocktails by a liquid scintillation counter. The inhibitors were screened at a concentration of 100 nmol/L and the IC50 of inhibitors were measured at more than seven suitable concentrations for at least three times. The IC50 values were calculated by nonlinear regression. Sildenafil citrate serves as the reference compound with an IC50 of 5.1 nmol/L for PDE5.

Pharmacokinetics analysis in vivo

Pharmacokinetic properties of L12 were analyzed by Medicilon Company, Shanghai, China. Six male SD rats with body weight of 230–260 g were purchased from Shanghai SIPPR-BK LAB Animal Ltd., Shanghai, China, and used for the pharmacokinetic analysis of L12. It was dissolved/suspended in 5% DMSO, 10% Solutol, and 85% water for intravenous administration (iv) and for oral administration (po). A final dosage of 2.5 and 5 mg/kg rat of the formulated compounds was administrated for the iv and po purposes, respectively, and the blood samples were taken at various time points in 24 h. The concentration of the compounds in 23 was analyzed by LC–MS/MS (Shimadzu liquid chromatographic system and API4000 mass spectrometer, Applied Biosystems, Ontario, Canada). All animal care and experimental protocols were in accordance with “Guide for the Care and Use of Laboratory Animals” (National Institutes of Health Publication, revised 1996, No. 86-23, Bethesda, MD, USA) and were approved by the Institutional Ethical Committee for Animal Research of Sun Yat-sen University (Guangzhou, China).

Pharmacodynamics effects of compound L12 against PAH in rats

Forty-eight Wister rats (8 weeks, 180–220 g), purchased from the Laboratory Animal Center of Southern Medical University (Guangzhou, China), were used to evaluate the pharmacodynamics effects of L12 on PAH. The rats were randomly divided into four groups: control, model, compound L12 (5.0 mg/kg), and positive (sildenafil citrate, 10 mg/kg). Rats were maintained on a 12 h light/dark cycle (light from 7:00 to 19:00) at 24 ± 1 °C and 60%–70% relative humidity. Sterile food and water were given according to the institutional guidelines. Prior to each experiment, the rats were fasted overnight and allowed free access to water. All the rats were administrated with MCT 60 mg/kg except group control. Then, the rats were orally treated with drug vehicle (control and model groups), compound L12 (5.0 mg/kg) and sildenafil citrate (10 mg/kg) for 3 weeks, respectively. Compound L12 and sildenafil citrate were dissolved in 5% DMSO/10% Solutol/85% water solution and orally administrated 0.4 mL per 100 g weight. The method of right cardiac catheter was applied to measure the pulmonary artery pressure and the mean pulmonary artery pressure (mPAP) was used to conduct statistics. Subsequently, the rats were killed and the hearts were dissected into right ventricle (RV) and left ventricle and interventricular septum (LV + S); the 2 parts of the hearts were weighed with electronic scales, the value of RV/(LV + S) was used to conduct statistics.

Acute toxicity of compound L12

The acute toxicity was tested following the similar protocols described in our previous study. Thirty KM mice (22 days, 18–20 g), purchased from the Laboratory Animal Center of Sun Yat-sen University (Guangzhou, China), were used to evaluate the acute toxicity of L12. Mice were randomly divided into three groups, each of which was given in single oral dose of 0, 1000, or 1500 mg/kg L12 on the first day of the experiment. Mice were maintained on a 12 h light/dark cycle (light from 7:00 to 19:00) at room temperature and 60%–70% relative humidity. Sterile food and water were given according the institutional guidelines. Prior to each experiment, mice were fasted overnight and allowed free access to water. Compound L12 was dissolved in 5% DMSO/10% Solutol/85% water solution and orally administrated. Mice were observed for any abnormal behavior and mortality and weighed at the fourth hour of L12 administration and then every 24 h for 14 days. Animals were sacrificed on the 14th day, and tissue samples of heart, liver, and kidney were macroscopically examined for possible damages.

Statistical analysis

All experiments were performed in triplicate and repeated at least twice; representative data were selected for generating Figs. The statistical difference between treatments and controls was analyzed using Student's t-test. P ≤ 0.05 was considered statistically significant.

Accession codes

The atomic coordinates and structure factors have been deposited into the RCSB Protein Data Bank with accession number 7FAQ and 7FAR. Authors will release the atomic coordinates and experimental data upon article publication.