An integrated approach of network-based systems biology, molecular docking, and molecular dynamics approach to unravel the role of existing antiviral molecules against AIDS-associated cancer

A serious challenge in cancer treatment is to reposition the activity of various already known drug candidates against cancer. There is a need to rewrite and systematically analyze the detailed mechanistic aspect of cellular networks to gain insight into the novel role played by various molecules. Most Human Immunodeficiency Virus infection-associated cancers are caused by oncogenic viruses like Human Papilloma Viruses and Epstein–Bar Virus. As the onset of AIDS-associated cancers marks the severity of AIDS, there might be possible interconnections between the targets and mechanism of both the diseases. We have explored the possibility of certain antiviral compounds to act against major AIDS-associated cancers: Kaposi’s Sarcoma, Non-Hodgkin Lymphoma, and Cervical Cancer with the help of systems pharmacology approach that includes screening for targets and molecules through the construction of a series of drug–target and drug–target–diseases network. Two molecules (Calanolide A and Chaetochromin B) and the target “HRAS” were finally screened with the help of molecular docking and molecular dynamics simulation. The results provide novel antiviral molecules against HRAS target to treat AIDS defining cancers and an insight for understanding the pharmacological, therapeutic aspects of similar unexplored molecules against various cancers.


Introduction
Some cancers are associated with AIDS. Cancer development in AIDS patients is generally due to the failed immune system. AIDS is a virus-associated opportunistic infection disease caused by the Human Immunodeficiency Virus (HIV). Patients with HIV infection have a high risk of developing certain types of cancer that is a vast group of various diseases, all involving uncontrolled cell division and growth. Certain types of cancer like Kaposi's sarcoma, Non-Hodgkin Lymphoma, and Cervical Cancer are the most common one and collectively termed as "AIDS defining Cancers" that confirm AIDS in HIV-infected person. As the onset of these AIDSdefining cancers alarm the severity of AIDS, there may be a possible connection between the targets of these diseases with some molecules (Akhter, Luthra, & Rajeswari, 2016;Board 10/2013;Casper, 2011;Cutolo, Basdevant, Bernadat, Bachelerie, & Ha-Duong, 2016;Karami & Jalili, 2015;Lu et al., 2016).
Many viruses like Epstein-Barr Virus, Hepatitis C Virus, Merkel Cell Polyoma Virus, Hepatitis B Virus, Human Papilloma Virus (HPV), Kaposi's Sarcoma Herpes Virus, Human T-lymphotropic Virus-1, etc. are known as oncogenic viruses (Morales-Sanchez & Fuentes-Panana, 2014). Thus, any infections with tumorassociated viruses can be linked to risk factors involved in developing cancers. Thus, we can presume existence of some common proteins having role in the development of AIDS and AIDS-defining cancers. As AIDS is a virus infection and linking viruses infection with the cancer progression is difficult, but certain facts like HPV infection is necessary for Cervical cancer, the existence of HIV transactivator (tat) protein with known oncogenic properties, the fact that 15% of total cancer worldwide is caused by viruses makes our presumption stronger (Bosch, Lorincz, Munoz, Meijer, & Shah, 2002;zur Hausen, 1991;Nunnari, Smith, & Daniel, 2008).
The biological molecules perform in collaboration with each other. Any outcome is the result of nature and characteristic of these associations. The identification of interactions between molecules and target proteins is the prime area in drug discovery. Network approaches in pharmacology are focused on arranging high-dimensional biological data-sets like chemical molecules and targets, to extract meaningful information and is the current trend in drug discovery. Docking is a very effective way to predict the possible binding of a ligand to its receptor, and there are a number of studies which have successfully implemented these docking studies (Ch, Jamil, & Raju, 2011;Jamil & Mustafa, 2012;Jayaraman & Jamil, 2014). Docking method can be used to screen a number of compounds library in a short time. These docking protocols can be accompanied by molecular dynamics (MD) simulation that is much more exhaustive and expensive but more accurate method for predicting and analyzing protein-ligand binding interactions.
In this article, systems pharmacology approach that includes target identification through a series of drugtarget (DT) and drug-target-disease (DTD) network and interaction analyzes in conjunction with molecular docking and MD simulation approach have been explored to estimate the possible correlations between the various targets responsible for AIDS-associated cancers with the known antiviral molecules that can act against AIDS and these cancers simultaneously.

Preparation of compound library
The data-set of antiviral compounds was prepared with the help of rigorous literature search; anti-HIV compounds were also included. Compound structures were downloaded from online chemical repositories like Pub-Chem (http://pubchem.ncbi.nlm.nih.gov/) and Drugbank (http://www.drugbank.ca/) in their two-dimensional Structure Data File (SDF) format. All the structures of these compounds were optimized by Merck Molecular Force Field (MMFF94) through Chemdraw v12 (Cousins, 2011;Law et al., 2014).

Preparation of target library
The identification of novel targets for the already known compounds is the urgent requirement for exploring novel functions. As traditional methods are time consuming and expensive, the computational approaches can be used to overcome these problems and to extract information for novel targets from the genome-wide scale (Kuruvilla, Shamji, Sternson, Hergenrother, & Schreiber, 2002;Li et al., 2012). A pharmacophore modeling technique was applied to extract the potential targets of each molecule. Based upon the pharmacophoric sites of the chemical potential targets were identified and ranked accordingly. The PharmMapper server was used for this task, which is supported by the information obtained from the database of PDTD, DrugBank, BindingDB, and TargetBanks (Liu et al., 2010). Combined algorithms of triangle hashing and genetic algorithm optimization were used sequentially to solve the ranking task of the molecule. The target was set to the human targets, while rest of the parameters were kept default (Liu et al., 2010). Each compound based on different conformations resulted in 300 targets, for each molecule top 50 targets were searched from PharmMapper server and out of them top 10 targets were screened after cross-validating their presence in human HIV targets (4431) present in NCBI's gene database.

Compound selection through DT network construction
The approach behind the modern drug discovery is to inhibit the primary target responsible for the disease; however, there are a number of targets responsible for a complex disease and merely altering an individual target may be insufficient to restore the healthy state of the system. Such cases require activity modulation of multiple targets to achieve the optimal therapeutic response (Li et al., 2012). Network biology can be a good approach to utilize single target or drug for multiple purposes (Kitano, 2007). The network approach will not only be helpful in exploring the complex mechanism of novel chemicals but will also unravel some novel interconnections between the already existing chemicals and targets.
A DT network with 132 molecules (D) and 138 targets (G) was constructed. In the network diagram, compounds and targets are represented by triangle shaped and circular shaped nodes, respectively. The chemical signals are the messenger to communicate between cells; cells grow, divide, and die according to the information encoded by the chemical. These informations are passed to each other through a process called signal transduction. This process signifies the association between two nodes of any signaling network, represented by edges in our network. Diseases are represented by diamond-shaped boxes. All networks were generated in Cytoscape v2.8.1 and v3.0, a software package for biological network visualization and analysis (Smoot, Ono, Ruscheinski, Wang, & Ideker, 2011).

Utilization of the network to prioritize targets (using network parameters)
All the targets of AIDS-associated cancers (110 targets for lymphomas; 561 targets of cervical cancer; 69 targets of Kaposi Sarcoma) were obtained from NCBI gene database (Accession date: 04-02-2015) and each target was individually matched with the 4431 targets fetched from NCBI's gene database. We enlisted all the targets common to HIV and Sarcoma; HIV and Lymphoma; HIV and Cervical cancer (Smoot et al., 2011).

Network analysis
The quantitative properties of these networks were analyzed by NetworkAnalyzer plugin and Centiscape plugin (Scardoni, Petterlini, & Laudanna, 2009). Analysis and visualization tools are needed to understand the functions of individual node masked by the overall network complexity. For each node in the network, NetworkAnalyzer estimates degree, clustering coefficient, and a number of other parameters. NetworkAnalyzer also calculates edge betweenness for each edge in the network. In undirected networks, the degree of a node corresponds to the number of edges linked to the node. The nodes with the highest degree and betweenness centrality score were screened.

Molecular docking
It is a robust approach which possesses a significant share in rational designing of drugs. The technique can not only predict the binding orientation of the molecules on the active site of their targets quickly and efficiently but can also estimate their binding affinity. To validate the target molecule link, docking was performed with the help of Autodock vina software. Its scoring function is based on advantages of knowledge-based potentials and empirical scoring functions. It is compatible with the PDBQT file format. The protein structure was directly downloaded from the RCSB protein data bank. Receptor file was prepared using AutoDock tools (ADT) (version 1.4.5). All the molecules in the library were in SDF format which were converted to PDBQT format using ADT software. Before proceeding for docking, the protein molecule was prepared with the help of ADT, hydrogen atoms were added, and water molecules were deleted. The auxiliary program Autogrid was utilized to generate the grid around the active site of the target protein. The docking area was defined by a 76 × 62 × 52 3D grid centered around the ligand-binding site with a default .375 Å grid space (O'Boyle et al., 2011;Trott & Olson, 2010).

MD simulation
In this work, we computed MD simulation studies on three protein-ligand complexes. MD simulations were performed using Gromacs 5.0.6 version on Linux platform (Bowers, Chow, et al., 2006;Pronk et al., 2013) with GROMOS 54A7 force field (Schmid et al., 2011). Prepared structures were processed through GROMACS and they were solvated in a Cubic box of SPC water molecules and neutralized using appropriate number (eight NA+) of counter ions. A distance of 1 nm was set between the box edges and protein complex to avoid direct interaction with its periodic image. The system was solvated and then subjected to energy minimization to a maximum of 5000 steps to remove steric clashes. The equilibrated systems were further carried to perform final MD simulations of 5,000,000 steps for 10,000 ps at the constant temperature of 300 K with a time step of 2 fs. The constant pressure of 1 atm was maintained by the Parrinello-Rahman algorithm. Trajectory snapshots were taken every 1 ps (Figure 1).

Results and discussion
Data-set preparation With reference to literature search 132 antiviral chemical constituents were searched and their structures were downloaded from various online repositories. Similarly the target data-set of 136 unique targets was prepared based on the pharmacophoric sites of the ligand molecules using Pharmmapper server. The list of fetched targets, molecules, and the symbols that are used to represent respective target and chemical molecule along with their identifier name are enlisted in the Table S1 (Supplementary material).

DT network construction
A DT network was constructed between 136 targets and 132 chemical molecules resulting into the network consisting of 268 nodes and 1293 edges.
Retaining targets common to AIDS and AIDSassociated cancers and refining the molecules Targets were further screened and only those targets common to HIV and Sarcoma (three targets); HIV and Lymphoma (four targets); HIV and Cervical cancer (30 targets) were retained. Now only those molecules associated with screened targets were kept. The molecules on the basis of their degree were filtered and 25 molecules showing four or more than four connections were finally screened. The screened molecules along with their number of associations and represented symbols are enlisted in Table 1.

Utilization of the network parameters to prioritize targets
The initial network was also analyzed with the help of NetworkAnalyser plugin to estimate the degree of targets. The nodes are represented with different colors depending upon the number of associations, the biggest and brightest colored nodes represent the most connected one (Figure 2). Various parameters like degree, betweenness centrality were used to screen the best target HRAS (G3). Figure 2 represents the Organic view of the DT network. The smaller and dark colored nodes represent less associated nodes. D84 and D113 nodes at lower left corner of the network show no associations as no targets were fetched from Pharmmapper server. The degree is one of the topological indexes that correspond to the number of nodes directly connected to a given node. Figure 2 also shows the plot between number of neighbors and betweenness centrality. Betweenness centrality score shows a score of around .15 highest for HRAS (G3).
The Figure 3(a) shows the DTD network with cart wheel structure, containing the 28 targets common to AIDS and associated cancers along with their linked molecules (25). After screening the targets common between the AIDS and its associated cancers, the most connected target was analyzed that results into same target HRAS (G3) shown in Figure 3(b).
The node with the highest degree score can be the regulatory hub and suggests the central and influential role of that node in signaling network. Another fundamental parameter, betweenness centrality was considered to assess the relevance of the position of the node in a network. The parameter is crucial to judge and maintain functionality and coherence of signaling mechanisms (Scardoni et al., 2009).

Molecular docking
Mutations occurring in the Ras gene are responsible for the uncontrolled cell growth and proliferation. Thirty percent of cancer is caused by mutation in Ras gene. Any lead molecule that prevents complex formation with Raf protein shall putatively inhibit the Ras-stimulated signal transduction pathway and thus can be a potential molecule to treat cancer (Jayakanthan et al., 2009). Ras protein initiates activation of c-Raf-1, MAP kinase, or ERK any mutation in the effector region, Tyr32-Tyr40 should inhibit its signal transducing activities (Sigal, Gibbs, D'Alonzo, & Scolnick, 1986;de Vries-Smits, Burgering, Leevers, Mar shall, & Bos, 1992). On the basis of current drug discovery principle "multiple target, multiple drug" the potential compounds interacting with the HRAS target were explored. The strength of interaction is the most effective parameter to rank different target-drug complexes. For this the final screened molecules (25) were docked with the HRAS target. The crystal structure of HRAS was downloaded from Protein Data Bank (PDB) with PDB ID: 3OIW (Ahmadian, Stege, Scheffzek, & Wittinghofer, 1997).  The Table 2 shows the binding affinity of the docked molecules. Molecules with score −7.0 kcal/mol or more were screened. Final three molecules were found with intact and stable binding associations with the target active site. Figure 4 clearly depicts the interactions between the respective molecules with the active site residues of HRAS protein. The Acteoside molecule is forming 10 hydrogen bonds; Chaetochromin B forms seven hydrogen bonds while Calanolide A is not forming any hydrogen bonds. As stated the main effector region of HRAS    MD simulation MD simulations were performed to obtain a better understanding of the dynamic behavior and stability of 3OIW and its ligand complex. Total of three MD simulations were performed for 10,000 ps each. The change in the structural integrity was analyzed by calculating the root-mean-square deviation (RMSD) as a function of period over backbone atoms. The RMSD of the Chaetochromin B complex was found to attain the stability around the 3 ns time period at around .9 nm ( Figure 5). The RMSD value of Calanolide A attained the stability after 1 ns timescale at around .4 nm while Acteoside RMSD values seem to fluctuate throughout the simulation period. It can be concluded that the ligands Chaetochromin B and Calanolide A seem to be stable inside the binding cavity.
Further, root-mean-square fluctuations (RMSF) were calculated to interpret the fluctuations of the protein structures at residue level. Figure 6 shows the residue fluctuations in 3OIW and its ligand complexes. It clearly depicts that Chaetochromin B protein complex has high residue fluctuations specifically in residues from 30 to 45 residues in comparison to other two complexes, which may have been resulted from the interaction between protein and ligand. Figure 7 shows the number of hydrogen-bond interactions between protein and ligand throughout the simulation period. The hydrogen bonds formed shows the strength of binding between the ligand and protein. The stability of the system can be accounted on the basis of RMSD values and H bonding results.
The potential energy values (Figure 8) show the stability of the complexes. Altogether, Chaetochromin B shows stable interaction, and possesses good binding potential toward HRAS target. The deviations in the results depict the adjustments made by the molecules while they remain bound within the active site.

Conclusion
In this study, a combined computational approach was applied to unravel the possibility of the selected antiviral compounds to act as anti-AIDS-associated cancers drug in alignment screening a number of AIDS-related targets as a target for these cancers. Ligands were screened in a stepwise manner at various stages. A number of criteria like affinity score, interactions, were used for screening ligands. Similarly, criteria like number of associations were used to screen targets. The network pharmacology approach was applied to estimate the number of associations between drug-targetsdiseases, to narrow down the search space of ligands as well as targets. MD simulation was applied to validate our findings; out of top three molecules, Chaetochromin B and Calanolide A formed stable complex with the screened target "HRAS" but only Chaetochromin B was able to bind specifically with the effector region. It was also followed that Chaetochromin B is present in natural sources, which also provide the possibility to explore other naturally derived drugs in similar reference or other. Conclusively, the molecules like Chaetochromin B may have the potency to treat AIDS-associated cancer: Kaposi's sarcoma, Cervical cancer, and Lymphoma, along with their already known anti-HIV activity; while HRAS protein might act as a putative target for these molecules. Obviously, experiments play a crucial role in validating simulation methodology and comparison of both the method will provide the exact scenario, but the present study is definitely a motivating aspect to work in the same direction.