Structure-Based Virtual Screening and Molecular Dynamic Simulation Studies to Identify Novel Cytochrome bc1 Inhibitors as Antimalarial Agents

Cytochrome bc1 (EC 1.10.2.2, bc1) is an essential component of the cellular respiratory chain, which catalyzes electron transfer from quinol to cytochrome c and concomitantly the translocation of protons across the membrane. It has been identified as a promising target in malaria parasites.The structure-based pharmacophore modelling and molecular dynamic simulation approach have been employed to identify novel inhibitors of cytochrome bc1. The best structure-based pharmacophore hypothesis (Hypo1) consists of one hydrogen bond acceptor (HBA), one general hydrophobic (HY), and two hydrophobic aromatic features (HYAr). Further, hydrogen interactions and hydrophobic interactions of known potent inhibitors with cytochrome bc1 were compared with Hypo1, which showed that the Hypo1 has good predictive ability. The validated Hypo1 was used to screen the chemical databases. The hits obtained were subsequently subjected to the molecular docking analysis to identify false-positive hits. Moreover, the molecular docking results were further validated by molecular dynamics simulations. Binding-free energy analysis using MM-GBSAmethod reveals that the van der Waals interactions and the electrostatic energy provide the basis for favorable absolute free energy of the complex. The five virtual hits were identified as possible candidates for the designing of potent cytochrome bc1 inhibitors.


Introduction
The cytochrome bc1 complex (EC 1.10.2.2, bc1) is a vital component of the cellular respiratory chain and the photosynthetic apparatus in photosynthetic bacteria [1].It is found on the plasma membrane of bacteria and in the inner mitochondrial membrane of eukaryotes [2].It consists of two heme groups, cytochrome b, iron-sulfur protein (ISP), with a Rieske-type Fe2S2 cluster and cytochrome c1 that undergoes reduction and oxidation during turnover of the enzyme [3].The role of cytochrome bc1 complex is to catalyze the electron transfer from quinol to a soluble cytochrome c (cyt c), and this electron transfer couples to the translocation of protons across the membrane [4].Due to its important role in the life cycle, inhibition of the bc1 complex has become an important target in the discovery of new antimalarial agent [5].
Atovaquone is a competitive inhibitor of the quinol oxidation site of cytochrome bc1.It has also shown potential role against pneumocystis pneumonia and toxoplasmosis in immunocompromised individuals [6][7][8].So far, two separate catalytic sites within the bc1 complex have been identified and confirmed by X-ray crystallographic studies: the quinol oxidation site (Q o site) and the quinone reduction site (Q i site) [9].Based on the specific binding interactions and conformational changes observed in both cyt b and the ISP domain into the bc1 complex, Q o inhibitors were further divided into two subgroups.subgroup I (Pf inhibitors) contains stigmatellin (SMA), famoxadone (FMX), and UHDBT, while subgroup II (Pm inhibitors) includes azoxystrobin (AZ) and methoxyacrylate-type inhibitors (MOA), KM, pyraclostrobin (PY), and many others [10].The ISP domain is still mobile after binding of subgroup II inhibitors but becomes fixed after the binding of subgroup I inhibitors [11].
In the present study, we aim to identify the structural features essential for inhibition of cytochrome bc1 and thereby design of potent Q o site inhibitors.Structure-based virtual screening and molecular dynamic simulations approaches have been applied to achieve this goal.The best structurebased qualitative pharmacophore hypothesis (Hypo1) was generated, validated, and used as a 3D query to screen Specs, NCI, and ChemDiv databases.The hits were subsequently filtered by docking study.Further, the hits obtained were validated by molecular dynamics simulations.Finally, five novel compounds with diverse scaffolds were identified as possible candidates for the designing of potent cytochrome bc1 inhibitors.

Data Sets.
The most important step in the pharmacophore modeling is the selection of suitable training set, responsible for determining the quality of the generated pharmacophore [12].In this study, structure-based pharmacophore hypothesis was generated using a potent inhibitor (FMX) and cytochrome bc1 complex.The known cytochrome bc1 inhibitors were collected from the literature and used as a test set for validation of developed pharmacophore hypotheses.

Pharmacophore Modelling.
Pharmacophore modeling studies were carried out using LigandScout3.0and Accelrys Discovery Studio2.5 (DS2.5)software package installed on IBM graphics workstation [13,14].From the available crystal structures of cytochrome bc1, one crystal structure was selected for generation of the structure-based pharmacophore models based on their IC 50 values of cocrystallized ligand, resolution, and source organism.LigandScout3.0,an automated tool for pharmacophore generation, was used to study the interaction between the inhibitors and amino acids in the active site of cytochrome bc1.It identifies protein ligand interactions such as hydrogen bond, charge transfer, and hydrophobic regions and also define excluded volume spheres based on the side chain atoms to characterize the inaccessible areas for any potential ligand.Initially, pharmacophore hypothesis was generated by Create Simplified Pharmacophore (Catalyst) module of LigandScout.The generated pharmacophore hypothesis was exported in .hypoeditformat and then converted into .chmformat using the Hypoedit tool in DS2.5 and used as a 3D query for the screening process [15].

Validation of Pharmacophore Model.
Validation of the developed pharmacophore model was done to determine its capability in differentiating active and the least active or inactive compounds and further performing virtual screening of databases.For the validation purpose two different methods, test set prediction and decoy test, were employed.In the test set predication, a set of 80 inhibitors of diverse chemical classes were prepared (see Supplementary   The test compounds were classified into the most active (IC 50 < 1000 nM) and the least active (IC 50 > 1000 nM) inhibitors based on their inhibitory activity.The molecules were sketched and minimized in Catalyst, and conformational analysis was carried out.For each molecule a maximum of 255 conformers (with an energy threshold of 4 kcal/mol) were generated and considered for model validation.Flexible method was employed for fitting the molecule to the pharmacophore hypothesis.This method ensures an exhaustive conformational mapping even for most complex molecules.Default values were used for all other parameters in the conformational analysis.All the test set inhibitors were mapped on the developed pharmacophore hypothesis, and FitValue was predicted for each inhibitor.Further, for evaluation purpose, test set inhibitors were divided into active and the least active based on the predicted FitValues.The hit rate of pharmacophore hypothesis was calculated to analyse efficiency of developed pharmacophore hypothesis in differentiating between active and the least active.
For further validation of generated hypothesis, decoy set was generated using DecoyFinder1.1 [20].Decoys are molecules that are supposed to be inactive against a target and used to validate the performance of the virtual screening workflow.To avoid bias in the enrichment factor calculation, decoys should resemble ligands physically, while still being chemically distinct from them.Decoys were selected if they are having similarity to that of the active ligands with respect to physical descriptors (molecular weight, number of rotational bonds, hydrogen bond donor count, hydrogen bond acceptor count, and the octanol-water partition coefficient) and deprived of the chemical descriptors to any of the active ligands [20].34 active cytochrome bc1 inhibitors were included in the database to calculate various statistical parameters such as accuracy, precision, sensitivity, specificity, goodness of hit score (GH), and enrichment factor (E value).Out of these, GH and E value are the two major parameters, playing significant role in identifying capability of the generated pharmacophore hypothesis [21].

Virtual Screening.
Virtual screening of chemical databases is a fast and accurate method, which helps to identify novel and potential leads suitable for further development.It has an advantage over any de novo design methods in which retrieved hits can be easily obtained for biological testing.Three commercially available databases of diverse chemical compounds were used in virtual screening.The databases were first screened for their drug-like properties using Lipinski's rule of five and subsequently submitted to DEREK for different toxicity filters.Fast/Flexible and Best/Flexible are the two databases searching options available in DS2.5.In our study, we performed virtual screening using Best/Flexible search option.The validated pharmacophore hypothesis Hypo1 was used as a 3D query in database screening.The hits obtained from pharmacophore based virtual screening were evaluated for drug-like properties using QED value.For calculation of QED value, physicochemical properties such as MW, ALOGP, number of HBDs, number of HBAs, molecular PSA, number of ROTBs, and the number of AROMs were calculated using the DS2.5.Finally, a substructure search was performed against each compound using a curate reference set of 94 functional moieties that are potentially mutagenic, are reactive, or have unfavorable pharmacokinetic properties.The number of matches for each compound was captured (ALERTS) [22].The unweighted and weighted QED values were calculated based on the previously mentioned molecular properties by using following formulae: where  is the derived desirability function corresponding to different molecular properties;  is the weight applied to each function;  is the number of molecular properties [23].
Hit compounds that passed all of these tests were taken for molecular docking analysis using Glide5.5.

Molecular Docking.
To investigate the detailed intermolecular interactions between the virtual hits and cytochrome bc1, an automated docking program Glide5.5 was used [24].Three-dimensional structure information about the target protein was taken from the protein data bank (PDB ID: 1L0L).Cocrystallized ligand (Famoxadone) was docked into cytochrome bc1 to validate the docking protocol.The protocol followed for docking studies of virtual hits included processing of the protein and ligand preparation.During protein preparation, ligand molecules were deleted, hydrogen atoms were added, solvent molecules were deleted; and bond orders for crystal protein were adjusted and minimized up to 0.30 Å RMSD [25,26].An active site of 10 Å was created around the cocrystallized ligand.The extraprecision (XP) mode and other default parameters of Glide software were used for the docking studies.

Molecular Dynamics Simulation.
The screened virtual hits were further validated, for its binding affinity with the active site residue of cytochrome bc1, by performing the molecular dynamics (MD) simulations.MD studies of the cytochrome bc1 complex with cocrystallized ligand (Famoxadone), and top two virtual hits were performed using AMBER11 [27].The input files for the MD simulations were prepared in tleap module of AMBER, in which the system was hydrated by TIP3P water box of 10 Å and neutralized by counter ions.The minimization of the initial solvated system was performed in two steps; in the first step, only water was permitted to move employing 500 cycles of the steepest descent followed by 500 cycles of conjugate gradient method.In the second step, the entire system was subjected for minimization using 500 cycles of the steepest descent followed by 500 cycles of conjugated gradient method.The heating of the system was carried out at 310 K using Langevin thermostat.SHAKE algorithm was implemented for hydrogen constraints.The force of 2 kcal/mol/ Å2 on the heavy atoms of the complex was restrained during density equilibration.Whole system equilibration was performed over a period of 1000 ps.Production phase of 5000 ps was carried out using NPT ensemble at 310 K and one atmospheric pressure.Nonbonded cutoff of 8 Å and step size of 2 fs was used for the simulations.Furthermore, the binding affinity of the cocrystallized ligand and top two virtual hits were calculated by the molecular mechanics generalized Born surface area (MM-GBSA) method integrated into AMBER11 [28].Binding-free energy (Δ Bind ) of inhibitors at the binding site was calculated by using ( 2) where  Complex ,  Protein , and  Ligand represent the free energies of complex, protein, and the ligand, respectively.Free energy () of each state was calculated using the following equations: where  MM represents the molecular mechanical energy;  GB and  SA are the contributions from polar and nonpolar terms of the free energy of the solvent continuum, and TS represents the entropic contribution of the solute. MM was obtained using (4), in which  ele is electrostatic energy,  vdW is van der Waals energy, and  int is internal energy. GB was calculated based on the generalized Born model. SA was computed by the molsurf module of AMBER according to (5), which is proportional to the solvent accessible surface area (SASA):

Results and Discussion
3.1.Generation of Pharmacophore Model.The overall workflow followed during identification of novel cytochrome bc1 inhibitors is as shown in Figure 1.The reported crystal structures of cytochrome bc1 from the protein data bank were analyzed for organism source, conformations, bound inhibitors, and their activity.Finally, one cocrystallized structure of cytochrome bc1 was selected for development of structurebased pharmacophore hypothesis.The binding site analysis was performed for better understanding of the specificity and pharmacophore requirement of inhibitors.The binding site was characterized by several direct interactions such as hydrogen bond interactions, - interaction, and lipophilic side chains complement, which reflects the hydrophobic nature of an active site.For structure-based pharmacophore studies, we have used LigandScout3.0,which provides the interactions between protein and ligand as well as some excluded volume spheres corresponding to their 3D structures of protein.In this study, the coordinates of cytochrome bc1 complex (PDB ID: 1L0L) and insight from molecular docking analysis of known cytochrome bc1 inhibitors were used for generation of structure-based pharmacophore hypothesis.In developed structure-based pharmacophore hypothesis, one HBA was pointed towards Glu271.Also, molecular docking analysis reveals that the oxygen atom of most potent inhibitors (IC 50 value < 100 nM) was accepting an H-bond from NH of Glu271.Hence, HBA was considered as an important chemical feature to identify novel cytochrome bc1 inhibitors.Finally, four-feature structure-based hypothesis (Hypo1) was generated, which consists of one HBA, one general hydrophobic (HY), and two hydrophobic aromatic features (HYAr).

Pharmacophore Model Validation.
Figure 2 shows the chemical features of a pharmacophore hypothesis (Hypo1) with their interfeature distance constraints.The test set prediction was performed to validate the developed pharmacophore hypothesis.The Hypo1 pharmacophore hypothesis was used to estimate the FitValue of test set compounds.The inhibitors showing FitValue higher than 2.38 considered as active and the rest were considered as the least active.The hit rate of pharmacophore hypothesis suggested that the developed pharmacophore hypothesis can differentiate between active and the least active.In detail, 50 of 57 highly active and 18 of 23 the least active compounds were predicted correctly by the Hypo1.Seven active compounds were underestimated as the least active, whereas five compounds are overestimated as active by Hypo1. Figure 3 shows the mapping of the most active and the least active compound from the test set over Hypo1.The highly active compound was mapped accurately to all features of Hypo1.While, in case of the least active compounds, one HBA was missing and two HYAr features were mapped partially.Finally, DecoyFinder1.1 was employed to generate small database () containing 1258 compounds, which includes 34  actives and 1224 decoys.This database was used to validate the generated hypothesis (Hypo1), whether it is capable of distinguishing the actives from decoys or not.The hypothesis Hypo1 was used as a 3D structural query to perform screening of decoy database and then accuracy, precision, sensitivity, and specificity were calculated.The hits were selected based on the FitValue.Furthermore, for the analysis of results,  Total number of actives in database () 3 4 3Total hits (  ) 4 2 4Active hits (TP) the enrichment factor ( value) and goodness of hit score (GH) were calculated using the following formulae: where , , Ht, and TP represent the total number of compounds of the database, total number of actives, total number of compounds screened by a pharmacophore model, and total number of active compounds screened, respectively [29].Hypo1 has shown an  value of 26.42.The calculated GH score was greater than 0.5, specifies that the quality of developed pharmacophore was significant (Table 1).From the overall validation results, we assure that the Hypo1 hypothesis could discriminate between the active and the least active compounds.Hence, we have used Hypo1 hypothesis for virtual screening to select or discriminate the suitable cytochrome bc1 inhibitors.

Sequential Virtual
Screening.Sequential virtual screening was performed as depicted in the flowchart (Figure 4).The NCI, ChemDiv, and Specs databases were screened using Lipinski's rule of five, which comprised 87374, 843113, and 276807 compounds, respectively.Further, the hits obtained were screened for several toxicity filters, such as carcinogenicity, chromosome damage, genotoxicity, hERG channel inhibition, hepatotoxicity, mutagenicity, and thyroid toxicity using DEREK.After applying toxicity filter 622196 compounds were obtained and subsequently screened by the validated pharmacophore hypothesis (Hypo1).381734 compounds were mapped to the Hypo1 pharmacophore hypothesis, which included some structurally similar compounds to that of existing cytochrome bc1 inhibitors, and some novel scaffolds were also identified.A set of 340 hit compounds from Hypo1 screening were selected based on the estimated FitValue > 2.50.Out of these, 236 compounds having QED value > 0.50 were selected and subjected to further analysis by molecular docking to avoid the false-positive hits from virtual screening.

Molecular Docking Studies.
Recently in molecular dynamic study, Zhao et al. have reported that the side chains of hydrophobic residues of the active site showed the conformational stability, except the phenyl group of Phe274, which exhibited the significant conformational flexibility in different complexes [1].The conformational flexibility of phenyl group of Phe274 allows optimizing the corresponding - interactions.The phenyl ring of the cyanophenoyl group is almost vertical to the pyrimidyl ring due to the presence of the linear cyano group in AZ.As a result, the pyrimidyl ring is not able to form ideal - interaction parallel with the phenyl ring of Phe274.On the contrary, myxothiazol is more potent than AZ because its thiazole ring parallels well with the phenyl group of Phe274 and from the ideal - interaction [1].These facts suggested that the - interaction between inhibitors and the phenyl of Phe274 is essential for inhibitory activity.Thus - interaction with Phe274 is considered as essential during molecular docking analysis.Docking studies were carried out to identify false-positive hits from virtual screening process and to explore the interaction between the virtual hits and cytochrome bc1.All the test set compounds along with virtual hits were docked into the active site of cytochrome bc1 (PDB ID: 1L0L) using Glide program.Glide score, which differentiates compounds based on their interacting pattern, is calculated for all molecules.The cocrystallized ligand famoxadone (FMX) was docked into the active site of cytochrome bc1 to validate the docking protocol with RMSD of 0.320 Å. FMX contains three aromatic rings, the central five-membered oxazolidinedione ring carries a phenyl amino group and a phenoxyphenyl group.In detail, it forms one strong hydrogen bond with an amide group of the conserved Glu271.The hydrophobic environment is formed by the Phe274, Tyr131, Tyr273, Phe91, and Tyr278 allowing formation of a network of - interactions with the two phenyl groups of famoxadone.Top virtual hits were showing similar interaction pattern to that of FMX along with some additional hydrogen bonds and hydrophobic interactions.The binding interaction pattern observed during docking studies was complementary with that of the hydrogen bond acceptor and hydrophobic features of pharmacophore hypotheses.On the basis of molecular docking analysis five virtual hits (Figure 5) with diverse scaffolds were identified as potential cytochrome bc1 inhibitors.Figure 6 shows 2D interaction pose of FMX and top virtual hit at the active site of cytochrome bc1.The Glide score, FitValue, and QED value of top five virtual hits are shown in Table 2.
Further, these five hit compounds were confirmed by the PubChem [30] and SciFinder [31] scholar search tools, as unreported for cytochrome bc1 inhibition.Hence, we suggest that these virtual hits with diverse scaffolds are novel as cytochrome bc1 inhibitors.

Molecular Dynamic Simulation.
In order to verify whether the results obtained by the molecular docking analysis were robust or a fortuitous, we have carried out molecular dynamics simulations with cocrystallized ligand (FMX) and top two virtual hits.The atomic RMSDs of the C, C, and N atoms of the protein and the ligands obtained during the trajectories and the initial structures were monitored during the simulations.The sharp rise was observed up to first 50 ps, and the function remained stable for the rest of the simulation.The system reached a constant temperature.The total energy was fluctuating around an average energy with the system stabilization at a given temperature of 310 K and 1 atmospheric pressure.The RMSD of the molecular dynamics simulation was stable after 1 ns of the equilibrium.Both protein-ligand trajectories Further, insights into the forces involved in the binding of inhibitors were obtained by analysis of the MM-GBSA-free energy contributions, as listed in Table 3.The calculated ΔG binding for the cytochrome bc1 with top two virtual hit complexes was slightly higher than the value of the cytochrome bc1/and cocrystallized ligand complexes, demonstrating that the potential activity of top virtual hit is as high as the known potent inhibitors.Both, the intermolecular van der Waals and the electrostatic interaction energies have shown significant contributions in the binding, whereas polar solvation terms are calculated by GB counteract binding.Nonpolar solvation terms that correspond to the burial of the solvent-accessible surface area contribute slightly in binding.Comparing the van der Waals contributions with the electrostatic contributions, we noticed that the association between cytochrome bc1 and inhibitors is mainly driven by more favorable van der Waals interaction in the complex than in solution.The van der Waals/nonpolar contributions of the cytochrome bc1/inhibitors reveal that inhibitor fits more tightly within the active-site cavity.Taking all energy contributions into account, virtual hits present slightly higher potential activity than cocrystallized ligand.Hence, inspired by the molecular dynamics studies, we conclude a possible strategy to design a new series of cytochrome bc1 inhibitors with high potential activity, which is to retain the important hydrogen bond and enhance the hydrophobic interactions with the hydrophobic residues in the pocket.In our study, top virtual hits exhibit the most favorable affinity with cytochrome bc1, because it obeys the previously mentioned strategy that it retains important hydrogen bonds with Glu271 and increases the hydrophobic interactions with receptor.

Conclusions
In this study, we have developed a structure-based pharmacophore model, which was validated to evaluate its predicting power over the diverse test set compounds.The best structure-based pharmacophore hypothesis consists of one HBA, one HY, and two HYAr features.This predictive hypothesis was further used in virtual screening to identify novel cytochrome bc1 inhibitors.Molecular docking analysis reveals that hydrogen bond interaction with Glu271 is essential for cytochrome bc1 inhibition.Finally, compounds were selected as positive hits based on the molecular docking analysis.Novelty of these hits was confirmed with PubChem and SciFinder scholar search tool.Molecular dynamics simulation of cytochrome bc1 with cocrystallized ligands and top two virtual hits was performed for getting deeper insights into the driving forces between the inhibitors and cytochrome bc1.Free energy calculations were done by MM-GBSA method.Based on the observation, a new strategy is proposed to design a new series of cytochrome bc1 inhibitors.It includes the retention of hydrogen bond interaction with Glu271 and increasing hydrophobic interaction.Combining all these results, five novel compounds were identified as possible lead candidates to be used as a novel and potent cytochrome bc1 inhibitors.Further, in vitro testing of the virtual hits will confirm the success rate of this study to optimize the hits thereafter.

Figure 1 :
Figure 1: The overall work flow followed during identification of novel cytochrome bc1 inhibitors.

Figure 2 :
Figure 2: Chemical features of the structure-based pharmacophore hypothesis with interfeature distance constraints.

Figure 3 :
Figure 3: The most active and the least active compound from the test set aligned with pharmacophore hypotheses.

Figure 4 :
Figure4: The sequential virtual screening followed during identification of novel cytochrome bc1 inhibitors.

Figure 5 :Figure 6 :
Figure 5: Chemical structures of top five virtual hits obtained through sequential virtual screening using Hypo1.

Table 1 :
The statistical parameters obtained from decoy test.

Table 2 :
The Glide Score, FitValue, and QED values of top five virtual hits.

Table 3 :
Binding-free energy and individual energy components calculated using MM-GBSA method for cocrystallized ligands and top two virtual hits.