Fragment-based design of potent non-nucleoside inhibitors targeting the thumb domain of HCV NS5B

2.1. Dataset preparation, QSAR model construction, and validation

In this study, 2D structures of 35 NNIs of HCV NS5B inhibitors, obtained from a previously published work [11], were converted to three-dimensional (3D) conformers, and the lowest-energy conformer of each molecule was retained following geometry optimization using the MMFF94 force field and the steepest descent algorithm for up to 1000 steps using RDKit [12]. Partial atomic charges were assigned using the MMFF94 protocol. Molecular descriptors were calculated using alvaDesc software (version 1.0.8) based on the optimized 3D structures to develop the QSAR model [13]. To address multicollinearity in the QSAR model, descriptors with constant values exceeding 95% were removed, and only one descriptor from pairs with a correlation coefficient greater than 0.9 was retained. The experimental half-maximal inhibitory concentration (EC₅₀) values were converted into negative logarithmic form (-log EC₅₀). The molecular structures and their corresponding EC₅₀ and pEC₅₀ data are provided in Table S1. The QSAR model was developed using the genetic algorithm-multiple linear regression (GA-MLR) approach, implemented in QSARINS software [14]. The validity of the model was comprehensively confirmed using a variety of validation techniques, including application domain evaluation (AD), internal and external validation, and Y-randomization tests [15,16]. The statistical reliability and robustness of the GA-MLR QSAR model were assessed using four validation steps: (a) Y-randomization tests; (b) external validation; (c) internal validation using the leave-one-out (LOO) cross-validation approach; and (d) adherence to preset threshold values for important statistical parameters. These parameters included a determination coefficient for the training set (R²_tr) and for the test set (R²_test) of at least 0.6, a cross-validated correlation coefficient squared from the LOO procedure (Q²_loo) of at least 0.5. Furthermore, the model’s predictive capability was evaluated using additional statistical metrics, such as Q²_Fn, the Fisher ratio (F), and the standard deviation of errors of prediction (s). The Q²_Fn metric assesses the model’s predictive performance on the external test set, while the Fisher ratio (F) evaluates the overall significance of the regression model. The standard deviation of errors of prediction (s), external concordance correlation coefficient (CCC_test), root-mean-square error (RMSE), and mean absolute error (MAE) were also analyzed. Lower values of these metrics indicate superior model performance. Models that failed to meet any of these validation criteria were eliminated to ensure the selection of the most reliable and predictive QSAR model.

2.2. Fragment-based database preparation

The FBDD approach was employed to optimize lead compounds activities targeting HCV NS5B polymerase [17]. Two ligands, L0 and C0 (lead compound, co-crystallized compound), were selected from the previous dataset [11]. L0 exhibits the highest biological activity, and C0 demonstrates notable selectivity and preclinical potential. The experimental 3D structure of HCV NS5B polymerase (PDB: 5CZB) was retrieved from the RCSB Protein Data Bank to guide fragment replacement strategies (Figure 1). Prior to fragment generation, missing residues in the 3D structure were modelled using the SWISS-MODEL web server [18]. This platform uses a template-searching approach based on BLAST and HHblits to identify suitable homologous structures from the SWISS-MODEL Template Library (SMTL). For this, chain B was selected for modeling due to its relatively fewer missing residues compared to chain A allowing more reliable structural reconstruction. The obtained complete structure ensures an accurate fragment generation, molecular docking, and MD simulations. Key regions for modification were identified based on their critical roles in target interaction and prior structure-activity relationship (SAR) studies indicating their potential to enhance activity (Figure 1). Fragment substitutions were executed using FragRep, a web server that evaluates geometric compatibility and local protein environmental alignment from a database of over 500,000 fragments [19]. Core modifications focus on optimizing hydrogen bonding, pi-pi stacking, and hydrophobic interactions to significantly enhance binding affinities, selectivity, and pharmacokinetic profiles. Generated compounds were initially represented in SMILES (simplified molecular input line entry system) notations and were converted to 3D structures using RDKit, following the protocol described in Section 2.1. The obtained 3D structures were subjected to further computational analyses and optimization within the drug discovery pipeline.

Close

Figure 1. Structural insights into NS5B binding site and fragment-based library design: (a) The left panel illustrates the binding site of the NS5B RNA-dependent RNA polymerase (PDB: 5CZB), (b) The right panel presents the core structures of selected compounds (L0 & C0) used to generate the FBDB, with highlighted regions indicating modified chemical moieties.

Export to PPT

2.3. QSAR-based screening, Lipinski’s rule of five, and synthetic accessibility assessment

The compounds generated through the FBDD approach were subjected to a rigorous screening process using the previously described machine learning-based QSAR model [20]. This filtration step prioritized candidates with strong predicted inhibitory activity against HCV NS5B polymerase based on key molecular descriptors and structural features influencing enzyme inhibition. With a predetermined cutoff based on the pEC₅₀ of the lead compound, compounds were ranked according to their predicted activity (pEC₅₀ values) and selectivity profiles. Only compounds that are above or at the cutoff value were selected for further study. After QSAR screening, Lipinski’s Rule of Five and synthetic accessibility evaluations were used for further filtering using DataWarrior software [21]. To assess the physicochemical characteristics necessary for oral bioavailability, Lipinski’s rule was used. The selected compounds were examined according to their molecular weight (MW <500 g/mol), lipophilicity (logP ≤5), hydrogen bond donors (HBD ≤5), and hydrogen bond acceptors (HBA ≤10). Although only one of these violations is usually tolerated during drug design, the main observation drawn from this study was that molecular weight tended to exceed 500 g/mol. In parallel, synthetic accessibility was assessed to determine the feasibility of chemical synthesis [22]. A cutoff score of 5 was applied for synthetic accessibility, and all compounds with scores exceeding this threshold were excluded from consideration.

2.4. Molecular docking

Before the molecular docking studies, the docking protocol was first validated by a redocking of the co-crystallized inhibitor C0 (IDX17119) into the NS5B polymerase structure. The docking grid was defined based on the coordinates of IDX17119, with grid box parameters set at x = 27.7, y = 1.9, and z = 47.88, and dimensions of 30 ų. For protein preparation, all water molecules and the co-crystallized ligand were removed, polar hydrogens were added, and Gasteiger charges were assigned with AutoDockTools [23]. The protein and ligand structures were converted to PDBQT format for compatibility with AutoDock Vina, and the docking was performed with an exhaustiveness of 20 to ensure comprehensive conformational sampling [24]. The docking results were analyzed using Discovery Studio Visualizer and PyMOL to evaluate the important ligand-protein interactions. Following validation, docking calculations were performed to investigate the binding interactions between the modeled NS5B protein and the designed compounds, along with the reference lead compound.

2.5. Molecular dynamics (MD) simulations

To thoroughly evaluate the stability and dynamic behavior of the NS5B protein and its interactions with selected compounds, independent MD simulations were performed in triplicate for each system using GROMACS 2021.3 [25]. Prior to MD simulation, the protonation states of the structure were generated considering the physiological pH (7.0) using the H + + server, and all missing hydrogen atoms were added. The CHARMM-GUI web server was employed to generate initial input parameters based on the CHARMM36m force field [26]. Each system was solvated in a TIP3P (transferable intermolecular potential 3 point) water model within a rectangular grid box, ensuring a 15 Å margin in all directions [27]. To maintain charge neutrality and to achieve a salt concentration of 0.15 M, counterions (Na⁺, Cl^-) were added via Monte Carlo ion displacement. Periodic boundary conditions (PBC) were applied, and non-bonded interactions were handled using Lennard-Jones and Coulomb potentials with a 12 Å cut-off distance, employing the Verlet cut-off scheme for efficient neighbor searching [28]. Long-range electrostatic interactions were treated using the particle-Mesh Ewald (PME) method, while the linear constraint solver (LINCS) algorithm was used to constrain all covalent bonds, including those involving hydrogen atoms [29,30]. Energy minimization was carried out using the steepest descent method with a maximum of 50,000 steps and a force threshold of 10.0 kJ/mol. System equilibration was conducted in two steps: a 500 ps NVT ensemble at 303.15 K using the Nose-Hoover thermostat, followed by a 500 ps NPT ensemble at 303.15 K and 1.01325 bar using the Parrinello-Rahman barostat [31,32]. All systems were simulated in triplicate for 200 ns to rigorously evaluate the reproducibility of results and the reliability of the estimates and associated errors. Structural analyses, including root mean square deviation (RMSD), radius of gyration (Rg), root mean square fluctuation (RMSF), and solvent-accessible surface area (SASA), were performed using built-in GROMACS modules. The resulting data were visualized using Xmgrace.

2.6. Molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) calculations

The calculation of the binding free energies of protein-ligand complexes was estimated by the MM-PBSA method using gmx_MMPBSA version 1.4.3, a GROMACS-compatible tool designed for MM-PB(GB)SA analysis [33]. To ensure statistical robustness, binding energy calculations were performed in line with the triplicate MD simulations of each system. The calculations were performed on 2000 snapshots extracted from a total of 20,000 frames, with an interval of 10. For these calculations, the Poisson–Boltzmann (PB) solvation model was chosen.

The binding free energy (ΔG_binding) was determined using the thermodynamic relationship (Eq.1):

where ΔG_complex denotes the total free energy of the protein-ligand complex, and ΔG_receptor and ΔG_ligand correspond to the free energies of the isolated protein and ligand, respectively.

Furthermore, the total binding energy can be decomposed into gas-phase and solvation contributions as follows (Eq.2):

where ΔE_gas represents the gas-phase interaction energy, which includes van der Waals (ΔE_vdw) and electrostatic (ΔE_ele) components. The solvation free energy (ΔG_solv) accounts for both polar (ΔG_polar) and nonpolar (ΔG_nonpolar) contributions.

2.7. ADMET evaluation

The ADMET analysis plays a critical role in drug discovery and development. Based on the ADMET results, we could predict the efficacy and safety of the potential hit compound in the early stages of drug development. Accordingly, the ADMET predictions for the newly designed compounds in comparison with the reference compounds (L0 and C0) were performed using ADMET-SAR and pkCSM web servers [34,35]. The absorption and distribution parameters, including blood-brain barrier permeability (BBB), human intestinal absorption (HIA), total polar surface area (TPSA), and solubility (logS), were first calculated to evaluate the potential of the compounds’ bioavailability properties. Then, metabolic interactions were predicted by analyzing the inhibition potential of five key cytochrome P450 isoforms (CYP1A2, CYP2C9, CYP2C19, CYP2D6, and CYP3A4). The excretion profiles, including clearance rates and half-life (T_1/2), were estimated. For a full assessment of the safety of the compounds, toxicity predictions assessed AMES toxicity, tumorigenic potential, irritant effects, reproductive toxicity, and carcinogenicity. This ADMET analysis provides valuable insights into their pharmacokinetic properties.

3.1. GA-MLR QSAR model

The GA-MLR model was constructed to predict the biological activity (pEC₅₀) of the studied compounds. The model was built based on four key molecular descriptors SHED_AL, P_VSA_charge_1, CATS2D_09_DL, and MaxaasC selected from a comprehensive set of descriptors using genetic algorithms. These descriptors were chosen for their significant contribution, explaining the relationship between the chemical structures and their corresponding pEC₅₀ values. Their definitions and relevance are provided in Table S2.

The generated model using the GA-MLR approach (Eq.3) and its corresponding statistical parameters are shown below:

$\begin{array}{l}{R}^2{}_{Tr}=0.908,Q^2{}_{LOO}=0.880,R^2{}_{test}=0.907,MA{E}_{test}=0.589,\\ RMS{E}_{test}=\mathrm{0.639.}\end{array}$

$\begin{array}{l}{Q}^2{}_{F1}=0.833,Q^2{}_{F2}=0.833,Q^2{}_{F3}=0.841,CC{C}_{ext}=0.924,\\ F=56.533,s=\mathrm{0.536.}\end{array}$

The assessment of the parameters in the constructed GA-MLR model meets the validation criteria specified by the OECD (Organisation for economic co-operation and development) guidelines. Additionally, the scatter plot in Figure 2(a)

Close

Figure 2. (a) Experimental (X-axis) vs. Predicted (Y-axis) pEC₅₀; (b) Williams plot: Leverage (X-axis) vs. Standardized residuals (Y-axis).

Export to PPT

3.2. Fragment-based drug design (FBDD)

The FBDD approach was used to generate a structurally diverse library of derivatives from both L0 and C0 ligands, resulting in 7,768 and 9,034 compounds, respectively. This approach aimed to enhance molecular interactions within the NS5B active site by integrating key structural features of the parent compounds with novel modifications designed to improve binding affinity and pharmacokinetic properties. Generated compounds were initially represented as SMILES notations and subsequently processed to remove duplicates by identifying and eliminating redundant entries. After this verification step, the unique SMILES were converted into their corresponding 3D structures using RDKit and saved in structure data file (SDF) format for further computational analyses.

3.3. Filtration process

The filtration of the generated compounds was conducted in a stepwise manner to refine the library and identify promising candidates for further analysis. Initially, the developed GA-MLR model was used to exclude compounds with predicted pEC₅₀ values lower than 9.53 (pEC₅₀ of L0), resulting in 1881 and 2558 modified compounds from the L0 and C0 compounds, respectively. These selected compounds (4439) were then screened according to Lipinski’s rule of five to assess their drug-likeness based on the key physicochemical properties such as MW, LogP, and hydrogen bond acceptors and donors. While MW is a crucial factor influencing absorption and permeability, its threshold is not absolute, as larger molecules can still exhibit drug-like properties through other compensatory mechanisms. Following this filtration step, 1,076 L0-derived compounds and 967 C0-derived compounds met the criteria and were retained. Finally, synthetic accessibility (SA) was evaluated to prioritize compounds that are easier to synthesize, applying a cutoff value of 5. Only those with SA scores below this threshold were selected, reducing the library to 208 L0-derived and 374 C0-crystallized-derived compounds (Supplementary materials, Spreadsheet). This systematic refinement ensured the selection of candidates with high predicted activity, favorable drug-like properties, and feasible synthetic accessibility, streamlining the subsequent phases of the drug discovery process.

3.4. Molecular docking

To validate the docking protocol, a redocking experiment was performed using the co-crystallized ligand C0 within the active site of NS5B. The predicted conformation obtained from AutoDock Vina closely matched the experimental conformation with an RMSD of 0.26 Å, confirming the accuracy and reliability of the docking methodology used in reproducing native poses (Figure S1). This validation ensured the robustness of the approach for subsequent screening of potential inhibitors. Following the validation step, all 208 compounds from the L0-designed library (L1-L208) and 374 compounds from the C0-designed library (C1-C374) were subjected to molecular docking to evaluate their binding potential. Among the 582 docked compounds, seven (C341, C345, C350, C360, L151, L168, and L171) exhibited the most promising binding affinities, ranging from -10.2 to -10.6 kcal/mol. Their chemical structures, binding affinity values, and interaction details are given in Figure 3, Table 1, and Table S3 respectively.

Close

Figure 3. Chemical structures of the top newly designed compounds.

Export to PPT

Figures 4(a-c) shows the 2D representation of the selected compounds within the predicted binding site pocket of NS5B. The binding affinities of the ligands were evaluated, with L171 and C341 showing the highest affinity (-10.6 kcal/mol), supported by an extensive interaction network. Ligand L171 formed hydrogen bonds with SER477 and PRO496, a carbon-hydrogen bond with SER477, and multiple π-sigma and π-π stacked interactions with LEU490 and TRP529. Hydrophobic interactions with ALA487, LEU490, ARG491, LEU420, ARG423, MET424, and LEU498, along with π-alkyl and π-sulfur interactions, further stabilized its binding in the active site. C341 also exhibited strong interactions, including hydrogen bonds with TYR478 and ALA487 and a carbon-hydrogen bond with TRP529. Additional π-sigma, π-π stacking, and hydrophobic interactions with residues like PRO496, ALA487, ARG491, LEU420, ARG423, and MET424 reinforced its stability. C360 (-10.5 kcal/mol) formed hydrogen bonds with TYR478 and PRO496, supported by π-sigma, π-π stacked, and π-sulfur interactions with TRP529 and MET424, along with hydrophobic contacts with PRO497, LEU420, ARG423, and ARG491. L168 (-10.5 kcal/mol) demonstrated strong binding through hydrogen bonds with ARG491 and PRO496, π-sigma interactions with LEU420, LEU490, and TRP529, and π-sulfur interactions with MET424. C345 (-10.4 kcal/mol) formed three hydrogen bonds with SER477, TYR478, and PRO496, while π-sigma and π-π stacking interactions with TRP529, along with extensive hydrophobic and π-alkyl contacts with ARG491, LEU498, LEU420, ARG423, and MET424, further stabilized the ligand. C350 (-10.2 kcal/mol) displayed a distinctive binding pattern with hydrogen bonds to VAL495 and PRO496, π-donor hydrogen bonding with LEU498 and SER477, π-π stacking interactions with TRP529, and extensive hydrophobic interactions with ARG491, LEU420, ARG423, and LEU498. L151 (-10.0 kcal/mol) was stabilized through hydrogen bonds with SER477 and TYR478, π-sigma and π-π stacking interactions with TRP529, and hydrophobic contacts with LEU420, ARG423, and MET424. These results suggest that the binding affinities are significantly influenced by the ability of the designed ligands to form key interactions with critical residues in the active site, such as LEU420, ARG423, MET424, ARG491, and TRP529, which are essential for stabilizing the ligand-protein complex and enhancing the overall potency.

Close

Figure 4. (a) 2D representation of the ligand-protein interactions for compounds C0 and L0 within the NS5B binding pocket. (b) 2D representation of the ligand-protein interactions for compounds C341, C345, C350, and C360 within the NS5B binding pocket. (c) 2D representation of the ligand-protein interactions for compounds L151, L171, and L168 within the NS5B binding pocket.

Export to PPT

The QSAR model (Eq. 3) provides a mechanistic explanation for the enhanced activity predicted for the newly designed compounds C341, C345, C360, C350, L151, L168, and L171. The positive coefficient associated with SHED_AL (0.5586) highlights the importance of increasing the spatial diversity between hydrogen-bond acceptor groups and hydrophobic regions within a molecule. The large bicyclic or polycyclic fragments introduced into these compounds play a central role in achieving this effect. In C341, C345, C360, and L151, the introduced bicyclic amine-containing scaffolds substantially enlarge the hydrophobic domain, while in L168 and L171, the more complex nitrogen-containing polycyclic systems provide both extended lipophilicity and additional acceptor/donor sites. These fragments are positioned at the opposite end of the scaffold from the acceptor-rich functionality, thereby increasing the entropy-based separation between polar and lipophilic features and producing SHED_AL values consistent with improved predicted potency. The QSAR model also indicates a positive contribution for P_VSA_charge_1 (0.0255), suggesting that localized positive electrostatic potential enhances activity. This aligns with the structural modifications in C341, C345, C360, and L151, where primary amines attached to the hydrophobic scaffold contribute to a larger positively charged surface area. Similarly, L168 and L171 contain protonatable heterocyclic nitrogen atoms within their polycyclic fragments, providing an additional source of positive polarization. These features increase the likelihood of favorable electrostatic interactions with electron-rich regions of the NS5B binding pocket. The negative coefficient of CATS2D_09_DL (-0.4384) further suggests that excessive separation between donor and lipophilic atoms is detrimental to activity. The newly introduced fragments across all compounds address this requirement by maintaining a compact spatial relationship between donor functionalities and hydrophobic carbons. The amine donors in C341, C345, C360, and L151, and the embedded donor sites within the polycyclic frameworks of L168 and L171, remain in close topological proximity to surrounding lipophilic atoms. This structural compactness naturally keeps the donor-lipophilic distances short, yielding lower values for the CATS2D_09_DL descriptor, which is favorable for potency. Finally, MaxaasC shows the strongest positive coefficient (3.8303), indicating that the presence of atoms with high charge polarization significantly enhances predicted activity. The amine-bearing fragments in C341, C345, C360, and L151 introduce strongly polarized nitrogen atoms, while L168 and L171 provide even greater charge diversity through multiple heteroatoms with varying degrees of protonation and electron density. These highly polarized atoms elevate the MaxaasC descriptor, reinforcing their contribution to inhibitory potency. This QSAR interpretation validates the design strategy, demonstrating that the new fragments integrated into C341, C345, C360, C350, L151, L168, and L171 are structurally aligned with the molecular features identified as essential for improved NS5B inhibition.

3.5. Molecular dynamics (MD) simulations

MD simulations were performed to assess the stability of the newly designed compounds when bound to the NS5B protein. Four key parameters were analyzed: RMSD, RMSF, Rg, and SASA, with data averaged over three independent simulations (Table 2 and Table S4). These metrics offer insights into the structural dynamics and binding stability of the protein-ligand complexes. RMSD, which measures the deviation of the protein backbone from its initial structure, serves as an indicator of complex stability. Lower RMSD values correspond to more stable interactions. After 100 ns, all the designed compounds complexed with NS5B reached a convergence point and showed greater stability compared to the apo protein (NS5B-apo) and the complexes of NS5B with reference compounds (NS5B-L0 and NS5B-C0) (Figure 5). Notably, NS5B-L151 and NS5B-L168 complexes exhibited the lowest RMSD values (0.133 ± 0.0124 nm and 0.134 ± 0.0141 nm, respectively), indicating strong binding with minimal structural change. NS5B-L171 and NS5B-C350 also showed low RMSD values (both around 0.136 nm), while NS5B-C360 maintained a stable interaction at 0.141 ± 0.0148 nm. In contrast, NS5B-C345 and NS5B-C341 displayed slightly higher RMSD values (0.149 ± 0.0202 nm and 0.152 ± 0.0198 nm), suggesting greater fluctuation and somewhat less stable binding. RMSF values provide information about the flexibility of the protein-ligand complexes; lower RMSF values indicate more rigid, stable interactions (Figure 6). NS5B-L151, NS5B-L168, and NS5B-L171 complexes showed RMSF values around 0.077-0.079 nm, reinforcing the observation from the RMSD results. Similarly, NS5B-C350 and NS5B-C360 demonstrated low RMSF values (approximately 0.078 and 0.083 nm, respectively). However, NS5B-C345 had a slightly elevated RMSF of 0.088 nm, and the RMSF value of NS5B-C341 (about 0.085 nm) further supports the idea of increased flexibility and less stable binding. The Rg metric indicates the compactness of the protein structure; a lower Rg value suggests a more compact and stable complex (Figure 7). Here, L151, L168, and L171 complexed with NS5B formed very compact complexes with Rg values in the range of 2.424-2.430 nm. The complexes of NS5B with C350, C345, and C360 also maintained compact structures with similar Rg values, while NS5B-C341 exhibited a slightly higher Rg, hinting at a marginally expanded structure that might correlate with less stable binding. SASA values reveal how much of the protein surface is exposed to solvent. Lower SASA values generally correlate with more stable, less solvent-exposed interactions (Figure 8). Complexes with L151, L168, and L171 had SASA values around 248-250 nm², indicating moderate exposure and stable binding. NS5B-C350, NS5B-C345, and NS5B-C360 showed comparable SASA values, whereas NS5B-C341 had the highest SASA at approximately 253 nm², suggesting a more extended structure with increased solvent exposure that could weaken its binding interactions. Overall, the MD simulation results indicate that complexes NS5B-L151, NS5B-L168, and NS5B-L171 provide robust binding with strong structural stability, as evidenced by their low RMSD, RMSF, and Rg values coupled with moderate SASA levels. While the complexes involving NS5B-C350, NS5B-C345, and NS5B-C360 also display stable interactions, they show slightly higher solvent accessibility. In contrast, NS5B-C341 appeared to be less stable, marked by higher flexibility and increased SASA, which may point to weaker binding interactions. These findings offer a comprehensive view of the structural dynamics and stability of the NS5B protein when complexed with these ligands.

Close

Figure 5. Time-dependent RMSD of carbon-α backbone of (a) NS5B-apo, NS5B-C0, NS5B-C341, NS5B-C345, NS5B-C350, NS5B-C360 and (b) NS5B-apo, NS5B-L0, NS5B-L151, NS5B-L168, NS5B-L171.

Export to PPT

Close

Figure 6. Plot of RMSF for c-α atoms of (a) NS5B-apo, NS5B-C0, NS5B-C341, NS5B-C345, NS5B-C350, NS5B-C360 and (b) NS5B-apo, NS5B-L0, NS5B-L151, NS5B-L168, NS5B-L171.

Export to PPT

Close

Figure 7. Plot of Rg vs. time for of c-α backbone of (a) NS5B-apo, NS5B-C0, NS5B-C341, NS5B-C345, NS5B-C350, NS5B-C360, and (b) NS5B-apo, NS5B-L0, NS5B-L151, NS5B-L168, NS5B-L171.

Export to PPT

Close

Figure 8. SASA plots vs. time of (a) NS5B-apo, NS5B-C0, NS5B-C341, NS5B-C345, NS5B-C350, NS5B-C360 and (b) NS5B-apo, NS5B-L0, NS5B-L151, NS5B-L168, NS5B-L171.

Export to PPT

3.6. MM-PBSA calculations

The MM-PBSA method is employed to evaluate the binding affinity of compounds to target proteins by calculating the binding free energy. A more negative binding free energy generally correlates with a stronger interaction between the ligand and the target protein. This approach considers contributions from gas-phase energy, solvation effects, and entropy. To enhance the accuracy and reliability of the binding energy estimate, the values are averaged over multiple snapshots from the molecular dynamics trajectory. All the results are presented in Table 3 and Table S5. Both NS5B-C345 and NS5B-L168, with ΔH total values of -33.78 ± 0.70 kcal/mol and -33.20 ± 0.28 kcal/mol, respectively, demonstrate strong binding to the NS5B protein. These compounds show favorable interactions, as reflected in their highly negative ΔH_total values, suggesting that they form stable ligand-protein complexes. The significant electrostatic interactions, particularly for NS5B-C345, with a ΔE_ele of -168.13 ± 4.01 kcal/mol, likely play a crucial role in their high binding affinity. Additionally, both compounds have relatively favorable solvation energies (ΔE_gb and ΔE_surf), contributing to the stability of the complexes in aqueous environments. However, consideration of the selected compounds dataset reveals that L171, C360, L151, and C341 also demonstrated markedly improved binding affinity relative to the reference inhibitors (C0 and L0). These compounds showed ΔH_total values ranging from -23.20 to -28.23 kcal/mol, forming stable complexes with NS5B and displaying well-balanced energetic contributions comparable to the top-ranked ligands. These results highlight C345 and L168 as top-performing ligands, while also demonstrating that several designed compounds possess favorable thermodynamic binding profiles, supporting their potential as viable NNIs for further investigation, as their binding affinity is likely to lead to potent inhibitory effects on the NS5B protein. Their strong binding, coupled with favorable solvation and enthalpic contributions, positions them as promising candidates for drug development targeting hepatitis C.

3.7. ADMET results

The pharmacokinetic and toxicity profiles of the new compounds were compared to the reference compounds (C0 and L0) as summarized in Table 4. In the absorption and distribution domain, nearly all compounds demonstrated excellent BBB characteristics (values between 0 and 0.3), indicating limited central nervous system exposure; however, C345 exhibited a medium BBB value (0.583), which may suggest an increased potential to cross the barrier. Regarding HIA, most compounds achieved an excellent rating (HIA < 0.3), with L168 (0.38) and L171 (0.323) classified as medium, thereby showing improved absorption compared to the references (C0: 0.01, L0: 0.013). Additionally, the TPSA values ranged from 114.54 to 149.66 Ų, with the reference compounds having the highest values, implying that some of the new compounds might offer enhanced membrane permeability. In terms of metabolism, only C350 was identified as an inhibitor of several CYP450 isoforms (2C19, 3A4, 2C9, and 1A2), potentially posing drug-drug interaction risks, while all other compounds, including the reference molecules, did not show significant CYP450 inhibition. The excretion profile revealed that C350 had an excellent clearance value (6.356 mL/min/kg, with ≥5 considered excellent), indicating rapid elimination; the remaining compounds, however, demonstrated poor clearance (<5). Notably, all compounds exhibited an excellent half-life (T_1/2 < 0.3), suggesting rapid elimination that might affect bioavailability. Regarding toxicity, apart from C350, which showed positive AMES toxicity and was flagged for carcinogenicity, none of the compounds displayed mutagenic, tumorigenic, or irritant effects. The reference compound C0 was also marked for carcinogenicity, adding an extra layer of caution. It is important to note that all designed molecules, except C350, preserve a terminal carboxylic acid group, a functional motif commonly associated with improved solubility, controlled permeability, and reduced metabolic liabilities in antiviral scaffolds. In contrast, C350 contains a bulkier hydrophobic cyclic substituent, which likely contributes to its higher lipophilicity, predicted CYP450 inhibition, and toxicity alerts. This structure-property difference aligns with reported ADMET trends, where increasing hydrophobicity may enhance potency but often compromises solubility, metabolism, and safety. Overall, the newly designed compounds demonstrate improved absorption characteristics and metabolic profiles that are comparable to the reference compounds. However, some of these compounds display distinct pharmacokinetic and safety features, such as variations in blood-brain barrier permeability, CYP450 interactions, and clearance, that merit additional investigation as they progress through further stages of drug development.