In silico design and molecular modeling of finasteride analog compounds with potential inhibitory activity on enzyme 5α-reductase type 2: A target in treating benign prostatic hyperplasia

2.1. In silico construction of molecules

Seventeen molecules were designed based on the molecular structure of finasteride, which forms an adduct with NADPH. These molecules were developed through structural modifications using classical bioisosterism in the ‘N-tert-butyl 1-carboxamide’ region of the drug. The compounds were constructed and visualized using the Gauss View 5.0 program [20] and subsequently optimized geometrically with the semi-empirical PM6 method [21] implemented in the Gaussian 16 software [22]. After optimization, the molecular formats were converted using the Open Babel program [23] and further visualized with AutoDock Tools to prepare them for molecular docking studies [24].

2.2. System preparation

The 3D structure of the enzyme 5αR2 in complex with Dihydrofinasteride-nicotinamide adenine dinucleotide phosphate positive (NADP+) (NDX) was obtained from the Protein Data Bank with code 7BW1 [13]. The complex was visualized and analyzed with the Chimera 1.17.3 program [25].

In preparing the enzyme structure for molecular docking, the missing residues Ser39, Leu40, Lys41, Pro42, and Ala43 were initially added using the CHARMM-GUI web platform [26], as these were absent in the enzyme’s 3D structure. Subsequently, hydrogens were added using the ‘noBonderOrder’ method in AutoDock Tools 1.5.6. Non-polar hydrogens were removed, and Kollman charges were assigned [27]. Finally, the structure was converted from PDB format to PDBQT format.

The co-crystallized ligand (NDX) was extracted from the studied complex, its chemical structure was reviewed, and it underwent geometric optimization using the semi-empirical PM6 method, following the same protocol applied to the previously designed compounds. Like the other optimized compounds, it was then saved in .pdb format and subsequently converted to .pdbqt format using the AutoDock Tools 1.5.6 program [28].

2.3. Validation of the docking protocol

The NDX ligand was re-docked into the enzyme binding site using the AutoDock 4.2 program. Calculations were performed in triplicate, with a grid box defined at the binding site (X = -29.743; Y = 15.318; Z = 37.151) and dimensions of 60 × 60 × 60 Å, with a grid point spacing of 0.375 Å. 500 poses were calculated for each run. At the same time, all other parameters were kept at their default settings [29].

2.4. Molecular docking

The molecular docking process was conducted between the ligands from the created database and the 5αR2 enzyme using the Autodock 4.2 program, following the parameters established in the validation of the docking protocol. The results of the molecular docking calculations were visualized with the molecular operating environment (MOE) program and LigandScout [30,31].

2.5. MD simulations

The four enzyme-ligand complexes with the best binding energies from molecular docking calculations were subjected to MD simulations of 500 ns using the Gromacs 2020.2 program [32] to evaluate their stability during the simulation period and corroborate the energetic calculations. The force fields used for the enzyme and the ligand were the chemistry at harvard molecular mechanics (CHARMM) force field [33] and the CHARMM General Force Field (CGenFF) [34], respectively. Simulations were also performed to benchmark the enzyme with its co-crystallized ligand.

All complexes were immersed in a cubic periodic box, where each complex was solvated with TIP3P water under periodic boundary conditions [35]. The systems were neutralized, and the medium’s ionic strength (0.1 mol L⁻^-1) was adjusted by adding Na⁺ and Cl⁻ ions while keeping the number of particles constant. Energy minimization was performed until convergence was achieved. Next, an equilibration phase was conducted under constant pressure and temperature conditions [normal volume and temperature (NVT) and normal pressure and temperature (NTP) ensembles] at 300 K and 1.0 bar, with equilibration periods of 1.0 ns. The production runs lasted 500 ns, and trajectories were saved every 0.01 ns. The data obtained from the MD simulations were used to calculate the root mean square deviation (RMSD), root mean square fluctuation (RMSF), and hydrogen bond (H-bond).

The conformations obtained from the simulations were used to compute the total binding free energy. The molecular mechanics poisson-boltzmann surface area (MMPBSA).py script from the AMBER 21 suite was used to estimate the binding free energies of the protein-ligand complexes. For this purpose, the topology, coordinate, and production files generated in Gromacs were converted into their corresponding formats in Amber. The interaction energy and solvation-free energy for the complex, receptor, and ligand were calculated, and the results were averaged to estimate the binding free energy using the MM-GBSA approach [36].

2.6. ADMET properties

The ADME-T properties and other properties, such as physicochemical, drug-likeness, and characteristics such as lipophilicity and water solubility, were evaluated. Using the SwissADME server, the following parameters were calculated: molecular weight (Mw), number of hydrogen bond acceptors (HBA), number of hydrogen bond donors (HBD), number of rotatable bonds (Nrot), and topological polar surface area (TPSA) [37]. Lipophilicity, water solubility, and drug-likeness criteria were assessed. For the latter, compliance with the rules of Lipinski, Ghose, Veber, Egan, and Muegge were analyzed [37].

The evaluated parameters for the pharmacokinetic property “Absorption” included gastrointestinal absorption, Caco-2 permeability, and P-glycoprotein substrate status. Parameters related to “Distribution” included plasma protein binding, volume of distribution, and blood-brain barrier (BBB) permeability. For “Metabolism,” the involvement of Cytochrome P450 3A4 (CYP3A4) as a substrate and inhibitor was considered, as this enzyme is predominant in hepatic metabolism according to DrugBank [38]. For “Elimination,” the half-life (t₁/₂) and clearance rate were calculated. Lastly, for “Toxicity,” factors such as human hepatotoxicity, AMES mutagenicity, skin sensitization, lethal dose 50 (LD₅₀), and drug-induced liver injury were predicted. The ADMET-lab 2.0 server calculated these properties [39].

3.1. Construction of molecules

The NDX ligand was divided into two regions: Region A corresponds to the molecular structure of Finasteride in its Dihydrofinasteride (DHF) form, and Region B corresponds to NADP (see Figure 1).

Close

Figure 1.

Chemical structure of the NDX ligand. Region A corresponds to the molecular structure of Finasteride in its DHF form, and Region B corresponds to NADP+.

Export to PPT

Table 1 presents the 17 compounds derived from the NDX ligand, constructed using the ChEMBL server (ChEMBL Database, ebi.ac.uk). This platform enabled us to illustrate the structural modifications at positions R₁ and R₂ within the “N-tert-butyl 1-carboxamide” region (Region A) of Finasteride through classical bioisosterism [40]. The Nitrogen and Oxygen atoms are part of the original structure of Finasteride, as well as the amino group and the N-tert-butyl substituent in the R₂ region. However, modifications were made to R2 by altering the substituents attached to the amino group, introducing benzene, trifluoromethyl, dimethyl, propyl, ethyl, diethyl, and others. In some cases, Sulfur or Nitrogen atoms replaced the amino group itself.

Table 1. Compounds derived from the NDX ligand (Finasteride forming a molecular adduct with NADP+) through classical bioisosterism in Region A, specifically targeting the “N-tert-butyl 1-carboxamide” moiety of Finasteride.

Compounds	Substituents		Compounds	Substituents
	R1	R2		R1	R2
1	O		10	S
2	O		11	O
3	O		12	O
4	H2C		13	O
5	HN		14	O
6	HN		15	O
7	HN		16	O
8	HN		17	O
9	S

3.2. Molecular docking

3.2.1. Validation of the docking protocol

The structure of the NDX ligand was re-docked into the binding site of the 5αR2 enzyme using the AutoDock 4.2 program. The binding energy for the best pose was -21.6 ± 2.4 kcal/mol, and the RMSD value was also considered, which was 1.26 ± 0.23 Å. This pose adopts an orientation and conformation like that of the native ligand in the co-crystallized complex (see Figure 2) despite NADP having a high conformational flexibility. Figure 2 shows the superimposition of the co-crystallized ligand and the docked ligand, along with the enzyme’s amino acid residues interacting with the docked ligand in conjunction with NADP via a molecular adduct.

Close

Figure 2.

Binding at the active site of the 5αR2 enzyme. A-1: Alignment between the docked ligand (magenta) and the ligand present in the crystal structure (blue), average affinity, and RMSD of the top 4 poses of the docked ligand. A-2: Interactions of the native ligand with multiple amino acids present in the 5αR-2 enzyme.

Export to PPT

3.2.2. Molecular docking of the analogues

The Glu57 residue of the 5αR2 enzyme is considered the key amino acid residue involved in the inhibition of DHT formation by Finasteride [13]. This residue forms hydrogen bonds with the C-3 carbonyl and the N-4 nitrogen of the amino group in the quinoline of region A, facilitating the covalent bond formation between the quinoline and NADP.

Table 2 presents the binding energies of the four designed molecules with the best binding affinities to the 5αR2 enzyme, along with the NDX ligand, as well as the amino acids involved, and the types of interactions established. Binding affinity scores for the 13 missing compounds and 5αR2 are provided in Table S1.

Table S1

Download

Table 2. Binding energy values of the native ligand and the top four designed compounds.

Compound	Binding energy (kcal/mol)	Amino acid residues of 5αR2	Types of interactions
NDX ligand	-18.7	Arg114, Glu57, Glu197, Ser220, Arg179, Arg171, Asn160, Asn193, Trp201, His231, Lys35, Arg105, Tyr235, Tyr33, Asp164, Arg227, Tyr178, Phe118, Phe219, Phe223 y Ala24.	Hydrogen bonds, salt bridges, π-stacking, hydrophobic interactions
9	-17.7	Arg114, Glu57, Glu197, Asn160, Asn193, Ser220, Trp201, HSD231, Arg105, Tyr33, Asp164, Arg171, Arg227, Tyr178, Leu167, Gly104, Phe118.	Hydrogen bonds, salt bridges, π-stacking, hydrophobic interactions
2	-17.03	Arg114, Glu57, Glu197, Ser220, Asn160, Asn193, Trp201, Tyr33, Asp164, Leu167, Ser33, Gly34 Arg227, Arg171, Tyr178, Phe219 y Phe118.	Hydrogen bonds, salt bridges, π-stacking, hydrophobic interactions
14	-16.58	Arg114, Glu57, Glu197, Ser220, Arg179, Asn160, Asn193, Trp201, Lys35, Tyr235, Asp164, Leu167, HSD231, Tyr178, Arg171, Arg227, Phe219 y Phe118.	Hydrogen bonds, salt bridges, π-stacking, hydrophobic interactions
11	-16.32	Arg114, Glu57, Glu197, Ser220, Asn193, Tyr33, Asp164, Tyr96, Tyr91, Cys119, Arg171, Arg227, Phe118 y Ala24.	Hydrogen bonds, salt bridges, hydrophobic interactions.

COMP 9 exhibited the best binding energy, with a value of -17.7 kcal/mol, indicating strong docking affinity with 5αR2. This performance is comparable to the NDX ligand, which exhibited a binding energy of -18.7 kcal/mol. This molecule reproduces hydrogen bond interactions with the Glu57 amino acid in region A, which is critical for the action, as well as with the Arg114 residue interacting with sulfur in the same region. The hydrophobic interaction in this region is mediated by the Phe118 residue, which interacts with the C-9^a carbon of the methyl group in quinoline.

Additional hydrogen bonds in this region are observed with the following residues, as depicted in Figure 3: Trp201 interacts with the oxygen in the carbonyl group of nicotinamide, Tyr33 with the nitrogen in the amino group of adenines, Arg105 with the nitrogen in the amino substituent at the C-6 carbon, and Leu167 with the hydroxyl group of the phosphate group. Amino acids exhibiting hydrophobic interactions include the Glu197 residue interacting with the C-3 carbon of nicotinamide, the Arg105 residue with the five-membered ring of adenine, and Gly104 with the six-membered aromatic ring of adenine. Salt bridges in this molecule involve the Arg227 residue with the O-1 oxygen of the phosphate group attached to the ribose of nicotinamide and Arg171 with the oxygen of the phosphate group attached to the ribose of adenine.

Close

Figure 3.

3D and 2D views of complex 9. A-1: 3D view of the ligand at the binding site with the 5αR2 enzyme. A-2: Zoomed-in 3D view of the ligand at the binding site with the 5αR-2 enzyme. A-3: 2D view of the ligand with residue interactions.

Export to PPT

Finally, the π-stacking interaction is observed between the Tyr178 residue and the phosphate group attached to the ribose of adenine.

In COMP2, the interactions of residues Arg114 and Glu57 in region A were preserved. It was observed that the Ser31 residue interacts with the nitrogen of the amino group in the N-tert-butyl-1-carboxamide region via a hydrogen bond (see Figure 4). Additionally, the residues Phe219 and Phe118 maintain hydrophobic interactions with the C-10 carbon of the quinoline and the C-9^a carbon of the methyl group in quinoline, respectively.

Close

Figure 4.

3D and 2D views of complex 2. A-1: 3D view of the ligand at the binding site with the 5αR-2 enzyme. A-2: Zoomed-in 3D view of the ligand at the binding site with the 5αR2 enzyme. A-3: 2D view of the ligand with residue interactions.

Export to PPT

For region B, the hydrogen bond interactions of residues Trp201, Ser220, Asp164, Asn160, Asn193, and Tyr33, described in molecule 9, were conserved. At the same time, the Leu167 residue interacts via a hydrogen bond with the nitrogen of the amino substituent belonging to adenine, and the amino acids Gly34 and Arg105 also form hydrogen bonds with the O-3 oxygen of the phosphate group and the hydroxyl group oxygen of ribose, respectively.

The hydrogen bond interactions of COMP 14 in region A with the enzyme were also preserved (Glu57 and Arg114), as well as the hydrophobic interactions of residues Phe118 and Phe219, which exhibit Van der Waals interactions with the C-9^a carbon of the methyl group in quinoline and the C-10 carbon of quinoline, respectively.

As in the two previous compounds, in region B of this molecule, hydrogen bond interactions were conserved in the residues Ser220, Trp201, Asn160, Asn193, and Asp164 (see Figure 5). Additionally, hydrogen bond interactions were observed in the residues Glu197, Leu167, Lys35, Tyr235, Arg179, and HSD231, which formed hydrogen bonds as follows: the Glu197 residue interacts with the oxygen of the hydroxyl group of the ribose attached to nicotinamide; the amino acids Leu167, Lys35, Arg179, and Tyr235 form hydrogen bonds with the oxygen of the hydroxyl group of the phosphate group attached to the ribose of adenine; and the HSD231 residue interacts with the oxygen of the phosphate group of NADP.

Close

Figure 5.

3D and 2D views of complex 14. A-1: 3D view of the ligand at the binding site with the 5αR2 enzyme. A-2: Zoomed-in 3D view of the ligand at the binding site with the 5αR-2 enzyme. A-3: 2D view of the ligand with residue interactions, using LigandScout 4.0 software. A-4: 2D view of the ligand with residue interactions, using MOE software.

Export to PPT

In the 5αR2–COMP11 complex, hydrogen bond interactions are reproduced in region A with the residues Arg114 and Glu57, as in the previous compounds. Another hydrogen bond interaction was observed between the Tyr91 residue and the oxygen of the oxo group in quinoline. Hydrophobic interactions in this region occurred between the amino acids Ala24, Phe223, and Phe118, which interacted with the C-2 carbon of the ethyl group, the C-10 carbon of quinoline, and the C-9^a carbon of the methyl group in quinoline, respectively.

On the other hand, interactions in region B involving the amino acids Ser220, Asn193, Lys35, Trp201, and Asp164 were preserved. The Glu197 residue, in addition to forming a hydrogen bond with the nitrogen of the amide group of nicotinamide, presents a Van der Waals interaction with the C-1 carbon of the ribose attached to it. Other hydrogen bonds formed in this region occurred with the residues Tyr91 and Cys119, which form a hydrogen bond with the 7-oxo carbonyl group of the compound; the Tyr98 residue, which interacts in the same way with the oxygen of the hydroxyl group belonging to the phosphate group; and the Tyr33 residue, which forms a hydrogen bond with the nitrogen of the amino substituent group in adenine (see Figure 6).

Close

Figure 6.

3D and 2D views of complex 11. A-1: 3D view of the ligand at the binding site with the 5αR2 enzyme. A-2: Zoomed-in 3D view of the ligand at the binding site with the 5αR-2 enzyme. A-3: 2D view of the ligand with residue interactions.

Export to PPT

3.3. MD

Since molecular docking provides only a static analysis of the systems under study [41], MD simulations were conducted to gain deeper insight into the dynamic behavior of the interactions between the ligands and the enzyme. For this purpose, the four most stable complexes from the molecular docking results, listed in Table 2, were selected and subjected to 100 ns and 500 ns MD simulations per complex to assess their stability and potential conformational changes [42], as well as to calculate the binding free energy for each. The 5αR2-NDX complex was also simulated under the same conditions to compare the conformational and energetic changes, using RMSD, RMSF, Hbond, and MM-GBSA calculations.

3.3.1. RMSD analysis

These values suggest that complexes 2, 9, 11, and 14 do not undergo abrupt conformational changes throughout the simulation and remain in equilibrium for the majority of the time. The variation of RMSD of the backbone over time is shown in Figure 7, where the fluctuations of each complex are compared to the 5αR2 enzyme alone (Black color in all figures) and to the enzyme with the NDX ligand (Red color in all figures).

Close

Figure 7.

RMSD of the backbone of the 5αR2-COMP complexes during the simulation as a function of time. The X-axis represents the time scale in nanoseconds (ns), and the Y-axis represents the RMSD in nanometers (nm).

Export to PPT

It can be observed that the 5αR2-NDX complex does not exhibit significant conformational changes during the simulation, showing minimal variations over time. It presents an average RMSD value of 0.28 ± 0.021 nm. The 5αR2 enzyme alone, on the other hand, displays fluctuations between 0.15 and 0.25 nm, during the first 325 ns, then it presents a fluctuation that goes from 330 to 350 ns, with RMSD values ​​of 0.3 to 0.37nm. Its maximum RMSD value is 0.45 nm between 461 ns and 463 ns of the simulation. The 5αR2-COMP2 complex presents a maximum RMSD value of 0.52 nm around 435 ns of simulation, with an average RMSD value of 0.35 ± 0.06 nm. The 5αR2-COMP9 does exhibit minimal conformational changes during the simulation (see Figure 7), with an average RMSD value of 0.36 ± 0.05 nm. The 5αR2-COMP11 complex presents an average RMSD value of 0.34 ± 0.06 nm and Its maximum RMSD value is 0.48 nm between 88 ns and 90 ns of the simulation. After 300 ns, the RMS value decreases with values ​​ranging between 0.3 and 0.32 ns.

The 5αR2-COMP14 complex presents a uniform behavior throughout the simulation with an average RMSD value of 0.33 ± 0.03 nm. These values suggest that complexes 2, 9, 11, and 14 do not undergo abrupt conformational changes throughout the simulation and remain in equilibrium for most of the time. The variation of RMSD of the backbone over time has been shown in Figure 7, where the fluctuations of each complex are compared to the 5αR2 enzyme alone (Black color in all figures) and to the enzyme with the NDX ligand (Red color in all figures).

On the other hand, the RMSD values of the ligands has been described below. The NDX ligand, COMP2, COMP9, COMP11, and COMP14 exhibited values of 0.36 ± 0.028 nm, 0.29 ± 0.026 nm, 0.28 ± 0.025 nm, 0.36 ± 0.019 nm, and 0.33 ± 0.02 nm, respectively. In this case, the RMSD comparison was performed between the co-crystallized ligand (colored red in all figures), which forms an adduct with NADP, and the designed compounds. As shown in Figure S1, all compounds are highly stable, with no significant variability during the simulation, maintaining their RMSD values close to 0.28nm y 0.36 nm, which is very similar to that of the co-crystallized ligand.

Figure S1

Download

3.3.2. RMSF analysis

The RMSF was calculated to analyze the flexibility of amino acid residues over time [43], identifying the regions with the highest fluctuations relative to the native ligand [44,45]. The average RMSF values for 5αR2 alone, 5αR2-NDX, 5αR2-COMP2, 5αR2-COMP9, 5αR2-COMP11, and 5αR2-COMP14 were 0.13 ± 0.090 nm, 0.10 ± 0.06 nm, 0.13 ± 0.10 nm, 0.13 ± 0.1 nm, 0.15 ± 0.08 nm, and 0.13 ± 0.080 nm, respectively. In this case, the comparison was made between the enzyme alone, the enzyme bound to the co-crystallized ligand, and each of the four compounds, as illustrated in Figure 8. The enzyme alone exhibits fluctuations, with the most significant occurring from residues Ala24 to Ala44, where the RMSF values range from 0.2 nm to 0.4 nm, and from residues Leu170, Arg171, Lys172, Pro173, Gly174, and Glu175, with RMSF values ranging from 0.2 nm to 0.53 nm.

Close

Figure 8.

RMSF of the 5αR2-COMP complexes over 500 ns of simulation. The X-axis represents the residue number, and the Y-axis represents the RMSF of the backbone atoms in nanometers (nm).

Export to PPT

3.3.3. Hydrogen bonds

The comparison was made between the enzyme bound to the native ligand via a molecular adduct with NADP and the constructed complexes. As shown in Figure 9, the 5αR2-NDX complex forms approximately 11 hydrogen bonds up to 200 nanoseconds. Subsequently, a slight decrease to 7–8 hydrogen bonds is observed between 220 ns and 280 ns. After 300 ns, the number of hydrogen bonds increases and stabilizes at 14-15. These interactions remain consistent throughout the simulation, corroborating the molecular docking results, which show that most residues interact with the co-crystallized ligand via hydrogen bonds. These findings are directly related to the “Electrostatic Energy” component obtained from the Binding Free Energy calculations performed using the MM-GBSA methodology.

Close

Figure 9.

Number of hydrogen bonds formed by 5αR2 complexes with the four ligands under study throughout the 500 ns simulation. The X-axis represents the simulation time in nanoseconds (ns), and the Y-axis indicates the number of hydrogen bonds.

Export to PPT

In the complex formed by 5αR2 and COMP2, 15 hydrogen bonds are maintained up to 15 ns. However, from approximately 20 ns to 100 ns, up to 12 hydrogen bonds are formed. Beyond 100 ns, a decrease in the number of hydrogen bonds is observed, reaching an average of 10 bonds by 420 ns, which remains constant. From 420 ns until the end of the simulation at 500 ns, the number of hydrogen bonds slightly decreases to approximately 8 bonds. For the complex formed by 5αR2 and COMP9, the number of hydrogen bonds remains at an average of 6 bonds up to 100 ns) Subsequently, there is an increase up to 280 ns, where 11-12 hydrogen bonds are maintained constant. From that point until 500 ns, the number of hydrogen bonds decreases slightly, fluctuating between 10 and 11 bonds.

COMP11 exhibits a variation in the number of hydrogen bonds up to 300 ns. From 0 ns to 70 ns, there is a decrease in hydrogen bonds, ranging from 10 bonds to approximately 7 bonds. Subsequently, the number of hydrogen bonds increases to an average of 13 bonds, which remains constant until 180 ns. From that point onward, the number of hydrogen bonds stabilizes at approximately 12 bonds, persisting until 500 ns.

COMP14 maintains an average of 7 hydrogen bonds up to 250 ns. From 250 ns to approximately 400 ns, it exhibits 6 hydrogen bonds. From 400 ns until the end of the simulation, there is an increase in the number of hydrogen bonds, fluctuating between 10 and 12 bonds.

3.3.4. Calculation of binding free energies using MM-GBSA

In the Figure 10 presents the binding free energy results for the NDX ligand and the analogous compounds. The findings indicate that all the designed molecules are thermodynamically stable. COMP11 exhibited the highest binding free energy (-89.93 kcal/mol), demonstrating the most favorable binding with 5αR2. COMP2 and COMP9 followed with very similar energy values (-84.37 kcal/mol and -78.48 kcal/mol, respectively), while COMP14 showed a binding free energy of -49.74 kcal/mol. A representation of the binding free energy terms derived from MM-GBSA calculations for each ligand’s binding is shown in Table S2.

Table S2

Download

Close

Figure 10.

Representation of the binding free energy terms derived from MM-GBSA calculations for ligand binding.

Export to PPT

3.4. ADME-T properties

The assessment of pharmacokinetic and toxicological properties of new compounds is crucial for ensuring safety and efficacy in drug development [46]. Table 3 presents the results obtained using the ADMETlab 2.0 platform for the absorption and distribution of region A of the NDX ligand and region A of the four designed compounds that demonstrated the best molecular binding with the enzyme. It can be observed that these compounds exhibit high absorption at the gastrointestinal level. Regarding Caco-2 permeability (PAPP), the web server suggests that the optimal value should be greater than -5.15 log units, or ideally around -4.70 to -4.80 [47]. As shown in Table 3, while all compounds have values greater than -5.15, compound 2 exhibited the best PAPP value. In contrast to this parameter, all compounds demonstrated an exceptionally low probability of being substrates for P-Glycoprotein (Gpp).

On the other hand, all the molecules showed a high degree of binding to plasma proteins (UPP), indicating that the free fraction available to exert the therapeutic effect would be small (see Table 3) [48]. Additionally, the volume of distribution of the compounds was very low, likely because of the high UPP; however, it should be noted that this term is considered more theoretical than physiological [49]. Furthermore, the molecules displayed a high probability of crossing the BBB. These results are consistent with an adaptation of Lipinski’s Rule of Five, which postulates that a lipophilicity range of 2.0 to 3.5 is a predictive factor for the ability to cross the BBB through passive diffusion [50], and the molecules fall within this range (see Table S3).

Table 3 indicates that the molecules have a low probability of being inhibitors of the Cytochrome P450 3A-4 enzyme, while they show a high probability of acting as substrates for this enzyme. These enzymes participate in metabolizing approximately 90% of drugs, which can either minimize their therapeutic effect or make the body more susceptible to toxicity, depending on the extent to which the enzyme is induced or inhibited [51]. Regarding the excretion property descriptors, the half-life for the four compounds was low, based on the ranges provided by ADMETlab 2.0, where T1/2 ≤ 3h is considered short. For the clearance rate, the web server predicted a low clearance rate for all molecules, in accordance with the proposed ranges (>15 mL/min/kg: high; 5 mL/min/kg < Cl < 15 mL/min/kg: moderate; <5 mL/min/kg: low) [52].

The four compounds studied were unlikely to function as human ether-a-go-go-related gene (hERG) blockers or skin sensitizers. However, COMP14 exhibited the highest probability in both descriptors, with a moderate likelihood of blocking these channels and a slightly lower probability of generating sensitivity in the skin (see Table 4). Additionally, a low probability of significantly inducing genetic mutations was predicted for all molecules using Salmonella Typhimurium strains [53] although COMP11 had the lowest probability [52].

According to the ranges provided by the web server for the lethal dose 50 (LD₅₀), the molecules displayed moderate toxicity [51]. The analogues showed a low probability of inducing liver injury but had a greater than 50% chance of inducing human hepatotoxicity.

The selected molecules showed positive results for ‘Druglikeness,’ as they did not violate any of the proposed rules (see Table S3). This suggests that they are likely to have no issues with absorption. Veber’s rules support this conclusion, as indicated in Table 4: the sum of HBDs and acceptors was not greater than 12, there were fewer than 10 rotatable bonds, and the TPSA was less than 140 Ų [54]. Furthermore, these molecules were shown to be soluble in water, except for compounds 2 and 9, which are classified by the SwissADME server as moderately soluble [37].

Although the binding pocket of the NDX ligand is predominantly hydrophobic, the electrostatic interactions between the Arg114 residue and the carbonyl group’s carbon, as well as between the Ser31 amino acid and the nitrogen in the “N-ethyl-1-carboxamide” region (depicted in Figure 2), play a crucial role. The Arg114 residue forms an ionic bond, requiring deprotonation to interact with the nitrogen. Due to its resonance structures, it can simultaneously form hydrogen bonds and salt bridges.

On the other hand, the carbonyl group, when in its double-bonded state, favors hydrogen bonding between Ser31 and nitrogen. However, when the carbonyl group undergoes delocalization, it promotes the formation of a salt bridge. Among these interactions, the ionic interaction of Arg114 is among the strongest, likely contributing the highest energy and ensuring greater stability. Structural modifications in R₁ and R₂ are expected to influence these interactions, affecting the acidity and basicity of the carbonyl group in R₁, as well as the substituent group (in this case, sulfur) in R₂. These changes would, in turn, influence the binding energy, as predicted by molecular docking and MM-GBSA.

In COMP9, the addition of sulfur reduces the delocalization of the salt bridge with the Arg114 amino acid, while favoring the formation of a hydrogen bond with Ser31. This is due to sulfur’s lower electronegativity compared to oxygen. In COMP11 and COMP14, modifications with dimethyl and propyl hydrocarbon chains, respectively, exert a greater inductive effect on the nitrogen, making it more basic. This favors the formation of hydrogen bonds, similar to COMP9, while weakening the salt bridge between Arg114 and the carbonyl group. The dimethyl group has a more moderate inductive effect than the propyl group, and it also provides greater conformational flexibility, as it introduces less steric strain. In contrast, COMP2, with the dichloro group, reduces the electron density of nitrogen. This is because the dichloro group is highly electronegative and does not tend to donate electrons through resonance. As a result, the basicity of nitrogen decreases, favoring the salt bridge between Arg114 and the carbonyl group’s oxygen.

The 5αR2 enzyme, when unbound, has regions where amino acid residues exhibit high flexibility due to the β and γ turns they contain, which increases conformational freedom and facilitates backbone movement. However, NADP was able to significantly stabilize the enzyme molecule by reducing the RMSD of the backbone. Additionally, the flexibility of amino acids, as measured by RMSF, was decreased, as several flexible residues maintained polar interactions with NADP.

The A region of the native ligand binds to a predominantly hydrophobic binding pocket, as previously mentioned, interacting with amino acids such as Phe118, Phe219, and Phe223, which are also present in the docking protocol process (Figure 2). However, it also interacts with polar amino acid residues such as Glu57 and Arg114, located at the ends of DHF. These two residues are particularly important because Arg114 interacts with the C-3 carbonyl group of Finasteride, while Glu57 forms a hydrogen bond with the nitrogen of the amino group (N-4). Together, these interactions facilitate the positioning of NADPH so that the 4-pro-(R)-hydride of NADPH is positioned near the C-2 carbon, allowing the H⁺ transfer to the 1,2 bond of the drug. Therefore, it is crucial that these interactions were accurately reproduced in the designed analogues.

The docking analysis of these compounds revealed that COMP9 exhibited the highest affinity and binding energy with the 5αR2 enzyme. This compound successfully reproduced most of the interactions with the enzyme’s amino acid residues, particularly intermolecular hydrogen bonds. For instance, Glu57 functioned as an acid acceptor for H⁺ from the N-4 nitrogen, and Arg114 served as a basic donor of H⁺ to the sulfur atom of the drug, alongside interactions with amino acids that engage with NADP. This is further supported by the RMSD value of the ligand, which was superior to that of the native ligand, showing no fluctuations throughout the simulation. Additionally, the RMSF results indicated minimal fluctuations, similar to those observed for 5αR2-NDX (see Figure 8), except for a slight oscillation with the Leu40 residue, which was not observed in the molecular docking results (see Figure 3). This suggests that Leu40 does not interact with the analogue.

Electrostatic energy is another energetic component that demonstrates the stability of COMP11, as it exceeds that of the NDX ligand, suggesting enhanced stability of the complex. Additionally, Van der Waals interactions and π-stacking are evident in proportions very similar to those observed with 5αR2-NDX, which are important because they represent the initial drug-receptor interactions before polar interactions occur. This is further complemented by the selectivity conferred by these interactions, as hydrophobic interactions predominate in the NDX ligand. Polar interactions are crucial for molecular recognition between the enzyme and the constructed analogues, and, together with hydrophobic interactions, they contribute to the stabilization of the formed complexes [55].

COMP2 demonstrates that, despite having lower binding energy compared to COMP9, it exhibits an RMSD of ligands with consistent behavior, similar to that of COMP 9 (see Figure 8), and an RMSD of the backbone with only minimal, insignificant fluctuations. Its high electrostatic energy component is evident from the numerous polar interactions observed in the molecular docking results and the hydrogen bonds it maintains throughout the simulation. This compound displayed more fluctuations in the RMSF, particularly with residues Gly34, Lys35, His36, Thr37, Glu38, Leu73, Ser74, and Leu75, indicating greater flexibility in these amino acids and a decrease in stability in that region. However, these residues have not been identified as key amino acids in the inhibitory activity of the NDX ligand.

The fluctuations in COMP2 are less pronounced compared to those observed in the amino acid residues of COMP9, inferring a greater interaction with NADP during the simulation. This could explain the variation in electrostatic energy values observed in the binding free energy calculations using MM-GBSA (see Figure 10). Similarly, the Van der Waals force component is affected, as some of the fluctuating amino acids maintain hydrophobic interactions, which may be slightly altered over time.

COMP11 was the analogue that maintained the most hydrogen bonds throughout the simulation, ranging from 7.5 to 12.5, and stood out for its high electrostatic energy. The RMSD of the ligands was similar to that of the native ligand. However, the RMSD of the backbone from 50 ns to approximately 300 ns showed the greatest fluctuation compared to the other compounds during the same simulation time. This variability could negatively impact the stability of the complex.

COMP14 exhibited the highest number of hydrogen bonds compared to the other compounds, but only during the 20–45 ns range. For the rest of the simulation, these polar interactions were maintained at a similar level to those of the native ligand. Although the RMSF of COMP14 was very similar to that of 5αR2-NDX, this compound displayed the fewest interactions with amino acid residues in the molecular docking results. This lack of interactions may explain why the RMSD of the backbone exceeded 0.2 nm compared to the NDX ligand, suggesting that NADP did not consistently interact with all amino acid residues of the enzyme. This conformational freedom could have introduced flexibility in the backbone atoms.

Regarding the toxicology section measured through ADME-T properties, these analogues showed a probability of over 50% of causing harmful effects in some of the evaluated aspects. COMP2, which contains a dichloro group, could potentially cause damage to biological tissues due to its ability to form free radicals, as it has an unpaired electron [56], similar to the sulfur present in COMP9. The aliphatic substituents in compounds 11 and 14 may accumulate in liver cells, with their length and branching determining their affinity for liver enzymes. This, in turn, could affect their metabolism, increasing the potential for liver toxicity. However, since these aliphatic groups are not highly branched, the probability of significant toxicity is relatively low [57].

The pharmacokinetic properties were generally positive, except for the distribution process, which showed a high probability of binding to plasma proteins—a property already observed in the native ligand. This behavior may be influenced by the position of the substituents and the electrostatic and hydrophobic interactions within the complexes, as well as the conformational flexibility of NADP, which enhances its likelihood of binding to plasma proteins. While the compounds performed well overall in ADME assessments, COMP2 and COMP9 yielded the best results.