In silico testing of flavonoids as potential inhibitors of protease and helicase domains of dengue and Zika viruses

Identity analysis for the DENV and ZIKV polyprotein

Sequence alignments of polyproteins from DENV and ZIKV, and NS3 proteins were carried out using Clustal Omega (https://www.ebi.ac.uk/Tools/msa/clustalo/) (Mount, 2009). Viral sequences used for all identity analyses and to calculate the similarity matrix were obtained from the NCBI “Virus Variation Resource” database (https://www.ncbi.nlm.nih.gov/genome/viruses/variation/). We aimed to include up to six viral genomes for ZIKV and six for each DENV serotype from each world region: Africa, Asia, North America, Oceania, South America, and Europe. The selection was made considering recent complete sequences known to infect humans. In the case of ZIKV there were six genome sequences for each region, but for DENV fewer sequences were available: sequences for DENV3 and 4 were lacking in Africa as well as for all the serotypes in Europe (Table S1). Complete NS2B sequences were also analyzed for identity but are not shown.

Structural similarity analysis of the NS3 protein in DENV and ZIKV

To compare the 3D structure of the NS3 protein for DENV and ZIKV, the structures for ZIKV (PDBIDs for NS3-pro: 5YOD and for NS3-hel: 5TGX) and for each of the four DENV serotypes (PDBIDs for NS3-pro: 3L6P, 2FOM, 3U1I, 5YVV; for NS3-hel: 2BMF and 2JLS, no structures available for DENV1 and 3 NS3-hel) were obtained from the Protein Data Bank (PDB) (https://www.rcsb.org/) (Table S2). Superposition analysis was performed using ProFit v3.3 (Sumathi et al., 2006), first on the complete NS3 proteins and afterwards limited to the active sites: the protease and the ssRNA binding site of the helicase. Root-mean-square deviation (RMSD) expressed in angstroms (Å) for N, C, C α and O atoms of the protein backbone was calculated.

Selection and preparation of ligands and target molecules

A total of 15 molecules previously identified by our group in extracts of Taraxacum officinale and Urtica dioica (Flores-Ocelotl et al., 2018) were used for the screen (Table 1). These were expanded to 20 since mass spectrometry cannot distinguish between isomers. Their 3D structures were obtained from the PubChem database in SFD format and converted to mol2 format using Avogadro (Hanwell et al., 2012). Ligands were prepared prior to molecular docking, by structural optimization using the general AMBER force field (GAFF) to optimize drug geometries (He et al., 2020). As target receptors, the best quality crystal structures of NS3 proteins were chosen, to represent each DENV serotype and ZIKV (Cozzini et al., 2008; Scapin, 2006). The receptors structures in PDB format were loaded to VEGA ZZ version 3.2.0 (obtained from https://www.ddl.unimi.it/cms/index.php?Software_projects:VEGA_ZZ:Download). VEGA ZZ was then used to clean the structures by removing solvent molecules, salts, and ligands. Then, the structures were optimized, and partial charge assignment was done prior to docking.

Two ligands were added as negative and positive controls, respectively: dibasic phosphate and Bz-NIle-KRR-AMC. Phosphate is often used for crystallization and may even be solved in the final structure due to the high concentrations used. Thus, it is expected to be an unspecific binder which can be used as a substrate to measure enzyme kinetics for the DENV protease and is expected to have significant affinity.

Molecular docking

For docking, ADFR 1.2 rc1 (Ravindranath et al., 2015) was used. Ligands and receptors were converted to pdbqt using the scripts included in ADFR. All receptors were aligned to have the same reference coordinates. The grid for docking (trg files) was prepared using agfrgui and the bound peptide from structure 3U1I, as a reference for the binding site location. For helicases, we selected the biggest binding site identified by agfrgui, coordinates −16.525, 44.133, 23.142 and length 31.5, 42.5, 40.5. Boxes for the proteases were created with coordinates −35.132, −28.372, −24.665 and length 31.250, 26.750, 31.250 (spacing for both was 0.375). For NS3-pro docking NS2B polypeptides were included. For each ligand, 500 runs with 2.5 million energy evaluations were carried out. This set up keeps the receptor rigid but ligands are flexible. Results were visualized on UCSF Chimera (https://www.cgl.ucsf.edu/chimera/download.html) (Pettersen et al., 2004). The results were selected based on largest cluster size first and then by highest affinity to finally graph them.

Single point MMGBSA after short molecular dynamics

After docking with ADFR, UCSF DOCK (http://dock.compbio.ucsf.edu/DOCK_6/index.htm) was used to sample different protein-ligand conformations in a set up that lets the receptor move. (Graves et al., 2008). The most populated ligand pose from ADFR was selected as the initial conformation for molecular mechanics. Briefly, ligand was energy minimized against a rigid receptor then both were parametrized with AMBER forcefield. Then, a short (i.e., 10000 steps) molecular dynamic was run allowing for the movement of all the atoms in the complex. This was followed by 500 steps of energy minimization and MMGBSA evaluation of the energetics.

Molecular dynamics and trajectory MMGBSA calculations

Single point MMGBSA is ran after a short molecular dynamic simulation, usually a few thousand steps. Thus, they are unable to catch ligand dissociation or any change that takes significant evolution time, i.e. hundreds of nanoseconds. To try and observe those events we used AMBER20 to run molecular dynamics and MMGBSA calculations over tens of nanoseconds. The amber19SB forcefield was used for the proteins, OPC model for water molecules, and gaff2 for ligands. Charges for ligands were derived using ab initio QM (mp2/6-311g(d)). Neutralizing ions were added to a concentration of 0.15 M, water box was built with a10 angstroms distance from the solute to the box edge. Production simulations were run for 100 nanoseconds (ns) with five repeats. MMGBSA calculations were run with the MMPBSA.py module, using the GB method 8 (mbondi3), and a salt concentration of 0.15 M. For negative control we used dibasic phosphate and for a positive control the inhibitor aprotinin as its affinity is known to be in the nanomolar range.

Sequence similarity between DENV, ZIKV polyproteins and NS3

It is unclear if the sequence similarity between DENV and ZIKV is enough to allow a single molecule to target both pathogens. To explore this, we conducted sequence analyses. The four DENV serotypes shared 68–78% identity throughout their polyprotein sequence (3,430 aa), while comparison with ZIKV showed lower identity of ∼55% (Table S3).

On NS3, DENV serotypes shared 75–86% identity, while the identity was again lower (64–68%) between DENVs and ZIKV (Table S4). NS3-hel (438 residues) was more conserved between DENVs and ZIKV (69–72% identity) (Table S5), than NS3-pro, (180 aa; 51–58% identity) (Table S6).

Sequence and structure similarity of the NS3-hel from DENV and ZIKV

The RNA binding cleft in the NS3-hel, where the RNA substrate makes temporary contact during unwinding for replication, is made up of two alpha helices from domain II and two alpha helices from domain III (Fig. 1A) (Xu et al., 2005). Between DENVs and ZIKV, helices α7” DIII, α3” DIII, α1′DII and α2′DII had 69–84%, 60–73%, 53–66% and 45–81% identity, respectively and showed higher similarity within DENV serotypes (Table S7).

Superposition of the two DENV and one ZIKV helicase structures available (Table S2) also suggested more similarity between DENV serotypes than to ZIKV. Besides, dengue structures displayed a more closed conformation. Main differences are located at two alpha helices (α1′DII and α2′DII) that make the RNA entrance to the cavity, with residues D290 and R538 (Fig. 1A, in sticks) in different orientation. The RMSD was 0.8 Åbetween DENV2 and DENV4; 1.1 Åbetween DENV2-ZIKV, and 1.25 Åbetween DENV4-ZIKV (Table S8). Given the similarity in secondary and tertiary structures, we considered NS3-hel as a relevant target for docking in the next part of the study.

Sequence and structure similarity of the NS3-pro from DENV and ZIKV

The protease carries out an important step in viral maturation: processing the polyprotein into active single proteins. Sequence similarity within DENV serotype was above 93% except for DENV2 that showed 84%. Across dengue viruses, identity was as low as 62% (DENV1 vs DENV4) and as high as 76% (DENV1 vs DENV3). When comparing DENVs to ZIKV, identities dropped to ∼55% (Table S6).

Next, we analyzed NS3-pro binding site sequence conservation. We selected each of the catalytic-triad residues (H51, D75 and S135, numbering from structure 2FOM) in a block of 11 aa that define their immediate environment: catalytic residue plus 5 adjacent residues towards the amino- and 5 towards the carboxy-terminus. The block surrounding H51 was 90.90–100% conserved among DENV1, DENV2 serotypes and ZIKV, with 80.81–90.90% when comparing to ZIKV. For the block around D75, identities between DENV serotypes were 54.54–90.90% and 45-63.63% with ZIKV; while for the region around S135, DENV serotypes were 81.81–100% identical and 54.54–72.72% with ZIKV. Thus, the region around H51 was the most conserved (Table S9). Overall, the catalytic site showed high variation in terms of sequence, except for the catalytic triad.

To quantify differences in three dimensions, we analyzed five NS3-pro structures: one for each DENV serotype plus one for ZIKV (Table S2). The catalytic triad did not display much orientation variation and the main difference was observed in the loop with residues Y161 and N152, representing different catalytically important structural conformations in the different proteases (Fig. 1B, sticks). Within DENV serotypes, the structures with the greatest overlap were DENV1-DENV3 with an RMSD of 0.43 Å, followed by DENV1-DENV2, DENV1-DENV4 and DENV2-DENV4 with 0.64 Å, 0.79 Åand 0.83 Å. ZIKV was more like DENV3 with RMSD 0.51 Å, followed by DENV2, DENV1 and DENV4, with 0.63 Å, 0.67 Åand 0.74 Å, respectively (Table S10). Taken together, these data suggest that NS3-pro are diverse in sequence yet conserved in structure, making plausible the goal of identifying a single lead for drug discovery that can target several flaviviruses.

Molecular docking of ligands against the NS3-hel and -pro domains of DENV and ZIKV

Binding energies for the 20 ligands listed in Table 1 were determined using ADFR and were worked with each of the structures used for 3D superposition. While it is common to use the best pose, we chose to first rank results using the largest cluster, to focus on the most frequent results and to avoid bias towards infrequent poses with high affinity (Fig. 2).

NS3-helicase

Despite the different openings of the RNA-binding sites, docking showed similar overall distribution of the 20 ligands in the three NS3-hel structures tested (DENV2, DENV4 and ZIKV) with a mean binding energy of ∼8.3 kcal/mol. Five ligands stood out with affinities <1 standard deviation from the mean, at around -10 kcal/mol: 3,5-dicaffeoylquinic acid (DCA01), quercetin 3-rutinoside (QNR05), Quercetin 3,7-diglucoside (QND10), Quercetin 7,4′-Diglucoside (QND12) and luteolin 3′, 4′-diglucoside (LND17). The first one (DCA01; in Fig. 2A) binds with similar affinity to all the helicases compared. The remaining 4 bind strongly to a different pair of receptors: QNR05 (asterisk) binds tightly to NS3-hel DENV4 and less so to ZIKV, QND10 (+plus sign) binds well to the DENV2 structure at −10.8 kcal/mol while for DENV4 and ZIKV was −8 and −8.7 kcal/mol respectively, QND12 (x cross) binds better to NS3-hel DENV4 and ZIKV, while LND17 (slanted square) binds well to NS3-hel DENV2 and DENV4. However, we observed that there was no binding specificity within the RNA binding site, rather they presented a favorable binding energy but at significantly different locations within the site. This can be observed for DCA01 (Figs. 3A–3C) for QND12 (Figs. 3D–3F) and the rest of the ligands.

After single point MMGBSA (Fig. S1) the results did not show an improvement: no single ligand was observed as a common high affinity binder located at a specific site in the studied NS3-helicases. This is compatible with the ligands docking all along the RNA-binding cleft in the NS3-helicases (Fig. 3), which is large even in DENV structures that tend to be in a closed conformation. The large NS3-hel cleft likely provides little substrate specificity for our small molecules. Despite the high affinities calculated by ADFR or single point MMGBSA for some of the ligands tested, it is unlikely that these ligands are specific to NS3-hel, due to this observed drawback, an exhaustive analysis of molecular dynamics was not carried out.

NS3-pro

NS3-pro docking (ADFR) results showed that some receptors bind all the ligands slightly better than others: DENV2 had the weakest mean binding energy (−6.2 kcal/mol), followed by ZIKV and DENV1 (−7.4 kcal/mol), while DENV3 and DENV4 showed the highest mean affinities (−7.8 and −8.3 kcal/mol, respectively). Only on DENV4 NS3-pro, two ligands had results better than −10 kcal/mol: quercetin 3-rutinoside (QNR05) and quercetin 3,7-diglucoside (QND10) (Figs. 2B; 4D and 4F). These results are compatible with the low identities between NS3-pro sequences (Table S6) and the observation that binding site residues in sub-pockets S1 and S3 show a different orientation between the structures analyzed (Fig. 1B). In the results ADFR DENV4 had the best affinities with the mentioned ligands, when considering the analysis with DENV3, DCA01 presented an affinity of −9.7 kcal/mol followed by DENV1 with an affinity of −9.1 kcal/mol; therefore, it was also included in MMGBSA analyzes.

After single point MMGBSA (Fig. S2), DCA01 was found distant from the catalytic residues both in DENV3 and DENV4 NS3-pro (Fig. 4A–4B), occupying part of sub-pocket S1 and S3; in summary the results did not improve for DCA01.

QNR05 and QND10 converged between Y161 and N152 and near S135 (Figs. 4C through 4F) filling up sub-pockets S1′, S1 and S3. In particular, the quercetin core (bold bonds in Figs. 5A and 5B) lodged in sub-pocket S1 in a similar orientation in both DENV3 and DEVN4 NS3-pro, while the saccharides anchored the ligands to the other sub-pockets. Sub-pocket S2 was not visited by any of the ligands. With this analysis we determine that QNR05 and QND10 was the best binder for both DENV3 and DENV4 NS3-pro with a binding energy of −34.41 and −37.05 kcal/mol for QNR05 and −29.86, −36.26 kcal/mol for QND10 and at the same time opens the opportunity to make improvements to these molecules

NS3-pro molecular dynamics and MMGBSA

The complexes obtained after docking were used as starting points for standard molecular dynamics simulation. Five production simulations were run for each receptor–ligand pair for a time of 100 ns each: a total of four microseconds (µs) simulation time. During the simulations, no ligand remained static relative to the receptor. By examining the evolution of binding energy in time a few trends became apparent. DCA01 is the ligand with the least negative binding energy (Fig. S3). During some simulations the ligand left the receptor completely as the energies drop to 0 (Figs. S3A and S3B).

From the point of view of the receptor and the RMSD of its backbone atoms: all ligands allow for a relatively low for DENV3 (Figs. S4A to S4C). For DENV4, RMSD in the presence of DCA01 reached up to 4 Å (Fig. S4A). For QNR05, the energies when bound to DENV3 are slightly better than those of DCA01. For DENV4, two trajectories displayed almost −40 kcal/mol but 3 were above −20 for a significant simulation time (Fig. S3C and Fig. S3D). RMSD for these complexes is low for both DENV3 and DENV4 (Fig. S4B).

QND10 displayed a more converged evolution: on DENV3 affinities had values around −20 kcal/mol whereas for DENV4 the average was around −40 kcal/mol (Figs. S3E and S3F). RMSD is also low for both DENVs.

Frequency analysis of the hydrogen bonds between NS3-pro and the different ligands (Fig. S6) shows that QND10 shows increased hbonds when compared to DCA01 or QNR05.

These results are consistent with the initial docking results yet also put into manifest that the binding site of NS3-pro is a difficult binding site to target. Our ligands bind to the site but display a large apparent k_off.

Quercetin 3- β-D-glucoside (CAS no. 482-35-9) is the common core shared by QNR05 and QND10; it anchors these ligands to three of the four sub-pockets. A possible strategy to improve affinity is to use Quercetin 3-β-D-glucoside as lead. A clear path to do so would be to add a functional group to this core to occupy sub-pocket S2. Our results, guided by the electrostatics of that binding site (Behnam & Klein, 2020) (Fig. S5) led us to add a positively charged motif (Fig. 5C). To test this hypothesis, we used the QND10 structure as a template and computationally added a substituent on carbon 6 of the benzene ring. A new molecule, called MOD10, showed better affinity for both receptor structures, with binding energies of −47.84 kcal/mol for DENV3 and −40.40 kcal/mol for DENV4 (Figs. S3G–S3H). The RMSD values for both targets showed that MOD10 reduced the changes to the structures of DENV3 and DENV4 at around 2 angstroms (Fig. S4D) in contrast with the other more diverse trajectories with the other ligands.

A summary of the MMGBSA energies and comparison to negative and positive controls is shown in Fig. S7.