In silico Design of Novel N-hydrosulfonylbenzamides inhibitors of dengue RNA-dependent RNA polymerase showing favorable predicted pharmacokinetic profiles

Melalie Keita

Email ID: keitamelalie@yahoo.fr;

Tel.: +225-07-07-82-93-40

Kouakou Kouakou Jean-Louis ¹, Melalie Keita ^1*, Akori Elvice Esmel ¹, Brice Dali ¹, Aubin N'Guessan ¹, Affiba Florance Kouassi ¹

¹ Laboratoire de Physique Fondamentale et Appliquée (LPFA), University of Abobo Adjamé (now Nangui Abrogoua), Abidjan, Côte d’Ivoire;

Author Email ID: louisetjojo@yahoo.fr, keitamelalie@yahoo.fr, elvicee@yahoo.fr,

dalibrice@yahoo.fr, nguessaubin@gmail.com, akouassi859@yahoo.com  

Kiani F(farzad.kiyani@gmail.com)

Ionescu MI(mionescu@umfcluj.ro)

Sarma PVGK(pvgksarmasvims@gmail.com)

Melalie Keita, Kouakou Kouakou Jean-Louis, Akori Elvice Esmel, Brice Dali, Aubin N'Guessan, Affiba Florance Kouassi, In silico Design of Novel N-hydrosulfonylbenzamides inhibitors of dengue RNA-dependent RNA polymerase showing favorable predicted pharmacokinetic profiles(2021) Journal of Computational Chemistry & Molecular Modeling 5(2) p:561-584

Background: In recent years, there has been a growing interest in Denv NS5 inhibition, with several reported RdRp inhibitors such as sulfonylbenzamides, non-nucleo-side inhibitors without any 3D-QSAR pharmacophore (PH4) available. In this context, we report here, in silico design and virtual evaluation of novel sulfonylbenzamides Denv RdRp inhibitors with favorable predicted pharmacokinetic profile.

Methods: By using in situ modifications of the crystal structure of 5-(5-(3-hydroxyprop-1-yn-1-yl)thiophen-2-yl)-4- methoxy-2-methyl-N-(methylsulfonyl) benzamide (EHB)-RdRp complex (PDB entry 5HMZ), 3D models of RdRp-EHB_x complexes were prepared for a training set of 18 EHBs with experimentally determined inhibitory potencies (half-maximal inhibitory concentrations IC₅₀^exp). In the search for active conformation of the EHB1-18, linear QSAR model was prepared, which correlated computed gas phase enthalpies of formation ∆∆H_MM of RdRp-EHB_x complexes with the IC₅₀^exp. Further, considering the solvent effect and entropy changes upon ligand binding resulted in a superior QSAR model correlating computed complexation Gibbs free energies (∆∆G_com). The successive pharmacophore model (PH4) generated from the active conformations of EHBs served as a virtual screening tool of novel analogs included in a virtual combinatorial library (VCL) of compounds with scaffolds restricted to phenyl. The VCL filtered by the Lipinski’s rule-of-five was screened by the PH4 model to identify new EHB analogs.

Results: Gas phase QSAR model: -log₁₀(IC₅₀^exp) = p IC₅₀^exp =-0.1403 x ∆∆H_MM _ 7.0879, R² = 0.73; superior aqueous phase QSAR model: p IC₅₀^exp = -0.2036 x ∆∆G_com + 7.4974, R² = 0.81 and PH4 pharmacophore model: p IC₅₀^exp = 1.0001 x p IC₅₀^pre -0.0017, R² = 0.97. The VCL of more than 30 million EHBs was filtered down to 125,915 analogs Lipinski’s rule. The five-point PH4 screening retained 329 new and potent EHBs with predicted inhibitory potencies p IC₅₀^pre up to 30 times lower than that of EHB1 (IC₅₀^exp = 23nM). Predicted pharmacokinetic profile of the new analogs showed enhanced cell membrane permeability and high human oral absorption compared to the alone drug to treat dengue virus.

Conclusions: Combined use of QSAR models, which considered binding of the EHBs to RdRp, pharmacophore model and ADME properties helped to recognize bound active conformation of the sulfonylbenzamide inhibitors, permitted in silico screening of VCL of compounds sharing sulfonylbenzamide scaffold and identify new analogs with predicted high inhibitory potencies and favorable pharmacokinetic profiles.

Keywords: ADME properties prediction, Dengue, 3-(5-ethynylthiophen-2-yl)-N-hydrosulfonylbenzamides, in silico screening, RNA-dependent RNA polymerase.

Emerging" and / or "re-emerging" diseases have been public health deep concern in recent decades. The incidence of dengue fever is currently increasing dramatically, and is now included among the so-called "re-emerging" diseases. The World Health Organization (WHO) estimates the number of annual cases at 50 million, including 500,000 cases of dengue haemorrhagic fever which are fatal in more than 20% of cases [1]. According to data received from the Early Warning System, part of the WHO Public Health Case Management System - Event Management System (EMS) - 52 public health cases have been reported to the WHO Regional Office for Africa between January and September 2014, of which 94% (49/52) were due to infectious diseases; dengue (11%) after cholera (33%) [1].

Currently no specific medicine against dengue fever and the only approved vaccine, Dengvaxia®, developed by Sanofi Pasteur [2] is contraindicated to children under 9 and in adults over 45 years old. According to the WHO Dengvaxia® Vaccine Report (September 2018), the live attenuated dengue vaccine CYD-TDV has been shown to be effective and safe in clinical trials in people who previously had an infection with the dengue virus (HIV positive people). However, it carries an increased risk of severe dengue fever in those who experience their first natural dengue infection after vaccination (those who were HIV negative at the time of vaccination).

Dengue fever causative agents are four dengue viruses (Denv 1, Denv 2, Denv 3 and Denv 4). Dengue virus contains an 11 kb positive-sense, single-stranded RNA genome. The genome consists of a single open reading frame which encodes three structural proteins (capsid C, pre-membrane/membrane (prM/M), and envelope (E) protein), and seven non-structural proteins (NS1, NS2A, NS2B, NS3, NS4A, NS4B, and NS5) [3]. The structural proteins form the viral particle and the non-structural proteins participate in the replication of the RNA genome, virion assembly and invasion of innate immune response [4]. NS5 is the most conserved protein of the dengue proteome as it shares a minimum of 67% amino acid sequence across all four dengue serotypes [5, 6]. NS5 is essential for RNA replication and performs enzymatic activities required for capping and synthesis of RNA genome of virus. It consists of two domains with distinct functions, the N-terminal methyl transferase (MTase) and the C-terminal RNA-dependent RNA polymerase (RdRp) catalytic domain [6,7]. The tertiary structure of RdRp consists of palm thumb and finger subdomains. The catalytic site contains conserved aspartic residues.

In recent years, there has been a growing interest in Denv NS5, with several groups reporting RdRp inhibitors [8, 9, 10, 11, 12, 13, 14, 15]. A series of NS5 RdRp inhibitors recently has been reported which has led to the identification of 5-(5-(3-Hydroxyprop-1-yn-1-yl)thiophen-2-yl)-2,4-dimethoxy-N-((3-methoxyphenyl) sulfonyl) benzamide (IC₅₀^exp = 170nM 27) and 5-(5-(3-Hydroxyprop-1-yn-1-yl)thiophen-2-yl)-4-methoxy-2-methyl-N-(quinolin-8-ylsulfonyl)benzamide (IC₅₀^exp = 23nM 29) as promising avenues for further optimization and development [16]. The 3D-QSAR pharmacophores (PH4) for RdRp inhibition are not available so far, to our knowledge for these inhibitors.

The main objective of this work was to design novel potent 3-(5-ethynylthiophen-2-yl)-N-hydrosulfonylbenzamides (EHBs) based on a series of 18 (training set) and 4 (validation set) nanomolar inhibitors with observed inhibitory potencies as high as IC₅₀^exp = 23nM [16].Starting through in situ modification of the crystal structure of RdRp-EHB5 complex (PDB: 5HMZ) we have elaborated a QSAR model which correlated Gibbs free energies of RdRp-EHBx complex formation with the potencies IC₅₀^exp and determined the active conformation of EHBs bound at the active site of RdRp of Denv (MM-PB complexation approach). Based on this active conformation we have formulated 3D QSAR pharmacophore of RdRp inhibition (PH4). Large virtual library of compounds sharing the EHB scaffold has been generated and in silico screened with the PH4. The screening yielded virtual hits that exhibited predicted inhibitory potencies IC₅₀^pre more than 30 times higher than the most active training set compound EHB1. Several of the identified putative inhibitors displayed favorable ADME profiles.

2.1. Training and validation sets

Chemical structures and biological activities  (IC₅₀^exp)  of training and validation sets of 3-(5-ethynylthiophen-2-yl)-N-hydrosulfonylbenzamides inhibitors of RdRp used in this study were taken from literature [16]. The potencies of these compounds cover a sufficiently broad range of half-maximal inhibitory concentrations ( 23 ≤ IC₅₀^exp ≤ 734,000nM) to allow construction of a QSAR model. The training set (TS) containing 18 EHB inhibitors and the validation set (VS) including 4 EHBs were taken from the ref. [16].

2.2. Model building

Three dimensional (3D) molecular models of enzyme-inhibitor (E-I) complexes RdRp-EHBx, free enzyme RdRp and free inhibitors EHBx were prepared from high-resolution (1.99 Å) crystal structure of a reference complex containing the training set compound 5-(5-(3-hydroxyprop-1-yn-1-yl)thiophen-2-yl)-4-methoxy-2-methyl-N-(methylsulfonyl)benzamide (EHB5, Table 1) bound to the RdRp (Protein Data Bank [29] entry code 5HMZ [16]) using Insight-II molecular modeling program [30].

The structures of RdRp and the E-I complexes were at pH of 7 with neutral N- and C-terminal residues and all protonizable and ionizable residues charged. No crystallographic water molecules were included into the model. The inhibitors were built into the reference structure 5HMZ [16] by in situ replacing of derivatized groups in the molecular scaffold of the template inhibitor EHB5. An exhaustive conformational search over all rotatable bonds of the replacing function groups coupled with a careful gradual energy-minimization of the modified inhibitor and active site residues of the RdRp located in the vicinity of the inhibitor (within 5Å distance), was employed to identify low-energy bound conformations of the modified inhibitor. The resulting low-energy structures of the E-I complexes were then carefully refined by minimization of the whole complex. This procedure has been successfully used for model building of viral, bacterial and protozoal enzyme-inhibitor complexes and design of peptidomimetic, hydroxynaphthoic, thymidine, triclosan, pyrrolidine carboxamide, nitriles and chalcone-based inhibitors [22,31,32,33,34,35,36,37,38,39,40].

2.3. Molecular mechanics

Modeling of inhibitors, RdRp and E-I complexes was carried out by molecular mechanics as described earlier [22].

2.4. Conformational search

Free inhibitor conformations were derived from their bound conformations in the E-I complexes by gradual relaxation to the nearest local energy minimum as described earlier [22].

2.5 Solvation Gibbs free energies

The electrostatic component of solvation Gibbs free energy (GFE) that includes also the effects of ionic strength via solving nonlinear Poisson–Boltzmann equation [41, 42] was computed by the Delphi module in Discovery Studio [20] as described earlier [22].

2.6 Calculation of Binding affinity and QSAR Model

The calculation of binding affinity expressed as complexation GFE has been described fully earlier [22].

2.7 Interaction Energy

The calculation of MM interaction energy (E_int) between enzyme residues and the inhibitor CFF91 force field [42] was performed as described earlier [22].

2.8. Pharmacophore Generation

Bound conformations of inhibitors taken from the models of E-I complexes were used for constructing of 3D-QSAR pharmacophore (PH4) by means of Catalyst HypoGen algorithm [43] implemented in Discovery Studio [20] as described earlier [22].

2.9. ADME Properties

The pharmacokinetics profile of EHBs were computed by the QikProp program [26] as described earlier [22].

2.10. Virtual Library Generation

The virtual library generation was performed as described earlier [22].

2.11. ADME-Based Library Searching

The drug-likeness selection criterion served to focus the initial virtual library as described earlier [22].

2.12. Pharmacophore-Based Library Searching

The pharmacophore model (PH4) described in Section 4.8 and derived from the bound conformations of EHBs at the active site of RdRp served as library searching tool as described earlier [22].

2.13. Inhibitory Potency Prediction

The conformer with the best mapping on the PH4 pharmacophore in each cluster of the focused library subset was used for ∆∆G_com calculation and p IC₅₀^pre estimation (virtual screening) by the complexation QSAR model as described earlier [22].

3.1 Training and Validation Sets

The training set of 18 EHBs and validation set of 4 analogs (Table 1) were selected from a series of NS5 RdRp inhibitors with known experimentally determined inhibitory activities. The whole series was obtained by substitution at five positions R₁,R₂ ,R₃,R₄  and  of the phenyl ring and R group as shown in Table 1. The experimental half-maximal inhibitory concentrations (23 ≤ IC₅₀^exp ≤ 734,000 nM) [16] cover a sufficiently wide concentration range for building of a reliable QSAR model.

Table 1. Set (EHB1-18) and validation set (EHV1-4) of NS5 RdRp inhibitors [16] used in the preparation of QSAR models of inhibitor binding. The R groups are numbered as #R ≡ group index.

3.2 QSAR Model

3.2.1 One Descriptor QSAR Models

Each of the 18 training sets (TS) and 4 validation sets (VS) RdRp-EHBx complexes (Table 1), was prepared by in situ modification of the refined template crystal structure (PDB entry code 5HMZ [16]) of the complex RdRp-EHB5 as described in the Methods section. Further, the relative Gibbs free energy of the RdRp-EHBx upon complex formation (∆∆G_com) was computed for each of the 22 optimized enzyme–inhibitor complexes. Table 2 lists computed values of ∆∆G_com  and its components for the TS and VS of sulfonylbenzamides [16]. The QSAR model explained variation in the EHBs experimental inhibitory potencies (p IC₅₀^exp = _-log₁₀(IC₅₀^exp)) [16] by correlating it with computed GFE ∆∆G_com through a linear regression. In addition, significant correlation obtained in this QSAR relationship permitted to determine the active bound conformation of the EHBs at the RdRp binding site and enabled generation of the Denv RdRp inhibition PH4 pharmacophore. In search for a better insight into the binding affinity of EHBs towards Denv RdRp, we have analyzed the enthalpy of complexation in gas phase ∆∆H_MM by correlating it with the p IC₅₀^exp . The validity of this linear correlation (for statistical data of the regression see Table 3, Equation A) allowed assessment of the significance of inhibitor-enzyme interactions (∆∆H_MM) when solvent effect and loss of entropy of the inhibitor upon binding to the enzyme were neglected. This in extremis and unexpected correlation due to the non-homogeneity of molecules explained about 73% of the pIC₅₀^exp data variation and underlined the role of the enthalpic contribution to the binding affinity of the ligand. More, the advanced descriptor, namely the GFE of the RdRp-EHBx complex formation including all components: ∆∆H_MM, ∆∆TS_vib, and ∆∆G_sol, has been assessed (for statistical data see Table 3, Equation B). Relatively high values of the regression coefficient R², leave-one-out cross-validated regression coefficient R²_xv and Fischer F-test of the correlation show the importance of the term entropic in the biological environment and suggest a good relationship between the 3D model of inhibitor binding and the observed inhibitory potencies of the EHBx [16] .Therefore, structural information derived from the 3D models of RdRp – EHBx complexes is expected to lead to reliable prediction of RdRp inhibitory potencies for novel EHBs analogs based on the QSAR model B, Table 3.

Table 2. Gibbs free energy (binding affinity) and its components for the training set of RdRp inhibitors EHB1-18 and validation set inhibitors EHV1-4 [16].

^a  for the chemical structures of the training set of inhibitors see Table 1.

^b  is the molar mass of inhibitors.

^c ∆∆H_MMis the relative enthalpic contribution to the Gibbs free energy change related to E:I complex formation derived by molecular mechanics (MM): ∆∆H_MM≅ [E_MM {E:I_x} - [E_MM{I_x}]-[E_MM{E:I_ref} - E_MM{I_ref}], I_ref is the reference inhibitor EHB1;

^d ∆∆G_sol is the relative solvation Gibbs free energy contribution to the Gibbs free energy change related to E:I complex formation: ∆∆G_sol =[G_sol{E:I_x} - G_sol{I_x}] - [G_sol{E:I_ref} - G_sol{I_ref}]

^e ∆∆TS_vib is the relative entropic contribution of the inhibitor to the Gibbs free energy related to E:I complex formation: ∆∆TS_vib = [∆∆TS_vib{I_x}_E - ∆∆TS_vib{I_x}] - [∆∆TS_vib{I_ref}_E - ∆∆TS_vib{I_ref}];

^f ∆∆G_com is the relative Gibbs free energy change related to E:I complex formation: ∆∆G_com ≅ ∆∆H_MM + ∆∆G_sol + ∆∆TS_vib

^g IC₅₀^exp is the experimental RdRp half maximal inhibition concentration obtained from reference [16].

^h Ratio of predicted and experimental half maximal inhibition concentrations pIC₅₀^pre /IC₅₀^exp . pIC₅₀^pre = _-log₁₀(IC₅₀^pre) was predicted from computed ∆∆G_com using the regression equation for RdRp shown in Table 3, B.

Table 3. Analysis of computed binding affinities ∆∆G_com, its enthalpic component ∆∆H_MM and experimental half-maximal inhibitory concentrations pIC₅₀^exp = _-log₁₀(IC₅₀^exp) of EHBs towards Denv RdRp [16].

The statistical data confirmed validity of the correlation Equations (A) and (B) plotted on Figure 1. The ratio pIC₅₀^pre/IC₅₀^exp ≅ 1 (the  values were estimated using correlation Equation B, Table 3) calculated for the validation set EHV1-4 documents the substantial predictive power of the complexation QSAR model from Table 2. Thus, the regression Equation B (Table 3) and computed GFE ∆∆G_com can be used for prediction of inhibitory potencies IC₅₀^pre against Denv RdRp for novel EHB analogs provided they share the same binding mode as the training set sulfonylbenzamides EHB1-18.

Figure 1. (A) Plot of correlation equation between IC₅₀^exp and relative enthalpic contribution to the GFE ∆∆H_MM [kcal·mol⁻¹]. (B) Similar plot for relative complexation Gibbs free energies of the RdRp-EHBx complex formation ∆∆G_com [kcal·mol⁻¹] of the training set [16]. The validation set data points are shown in red color.

3.2.2 Binding Mode of EHBs

The 3D crystal structure of Denv RdRp adopts a classical polymerase hand shape with fingers, palm, and thumb subdomains [17]. This last part contains the initiation loop which triggers the polymerization process of the viral RNA close to the palm with the catalytic residues ASP 663 and ASP 664 [6,18,19]. The recently reported x-rays complex of Denv RdRp and non-nucleoside inhibitors reveals that the propargyl alcohol projected into the narrow cavity and formed two Hbond interactions with His 800 and Glu 802. The sulfonylbenzamide is involved in three Hbond contacts with the side chains of Thr 794 and Arg 729 and the backbone of Trp 795. Changes made to this methyl have improved the activity but limited still to the two digits nanomolar range IC₅₀^exp =23nM)[16]

Figure 2. (A) 3D structure of the Denv RdRp active site with bound inhibitor EHB5 (5HMZ.pdb). (B) 2D schematic interaction diagram of the inhibitor EHB5 [16] at the active site of Denv RdRp. (C) Hydrophobic surface of the active site of Denv RdRp with the most potent inhibitor EHB1 [16]. Surface coloring legend: red = hydrophobic, blue = hydrophilic and white = intermediate.

3.3. Interaction Energy

Other key structural information was provided by the interaction energy (IE, ) diagram obtained for each training set inhibitor. IE breakdown to contributions from Denv RdRp active site residue helpfully directs the choice of relevant R-groups able to improve the binding affinity of EHB analogs to the Denv RdRp and subsequently enhance the inhibitory potency. A comparative analysis of computed IE for training set EHBs (Figure 4) divided into three classes (highest: 23 – 340 nM, moderate: 2400 – 7500 nM, and lowest activity: 26,000 – 734,000 nM) has been carried out to identify the residues for which the contribution to binding affinity could be increased. The comparative analysis showed IE contributions of active site residues for the three classes of inhibitors that should be retained or even improved such as those of Leu 511, His 711, Arg 737 and Thr 794. However, interactions with residues such as Met 340, Glu 733 Met 765 and Gln 802 are accentuated from the low activity class to the high activity class via that of medium activity. It should be noted that these residues Met 340, Glu 733 and Met 765 had not been listed as belonging to the active site (PDB:5HMZ) [16]. Since specific substitutions could not be proposed, we have adopted a combinatorial approach to novel EHB analogs design and in silico screened a virtual library of EHB analogs with help of the PH4 pharmacophore of Denv RdRp inhibition derived from the complexation QSAR model.

The statistical data confirmed validity of the correlation Equation plotted on Figure 3. This correlation of 88 % shows that there are interactions to be made to maintain or improve the activity of the new analogs.

Table 4. Analysis of computed binding affinities ∆∆E_int and experimental half-maximal inhibitory concentrations pIC₅₀^exp = -log₁₀(IC₅₀^exp) of EHBs towards Denv RdRp [16].

Figure 3. Plot of correlation equation between pIC₅₀^exp and relative interaction energies on the active site of the RdRp-EHBs complex.

Figure 4. Mechanics intermolecular interaction energy E_int breakdown to residue contributions in [kcal.mol^-1]:(A) the most active inhibitors EHB1-5, (B) moderately active inhibitors EHB6-9, (C) a few of less active inhibitors EHB10-18, Table 2 [16].

3.4. 3D-QSAR Pharmacophore Model

RdRp inhibition 3D-QSAR pharmacophore was generated from the active conformation of 18 TS EHB1-18 and evaluated by 4 VS EHV1-4 covering a large range of experimental activity (23 - 734000 nM) spanning more than two orders of magnitude. The generation process is divided into three main steps: (i) the constructive step, (ii) the subtractive step, and (iii) the optimization step [20]. During the constructive phase, EHB1 alone was retained as the lead (since only the activity of EHB1 fulfilled the threshold criterion: IC₅₀^exp≤ 1.25 x 23 nM) and used to generate the starting PH4 features. In the subtractive phase, compounds for which: IC₅₀^exp >23×10^3.5 nM=72,732 nM were considered inactive. As a result, EHB14, EHB15, EHB16, EHB17 and EHB18 from the EHBx training set were inactive. Finally, during the optimization phase, the score for pharmacophoric hypotheses was improved. The assumptions were scored based on errors in the regression activity and complexity estimates via a simulated annealing approach. At the end of the optimization, the 10 highest-rated unique pharmacophore hypotheses were retained, all showing five-point characteristics. The cost values, correlation coefficients, root mean square deviation (RMSD) values, pharmacophore characteristics and max-fit value of the first 10 ranked hypotheses (Hypo1 - Hypo10) are listed in Table 5. They were selected based on statistically significant parameters, such as high correlation coefficient, low total cost and low RMSD.

Table 5. Parameters of 10 generated PH4 pharmacophoric hypotheses for RdRp inhibitors after CatScramble validation procedure (49 scrambled runs for each hypothesis at the selected level of confidence of 98%).

^a Root Mean Square Deviation; ^b squared correlation coefficient; ^c overall cost parameter of the PH4 pharmacophore; ^d cost difference between Null cost and hypothesis total cost; ^e lowest cost from 49 scrambled runs at a selected level of confidence of 98%. The Fixed Cost = 44.55 with RMSD = 0, the Null Cost = 1583.2 with RMSD = 13.129 and the Configuration cost = 11.57.

The generated pharmacophore models were then assessed for their reliability based on the calculated cost parameters ranging from 93.16 (Hypo1) to 143.55 (Hypo10). The relatively small gap between the highest and lowest cost parameter corresponds well with the homogeneity of the generated hypotheses and consistency of the TS of EHBx. For this PH4 model, the fixed cost (44.55) is lower than the null cost (1583.2) by a difference ∆ = 1 538.65. This difference is a major quality indicator of the PH4 predictability (∆ > 70 corresponds to an excellent chance or a probability higher than 90% that the model represents a true correlation [20]). To be statistically significant, a hypothesis must be as close as possible to the fixed cost and as far as possible from the null cost. For the set of 10 hypotheses, the difference ∆ ≥ 1 439.65, which attests to the high quality of the pharmacophore model. The standard indicators such as the RMSD between the hypotheses ranged from 2.31 to 3.30, and the squared correlation coefficient (R²) falls to an interval from 0.984 to 0.968. The first PH4 hypothesis with the closest cost (93.16) to the fixed one (44.55) and best RMSD and  was retained for further analysis. The statistical data for the set of hypotheses (costs, RMSD, ) are listed in Table 5. The configuration cost (11.57 for all hypotheses) far below 17 confirms this pharmacophore as a reasonable one.

The link between the 98% significance and the number 49 scrambled runs of each hypothesis is based on the formula S=[1-((1+X))⁄Y]×100, with  the total number of hypotheses having a total cost lower than the original hypothesis (Hypo 1) and  the total number of HypoGen runs (initial + random runs):X=0 and Y=(1+49)  and , hence 98% = {1-[(1+0)/49+1]} x  100.

The evaluation of Hypo 1 was performed first through Fischer’s randomization cross-validation test. The CatScramble program was used to randomize the experimental activities of the training set. At 98% confidence level, each of the 49 scramble runs created ten valid hypotheses, using the same features and parameters as in the generation of the original 10 pharmacophore hypotheses. Among them, the cost value of Hypo1 is the lowest compared with those of the 49 randomly generated hypotheses, as we can see in Table 5 where the lowest cost of the 49 random runs is listed for each original hypothesis, and none of them was as predictive as the original hypotheses generated shown in Table 5. Thus, there is a 98% probability that the best selected hypothesis Hypo1 represents a pharmacophore model for inhibitory activity of RdRp with a similar level of predictive power as the complexation QSAR model, which relies on the EHBx active conformation from 3D structures of the RdRp-EHBx complexes and computed GFE of enzyme–inhibitor binding ∆∆G_com. Another evaluation of Hypo 1 is the mapping of the best active training set EHB1 (Figure 5) displaying the geometry of the Hypo1 pharmacophore of NS5RdRp inhibition. The regression equation for pIC₅₀^exp vs pIC₅₀^pre estimated from Hypo1: pIC₅₀^exp = 1.0001×pIC₅₀^pre - 0.0017(n = 18, R² = 0.969,R²_xv =  0.967, F-test = 498.83, σ =  0,238, α > 98 %) is also plotted on Figure 5.

Figure 5. Features (A) coordinates of centers, (B) angles between centers of pharmacophoric features, (C) distances between centers, (D) mapping of pharmacophore of RdRp inhibitor with the most potent molecule EHB1. Features legend: HBD = Hydrogen bond Donor (magenta), HYDAr = Hydrophobic Aromatic (cyan), HBA = Hydrogen bond Acceptor (green). (E) Correlation plot of experimental vs. predicted inhibitory activity (open circles correspond to TS, red dots to VS).

We can carry out computational design and selection of new EHB analogs with elevated inhibitory potencies against Denv RdRp, based on a strategy using the noticeable presence of the hydrophobic features included in the best pharmacophore model at the position of  coupled with mapping of  to the HBD feature and the appropriate substitution to the others hydrophobic features in Hypo1 (Figure 5).

3.5.  Virtual Screening

In silico screening of a virtual (combinatorial) library can lead to hit identification as it was shown in our previous works on inhibitors design [20,21,22].

3.5.1  Virtual Library

An initial virtual library (VL) was generated by substitutions at positions for R₁,R₂,R₃,R₄ and R₅ (see Table 6) on the phenyl ring scaffold. During the virtual library enumeration, the 221 R-groups listed in Table 6 were attached on following way: 1-7, 11, 18 to position ; 4, 8-20 to positions R₂ and R₃ ; 2, 4, 8-90 to position R₄ then all 221 R-groups to positions R₅ of the phenyl ring. The combinatorial library size is R₁ x R₂ x R₃ x R₄ x R₅ =9×14×14×85×221=33,136,740 analogs. To design a more focused library of a reduced size and increased content of drug-like molecules, we have introduced a set of filters and penalties such as the Lipinski rule-of-five [23], which helped to select a smaller number of suitable EHBs that could be submitted to in silico screening. This focusing has reduced the size of the initial library to 125,915 analogs.

Table 6.  R₁ to R₅-groups (fragments, building blocks, substituents) used in the design of the initial diversity virtual combinatorial library of sulfonylbenzamides.

R₁-groups : fragments 1-7, 11, 18; R₂ and R₃-groups : fragments 4 ; 8 - 20. R₄-groups: fragments 2, 4, 8-20 ; R₅-groups: fragments 1-221.

3.5.2. In Silico Screening of Library of EHBs

The focused library of 125,915 analogs was further screened for molecular structures matching the 3D-QSAR PH4 pharmacophore model Hypo1 of RdRp inhibition. 329 EHBs mapped to at least 4 features of the pharmacophore. These best fitting analogs (PH4 hits) then underwent complexation QSAR model screening. The computed GFE of RdRp-EHBx complex formation, their components, and predicted half-maximal inhibitory concentrations IC₅₀^exp calculated from the correlation Equation B (Table 3) are listed in Table 7.

Table 7. GFE and their components for the top scoring 304 virtual EHB analogs. The analog numbering concatenates the index of each substituent  to  with the substituent numbers taken from Table 6.

^a Mw is molar mass of inhibitor.^b ∆∆H_MM is the relative enthalpic contribution to the GFE change of the NS5RdRp-EHB complex formation ∆∆G_com (for details see footnote of Table 2); ^c ∆∆G_sol is the relative solvation GFE contribution to ∆∆G_com;  ^d ∆∆TS_vib is the relative (vibrational) entropic contribution to ∆∆G_com; ^e ∆∆G_com is the relative Gibbs free energy (GFE) change related to the enzyme–inhibitor NS5RdRp-EHB complex formation ∆∆G_com = ∆∆H_MM + ∆∆G_sol - ∆∆TS_vib; ^f  IC₅₀^pre is the predicted inhibition potency towards NS5RdRp; ^g Experimental value IC₅₀^exp is given for the reference inhibitor EHB1 instead of the predicted value.

3.6 Novel EHB Analogs

The design of virtual library of novel analogs was guided by structural information retrieved from the EHBx active conformation and was used for the selection of appropriate substituents (R₁ - to R₂-groups). In order to identify which substituents lead to new inhibitor candidates with the highest predicted potencies towards the RdRp of Denv, we have prepared histograms of the frequency of occurrence of R₁- to R₅-groups among the 329 best fit PH4 hits (Figure 6). The histograms show that the R₁-groups 2, 1 and 7 were represented with the highest frequencies of occurrence (93), (86) and (59) among the 329 EHB hits. The R2-groups and R₃-groups most frequently represented in this subset are 17 (99) and (110); 9 (45) and (47); 11 (36). As for the R₄-groups the highest frequencies concern 17 (54); 9 (35); 11 and 18 (32). The R₅-groups are dominated by substituents 113 (54); 221 (33). The top five scoring virtual hits namely analogs are: 2-8-9-8-218 (IC₅₀^pre = 0.21 nM) , 2-8-9-8-219 (IC₅₀^pre = 0.31 nM), 2-8-10-18-108  (IC₅₀^pre = 0.58 nM), 2-8-9-8-121 (IC₅₀^pre = 0.65 nM)  and 2-8-9-8-220 (IC₅₀^pre = 0.72 nM). They include the following substituents at  position: 2: 5-(3-aminoprop-1-yn-1-yl)thiophen-2-yl (5), at  position: 9 Me (4). Despite the hydrophobicity of the pocket exploited by the  -groups, their orientation to invest it seems to be dictated by other R-groups.

Figure 6. Histograms of frequency of occurrence of individual R-groups in the 329 best selected analogs mapping to four features of the PH4 pharmacophore hypothesis Hypo1 (for the structures of the fragments see Table 6); R₁ = 5-(3-aminoprop-1-yn-1-yl)thiophen-2-yl (2); R₂= MeO (8); R₃=Me (9), Cl (10); R₄=MeO (8), NH2(18) and ((3,3-diMebutyl)sulfonyl)carbamoyl (108), ((5-(Me-amino)pentyl)sulfonyl)carbamoyl (121),(R)-((5-(Etamino)hexyl)sulfonyl)carbamoyl(218),(R)-((5-(isopropylamino)hexyl)sulfonyl)carbamoyl(219),(3R,5R)-((5-(Etamino)-3-Mehexyl)sulfonyl)carbamoyl (220).

The substitutions in  to  positions of EHBs led to an overall increase of affinity of RdRp binding as exemplified by the inhibitory potencies of majority of new designed analogs. The best designed sulfonylbenzamide EHB 2-8-9-8-218 displays predicted half-minimal inhibitory concentration of IC₅₀^pre =0.21 nM that is more than 100-times lower than that of the most active compound of the TS, namely the EHB1 with IC₅₀^exp = 23 nM , Figure 7,8.

Figure 7. (A)-Close up of virtual hit 2-8-9-8-218, the most active designed EHB analog (IC₅₀^pre = 0.21nM) at the active site of RdRp. Interacting residues are coloured by element. (B) - mapping of the EHB 2-8-9-8-218 to RdRp inhibition pharmacophore. (C) - 2D schematic interaction diagram of the EHB 2-8-9-8-218 at the active site of Denv RdRp. (D) - 2D schematic interaction diagram of the analog EHB2-8-9-8-219 (IC₅₀^pre = 0.31nM) at the active site of Denv RdRp. (E) - 2D schematic interaction diagram of the ligand EHB1 at the active site of Denv RdRp.

Figure 8. Surface of the active site of Denv RdRp with bound 5 best active designed EHB analogs. The binding site surface is colored according to residue hydrophobicity: red = hydrophobic, blue = hydrophilic and white = intermediate.

3.7 Pharmacokinetic Profile of Novel EHB Analogs

Many antivirals used in the treatment of viral infections exhibit a pharmacokinetic profile unsatisfactory with poor bioavailability and/or short half-life increasing the risk of ineffectiveness. Improving the pharmacokinetic characteristics of a drug requires improving its absorption profile. Improving the pharmacokinetic profile, i.e. increasing systemic exposure, increases the effectiveness of antiviral drugs, but also improves the quality of life of treated patients [24]. Although we do not have a specific antiviral against dengue, it seems advisable to compare, in Table 8, the properties linked to the ADME of our molecules with that of a known drug and sometimes used in the treatment of the dengue. Octanol-water partition coefficient, aqueous solubility, blood brain partition coefficient, Caco-2 cell permeability, serum protein binding, number of probable metabolic reactions and eighteen other descriptors related to absorption, the distribution, metabolism and excretion (ADME) of the new analogues were calculated by the QikProp program [25] based on the Jorgensen method [26,27]. Experimental data from over 710 compounds were used to produce regression equations correlating experimental and calculated descriptors resulting in an accurate prediction of the pharmacokinetic properties of the molecules. Since a value of less than 25% is considered poor, on the contrary, the 10 best predicted analogues show human oral absorption from the gastrointestinal tract (HOA) ranging from 48% to 83%, well above that of Ribavirin [28]. Drug likeness (#stars) - the number of property descriptors that fall outside the range of optimal values determined for 95%. The values of the best designed active EHBs are compared to those calculated for Ribavirin, Table 8. Most of our best designed analogs have #stars less than or equal to one. Thus, the designed EHBs exhibit a favorable pharmacokinetic profile.

Table 8. ADME-related properties of the best designed EHB analogs and known medicate either in clinical use to treat dengue or currently undergoing clinical testing computed by QikProp [26].