MOLECULAR MODELLING AND VIRTUAL SCREENING APPLICATION TO THE COMPUTER-AIDED DESIGN OF ANTICANCER INHIBITORS WITH A FAVORABLE PHARMACOKINETIC PROFILE AGAINST E6 PAPILLOMAVIRUS TYPE 16

Ingrid Clémence Bléhoué¹, Mawa Koné^2,3, Elvice Akori Esmel^l, Issouf Fofana¹, 

Mélalie Kéïta¹, Eugene Megnassan^*1,2,4

¹Laboratory of Fundamental and Applied Physics, University of Abobo-Adjamé (Now Nangui Abrogoua), Côte d’Ivoire.

²Laboratory of Constitution and Reaction of Matter, University of Cocody (Now Felix Houphouët Boigny), Côte d’Ivoire.

³National Laboratory for Quality Testing, Metrology and Analysis, BP V 174 Abidjan, Côte d’Ivoire.

⁴International Centre for Theoretical Physics (ICTP), Strada Costiera, Trieste, Italy.

Abstract

Background: Worldwide, cervical cancer (CC) is the fourth most common women cancer. It is crucial to develop more effective treatments for this disease. We aim at designing new anticancer compounds with a favorable pharmacokinetic profile, targeting the E6 oncoprotein of human papillomavirus type 16 (HPV16 E6).

Methods: Computer-Aided Molecular Design (CAMD) has been carried out by the elaboration of a Quantitative Structure-Activity Relationship (QSAR) model of molecular complexation, through the correlation between the relative Gibbs free energy (rGFE) and the observed biological activities of a series of flavonoids inhibitors separated in training (TS) and validation sets (VS). Starting from the 3D crystallographic structure of the oncoprotein HPV16 E6 (Protein Data Banck (PDB) input code: 4GIZ), enzyme – inhibitor complexes were built by in situ modification of the native ligand (FLAV1, IC₅₀^exp= 850 nM), replacing substitution R-groups at nine different positions R₁ – R₉.

Results: The complexation QSAR model equation (pIC₅₀^exp=-0.5494 × DDG_com+5.9983 (1); n=16; R²=0.98) explains 98% of the variation in biological activity by that of rGFE of the analogues used. The subsequent 3D pharmacophore model (PH4) generated from the active conformations of FLAVS (pIC₅₀^exp= 1.0177× pIC₅₀^pre – 0.0927(2); n=16; R²=0.90) was used as a virtual selection tool to identify new analogues from a virtual library (VL) reaching two digits nanomolar range predicted activity. 

Conclusions: The combination of molecular modelling and in silico screening of VL using the PH4 pharmacophore has led to the discovery of new promising anticancer candidates, with favorable pharmacokinetic profiles against HPV16 E6. Among them, the top two predicted respective inhibitory powers IC₅₀ ^pre (50 nM and 61 nM).

Keywords: Cervical cancer, flavonoids, HPV16 E6 oncoprotein, prediction of ADME properties, Quantitative Structure-Activity Relationship (QSAR).

INTRODUCTION

In 2022, cervical cancer was the fourth most common cancer among women worldwide, with approximately 660.000 new cases and around 350.000 deaths¹. The main cause of cervical cancer is persistent infection with human papillomaviruses (HPV), particularly the high-risk (HR) genotypes 16 and 18, responsible for about 50% and 20% of cases respectively^2,3. High-risk HPV infections are a major global medical issue, accounting for 99% of cervical cancer and contributing to approximately 5% of all cancers worldwide^4,5. HPV-induced keratinocyte transformation is a complex process, marked primarily by the intracellular accumulation of two key viral oncoproteins, E6 and E7, playing a key role in disrupting cellular control mechanisms, thus promoting progression to malignancy³. The oncoprotein E6 of papillomavirus of type 16 (HPV16 E6) plays a crucial role in cervical oncogenesis^4,5. It promotes the degradation of the p53 protein, a key regulator of the cell cycle and apoptosis, through the E6 ubiquitin ligase (E6AP)⁶. The inhibition of this destructive interaction is therefore a promising strategy for the treatment of cervical cancer. HPV vaccines have no therapeutic potential with limited access worldwide making design of new treatments a crucial need⁷. Recently, several innovative approaches have been explored E6 oncoprotein inhibition including the use of RNA interference⁸ to reduce viral protein expression. On the other hand, inhibitoric peptides⁹ have been developed against E6 and p53. Finally, small molecules such as luteolin¹⁰ were identified as potential candidates, due to their ability to restore the tumor suppressive activity of p53.The availability of the E6/E6AP crystal structure in the Protein Data Bank (PDB)¹¹ opened the gate to precisely identify amino acid residues involved in critical interactions with E6. Among them, residues in the LxxLL binding pocket are particularly important because their key role in the recognition of E6 protein binding partners. Natural products (NP) are drawing interest as anti-cancer agents, particularly due to their accessibility, affordability, and generally lower toxicity compared to other synthetic molecules. Among them, flavonoids stand out particularly, having demonstrated various biological effects, such as antibacterial, antiviral and antioxidant ones. Furthermore, these compounds are known to be able to induce apoptosis in cancer cells, making them promising candidates for the development of anti-cancer treatments based on NP¹².

In this context, Srikanth Kolluru and al.¹³, reported flavonols, with the most active, FLAV1, exhibits anIC₅₀^exp of 850 nM, and the x-rays crystallography structure of “HPV16 E6 – FLAV1”complex (PDB entry code 4GIZ)¹³. Analysis of the interactions from this seminal study and literature exploration led to a series of enough compounds with the same scaffold in order to build QSAR model of HPV16 E6 inhibition.

Therefore, CAMD is suitable to release novel FLAVs analogs with better inhibitory potency. To achieve our objective, we first developed a QSAR model that establishes a correlation between the free energy of Gibbs released during the formation of “HPV16 E6-FLAVx” complexes and their corresponding experimental inhibition activities. Next, we designed a 3D-QSAR pharmacophore model (PH4) for HPV16 E6 inhibition, based on the active conformations obtained through the complexation method. Subsequently, we generated and screened a virtual library using the PH4 model. Finally, we evaluated the predicted activity of the most promising PH4-identified analogs and calculated their ADME profile. In a last step the stability of the best analogs has been checked by molecular dynamics (MD) simulations. The expected results could provide new therapeutic leads in the design of novel small molecule inhibitors for the treatment of E6-targeted cancers.

METHODS

In this work, the chemical structures and biological activities studied belong to the flavonoid family. These structures were divided into two sets: one for training set (TS) and the other for validation set (VS).

Training and validation sets

Flavonoids belong to the large family of polyphenols. They share a common basic structure consisting of two benzene rings connected by a linear three-carbon chain, forming an oxygenated heterocycle¹⁴.

There are a total of 14 flavones and 5 flavonols, with experimental inhibition concentrations taken from the literature respectively^10,13. There are a total of nineteen (19) molecules with an experimental concentration of HPV16 E6 oncoprotein inhibition ranging from 850 to 671000 nM respectively, wide sufficiently to allow the construction of a reliable QSAR model. They were distributed as follows:

-  16 compounds for the trainingset (TS) and;

-  3 compounds for model validation (VS)set.

This distribution was done using a protocol called “Generate Training and Test Data” in the Discovery Studio 2.5 software¹⁵.

Model building

The construction of the QSAR model begins with the download of the crystallographic structure of the HPV16 E6–FLAV1 complex (PDB code:4GIZ). Indeed, the Protein – Inhibitor (IP) complexes were made from the crystallographic structure of HPV16 E6 –FLAV1 using the Molecular Modelling program Insight II¹⁶. No water molecule is included in the model. The “HPV16 E6 – FLAVx” complexes were built by in situ modifications on the reference structure of “HPV16 E6 – FLAV1” complex with FLAV1 inside the receptor active site by replacing atom by atom , the appropriate fragments of the molecular structure of FLAV1. During these in situ modifications, the addition of an atom requires a local minimization within a radius of 5 Å around it. When adding a group of atoms, a systematic conformational search around the torsion angles (torsion drive), is necessary, followed by a global minimization of the protein-ligand complex. This procedure has a double advantage to avoid explosive steric contacts and take into account the reality (flexibility of the lateral chains of the receptor active site residues).

The complete description of the computation of relative ligand binding affinity (ΔΔG_com) has been reported earlier¹⁷.

ΔΔG_com= ΔG_com(I) − ΔG_com(I_ref)=ΔΔH_MM– ΔΔTS_vib+ ΔΔG_sol                                            (1).

ΔΔH_MM describes the relative enthalpic contribution to the GFE change related to the intermolecular interactions in E:I complex derived by molecular mechanics (MM), ΔΔG_sol and ΔΔTS_vib represent, respectively, the relative solvation rGFE and simplified relative vibrational entropy.

Molecular mechanics

The FLAVs, HPV16 E6 –FLAVx complexes and HPV16 E6 have been refined by Molecular Mechanics (MM) using methodology previously described earlier¹⁸.

Solvation Gibbs free energies

The electrostatic component of solvation-free energy (GFE), which incorporates the effects of ionic force through the resolution of the nonlinear Poisson-Boltzmann (PB) equation¹⁹, was calculated by the Delphi module in the Discovery Studio¹⁵ suite as previously detailed¹⁸.

Molecular complexation QSAR

First, the thermodynamic quantities ∆∆H_MM, ∆∆G_sol, T∆∆S_vib, and ∆∆G_com were calculated for each ligand in the test set. The ligand FLAV1 being the most active of the set (IC₅₀^exp= 850nM) was chosen as the reference ligand. The values of the variation of the relative free enthalpy of complexation are obtained by the equation: ∆∆G_com= ∆G_com (I) - ∆G_com (I_ref) 

              =∆∆H_MM − T∆∆S_vib + ∆∆G_sol (3) with:

- ∆∆H_MM: relative variation of the complexation enthalpy for each inhibitor;

-               T∆∆S_vib:relative variation of the complexation entropy;

-               ∆∆G_solv: relative variation of complexation solvation enthalpy.

Then, the experimental biological activities IC₅₀^exp were correlated to the variation of the Gcom complexation enthalpy (model in solvent) relative to the reference ligand. Finally, to validate this 3D-QSAR model of molecular complexation, the statistical indicators of reliability of linear regression were calculated.

Interaction energy

The molecular mechanics (E_int) interaction energy calculation protocol available in Discovery Studio 2.5¹⁵ was used to determine unrelated interactions, including van der Waals and electrostatic interatomic potential terms, between two groups of atoms belonging one to P and the other to I in the P:I complexes to reach E_int breakdown into contributions of residues from individual as fully described earlier¹⁸.

Pharmacophore generation

The linked conformations of inhibitors, extracted from P:I complex models were used to build the 3D QSAR pharmacophore using the HypoGen algorithm of Catalyst²⁰, integrated in Discovery Studio¹⁵. The generation process is carried out in three main stages using the set of most active inhibitors: (i) the construction phase, (ii) the subtraction phase and (iii) the optimization phase, as described in the work of Kouassi et al.¹⁷. Inactive molecules were used to define the volume excluded. Five HypoGen features were selected: hydrogen binding acceptor (HBA), hydrophobic, negative ionizable (Neg-Ionizable), heavy hydrogen binding acceptor (HBA heavy) and projection hydrogen bonding acceptor (HBA projection). The majority of protocol parameters are maintained at their default value except of the biological activity uncertainty set to 1.5. The uncertainty is a value between 1 and 3. This adjustment has reduced the uncertainty interval of experimental biological activity from a wide range to a relatively narrow range  due to the accuracy and homogeneity of the values of the available experimental biological activities^10,13. The top ten pharmacophores were generated with a missing number of functionalities to the PH4 hypothesis set at 0. No new conformers are created when developing the PH4 model. The active site conformation is preserved during PH4 creation and superposition. Finally, the best pharmacophore model is selected.

Virtual library generation

Due to the high number of substituents (R₁to R₉) on the scaffold, the generation of the virtual library (VL) was done using the molecular modeling and simulation software MOE 2015²¹ (Molecular Operating Environment). This would require first to install and configure MOE 2015²² so that all the functionalities are accessible and the necessary licenses are activated. Then, prepare the molecular structures, which consists of importing structures (Load your molecules via files of type SDF, MOL, PDB or other supported formats). Molecules can be imported from an external database or by drawing them directly in MOE²². In this work, we imported the common basic structure of inhibitors in SD format that we converted to a SDF file type, taken into account by MOE and then we redefined the order of appearance of the different R_Groups on the scaffold, we also specified their attachment point using the MOE²² “builder” menu and then draw the fragments. Finally, we used the option “Databases>Build Database” to generate a database from the imported scaffold and fragments drawn.

ADME properties

The ADME properties of inhibitors were calculated using the QikProp^[i]program based on Jorgensen’s methods^2,3,4 as detailed previously¹⁸.

Pharmacophore-based library searching

The pharmacophore model (PH4) mentioned above was derived from the conformations of the FLAVs linked to the active site of HPV16 E6. The virtual library was then scanned using the "pharmacophore" mapping protocol available in Discovery Studio 2.5¹⁵. as described earlier¹⁸.

Inhibitory potency prediction

The conformer with the best match for PH4 pharmacophore in each group of the target library subset was selected for in silico screening by complexation QSAR model. The relative GFE of the formation of the P:I complex in water ( DDG_com) was calculated for each new analogue selected and then used to predict their inhibitory power of HPV16 E6.

Molecular Dynamics Simulations

In order to check 3D P:I complexes stability, molecular dynamics (MD) simulations to evaluate the stability of complexes formed between HPV16 E6 and flavonoid inhibitors (FLAVs), as well as the conformation flexibility of the active ligand FLAV1 and its three new potentials analogues in silico at the active site of HPV16 E6. The complexes obtained by in situ modification of the reference inhibitor FLAV1, followed by refinement through molecular mechanics (MM) methods, were used as starting geometries for the molecular dynamics simulations performed with the software Desmond²⁸_. Table 8 shows the averages of total energy (E_tot) and potential energy (E_pot) of the systems studied during the 2000time step simulation. Figure 9 shows the temporal evolution of the properties of the linked inhibitor, such as mean quadratic deviation (RMSD) from initial conformation, the radius of gyration (rGyr), the number of intramolecular hydrogen bonds (intraHB), the molecular surface (molSA), the solvent-accessible surface (SASA)and the polar surface (PSA). The interactions between the protein and ligand were examined throughout the simulation trajectory to identify specific interactions that were maintained during the calculation (Figure 10). Finally, interactions that occur more than 20% of the time during the simulation (from 0 to 200 ns) are illustrated in a detailed 2D diagram (Figure 11).

Furthermore, we superimposed the conformations of the ligands obtained after minimization of the complexes from the molecular dynamics simulations with those derived from the in situ modification and refinement of FLAV1 by molecular mechanics. 

RESULTS

Training and validation sets

All sixteen (16) training set compounds and all three (3) validation set compounds, including their experimental inhibitory activities, were retrieved from the literature and they are listed in Table 1^10,13. 

Their IC₅₀^exp cover a relatively wide range (about 10^2.9) 850 nM ≤IC₅₀^exp≤671000 nM sufficient to allow the construction of a valid QSAR model.

One-descriptor QSAR model

The relative Gibbs free energy (rGFE) of the P:I complex formation (ΔΔG_com) was calculated for the complexes HPV16 E6: FLAVs. The 16 formation complexes (TS) and V complexes (VS) were prepared respectively by in situ modification of the crystal structure of the refined model (input code PDB 4GIZ¹³ of HPV16/E6APin complex with FLAV1) as described in the Methods section.  Table 2 lists the values of the rGFE (ΔΔG_com). The QSAR model obtained highlights the correlation between the experimental inhibition potency of flavonoids and the variations in relation to a reference of the rGFE (ΔΔG_com) calculated by linear regression. The statistical data of linear regression are given in Table3. ΔΔG_com reflects the mutual binding affinity between the protein and the inhibitor. The resulting QSAR model of molecular complexation explains about 98% of the variation in biological activity by ΔΔG_com (Figure 3).

Binding mode of FLAVs

The analysis of interactions of HPV16 E6 − FLAV1 and HPV16 E6 − FLAV16 complexes corroborates the experimental inhibitory activities obtained by Srikanth Kolluru et al¹³. Indeed, Figure 3 below displays the binding mode of the most active ligand (FLAV1) and least active (FLAV16) in the active site of HPV16 E6. We can see Arg131 in two hydrogen bonds (HB) contact with FLAV1 instead of one with FLAV16. The oxygen atom binding domain at R₈ of the FLAV1 is in HB with Ala61and Tyr60 against only one with Tyr60 for FLAV16. FLAV1 is in one HB with the catalytic residue Cys51 whereas its loss with FLAV16. 

These specific HB contacts reported by Kolluru et al.¹³, were observed for the new class of HPV16 E6 inhibitors, except the interaction between the inhibitor and Tyr 60 which was not reported prior to the work of Kolluru et al.¹³, and is observed in this work. The superiority of interactions, notably of FLAV1’s hydrogen with residues from the active site of HPV16 E6 as opposed to FLAV16, could explain the important stability of FLAV1 within the active site of the E6 protein compared to other ligands in the series.

Interaction Energy

The distribution of contributions from HPV16 E6 active site residues is useful to select substituents (R groups) able to improve the binding affinity of FLAV analogues with HPV16 E6, and thus increase their inhibitory power. Indeed, in the SAR (Structure-Activity Relationship) study by Jonathan J., Cherry et al.¹⁰. On the training setflavones it was reported that the lack of the hydroxyl on the scaffold ring A and substitution on benzene or a heterocycle in ortho position restored the activity which jumped from 23 000nM (FLAV9) to 1100 nM (FLAV2).

Flavones with those structural modifications are among the most active (1100 nM ≤IC₅₀^exp≤12500 nM) of the TS compounds; others except for FLAV1 (more active ligand) are largely part of the moderately active ligands 23000 nM ≤IC₅₀^exp≤48000 nM and less active 382000nM ≤IC₅₀^exp≤671000nM. For more information on the different interactions between FLAVs and HPV16 E6 responsible for their inhibition potency, we first performed aQSAR model of molecular complexation allowing us to identify the active conformations of the FLAVs and derive the contributions of individual residues filling hydrophobic and hydrophilic pockets in view of the nature of the substituents of the FLAV1, FLAV2, FLAV9 and FLAV16. 

The interaction energy analysis is displayed on Figure 4. The comparative analysis of the per residue energetic contribution shows that the individual IE contributions of the catalytic residues of the HPV16 E6 active site are similar for the three classes of inhibitors, except for Arg102, which shows negligible energetic contributions for FLAV9 and FLAV16 as observed on Figure 3.

In lack of structural clues from IE breakdown to enzyme active site residues contribution, we retained a combinatorial approach to the design of new analogues virtual library and screen it with the help of HPV16 E6 inhibition pharmacophore (PH4) derived from the complexation model QSAR. The best hits are evaluated subsequently by our QSAR model (pIC₅₀^pre₌ - 0. 5494 ×∆∆G_com+5.9983). ∆∆G_com is the relative variation of the relative Gibbs free energy of each best hit formation.

Ligand-Based 3D-QSAR Pharmacophore Model of Inhibitory Activity

The 3D-QSAR pharmacophore protocol of the Discovery Studio molecular modeling¹⁵ program provides the 3D-QSAR pharmacophore for HPV16 E6 inhibition which was generated from the active conformations of the 16 ligands in the training set (FLAV1-16) covering a range of experimental activities 850 nM ≤IC₅₀^exp≤671 000 nM and evaluated by the other 3 in the validation set[FLAV17 (IC₅₀^exp₌40 000 nM), FLAV18 (IC₅₀^exp₌62 200 nM) and FLAV19 (IC₅₀^exp₌77 000 nM)]. The PH4 3D-QSAR generation was carried out in three stages: constructive, subtractive and optimization (Mehodes section). During the construction phase, FLAV1 was selected as a priority because it is the unique molecule that meets the threshold criterion IC₅₀^exp≤Uncert × min (IC₅₀^exp)  × 850.

It was then used to generate the characteristics of the starting pharmacophore. The uncertainty is  because the molecules are not from the same laboratory, so the experimental inhibition activities will not be measured under the same conditions. “Uncert” is a value between 1 and 3. During the subtraction phase, HypoGen eliminates the characteristics present in more than half of the molecules it considers as inactive, that is those with IC₅₀^exp> 850 × 10^3.5=2 687 936 nM. Therefore, none of the ligands in the formation set were found to be inactive, and no functionality of the initial pharmacophore was removed. Finally, during the optimisation phase, the pharmacophore hypothesis scores were improved. The hypotheses selected are validated on the basis of these scores, which assess the differences in the activity estimates obtained from the regression, as well as the complexity of the pharmacophore, using a simulated annealing approach¹⁶. This algorithm applies in fact perturbations to the pharmacophore developed during the construction and subtraction phases in order to optimise the score. After optimisation, the ten best pharmacophore hypotheses were retained. The cost values, correlation coefficients, RMSDs, pharmacophore functionalities and maximum goodness-of-fit values for the ten highest-ranked hypotheses (Hypothesis 1 to Hypothesis 10) are presented in Table 3. These hypotheses were selected on the basis of important statistical criteria, such as a high correlation coefficient, low total cost and a reduced RMSD value.

The reliability of the generated pharmacophore models was assessed as a function of the cost parameters, which ranged from 73.3 (hypothesis 1) to 117.8 (hypothesis 16). The small variation between the cost values reflects the homogeneity of the hypotheses generated and the consistency of the training set. For this pharmacophore model, the fixed cost (53.9) is lower than the zero cost (223.9), with a difference of ∆=169.9. This difference is a key indicator of the quality and predictability of the model (∆>70 indicates a high probability of more than 90% that the model represents a true correlation). 

A hypothesis is considered statistically significant if its total cost is close to the fixed cost and far from zero cost. For the 10 hypotheses, the minimum cost difference is 106.1, which indicates the high quality of the pharmacological model. In addition, standard indicators such as the root mean square deviation (RMSD) range from 1.610 (hypothesis 1) to 3.809 (hypothesis 10) and the correlation coefficient (R²) is between 0.82 and 0.95. The configuration cost, which is 14.68 for all hypotheses, well below 17, confirms the relevance of the model. Thus, the first hypothesis, with a cost of 73.3, close to the fixed cost (53.9) and with the best RMSD and R² correlation coefficient values (Table 3), was selected for the in silico screening of the virtual library of FLAV analogues. The additional statistical parameters for the regression of pIC₅₀^exp against pIC₅₀^pre calculated from Hypothesis 1 for the training setare: pIC₅₀^exp=1.0177 × pIC₅₀^pre - 0.0927 (n=16, R²=0.90, R²xv=0.89, F-test=129.07, σ=0.284, α> 98%). These parameters are in accordance with the OECD QSAR guidelines²⁸.

To assess the predictive power of the pharmacophore model, we calculated the ratio between the activities predicted by the PH4 model and those observed experimentally (pIC₅₀^pre/pIC₅₀^exp) for the validation set (FLAV17-19). The following ratios were obtained FLAV17: 1.021, FLAV18: 1.022 and FLAV19: 1.022. All ratios were close to one, demonstrating the strong predictive power of this regression for the optimal PH4 model. Another evaluation of hypothesis 1 is the mapping of the PH4 binding mode in the 3D QSAR (see Figure 5) of the most active flavonols FLAV1 and FLAV2. Jonathan J., Cherry et al.¹⁰, reported that the new flavones synthesised on the basis of the FLAV9 structure, notably, FLAV2, are all in a hydrophobic pocketat the interface between HPV16 E6 and ubiquitin ligase (E6AP). Figure 5 illustrates the correlation and the detailed geometry as well as the position of the hypothesis 1 features for HPV16 E6 inhibition by FLAV1 and FLAV2. 

In silico screening of FLAVs library

We have constructed a targeted virtual library (TVL) of novel flavonoid compounds with the aim of contributing to a more potent oral small-molecule targeted therapy against HPV16 E6. To do this, we substituted small fragments (R₁ to R₉) for the different aromatic rings of the FLAV scaffold based on structural information gathered¹⁰, including the absence of hydroxyls on ring A of the FLAV9 scaffold, which if replaced by a substituted benzene or heterocycle in the ortho position(R₂) would restore activity. With regard to the flavone analogues used¹⁰, we note that the presence of fragments (OCH₃, -CF₃) in the ortho position would contribute to improve biological activity 1100 nM ≤IC₅₀^exp≤5200 nM.

It was also proposed other substituents of a hydrophilic nature that would facilitate hydrogen bond contacts with the residues of the HPV16 E6 active site, which could improve the activity of the analogues proposed by Srikanth Kolluru et al¹³. Consequently, depending on the position of R-groups on the FLAVs scaffold, we have specified the type of suitable fragments. Table 5 shows the different fragments associated with the scaffold R-groups. Considering the nature of the substituents proposed on the flavonoid scaffold according to the position of the R groups, a combinatorial library of the following size was constructed: R₁ (7) × R₂ (5) × R ₃ (7) × R ₄ (7) × R ₅ (2) × R ₆ (7) × R ₇ (8) × R ₈ (7) × R ₉ (7)=9,411,920 new analogues. This library was then filtered without violating Lipinski's rules (molecular weight greater than 500 daltons)²⁹³⁰ and Veber's rules, which take into account 5 criteria (validity, efficacy, biocompatibility, economy and reproducibility) ³⁰. This led us to a virtual library (VL) of 304,628 molecules, screened using the PH4 3D QSAR Hypo1 model of HPV16 E6 inhibition. The best fitting analogues (PH4 best hits) were then screened using the QSAR complexation model. The relative Gibbs free energies (∆∆G_com) calculated for HPV16 E6-FLAVs complex formationwith their components and predicted inhibitory concentrations obtained from the QSAR model complexation correlation equation (Table 3 and Table 6).

Analysis of New Inhibitors

The design of the virtual library of new analogues was guided by structural information derived from the active conformation of FLAVs, used to select appropriate substituents at positions R₁, R₂, R₃, R₄, R₅, R₆, R₇, R₈ and R₉. In order to identify substituents that could generate new inhibitor candidates with enhanced predictive biological activities against HPV16 E6, we constructed frequency histograms for substituents at positions R₁, R₂, R₃, R₄, R₅, R₆, R₇, R₈ and R₉ from the 24 top hits from PH4 screening (Figure 6).

These histograms show that substituents (1) and (3), including hydrogen atom (H) and methoxyl (OCH₃), are the most frequent in position R₁, with frequencies of appearance of 11and 12 respectively.

The fragments in position R₂, in particular fragment (5), acetamidin is the most repetitive and is followed by the substituent (13) pyrrole with frequencies of appearance of 19and 4 respectively.

As for the substituents in position R₃, the substituents with the highest frequency of occurrence are substituents 1(H) and4 (COOH), with frequencies of 6 and 5 respectively. 

Concerning the fragments in position R₄, the relative histogram reveals that substituent (1) is the most important with an appearance rate of 17.

As for the frequency of appearance of the fragments in position R₆, the corresponding histogram shows that substituent 1(H) appeared most frequently in this position with a frequency of 11, followed by fragment 4 (COOH) which appeared 6 times. As for the histogram related to the substituents in position R₇, it is observed that fragment 9 (-CF₃) is in first place in this position with a significant number of occurrences of 16, followed by substituent (1) with an occurrence value of 8. 

The histogram showing the frequency of appearance of substituents in position R₈ shows that fragment (3), methoxyl (OCH₃) is the most repeated in this position, with a maximum appearance score of 23. 

As for the substituents in position R₉, according to the corresponding histogram, we note that the substituent with the most occurrences is hydrogen atom (1) with a frequency of occurrence of 12, and then we have hydroxyl (2) with a score of 7.

It is important to specify that the substituents at position R₅ are substituent 1 (H) or 2 (OH) with frequencies of appearance of 15 and 9, respectively. 

In view of the above, the proposed 24 new FLAV analogues are mainly flavones, since the hydrogen atom containing fragment (1) appeared most frequently at position R₅ with an appearance score of 15. The hydroxyl (OH) appeared only 9 times. In addition, the fragment 5 in position R₂, the heterocycle acetamidin has a strong appearance with a maximum score of 19 for 19 analogues out of 5 remaining analogues divided into four appearances of the fragment (13) namely pyrrole and one appearance of benzoic acid (fragment15).

ADME Profiles of Designed FLAVs

The properties related to ADME such as Caco-2 cell permeability, blood-brain partition coefficient, octanol- water partition coefficient, aqueous solubility, number of likely metabolic reactions, serum protein binding and another eighteen descriptors related to absorption, distribution, metabolism, and excretion (ADME) were calculated by the QikProp program²⁴ for the new best FLAV analogues (Table 7). The method of Jorgensen is used by this program. Empirical data from more than 710 compounds including about 500 drugs and related hetero-cycles were used to produce regression equations correlating experimental and computed descriptors resulting in an accurate prediction of pharmacokinetic properties of molecule. The pharmacokinetic profile of the best designed FLAV analogues is favourable compared with current anticancer compounds in use or in clinical trials.

DISCUSSION

Binding Mode of FLAV

The measurements of HPV16 E6 inhibition by small molecules belonging to the flavonoid family were reported by Srikanth Kolluru et al¹³. The Intermolecular stacking interactions between FLAV1(IC₅₀^exp=850 nM)and HPV16 E6 of hydrophobic nature (pi-alkyl, pi- sulphur) and hydrogen bonds (see Figure 3) were the main determinants of the better affinity with the target. Using the SAR study of new flavones synthesised on the basis of the chemical structure of FLAV9 against HPV16 E6 by Jonathan J., Cherry et al.¹⁰, we were able to extend our data set with known inhibition concentrations for the implementation of the QSAR molecular complexation model. This allowed us to identify the active conformations of the FLAV analogues used.

The specificity of the newly synthesised flavones, resulting from the absence of hydroxyls on the A-ring of FLAV9 in positions R₁ and R₃, replaced by a heterocycle or benzene substituted in the ortho position, in particular R_2, allowed the biological activity to be restored from 23,000 nM (FLAV9) to 1,100 nM (FLAV2), the most active flavone in this series of inhibitors. These molecules bind in a hydrophobic pocket between HPV16 E6 and the ubiquitin ligase (E6 AP). Using the known molecular specificity of FLAVs, we were able to build a virtual library of specific R₁, R₂, R₃, R₄, R₅, R_6, R₇, R₈ and R₉ fragments that can be accommodated inhydrophobic or hydrophilic pockets, such as hydrogen bonds.

The virtual screening of new FLAV analogues by the pharmacophore whose hypothesis1 allowed us to identify the best hits. These, evaluated by the QSAR model (pIC₅₀^pre=-0.5494×∆∆G_com+5.9983).Among them is a potential new flavone, FLAV 1-5-1-1-1-4-9-3-2(∆∆G_com= - 4.2; IC₅₀^pre= 5 nM), which is substituted in R₅ by the substituent 1(H) and respects the structural modification reported by Jonathan J., Cherry et al. In fact, the fragment in ortho position such as R₂ is a substituted heterocycle, in particular acetamidin (fragment 5), and the substituents in positions R₁ and R₃ are all hydrogen atoms (fragment 1).

In addition to this new potential flavone analogue developed in silico corroborating with the work of Jonathan J., Cherry et al.¹⁰, our study allowed us to discover important new flavones although the cycle A, substituted at position R₂ by a heterocycle, also has substituents in the R₃ or R₄ position, other than fragment 1 (H). Among these we can cite:

• FLAV 1-13-7-1-1-4-9-3-1(∆∆G_com =- 2.22; IC₅₀^pre= 61 nM), substituted in position R₃ by the fragment 7(-CH₂NH₂), which is 14 times more active than FLAV1;
• and FLAV 2-15-6-1-1-4-1-1-1 (∆∆G_com =- 2.38; IC₅₀^pre= 50 nM) where A-cycle is substituted in the R₁ and R₃ positions respectively by fragment 2 (OH) and fragment 6 (-NH₂) and which is 17 times more active than FLAV1. It is important to remember that flavonol FLAV1 is the most active ligand of the training set (IC₅₀^exp= 850 nM), which is followed by the flavone FLAV2 with IC₅₀^exp= 1100 nM.

The observed intermolecular interactions (see Figure 7) between HPV16 E6 oncoprotein and new flavones are responsible for the stability of the newly developed flavones, which is reflected in the different ∆∆G_com values used to account for it.

In fact, these interactions are hydrogen bonds, van der Waals, hydrophobic and halogen. Referring to Figure 7, it is noted that the interactions of hydrophobic nature are of the type:

• Pi-sulphur, for which the A-cycle FLAV 2-15-6-1-1-4-1-1-1 (IC₅₀^pre= 50 nM) formed a Pi-sulfur interaction with the residue Cys51;
• Pi-cation, where only the fragment 13 (pyrrole) in R₂ position of cycle A of FLAV 1-13-7-1-1-4-9-3-1(IC₅₀^pre= 61 nM) formed an interaction pi- cation with Tyr 70;
• Pi-Pi stacked, for which only the C cycle of the new FLAV 1-5-1-1-1-4-9-3-2 (IC₅₀^pre= 5 nM) and FLAV2-15-6-1-1-4-1-1-1 (IC₅₀^pre= 50 nM) maintained this interaction with Tyr 70;
• Pi -Alkyl, where virtually A and B cycles of our threenew potential flavones analogues formed a Pi -Alkyl interaction with the Leu50 and their respective C cycle also formed this interaction with the residue Val 62.

The residues of amino acid Leu50, Cys51, Val 62 and Tyr 70 also make from residues belonging to the hydrophobic pocket of HPV16 E6 for different FLAVs in complex such as FLAV9, taxifolin, alizarin³⁰.

It is important to specify the presence of a fluorine halogen interaction observed in the active conformations of FLAV 1-5-1-1-1-4-9-3-2 (IC₅₀^pre= 5 nM) and FLAV 1-13-7-1-1-4-9-3-1(IC₅₀^pre=61 nM) within the active site of HPV16 E6, in which respectively two atoms of fluorine substituted to the carbon atom (-CF₃ ) whose fragment 9 at position R₇ on the C cycle of the scaffold form respectively a halogen interaction with an electron donor atom of Ala 61 contrary to FLAV 1-13-7-1-1-4-9-3-1(IC₅₀^pre= 61 nM), where a fluoride atom forms a double halogen interaction with the catalytic residue Cys51 and Ala61 and the other two fluorine atoms form this type of interaction respectively with Cys51 and Ala61.

Thus, this new halogen interaction observed in the active conformations of these best new analogues proposed, marked by the presence of a halogen fragment at R₇ would be promising in the search of more potent oncoprotein HPV16 E6 inhibitors.

Considering the diversity and particularity of the structures of our three new potential analogues designed with 4.9 nM ≤IC₅₀^pre≤60.5 n Min relation to their in silico inhibition potency, these molecules are worth being synthesised and evaluated biologically in the development of small therapeutic molecules for the treatment of HPV infection and cervical cancer.

Molecular dynamics simulations

Small differences between the superposition of active FLAV inhibitor conformations modelled by in situ modification of FLAV1 using MM and those obtained via MD (Figure 12), suggest that the best modelled HPV16 E6 -FLAVs complexes are stable for FLAVs2-15-6-1-1-4-1-1-1 (IC₅₀^pre= 49.7 nM) and 1-13-7-1-1-4-9-3-1(IC₅₀^pre= 60.5 nM) respectively. These complexes respectively have RMSD values below 3Å (see Figure12), which generally indicates good stability of the complexes during simulation³². The FLAV 1-5-1-1-1-4-9-3-2 although by virtual screening and its evaluation by correlation equation of the QSAR model is potentially inhibiting with IC₅₀^pre= 5 nM, this one, is less stable referring to the molecular dynamics calculation (RMSD=9.01).

Consequently, the FLAVs 2-15-6-1-1-4-1-1-1 (IC₅₀^pre= 50 nM) and 1-13-7-1-1-4-9-3-1(IC₅₀^pre= 61 nM) newly designed flavones with predicted inhibitory potency 17 and 14 times that of the most active training set FLAV1 (850 nM) can be proposed for synthesis and subjected to biological evaluation³².

Limitations of the study

The main limitation of our MM-PB study is the fact that the inhibitors and their experimental inhibition concentrations used for the QSAR model of molecular complexation do not come from the same laboratory, although the relative variation in the free enthalpy of complexation ∆∆G_com allowed us to explain the variation in the biological activity of the inhibitors used for the QSAR model of molecular complexation.

CONCLUSIONS AND RECOMMENDATIONS      

The study of identification of potential flavonol pocket (FLAV1, FLAV3, FLAV10, FLAV11 and FLAV17) on HPV16 E6 and the structural SAR study of luteolin derivatives (FLAV9) as new target inhibitors, guided us in the preparation of a QSAR model for reliable complexation of HPV16 E6 inhibition which is correlated with the calculated relative Gibbs free energies to form a complex with observed HPV16 E6 inhibition potencies.

In addition, we developed a PH4 3D-QSAR model from active flavonoid conformation (FLAV) using a 16 FLAVs training set and a 3 FLAVs validation set with known experimental inhibition activities. A careful analysis of the interactions between the residues from the active site of HPV16 E6 and the FLAV has directed us to design a first virtual combinatorial library of new analogues of FLAV with multiple substitutions in position R₂ and R₇ for hydrophobic groups and hydrophilic groups in position R₁, R₃, R₄, R₅, R₆, R₈ and R₉. The library screening by matching analogues to the PH4 pharmacophore allowed for selection of a subset of FLAV in the library. This subset of the 24 best virtual results was submitted to the calculation of inhibitory potencies predicted by the QSAR complexation modelin the two digits nanomolar range of concentration.

By molecular dynamics simulations, we verified the stability of new potential analogues and retained the following FLAVs: FLAV 2-15-6-1-1-4-1-1-1 (IC₅₀^pre= 50 nM) and FLAV 1-13-7-1-1-4-9-3-1(IC₅₀^pre=61 nM) generated by computer-aided molecular design (CAMD) are recommended for synthesis and biological evaluation.

ACKNOWLEDGEMENTS

The authors would like to thank the Laboratory of Fundamental and Applied Physics at NANGUI ABROGOUA University, in Côte d’Ivoire, for providing the facilities necessary for this work.

The authors would like to thank the Laboratory of Constitution and Reaction of Matter, University of Cocody (Now Felix Houphouët-Boigny), Côte d’Ivoire, for this collaboration.

AUTHOR’SCONTRIBUTION

Bléhoué I: writing original draft, methodology, investigation. Koné M: editing, review. Kéita M: editing, review. Fofana I: formal analysis. Esmel A: writing, review, and editing. Megnassan E: writing, review, and editing, data curation. All authors read and approved the final manuscript for publication.

DATA AVAILABILITY   

Data will be available on request to anyone from the correspondence author.

CONFLICT OFINTERESTS          

None to declare.

REFERENCES

1. World Health Organization. Cervical cancer.https://www.who.int/news-room/fact-sheands/dandail/cervical-cancer

2. de Sanjosé S, Brotons M, Pavón MA. (2018). The natural history of human papillomavirus infection. Best Practice & Research Clinical Obstetrics Gynaecol 2018; 47: 2-13. https://doi.org/10.1016/j.bpobgyn.2017.08.015
3. Tommasino M. The human papillomavirus family and its role in carcinogenesis, Semin. Cancer Biol 26 (2014) 13-21.https://doi.org/10.1016/j.semcancer.2013.11.002
4. de Martel, Ferlay J, Franceschi S,  et al. Global burden 6 of cancers attributable to infections in 2008: a review and synthetic analysis, Lancet Oncol. 13 7 (2012) 607-615 (2012) 607-615.
5. de Martel, Plummer M, Vignat J, Franceschi S. World-wide burden of cancer attributable to 9 HPV by site, country and HPV type. Int J Cancer 114 (2017) 664- 670.
6. Martinez-Zapien D, Ruiz FX, Poirson J, Mitschler A, et al. Structure of the E6/E6AP/p53 complex.
7. Shrestha AD, Neupane D, Vedsted P, Kallestrup P. Cervical cancer prevalence, incidence and mortality in low and middle income countries: A systematic review. Asian Pac J Cancer Prev 2018; 19:319-324.
8. Parisa Shiri Aghbash et al. siRNA-E6 sensitizes HPV-16-related cervical cancer through Oxaliplatin: an in vitro study on anti-cancer combination therapy. European J Med Res (2023) 28 :42, 13. https://doi.org/10.1186/s40001-023-01014-9

9. Celegato, L. Messa, L. Goracci, et al. A novel small-molecule inhibitor of the human papillomavirus E6-p53 interaction that reactivates p53 function and blocks cancer cells growth, Cancer Landters.
10. Jonathan J, Cherry, et al., (2013) Structure based identification and characterization of flavonoids that disrupt human Papillomavirus-16 E6 Function 2013; 8:12.
11. Zanier K, Charbonnier S, Sidi AOMHO, et al. Structural basis for hijacking of cellular LxxLL motifs by papillomavirus E6 Oncoproteins. Science 2013, 339, 694.
12. Abotaleb M, Samuel SM, Varghese E, et al. Flavonoids in Cancer and Apoptosis. Cancers 2019, 11, 28.
13. Srikanth Kolluru et al. Identification of potential binding pocket on viral oncoprotein HPV16 E6: a promising anti-cancer target for small molecule drug discovery. BMC Molecular and Cell Biology.2019.https://doi.org/10.1186/s12860-019-0214-3

14. Elise Emeraux. Biological properties of flavonoids: bibliographic study and evaluation of antioxidant activity. Sciences Pharmaceut 2019 https://hal.univ-lorraine.fr/hal-03297878

15. Discovery Studio molecular modeling and simulation program, version 2.5, Accelrys, Inc., San Diego, CA, California. USA. 2009
16. Insight II, Version 2005, Molecular Modeling package and Discover Version 2.98. Simulation package, 200x, Accelrys, Inc, San Diego, California, USA. 2005.
17. Kouassi AF, Kone M, Keita M, et al. Computer-aided design of orally bioavailable pyrrolidine carboxamide inhibitors of enoyl-acyl carrier protein reductase of mycobacterium tuberculosis with favorable pharmacokinetic profiles. Int J Mol Sci 16 (12): 29744–71, 2015.
18. Kouassi AF, Kone M, Keita M, et al. Computer-aided design of orally bioavailable pyrrolidine carboxamide inhibitors of enoyl-acyl carrier protein reductase of mycobacterium tuberculosis with favorable pharmacokinetic profiles. Int J Mol Sci 2015; 16: 29744-29771. https://doi.org/10.3390/ijms161226196
19. Fogolari F, Brigo A, Molinari H. (2002). The poisson-boltzmann equation for biomolecular electrostatics: A tool for structural biology. J Mol Recog 15: 377–92.
20. Li H, Sutter J, Hoffmann R. 2000 Hypo Gen: An automated system for generating 3D predictive pharmacophore models. Güner 2000; 171–189.
21. Chemical Computing Group. MOE (Molecular Operating Environment) User Guide 2015 Edition.
22. Leach AR. Molecular modelling: Principles and applications. Pearson Education, 2016.
23. Qik Prop Version 3.7, Release 14. 2014 NY: X Schrodinger, LLC, New York.
24. Jorgensen WL, Duffy EM. Prediction of drug solubility from Monte Carlo simulations. Bioorg Med Chem Lett 2000; 10: 1155–8. https://doi.org/10.1016/s0960-894x(00)00172-4

25. Jorgensen WL, Duffy EM. 2002 Prediction of drug solubility from structure. Adv Drug Deliv Rev 2002; 54:355–66. https://doi.org/1016/s0169-409x(02)00008-x

26. Duffy, Jorgensen W. 2000 Prediction of properties from simulations: free energies of solvation in hexadecane, octanol, and water. J Am Chem Soc 122, 78-88.https://doi.org/10.1021/ja993663t

27. Maestro-desmond interoperability tools. In Desmond Molecular Dynamics System; Version, 3.6; Maestro-Desmond Interoperability Tools, Schrödinger; D. E. Shaw Research: New York, NY, USA, 2021.
28. Guidance Document on the Validation of (Quantitative) Structure-Activity Relationship [(Q)SAR] Models, OECD series on testing and assessment 2014; No. 69, OECD Publishing, Paris. https://doi.org/10.1787/9789264085442-en

29. Lipinski CA, Lombardo F, Dominy BW, Feene, PJ. Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings. Adv Drug Del Rev 2001, 46, 3–26.
30. Veber, DF, Johnson, SR, Cheng, HY, Smith, BR, Ward, KW, and Kopple, KD. Molecular properties that influence the oral bioavailability of drug candidates. J Med Chem 2002 Jun 6;45(12):2615-23.https:// doi.org/10.1021/jm020017n    

31. Gomes D, Yaduvanshi S, Silvestre S, et al. Taxifolin and Lucidin as Potential E6 Protein Inhibitors: p53 Function Re-Establishment and Apoptosis Induction in Cervical Cancer Cells. Cancers 2022; 14: 2834.https:// doi.org/10.3390/cancers14122834

32. Lemkul JA. From proteins to perturbed Hamiltonians: A suite of tutorials for the GROMACS-2018 molecular simulation package. Biophys J 2019; 117(3), 354-361.