In silico analysis of natural compounds targeting structural and nonstructural proteins of chikungunya virus

Target identification and homology modeling

Structures of CHIKV proteins E1, E2, E3, nsP2 and nsP3 were downloaded directly from RCSB Protein Data Bank (PDB). For the rest of the CHIKV proteins, CP, 6K, nsP1 and nsP4, whose structures are unavailable, CHIKV sequences present in NCBI, belonging to ECSA (East/Central/South Africa) genotype were downloaded. These sequences were utilized to form a consensus sequence with MEGA 6²⁸ using clustalW pairwise multiple alignment algorithm with all parameters set at default. Using these consensus sequences, homologous proteins from the PDB were identified using Protein BLAST²⁹ where the algorithm parameters were as follows: Max target sequences=100, Expect threshold=10 using BLOSUM62 scoring parameters, Gap cost=Existence:11 & Extension:1 with conditional compositional score matrix adjustment. The suitable templates for nsP1 and CP with highest query coverage, sequence identity and lowest E-value were selected for homology modeling. For proteins 6K and nsP4, no templates were available, and thus these structures were created using threading and looping method (see next section).

The template and target sequences of all CHIKV proteins were then aligned using CLUSTALW³⁰. MODELLER (version 9.16) was used to generate homology models³¹. Further, the homology model having the lowest MODELLER objective function (molpdf) or DOPE or SOAP assessment scores, or the one having highest GA341 score was selected as the best model structures and were further utilized for model validation. Nonstructural protein, nsP4, and the small accessory peptide of structural protein 6K did not have any template in PDB; therefore a threading and looping approach was implemented for them using LOMETS (Local Meta Threading Server)³². Both online server and standalone program present as a module of I-TASSER Suite version 5.1, which provides 3D models by combining alignment scores of template to target of 9 different threading programs (FFAS-3D, HHsearch, MUSTER, pGenTHREADER, PPAS, PRC, PROSPECT2, SP3, and SPARKS-X). All parameters were set as default. All structural and nonstructural CHIKV protein sequences were selected as potential drug targets.

Validation of homology modeled structures

Generated models were validated using MolProbity-(v4.3.1)³³. Ramachandran plot analysis was performed for the best protein models by analyzing the phi (Φ) and psi (Ѱ) torsion angles. To check reliability of the modeled structures, the root mean square deviation (RMSD) was calculated by superimposing it on template protein structure using PyMOL (v1.7.0.0) visualization software³⁴. Consistency between templates and the modeled structures were assessed by ProSA-web³⁵ (online server), a statistical analysis tool of all the proteins structures available at RCSB PDB. Here, a statistical average is obtained over the known structures with the help of combined potentials of mean force from the PDB database.

Molecular dynamic simulations

Stability of the domain regions of CHIKV protein structures was examined by molecular dynamics (MD) simulation using GROMACS (version 5.0) software package³⁶. Optimized Potential for Liquid Simulations All-Atom³⁷ force field was used to energy minimize the structures. Through this energy minimization, the high-energy intramolecular interactions were discarded. In order to avoid the steric clashes, overall geometry and atomic charges were also optimized. The proteins were kept at the center of the rectangular box, which was filled with SPC water model system to create the same environmental behavior of the molecules. All the atoms of the protein and the boundary of the rectangular box were separated by a minimum distance of 10 Å. 0.01M NaCl was used as a solvent exposure.

The system was further energy minimized without any restraints for 50,000-time steps; the steepest descent having step size of 0.01 ps. Then the system was equilibrated to reach a stable temperature by conducting NVT ensemble. Pressure was further equilibrated by NPT ensemble performance. The long-range electrostatic interactions were calculated by using particle mesh Ewald³⁸ method with a cut-off of 0.9 nm for Vander Waals interactions. All the bonds were constrained by LINCS³⁹, where only the water molecule moves to equilibrate with respect to protein structure keeping protein molecule as static. To couple the system Berendsen thermostat (V-rescale) and Parrinello-Rahman barostat were utilized to maintain the constant temperature (300 K) and pressure (1 bar). Further MD analysis was performed to observe structural changes and dynamic behavior of the protein by calculating RMSD, radius of gyration and root mean square fluctuation (RMSF) along with changes in temperature, pressure, density and total energy.

Virtual screening and molecular docking

Simulated computational models of CHIKV proteins were prepared and their binding sites were predicted using SiteMap (Version 2.3, 2009, Schrödinger, LLC, New York, USA). These were then used to perform molecular docking. The protein preparation wizard was used to prepare CHIKV proteins and a natural remedies library from MolBase database was prepared using the LigPrep module⁴⁰. Virtual screening of modeled proteins against a natural remedy library from MolBase was done by using GLIDE module in an Extra Precision (XP) mode (Version 5.5. 2009, Schrödinger, LLC). It produces the minimal ranks of inappropriate poses and determines the appropriate binding energy of the three dimensional (3D) structure of the protein along with a ligand^41,42.

Analysis and output visualization of drug target and protein

After the completion of molecular docking, the docked poses were listed depending upon the respective docking scores. Glide Score (obtained using GLIDE Module of Schrödinger Software Suite 9.0) was used as an empirical scoring function to predict free energy for ligands binding to the receptor. The structure showing minimum binding energy was filtered and subjected for further analysis. The 3D conformation ligand receptor was analyzed using PyMOL³⁴ and Chimera⁴³ v1.10.1 visualization software.

Absorption, distribution, metabolism and excretion (ADME) screening and toxicity analysis

Pharmacokinetics properties and percent human oral absorption values were further predicted for the potential lead molecules using QikProp module (Version 3.2, 2009, Schrödinger, LLC)⁴⁴. Both the physically remarkable descriptors and pharmaceutically admissible properties were predicted for neutralized ligands by QikProp. The program predicts 44 different properties, including log P (octanol/water), % human oral absorption in intestine (QP%) and predicted IC₅₀ value for blockage of HERG K+ channels (log HERG). The Lipinski’s rule of five⁴⁵, an important criteria for oral absorption, was evaluated for the acceptability of the compounds. In addition, admetSAR⁴⁶ v1.0 was used to calculated various attributes of the drugs, including the blood brain barrier (BBB), human intestinal absorption, Caco-2 permeable, aqueous solubility, P-gp substrate and inhibitor, CYP450 substrate and inhibitor, CYP IP, ROCT, HERG inhibition, and toxicity parameters. For Lipinski score calculations, the ligand in SMILE format was uploaded to QikProp. The physicochemical properties and Lipinski Rule of Five were also analyzed by PERL script, “CalculatePhysicochemicalProperties.pl” of MayaChemTools⁴⁷.

Ligand-protein interaction studies

The protein-ligand complex interaction at the atomic level was analyzed using Maestro 11.0 (LLC Schrodinger 2016)⁴⁸ and LigPlot+⁴⁹ v1.4.5. The protein and the docked ligand were merged together and uploaded to Maestro Suite vMaestro 11. Further, the “Ligand Interaction Diagram” option was selected to draw the protein-ligand binding interactions in the 2D visualization workspace.

In silico protein preparation, homology modeling and validation

CHIKV consists of four nonstructural proteins (nsP1-nsP4), three structural proteins (E1-E3), along with two sub-pro regions 6K and CP, which makes a part of structural protein unit (Figure 1). Structures of CHIKV proteins E1, E2, E3, nsP2 and nsP3 were downloaded directly from PDB (Figure 2a), and for other CHIKV proteins (CP, 6K, nsP1 and nspP4) homology modeling and threading and looping methods were utilized to predict their structures. For proteins with templates available, homology modeling was done with five models for every protein created based on sequence similarity using different model generation tools (MODELLER and LOMET), and validated by their internal scoring functions (molpdf, DOPE, SOAP and GA341 scores). Further, ProSA Z-score for all modeled structures were calculated to analyze the quality of models based on the Cα positions. Individual validation and ProSA Z scores for top ranked models are given in Table 1 and their structures are given in Figure 2b. The top ranked models were also analyzed by Ramachandran plot (Figure 3). The Ramachandran plot shows the distribution of phi (ϕ) and psi (ψ) angles for each amino acid residues of the modeled structures. The respective percentages of the favored and allowed regions for all the residues of all those validated are also shown in Table 1.

Figure 1. Organization of chikungunya virus genome.

The genome consists of two open reading frames (ORFs) separated by an untranslated junction (J). The first ORF encodes for a polyprotein and acts as a precursor of the non-structural proteins (nsP1, nsP2, nsP3 and nsP4). The second ORF encodes the structural proteins (Capsid, E3, E2, 6K and E1). The genome has 5` cap and 3` poly A tail.

Figure 2. Structures of chikungunya virus proteins.

(A) X-Ray structures of nsP2, E3, nsP3, E1 and E2. (B) Homology modelled structures of nsP1, Capsid, nsP4 and 6K.

Figure 3. Ramachandran plot of chikungunya virus proteins obtained from MolProbity.

(A) nsP1 and Capsid (homology modeling); (B) nsP4 and 6K (Threading/Looping).

Molecular dynamic simulation and analysis

Molecular dynamic simulations were employed to analyze the protein structure-function complexities, such as structural stability, conformational flexibility and folding. Domain regions of the structures (Table 2) were simulated for 20 ns. Moreover, various parameters, such as temperature, pressure, density and total energy, were calculated to check the stability of these structures along with steric properties. Further, RMSD values for the backbone atoms of proteins were plotted against time of MD simulations. Average RMSD during the simulations was 22.93. Radius of gyration on the other hand also supports the stability and compactness of the proteins. The RMSF with respect to each residue depicts the flexibility of the proteins. Average RMSF during the simulations was 1.45. The RMSD, radius of gyration and RMSF plots for all CHIKV proteins are shown in Figure 4A–C. The resulting graphs contributed to protein modeling, as they show a constant RMSD deviation throughout the 20ns simulation except for a small deviation in E2 after 14ns. Depending upon these simulation parameters, the proteins showed conformational stability.

Figure 4. Molecular dynamics profiles of the chikungunya virus (CHIKV) proteins tertiary domain structure optimization.

(A) root mean square deviation (RMSD), (B) Radius of Gyration, and (C) root mean square fluctuation (RMSF). A–C graphs are vs Time, and F vs Atoms. Each set shows the graph for both Non-structural (upper) and Structural (lower) CHIKV proteins. Non-structural Protein: nsP1 (green), nsP2 (blue), nsP3 (yellow) and nsP4 (red); Structural Proteins: Capsid (orange), E3 (mustard), E2 (purple), 6K (cyan) and E1 (pink).

Molecular docking

Drug discovery relies heavily on molecular docking to understand the interactions between ligand/inhibitor and target protein⁵⁰. In this study, we resorted to the docking of available protein structures (wherever applicable), as well validated, refined and simulated modeled proteins to screen against a natural remedy library from MolBase. The binding sites of all protein structures were predicted by SiteMap. The predicted binding pockets were further validated using Glide in XP mode. Top ten ligand/compounds having docking score (Glide Score) above -4, glide energy of -20 kcal mol^-1 and potential energy of a considerable range were considered for the next level of screening. The combined results of all the docked ligand along with the glide energy and potential energy have been provided in Table 3. Of these, two ligands, (1,3,6, -Trigalloyl-β-D-Glucose and Quercetin-3-rutinoside (Compound ID 164 and 153) were found to interact with all the proteins and were discarded from further analysis.

For the non-structural proteins, the top ligands included Rebaudioside A and Withanoside IV (Compound ID 149 and 179) for nsP1; Stevioside, Bacopaside II and Jujubogenin isomer of bacopasaponin C (Compound ID161, 26 and 113) for nsP2; Chebulinic acid and Corilagin (Compound ID 44 and 47) for nsP3; Rebaudioside A and Stevioside (Compound ID 149 and 161) for nsP4. For structural proteins, Catechin-5-O-gallate, Rosmarinic acid and Agnuside (Compound ID 42, 151 and 18) for CP; Bacopaside II, Mangiferin and Arjungenin, (Compound ID 26, 122 and 12) for E3; (Rebaudioside A, Tribulosin and Asiaticoside (Compound ID 149, 165 and 17) for E2; Arjunetin and Stevioside (Compound ID 10 and 161) for 6K; Chebulinic acid, Stevioside and Asiaticoside (Compound ID 44, 161 and 17) for E1. Top four docked poses of the modeled proteins and the small molecules having lowest docking score are shown in Figure 5 (ligand wise).

Figure 5. Binding interaction with the potential lead compounds and their representative binding pocket conformation for the top four docked poses of all chikungunya virus proteins.

Ligands are cyan sticks and receptors as pink ribbon/surface.

Unique ligand-protein partners were taken forward for ADME and toxicity analysis. In case of nonstructural proteins, these ligand-protein pairs were Withanoside IV (Compound ID 179)-nsP1, Jujubogenin isomer of bacopasaponin C (Compound ID 113)-nsP2 and Corilagin (Compound ID 47)-nsP3. In case of structural proteins, these pairs were Catechin-5-O-gallate (Compound ID 42), Rosmarinic acid (Compound ID 151) and Agnuside (Compound ID 18) against CP, Mangiferin (Compound ID 122) and Arjungenin (Compound ID 12) against E3, Tribulosin (Compound ID 165) against E2, and Arjunetin (Compound ID 10) against 6K.

ADME analysis of all potential leads

ADME screening was performed for all the top hits. Here, 44 various physically remarkable descriptors⁵¹ and pharmaceutically admissible properties of the top four lead compounds for every CHIKV protein were calculated using QikPro-P (Table 4). The Lipinski’s rule of five was further employed to evaluate oral absorption along with ADME. Compounds violating more than 2 Lipinski’s rule of 5 were discarded from further analysis.

Compounds Catechin-5-O-gallate (Compound ID 42), Rosmarinic acid (Compound ID 151) and Agnuside (Compound ID 18) against CP; Mangiferin (Compound ID 122) and Arjungenin (Compound ID 12) against E3; and Arjunetin (Compound ID 10) against 6K were studied further in greater detail for their toxicity.

Toxicity analysis

The efficacy and unexpected toxicity of a drug to penetrate biological barriers, such as the intestinal wall or BBB, were considered as a prime determinant of the compounds taken forwards for toxicity tests. CHIKV is an old world virus, but is now seen to affect the CNS as well; therefore, compounds that were predicted to cross the BBB were also considered as potential antivirals. Of all the compounds considered for toxicity analysis using AdmetSAR, Arjunetin (Compound ID 10) was considered ineffective for oral consumption and is also carcinogenic. Also, Agnuside (Compound ID 18) and Mangiferin (Compound ID 122) were not considered as potential antivirals as they are predicted to have positive AMES toxicity (Table 5).

The compounds that were judged to be potential antivirals were Catechin-5-O-gallate (Compound ID 42) and Rosmarinic acid (Compound ID 151) against CP and Arjungenin (Compound ID 12) against E3 structural protein of CHIKV. Thus, the ligand/drug-protein interaction was studied for these three compounds to understand their interaction pattern and strength of interaction with the protein for their role as potential antivirals against CHIKV (Table 5).

Ligand protein interaction

A ligand protein interaction study was done for validated protein structures as discussed earlier. CP residues (Peptidase S3 domain) were predicted to bind to Catechin-5-O-gallate and Rosmarinic acid (Compound IDs 42 and 151, respectively) and E3 residues (Endopeptidase domain) bind to Arjungenin (Compound ID 12). The top docking conformation of Catechin-5-O-gallate showed a predicted binding energy of -6.26 kcal mol-1, whereas Rosmarinic acid and Arjungenin showed similar binding energy of -6.11 kcal mol-1 and -6.01 kcal mol-1, respectively. The binding energy (Glide Score) and the interaction energy (Potential, Vander Waals and Electrostatic) are shown in Table 3. The intermolecular hydrogen bonds and hydrophobic residues showing close contact between receptor proteins (CP and E3) and ligand (Compound ID 42, 151 and 12) are shown in Table 6 and Figure 6A–C, respectively.

Figure 6. Hydrogen bonding interactions between ligand and chikungunya virus proteins.

(A) Hydrogen bonding interaction between Catechin-5-O-gallate [CompID - 42] and capsid, binding affinity of - 6.26 kcal/mol was obtained. The zoomed region shows the receptor-binding pocket. Residues that form hydrogen bond interaction are Glu260 (Distance - 2.57 Å) and Lys 141 (Distance - 2.93 Å); His139, Val140, Asp161, Glu259 and Trp261 forms hydrophobic interaction. (B) Hydrogen bonding interaction between Rosmarinic acid [CompID - 151] and capsid, binding affinity of - 6.11 kcal/mol was obtained. The zoomed region shows the receptor-binding pocket. Residues that form hydrogen bond interaction are Glu260 (Distance - 1.71 and 1.66 Å), Trp261 (Distance - 2.04 and 1.93 Å) and Lys 141 (Distance - 2.37 and 1.99 Å); His139, Val140, Asp161 and Glu259 forms hydrophobic interaction. (C) Hydrogen bonding interaction between Arjungenin [CompID - 12] and E3, binding affinity of - 6.01 kcal/mol was obtained. The zoomed region shows the receptor-binding pocket. Residues that form hydrogen bond interaction are Gln52 (Distance - 3.11 Å) and Arg63 (Distance - 2.72 and 2.71 Å); Pro5, Ser18, Gln19, Gln49, Ala53, Ser58 and His60 forms hydrophobic interaction.

The interaction result showed that most of the hydrogen bond donors are from the protein that bind to the acceptors on the respective ligands. The compound Catechin-5-O-gallate (Compound ID 42) binds to Glu260 and Lys141 residues (HBond distance of 2.57 and 2.93 Å) of the CP protein and forms hydrophobic interactions with Asp161, His139, Val140 and Trp261 residues (Figure 7a). Further 2-D workspace revealed that when the ligand-protein interactions were observed both in the presence and absence of solvent the compound Catechin-5-O-gallate binds to the CP protein, HIS139 forms the hydrogen backbone; GLU259, GLU260, ASP161 form the hydrogen side chain. The ligand forms hydrophobic interactions with TRP261, VAL140, LYS141 (Figure 7b). We were unable to acquire the Ligplot for the interaction of Rosmarinic acid (Compound ID 151) with CP protein as the coordinates were undetectable; however, using 2-D workspace, we identified that Rosmarinic acid binds to the CP protein, TRP261 forms the hydrogen backbone; GLU260, LYS141 form the hydrogen side chain. The ligand forms hydrophobic interactions with HIS139, VAL140, GLU259, ASP161 (Figure 7b). The third compound Arjungenin (Compound ID 12), binds with Arg63 and Gln52 residues (HBond distance of 2.72 and 3.11 Å) of the E3 protein and Ser18, His60 and Ser58 residues are involved in hydrophobic interactions (Figure 7a). Its 2-D workspace revealed that SER18, GLN49 form the hydrogen backbone; GLN52, SER58, ARG63 form the hydrogen side chain. The ligand forms hydrophobic interactions with PRO5, GLN19, ALA53, HIS60 (Figure 7b). Overall docking and interaction results for the best three natural compounds have been compiled in Table 7.

Figure 7. Intermolecular hydrogen bonding in 2D view.

(A) LigPlot of Comp 42 (Capsid) and Comp 12 (E3). (B) Maestro ligand interaction diagram of Comp 42 and 151 (Capsid) and Comp 12 (E3).