Identification of core therapeutic targets for Monkeypox virus and repurposing potential of drugs against them: An in silico approach

2.1. Global dataset of publicly available mpox virus genomes and quality assessment

All publicly available mpox virus genome assemblies were downloaded from the National Centre for Biotechnology Information (NCBI) on 3^rd June 2022.

2.2. Annotation of mpox virus genomes

To avoid artificial differences resulting from different genome annotation pipelines, the genome assemblies downloaded from GenBank were reannotated with the Prokka v 1.14.6 [43] with kingdom = Virus, genus = orthopox virus and species = ‘Monkeypox virus’.

2.3. Identification and re-annotation of mpox virus core genome

Clustering orthologs in viruses can ease the screening for drug candidates [44,45]. Therefore, to obtain the core proteome of mpox virus, the proteomes from all the genomes were clustered into groups of orthologs using the OrthoFinder v2.5.4 [46,47] with default parameters. The OrthoFinder surmises homologous regions and determines the orthogroups (OG's) by using a combination of the BLAST search and the Markov Cluster Algorithm (MCA) [47]. For this study, we defined core orthologs as ORFs present in all the strains of mpox virus. Therefore, only the ORFs shared by all (100%) the strains of mpox virus were considered as core orthologs. Further, to improve the functional annotations of all core proteins, we searched the consensus sequence of these proteins against the NCBI, UniProt and KEGG databases.

2.4. Subtractive proteomics and identification of druggable core proteome of mpox virus

We devised a suitable subtractive proteomics strategy [48] to pinpoint the highly druggable core proteins of mpox virus. We did so by screening out undesirable proteins from the core proteomic dataset of mpox virus by applying a set of rationally chosen parameters in the order listed below.

2.5. Physicochemical characterization of druggable core proteins

Physicochemical properties of the screened druggable enzymes were further assessed as these features have a critical role in the identification of potential druggable targets in the early stages of drug development. Physicochemical properties like isoelectric point (pI), molecular weight, aliphatic index, instability index, extinction co-efficient and residues accessibility were computed with the ProtParam tools (https://www.expasy.org/resources/protparam).

2.6. Collection of FDA approved/investigational drugs

For high-throughput virtual screening, we prepared a non-redundant dataset of FDA approved/investigational drugs. The Probes & Drugs portal (P&D) is a collection of high-quality bioactive compounds including drugs and probes with their targets and experimental values collected from different sources and updated on a daily basis [66]. We downloaded the approved/investigational drugs from the Probes & Drugs portal (P&D) on 13^th July 2022. This dataset included molecules approved for clinical use, clinical testing, veterinary use or neutraceutical applications.

2.7. In silico REOS and PAINS filtering

The SMILES of all the molecules were subjected to rapid elimination of swill (REOS) and pan-assay interferences (PAINS) filter using the RDKit [67], a python-based software (https://rdkit.org). Further, we implemented our in-house workflows using the Konstanz Information Miner (KNIME, v4.2.3; http://knime.org) followed by removal of redundancy.

2.8. Ligand preparation

All the non-redundant ligands were pre-processed using the OPLS_2005 force field in the LigPrep v5.3 module to obtain the accurate 3D energy minimized Lewis structures. Furthermore, EpiK v5.3 was used to generate the subsequent broad chemical, structurally diversified stereoisomers at pH 7.0 ± 2.0.

2.9. Virtual screening of FDA approved/investigational small molecules in the highly druggable binding sites

The Glide program and its virtual screening workflow were applied to identify hits against the highly druggable target proteins by incorporating three docking protocols; high throughput virtual screening (HTVS), Standard Precision (SP) mode and Extra Precision (XP) mode [68,69]. A grid box of size 12 Å was assigned around the centroid of the binding sites predicted for each target structure. The van der Waals radii (1.00 Å) of the receptor along with partial atomic charge (0.25) was kept unchanged. Default docking parameters were used, and no constraints were included. For screening curated FDA library against all the targets, we used the xglide module of the Schrodinger suite at every step of the glide docking protocol. All the docked poses of each drug were transferred from HTVS to SP to filter out the false-positive results. The SP docked poses were then piped to XP docking where the false positives were further eliminated based upon ligand-receptor shape complementarity. These XP docked poses were then ranked according to a more stringent scoring function i.e., XP glide docking score (XP GScore).

2.10. Molecular dynamics simulations: system preparation and data generation

To evaluate the binding stability of the selected common top hits at the predicted highly druggable active sites of the non-host homologous core proteins (enzymes), we employed the molecular dynamics (MD) simulations of the protein-ligand complexes in the water environment. Firstly, we prepared the conformational ensembles of all four apo proteins over the course of 400 ns MD simulation at a temperature of 300 K followed by simulation of the protein-ligand complexes over the course of 150 ns at the same temperature. The GROMACS suite v2021.1 [70] with charmm36-jul2021 force field was used to perform the MD simulation of apo proteins for the generation of ensembles and evaluation of the stability of bound ligands. The topology of proteins was generated by the in-built module of GROMACS and the topology of ligands was generated using the CHARMM General Force Field (CGenFF) web server (https://cgenff.umaryland.edu). Each system was solvated in a cubic periodic box with Tip3P charm-modified water model by keeping a 12 Å distance between the system and the water box. All the systems were neutralized by adding NaCl salt of concentration 0.15 M.

The structure of each system was minimized through a maximum of 50,000 steps using the steepest descent method of energy minimization. The equilibrations of each system were performed in two phases i.e. NVT (constant number of particles, volume and temperature) and NPT (constant number of particles, pressure and temperature) ensembles over the course of 5 ns each followed by MD production at contant temperature and pressure (1 bar and 300K). For ensemble simulation and MD production, each whole system was divided into two temperature-coupling groups of protein/protein-ligand and water-ions respectively. LINear Constraint Solver (LINCS) algorithm [71] was used to constraint all the bonds and isolated angles. Similarly, the velocity Verlet algorithm was used to integrate Newton's equations of motion. The temperature and pressure coupling on each ensemble was performed by V-rescale (modified Berendsen thermostat), Berendsen (NPT) and Parrinello-Rahman methods. However, no pressure coupling was used in NVT ensembles. For calculating the long-range (Electrostatic) interactions, the Particle Mesh Ewald (PME) method with a 1.2 nm cutoff was used. The output trajectories were recorderd every 10 ps for subsequent analysis.

2.11. Molecular dynamics simulations: stability, clustering, essential dynamics and binding free energy analysis

Following MD data productions, the trajectories of all the systems was compared by analysis of root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyrations (Rg), solvent accessible surface area (SASA) and amount of hydrogens bonds (HB) using the in-built module of the GROMACS. The analysis of positional changes of secondary structure elements (SSE) was accomplished by in-built module of biotite python library [72]. The dominant conformations of the trajectory data was inferenced by clustering Cα atoms RMSD of all the systems using TTClust [73]. The optimum cluster numbers was identified by elbow method. To understand the structural changes at Cα atoms level over 150 ns time of MD simulation, MODE-TASK tool were used the study the principal component analysis (PCA) using singular value decomposition (SVD) of the cartesian coordinates as well as internal coordinates [74,75]. The first two eigenvectors (principal component 1 and 2) of PCA was plotted using the ggplot2 R package [76]. To compare the collective motions of Cα atoms over the simulation time, the first essential principal component mode of PCA was represented as porcupine plots using Normal Mode Wizard (NMWiz), a ProDy python package based plugin in VMD [77,78]. The binding free energies and energies decomposition of key-interacting residues for all the complexes were calculated by AMBER's MMPBSA.py engine based gmx_MMPBSA tool [79,80] using molecular mechanics/Poisson-Boltzman surface area (MM/PBSA) model. The analysis module of the gmx_MMPBSA tool were used to infer and visualize the energies data. The average distance of bound ligands from the key-interacting residues, obtained from MM/PBSA analysis, were then calculated using the in-built gmx distance dismodule of the GROMACS.

3.1. The meaning of the open reading frame (ORF) and core proteins in this study

A clarification is necessary on what we mean by an ORF. This is because there is a slight discrepancy in how different authors like to define the term “Open Reading Frame” or “ORF”. In the bioinformatics community it is a standard practice to identify an ORF such that the evidence of translation is not required, a definition that virological community may be less acquainted with. For the purpose of this study, we defined an ORF having start codon, stop codon and no internal stop codon as a “protein-coding gene” which translates into a functional protein that contributes to viral replication, transmission, immune evasion or overall fitness. However, we also acknowledge that translation itself may not be a criterion significant enough for an ORF to be a protein-coding gene. We acknowledge that an act of translation may have a function even if the peptide produced is non-functional as in the case of regulatory uORFs [81]. Therefore, each predicted ORF can be important even if it is annotated as a hypothetical protein.

Likewise, the term “core proteins” used here frequently has been used to describe those proteins that are highly conserved across all the strains of mpox virus. This clarification is necessary as some members of virology community may be more accustomed to the term “core protein” as the proteins associated directly with the nucleocapsid.

3.2. A compendium of 125 mpox virus genomes was developed and explored

All the genome assemblies of mpox virus available as of 3^rd June 2022, were retrieved from the NCBI. Our dataset included 24 newly assembled mpox virus genomes from the latest mpox outbreak described in non-endemic countries. In total, our final dataset comprised of 125 complete genome assemblies of mpox virus. The mean genome size of mpox virus was 196318 bp (∼0.2 Mb). The highest number (275) of ORFs was identified in the strain mpox_FRA_2022_TLS67 (ON602722) whereas the strain mpox_Nig_2017_298464 (MG693724) harboured the least number of ORFs (193). The mean ORFs per genome was ∼213. Due to lack of pre-existing annotations, approximately half of the total ORFs could be successfully annotated whereas the other remaining half remained unannotated or hypothetical. The detailed metadata regarding the sample isolation information and genome statistics like genome size, GC content (%), gene count and completeness etc have been summarised in Supplementary Table S1.

3.3. Core proteome assessment identified sixty-nine core proteins from mpox virus genomes

Core proteins are major contributors to the replication and survival of a virus. Hence, developing novel therapeutics against such proteins could stop infectious pathways. Moreover, screening core proteins of an organism are advantageous as it helps to reduce the dependency on genome completeness [82]. In our core proteome assessment of 125 mpox virus strains, OrthoFinder assigned 26293 genes (99.4% of the total) to 294 ortho groups with 166 unassigned genes (0.6% of genes). Fifty percent of all genes were in orthogroups with 124 or more genes (G50 = 124) and were contained in the largest 101 ortho groups (O50 = 101). There were 69 ortho groups (core ortho groups) with all species present and 58 of these consisted entirely of single-copy genes. Two ortho groups were species-specific ortho groups with only 4 genes. The gene count in every orthogroup for each strain and strain-wise statistics is presented in Supplementary Tables S2 and S3. Every set of the 69 core proteins was then aligned to create a consensus protein sequence of each ortholog. Thus, in doing so, we created a dataset comprising the core proteome of all the 125 mpox virus genomes. Unfortunately, many predicted core proteins remained unannotated for the reasons mentioned previously. To overcome this problem, we reannotated every protein of our core proteomic dataset through a series of extensive manual curations. The detailed annotation of the core orthogroups dataset has been presented in Table 1 .

3.4. Subtractive proteomics successfully identifies four highly druggable, non-host homologous core proteins (enzymes)

The reannotated 69 core proteins of mpox virus were subjected to a set of multi-stage screenings to identify the druggable core proteins of mpox virus as shown in Table 2 . In the first stage, 56 globular proteins were identified based on their sub-cellular localization. These are most likely to be comprised of the proteins associated with the nucleocapsid, transcription factors and enzymes that may be critical to the virus lifecycle. In the second stage, we further screened out enzymes from the set of 56 globular core proteins. From the in silico assessment, we identified 23 core globular proteins which are most likely enzymes. Subsequently, in the third stage, the shared homology with human host proteins was assessed to limit the off-target issues. Eight enzymes namely, DNA-dependent RNA polymerase 132 kDa subunit (Rpo132), Phospholipase-D-like protein K4 (K4L), Ribonucleoside-diphosphate reductase small chain (OPG048), dual specificity protein phosphatase H1 (H1L), 3-beta-hydroxy-Delta(5)-steroid dehydrogenase (SALF7L), Cu–Zn superoxide dismutase-like protein (A46R), thymidylate kinase (TMK) and pseudokinase B12 (B12) shared 24.95%, 47.23%, 81.25%, 28.78%, 41.24%, 30%, 41.55% and 32.88% sequence identity with human proteins DNA-directed RNA polymerase II 140 kDa polypeptide (3J0K_B), 5′-3′ exonuclease PLD3 (Q8IV08), Human ribonucleotide reductase subunit R2 (D6W4Z6), Dual specificity protein phosphatase 22 (6LVQ_A), 3 beta-hydroxysteroid dehydrogenase type 7 (Q9H2F3), Superoxide dismutase [Cu–Zn] (3GTV_A), Thymidylate Kinase (1E98_A) and Vaccinia-related kinase 1 (2LAV_A) respectively. These proteins were subsequently discarded as they exceeded the permitted threshold (e-value = 0.0001). The 15 remaining non-host homologous enzymes were processed for druggability studies. 12 out of these 15 enzymes shared excellent coverages and sequence similarities with the solved crystal structures of proteins from the PDB database. On the contrary, DNA polymerase processivity factor component A20 (A20R) shared 96.75% similarity with 6ZXP_A at very low coverage (28%) whereas crystal structures even remotely similar to core cysteine protease (I7L) and metalloendopeptidase G1 (G1L) were not available (Table 1). In the absence of experimentally derived high-resolution structures, modeled proteins present suitable alternatives to examine structural perturbations in proteins. Therefore, to overcome this problem we applied the hybrid, a top performer of CASP-5 and highly accurate deep learning-based modeling method the RoseTTAFold to model the 3D structures of all the 15 non-host homologous enzymes. The RoseTTAFold [55] employs a 3-track neural network to fold a protein sequence into a model. The accuracy and validation assessment of all the predicted structures were assessed through the MolProbity web server [56]. Out of these 15, the 7 enzymes successfully modeled by the RoseTTAFold were subjected to druggability assessments for the time being. Of these 7 core non-host homologous enzymes, only 4 enzymes namely the I7L, DNA topoisomerase 1B (Top1B), early transcription factor 70 kDa small subunit (VETFS) and A20R complied with the thresholds of druggability prediction tools. The druggability scores (DScore), simple score (SScore), binding affinity (pKd), volume and surface area associated with the predicted highly druggable binding site of the four proteins have been summarised in Supplementary Table S4. The model quality statistics of these four druggable proteins are presented in Table 3 . Unfortunately, this study was impelled to abandon the druggability assessment of 8 out of these 15 enzymes (Group 5, Table 2). These 8 proteins will be assessed for their druggability in our future endeavors. Hence, by combining genomics and an extensive subtractive proteomic screening on the core proteomic dataset of mpox virus we finally identified 4 highly druggable, globular, non-host homologous core proteins (Table 2). The modeled proteins and their predicted druggable sites have been illustrated in Fig. 2 .

3.5. The druggable proteins are soluble and highly stable over a wide range of temperatures

The physicochemical properties of the 4 highly druggable proteins were predicted (Table 3). Proteins with instability index above 40 are considered unstable whereas proteins with an instability index below 40 are considered stable. The instability value of Top1B was slightly above 40. Therefore, this enzyme could be slightly unstable. The instability values of the remaining 3 proteins were below 40 and reflect their highly stable nature. This prediction was also supported by the high aliphatic index of these proteins which suggests their stability over a wide range of temperatures. Similarly, the hydrophilic features and solubility of a protein are reflected in its grand averages of hydropathy (GRAVY) score. The GRAVY scores of all 4 enzymes were below 0 which reflects upon their globular and soluble nature.

3.6. A curated library of 5893 FDA approved/investigational drugs was prepared

For high-throughput virtual screening (HTVS), we collected ∼19 thousand FDA-approved/investigational drugs of approved, investigational and experimental categories from the Probes and Drugs Portal, a highly curated chemical space portal comprising data from ∼50 different sources [66]. Discontinued FDA-approved/investigational drugs were excluded from our selection. For filtering problematic functional groups and false positives, we implemented the REOS and PAINS filters on the collected library. The PAINS and REOS filter are both based on the RDKit substructure counter and compare the substructures present in the input database with a list of problematic functional groups. After filtering, ∼32% (5893) of drugs remained in the dataset and were subjected to the ligand preparation stage. Finally, a set of 12,048 stereoisomers corresponding to 5893 drugs were obtained following library preparation with the LigPrep module of the Schrodinger suite.

3.7. Burixafor, batefenterol and eluxadoline are the top hits common to the four highly druggable enzymes of mpox virus

We performed the multi-step (HTVS, SP and XP) docking of 12,048 stereoisomers corresponding to 5893 drugs against all the identified target proteins by using the Schrodinger suite's Glide module. Numerous binding modes of all the drugs across the 4 targets were predicted. The docking scores varied between −13.36 kcal/mol to −0.50 kcal/mol for core protease I7L, −10.74 kcal/mol to 2.19 kcal/mol for Top1B, −11.55 kcal/mol to 0.75 kcal/mol for VETFS and −12.30 kcal/mol to 0.35 kcal/mol for A20R as shown in Fig. 3 . Interestingly, the drug batefenterol (BAT) was among the top 30 hits for all the proteins namely; A20R, I7L, Top1B and VETFS whereas burixafor (BUR) and eluxadoline (ELU) were amongst the top 30 hits for I7L, Top1B and VETFS. The binding affinities of BAT with A20R, I7L, VETFS and Top1B ranged between −12.07 kcal/mol to −8.89 kcal/mol. The binding affinities of BUR with I7L, VETFS and Top1B ranged between −11.84 kcal/mol to −8.29 kcal/mol. Similarly, the binding affinities of ELU with I7L, VETFS and Top1B ranged between −12.73 kcal/mol to −8.14 kcal/mol respectively. Likewise, tobramycin (PD001728), dibekacin (PD001070) and GLPG-0187 (PD058153) were identified as hits common only to I7L and VETFS and were among their top 30 hits. The XP GScore, XP HBond energy and ligand efficiency of the top 30 hits for each target is presented in Supplementary Tables S5–S8. The primary indications of the top five hits against all the targets is summarised in Table 4 .

Table 4.

Mode of actions, application and 2D structure of shortlisted top five docked drugs with high binding affinity for each target. The mode of actions, applications was collected from DrugBank, Inxight Drugs portal and PubChem database.

Target	Drug	PDID	Drug Status	Primary Indications (1. Mode of Action; 2. Treatment)	2D Structure
I7L	Eluxadoline	PD008978	Approved	1. An agonist of mixed mu-opioid receptor.
				2. Used for the treatment of irritable bowel syndrome with diarrhea.
I7L, VETFS	Burixafor	PD058763	Investigational	1.Inhibitor of CXC chemokine receptor 4 (CXCR4)
				2. Used in trials studying treatment of Hodgkin's disease, Non-hodgkin's Lymphoma and Multiple Myeloma.
I7L	Cefotiam	PD010157	Approved	1.Inhibits the bacterial cell wall biosynthesis.
				2. Used for the treating a myriad of bacterial infections.
I7L	Tobramycin	PD001728	Approved	1.Inhibits the synthesis of protein by binding to ribosome 30S subunit in bacterial cells.
				2. Used for the treating a myriad of bacterial infections.
I7L	Danegaptide	PD058923	Investigational	1.A selective 2nd generation gap junction modifier in adjacent cardiomyocytes.
				2. Used in studying treatment of chronic atrial fibrillation (AF) and postoperative AF in large animal models.
Top1B	Carbaphosphonate	PD059915	Investigational	1.3-dehydroquinate synthase inhibitor (Staphylococcus aureus)
Top1B	α-d-galacturonic acid	PD006851	Investigational	1. Backbone of pectin, cellular binders in the peel of many different fruits and vegetables.
				2. Used as: antidiarrheal drug, food emulsifier, food stabilizer, food thickening agent and food gelling agent.
Top1B	Adenosine phosphonoacetic acid	PD059913	Investigational	1. ndp kinase (human) inhibitor.
Top1B	Guanosine-2′,3′-O-methylidenephosphonate	PD006419	Investigational	1. catalyze the phosphorolytic breakdown of the N-glycosidic bond in the beta-(deoxy) ribonucleoside molecules
Top1B	Deferitazole	PD058219	Investigational	1. An iron cheater
				2. Used in trials studying the treatment and basic science of Beta-thalassemia.
VETFS	DB02785	PD059991	Investigational	1. Gag-Pol polyprotein inhibitor
				2. Used as an antiviral (antiHIV-1)
VETFS	Radezolid	PD012727	Investigational	1. A 2nd generation's oxazolidinone antibiotic, inhibits the RRBP1.
				2. Used in trials for studying the treatment of Abscess, Infectious Skin Diseases.
VETFS	Rapastinel	PD070040	Investigational	1. NMDA receptor modulator with glycine-site partial agonist properties.
				2. Used in trials studying the treatment of Major and Obsessive-Compulsive Disorder (OCD).
I7L, A20R	Batefenterol	PD058591	Investigational	1. A bifunctional muscarinic (M2 and M3 receptors) antagonist and β2-agonist
				2. Used in trials for studying the treatment of Chronic Obstructive Pulmonary Disorder (COPD).
A20R	Nebivolol	PD009358	Approved	1. β-1 adrenergic receptor antagonist (beta blockers).
				2. Used to treat high blood pressure either alone or in combination with other medications.
A20R	Pimozide	PD001879	Approved	1. Dopamine type 2 receptors blocker
				2. Used to control motor or verbal tics caused by Tourette's disorder.
A20R	BAZ2-ICR	PD046767	Investigational	1. Epithelial sodium channel blocker
				2. Used in trials for studying cystic fibrosis and chronic bronchitis.
A20R	PF-610355	PD058896	Investigational	1. A novel Ultra-Long-Acting β2-Adrenoreceptor agonist.
				2. Used for studying the treatment of Asthma and Chronic Obstructive Pulmonary Disease.

Open in a new tab

3.8. RMSD, RMSF, Rg and SASA calculations suggests stability of the ligand-bound complexes

Two independent sets of MD simulations were undertaken in this study by using the GROMACS package. The first set involved the model refinement of four proteins A20R, I7L, Top1B and VETFS to reduce stearic clashes among their residues. This was achieved by simulating every protein for 400 ns. The trajectories were subjected to clustering using the TTClust. The centroid of the cluster with the least mean rmsd between the frames were selected for protein-ligand complex MD simulation. The summary of cluster analysis of the trajectories after 400 ns of simulation has been presented in Supplementary Fig. S1. Thereafter, to access the conformational stabilities of the 12 protein-ligand complexes during the 150 ns of MD production, various quality control parameters like RMSD, Rg, SASA, RMSF, H-Bonds, secondary structure, essential dynamics, MM/GBSA binding free energies and distance of ligand from active site were examined.

RMSD characterises a protein-ligand complex's conformational stability in its dynamic state during simulation. A low difference in RMSD indicates low and consistent fluctuation between the ensembles which suggests a protein-ligand complex is stable. To ensure the sampling method's reliability, the RMSDs of the proteins' alpha carbon atoms (C-α) were analysed by plotting them against the time scale of 150 ns from the starting structure (Fig. 4 ). None of the A20R ligand associated complexes could achieve equilibrium after 150 ns of simulation (Fig. 4A). The RMSD plot of the I7L-BAT complex (Fig. 4B) depicted an increasing trend for the first 45 ns. However, the ensembles that followed gradually traced a constant trajectory for the next 105 ns i.e., from 45 ns to 150 ns with consistent minor fluctuations within a permissible window of 3 Å to 4.5 Å and an average RMSD of 0.63 ± 0.26 Å. The trajectory of the I7L-BUR complex (Fig. 4B) followed a stable trend for the first 55 ns while fluctuating between a small window of 2 Å to 3 Å. The complex then encountered a minor fluctuation before finally settling into a relatively flat trajectory with an average RMSD of 0.9 ± 0.3 Å for the next 95 ns i.e., from 55 ns to 150 ns within a permissible window of 2.5 Å to 4.5 Å. Similarly, the trajectory of I7L-ELU complex (Fig. 4B) projected an increasing trend for the first 35 ns. Thereafter, the ensembles adopted a relatively flat trajectory with an average RMSD of 0.72 ± 0.28 Å for the last 115 ns i.e., from 35 ns to 150 ns within a small window of 3 Å to 5 Å. These results indicate the protein I7L has successfully achieved conformational stability with the three ligands. The ensembles of Top1B-BAT complex (Fig. 4C) traced out a relatively flat trajectory with consistent minor fluctuations for the first 85 ns within a confined permissible space between 2 Å to 4 Å. Minor blips within very narrow time frames of (25–30) ns and (45–48) ns were also encountered. This was followed by high fluctuations for the next 30 ns i.e., from (85–115) ns following which the complex again reverted back to its original 2 Å to 4 Å window for the last 35 ns. All the stable ensembles of Top1B-BAT complex that fluctuated within the 2 Å to 4 Å permissible window had an average RMSD of 1.16 ± 0.85 Å. Strangely, this pattern repeated itself again in Top1B-BUR and Top1B-ELU complexes respectively (Fig. 4C). The Top1B-BUR complex retained a consistent trajectory for the first 40 ns within an allowable window of 2 Å to 4 Å. This complex then navigated through a similar violent RMSD fluctuations for the next 55 ns i.e., from (40–95) ns before finally returning back to its initial 2 Å to 4 Å window for the last 55 ns with consistent minor fluctuations. All the stable ensembles of Top1B-BUR complex that fluctuated within the 2 Å to 4 Å permissible window had an average RMSD of 0.88 ± 0.59 Å. The ensembles of Top1B-ELU complex followed an increasing trend for the first 20 ns. Thereafter, the ensembles retained a stable trajectory throughout the entire 150 ns between a 2 Å to 4 Å window with an average RMSD of 1.17 ± 0.77 Å with the exception of two small time frames between 75 ns to 85 ns and 115 ns–120 ns respectively where major fluctuations in RMSD were clearly distinguishable (Fig. 4C). The ensembles of VETFS-BAT complex pursued an upward trajectory for the first 10 ns (Fig. 4D). Thereafter, the ensembles retained a stable trajectory for the next 140 ns i.e., from 10 ns to 150 ns with constant but minor fluctuations within an allowable window between 3 Å to 5 Å with an average RMSD of 0.57 ± 0.36 Å. The ensembles of VETFS-BUR complex traced an ascending step ladder trajectory for the first 70 ns before gradually settling into a stable trajectory for the remaining 80 ns i.e., from 70 ns to 150 ns (Fig. 4D). The stable ensembles during this time frame constantly fluctuated within a narrow and acceptable window between 4 Å to 6 Å with an average RMSD of 0.88 ± 0.39 Å. A rising trend was observed in the trajectory of VETFS-ELU complex for the first 40 ns (Fig. 4D). The trajectory thereafter declined which was followed by a steady rise in the next 45 ns before finally achieving stability for the last 55 ns i.e., from 85 ns to 150 ns. The stable ensembles during this time frame fluctuated within a 6 Å to 8 Å window with an average RMSD of 0.88 ± 0.44 Å. It is also interesting to note that the VETFS-BAT system achieved equilibrium within a very short duration in comparison to its BUR and ELU bound counterparts. This suggests, in comparison to VETFS-BAT complex, VETFS-BUR and VETFS-ELU complexes had to navigate through several conformational transitions to accommodate their respective ligands before finally achieving stable conformations. Overall, these results indicate that the enzyme VETFS has successfully acquired conformational stability with all the three ligands. Taken together, only nine complexes that showed conformational stability were considered for further examination. The density distribution of RMSDs of these nine stable complexes have also been provided in Supplementary Fig. S2.

Subsequently, the radius of gyration (Rg) of the 9 stable protein-ligand complexes were also assessed (Fig. 5 ). Rg examines a protein's compactness where a constant Rg suggests that the system has achieved a relatively constant shape and size during the entire simulation. The average radius of gyration of I7L-BAT, I7L-BUR and I7L-ELU complexes were 23.16 ± 0.19 Å, 22.74 ± 0.14 Å and 23.16 ± 0.24 Å respectively. These complexes successfully retained their compactness without undergoing any major structural re-arrangements during the entire 150 ns simulation as demonstrated by the consistency of their time dependent Rg plot (Fig. 5A). Similarly, the Rg of the Top1B-BUR and Top1B-ELU complexes also remained fairly consistent during the entire simulation with average Rg's of 25.51 ± 0.51 Å and 25.34 ± 0.52 Å respectively (Fig. 5B). The average Rg of Top1B-BAT complex was 24.71 ± 0.69 Å during the entire simulation with some minor fluctuations. The Rg of Top1B-BAT complex remained consistent for the first 75 ns following which it experienced a sharp decline and a subsequent upturn in the next 55 ns i.e., between 75 ns and 130 ns before retracing its original trajectory in the last 20 ns (Fig. 5B). Overall, the consistent Rg of the three Top1B-ligand bound complexes suggest they have retained their structural integrity and are relatively compact after 150 ns of simulation. The VETFS-BAT complex also remained compact with an average Rg of 30.05 ± 0.27 Å (Fig. 5C). This was evident from its time-evolution Rg plot which remained fairly consistent during the entire simulation. The Rg of VETFS-BUR and VETFS-ELU complexes remained steady for the first 70 ns (Fig. 5C). Thereafter, the Rg of VETFS-BUR declined linearly for the next 30 ns before retaining its steadiness for the last 50 ns. On the contrary, the Rg of VETFS-ELU complex witnessed a linear increment for the next 10 ns before becoming relatively constant for the last 70 ns. This suggests that the presence of BUR in the active site of VETFS has enhanced its compactness. It also suggests that BUR is probably better adapted to VETFS in comparison to ELU. The Rg of VETFS-BUR and VETFS-ELU complex was 29.59 ± 0.32 Å and 30.64 ± 0.39 Å respectively. Taken together trends of Rg of all the nine complexes are in congruence with their corresponding RMSD trends that reflect the relative stabilities and compactness of the 9 ligand bound enzyme complexes.

SASA is another frequently studied quality control parameter that is examined to scrutinize a protein-ligand complex's conformational stability. The solvent environment around the protein plays a crucial role in retaining a protein's folding and governs the protein–ligand interaction processes, orientation, and stability. The time evolution SASA plots clearly indicate that all the nine complexes are relatively well equilibrated during the entire 150 ns simulation (Fig. 5D, E and 5F). The SASA values varied within a narrow range between 233.4 ± 4.43 nm² - 242.11 ± 6.67 nm² for I7L-ligand complexes, 191.76 ± 3.32 nm² - 192.50 ± 3.04 nm² for Top1B-ligand complexes and 339.96 ± 4.35 nm²–348.41 ± 5.27 nm² for VETFS-ligand bound complexes respectively. These evidences are in agreement with their corresponding RMSDs and suggest that all the proteins have successfully attained conformational stabilities with their respective ligands without undergoing any major alterations in their available surface area.

Root Mean Square Fluctuation (RMSF) was computed to examine the positional fluctuation of each amino acid residue around its mean position. The Cα atoms of all the I7L complexes shared very similar RMSF trends (Fig. 6 A). This case was also similar for the Cα atoms of Top1B (Fig. 6B) and VETFS (Fig. 6C) ligand bound complexes. The Cα atoms associated with the loops experienced high fluctuations whereas the Cα atoms associated with α-helix and β-sheets were comparatively rigid in all the complexes. High fluctuations were were also associated with the C terminal residues in all the nine complexes. Amongst the I7L complexes, the highest average atomic fluctuation of 1.83 ± 0.97 nm was recorded in the residues of the ELU bound complex. Similarly, amongst the Top1B and VETFS complexes, highest average atomic fluctuations of 3.19 ± 1.6 nm and 2.71 ± 1.27 nm was recorded in the residues of BUR and ELU bound complexes respectively. From Fig. 6B it is quite evident that high atomic fluctuations of all the Top1B complexes are mostly restricted to the residues from 1 to 70 which encompass the protein's N-terminal domain. The fluctuations in the residues that encompass the ligand bound C-terminal DNA binding domain (76–314) are relative restricted which hints at the conformational rigidity of this domain. This suggests that the C-terminal DNA binding domain of Top1B has been successfully stabilized by the three ligands and the wild fluctuations in Top1B's RMSD trajectories are most probably related to the flexibility of its N-terminal domain.

The molecular dialogues between a ligand and the flexible active site amino acid residues are largely driven by H-bond interactions. We further examined the time evolution plots of H-bonds involved in the molecular interaction of the three ligands BAT, BUR and ELU with each of the three proteins namely, I7L, Top1B and VETFS. The nine protein-ligand complexes depicted differential intermolecular H-bonding patterns during the 150 ns of MD run (Fig. 6D, E and 6F). The equilibrated complexes of protein I7L formed an average of 2.31 ± 1.36, 2.1 ± 1.19 and 1.5 ± 1.03 H-bonds per frame with the ligands BAT, BUR and ELU respectively. The equilibrated complexes of Top1B formed an average of 1.28 ± 1.03, 1.79 ± 1.17 and 1.55 ± 1.06 H-bonds per frame with the ligands BAT, BUR and ELU respectively. Similarly, the fully equilibrated complexes of VETFS formed an average of 1.14 ± 0.86, 1.705 ± 1.22 and 1.59 ± 0.65 H-bonds per frame with the ligands BAT, BUR and ELU respectively. H-bond interactions govern the correct orientation and stability of a ligand in a protein's active site which subsequently drives the protein's overall compactness and structural conformity. Therefore, the compactness and structural rigidity of I7L, Top1B and VETFS as inferred from their stable time evolution RMSD, Rg and SASA plots are most probably driven by the these H-bond interactions.

3.9. The stabilities of the protein-ligand complexes are largely driven by switching between secondary structural elements

Changes in secondary structure as a function of time were also mapped for every complex using the Biolite, a Python-based molecular biology suite. Helices and β-sheets are stable amongst the secondary structure elements whereas coils, loops and turns are highly flexible. Constant switchings between coils, bends and turns were clearly visible in all the complexes during their simulations. These switchings probably contribute in parts to the molecular motions of all the complexes. The integrity of all the major α-helices and β-sheets remained intact in all the complexes. These α-helices and β-sheets are most likely responsible for the global conformational rigidity of the proteins. Analysis of the time evolution map of the secondary structure elements resulted in some interesting observations. The flexible turn between the residues 40 and 50 in the I7L-ELU complex gradually transformed into a more rigid α-helix after 10 ns and remained constant for the remaining duration of the simulation (Supplementary Fig. S3). This correlates well with its RMSD plot where the stability of the complex was achieved immediately after 10 ns. This observation suggests that ELU mediated stabilization of the enzyme I7L was most likely driven by the metamorphosis of a turn between the residues 40 and 50 into an α-helix. The turn at the residues between 70 and 75 of Top1B-BUR complex was replaced by a more flexible coil between 50 ns and 95 ns which probably explains the corresponding sharp RMSD and RMSF fluctuations (Supplementary Fig. S4). The residues between 420 and 430 in VETFS-BUR complex continued to switch between a turn and α-helix for the first 85 ns before smoothly transitioning into an α-helical structure for the remaineder of the simulation (Supplementary Fig. S5). This correlates well with the RMSD trajectory of this complex and suggests that ELU driven stabilization of VETFS is regulated by the formation of an α-helix. Overall, these results suggest that the transition of unstable protein-ligand complexes into their equilibrated states and vice versa are largely driven by switching between secondary structural elements.

3.10. The equilibrated ensembles of the protein-ligand complexes are grouped into many clusters

The ensembles of I7L, Top1B and VETFS ligand bound complexes generated after 150 ns of the simulation were clustered with the TTClust to identify the clusters with stable ensembles and the corresponding centroid structures that represent these clusters. The TTClust analyses the trajectories to generate graphical imagery like distribution within the clusters, the relative distance between clusters and hierarchical cluster dendrograms. Fig. 7 and Supplementary Fig. S6 depict the results of the cluster analysis of I7L complexes during the 150 ns of MD simulation. The equilibrated ensembles of I7L-BAT complex was subgrouped into 3 clusters (C2, C3, C4), whereas those of I7L-BUR and I7L-ELU complex were further subgrouped into 2 clusters (C2, C3) each. The equilibrated ensembles of Top1B-BAT, Top1B-BUR and Top1B-ELU complexes were also subgrouped into 2 clusters (C1, C2) each (Fig. 8 and Supplementary Fig. S6). Similarly, the equilibrated ensembles of VETFS-BAT, VETFS-BUR and VETFS-ELU complexes were subgrouped into 3 (C1, C2, C3), 2 (C3, C4) and 1 (C3) clusters respectively (Fig. 9 and Supplementary Fig. S6). The distribution of these equilibrated ensembles into different clusters suggests that they probably sample small conformational subspaces with minor structural differences. The protein-ligand interactions experienced by the centroid frames of these stable sub-clusters are summarised in Table 5 .

3.11. PCA suggests the motions of the protein ligand complexes are relatively constricted

Protein function is controlled by shifting between various conformations. This ability to flexibly switch between multiple conformations is regulated by a set of their collective motions. These motions are critical to different biological processes and have a key contribution in the transmission of biological signals. A fair amount of flexibility as well as rigidity is necessary, especially of the binding site residues for any protein to be functional. A tighter ligand interaction would essentially restrict a protein's motion, thereby preventing it from sampling the critical conformations necessary for its functionality. Therefore, essential dynamics (ED) analysis was attempted to better understand the collective motions of all the nine ligand bound protein complexes in the 2D conformational phase space during the 150 ns simulation. As a multi-variate statistical method, ED is reliant upon diagonal covariance matrix of a protein's Cα atom to trace its global motion through Eigen vectors (EVs) popularly called principal components (PCs) and eigenvalues [83]. The EVs describe the global direction of motion of the atoms and the associated eigenvalues depict the atomic contribution of motion in MD trajectories of a ligand bound system which is essentially controlled by a proteins' secondary structure. The collective dynamic motions of the proteins captured through the projections of EVs (PC1 and PC2) are shown in Fig. 10 . The PCA statistics and the percentage of variance covered by the first three eigenvectors of I7L, Top1B and VETFS complexes are summarised in Supplementary Table S9. As visible from the figure, the first two PCs accounted for the bulk of all the internal motions. This indicates that the vectors mapped within this conformational subspace elucidate the essential subspace of the system. The trace values calculated from the covariance matrices of I7L-BAT, I7L-BUR, I7L-ELU, Top1B-BAT, Top1B-BUR, Top1B -ELU, VETFS-BAT, VETFS-BUR and VETFS-ELU complexes were 17.09 nm², 16.99 nm², 18.22 nm², 30.52 nm², 40.01 nm², 23.69 nm², 44.98 nm², 42.12 nm² and 57.12 nm² respectively. The trace values of Top1B-BUR, VETFS-BAT, VETFS-BUR and VETFS-ELU complexes were relatively higher compared to the other five complexes. These results are consistent with higher occupation of the conformational subspace by these complexes which is clearly visible in the form of dispersed spectral dots and further hints at greater flexibility of these proteins. Therefore, as a consequence of higher fluctuations these complexes had to sample a much wider conformational space most probably to accommodate the ligands before accomplishing the ensembles in their dynamically equilibrated states. It is interesting to note that the molecular motions of the 3 Top1B complexes are relatively constricted (as supported by their low phase space occupancy). This further strengthen the idea that ensembles of all the Top1B complexes have achieved as state of dynamic equilibrium with their respective ligands. Overall, ED analysis suggests the collective molecular motions of all the 9 complexes are restricted to a small and localised conformational space. These results are consistent with their corresponding RMSD, SASA and Rg plots which further strengthen the notion that the three proteins namely, I7L, Top1B and VETFS have been stabilized by the ligands BAT, BUR and ELU.

3.12. Porcupine plot assessment suggest stability of Top1B's ligand bound DNA binding domain

Porcupine plots were developed to gain quantitative insights into the molecular motion of the 9 protein-ligand complexes. To picture the molecular motions, a cone was drawn along the C-alpha atom of each residue in the direction traced by its extreme projection on PC1 (Fig. 11 ) where the length of the cone signifies the amplitude of motion of the Cα atom. It shows that a vast majority of all the fluctuations are closely associated with the regions representing the loops/coils in all the protein ligand complexes. These results are completely in tune with their corresponding RMSF plots. Similarly, cones were also observed on the residues neighbouring the ligands in all the complexes. These fluctuations have most likely occurred to accommodate the ligands in the active site cavity. Interestingly, large protruding cones were clearly visible in the N-terminal domains of all the Top1B complexes which was not visible in their ligand bound C-terminal DNA binding domain. This clearly suggests the large fluctuations in Top1B complexes are particularly associated with the N-terminal domain. This is perfectly in sync with their RMSF plots where large fluctuations were specifically confined to the residues of the N-terminal domain. Combined together, the rmsf plots and porcupine plots of the Top1B complexes evidently explain the large fluctuations in RMSDs of all the Top1B complexes. These correlations beyond doubt clarify that the target C-terminal DNA binding domain of Top1B has successfully established an equilibrium with the ligands BAT, BUR and ELU.

3.13. MM/PBSA assessment suggests stable binding of the ligands BAT, BUR and ELU with the three proteins I7L, Top1B and VETFS

The assessment of RMSD trajectories suggest all the nine protein-ligand systems stabilized in the last 50 ns. Therefore, MM/PBSA was used to estimate the binding free energy and understand the associated molecular interactions between the proteins I7L, Top1B, VETFS and the ligands BAT, BUR and ELU for the last 50 ns of MD run. Supplementary Fig. S7 suggests the trajectories of binding free energy of all the nine protein-ligand complexes remained almost constant in the last 50 ns. Fig. 12 suggests that BAT, BUR and ELU bind favourably with all the three proteins. The three proteins, I7L, Top1B, VETFS, had highest binding affinity (ΔG_mmpbsa) of −29.63 ± 7.53 kcal/mol, -32.72 ± 5.17 kcal/mol and -24.96 ± 4.48 kcal/mol respectively with BAT suggesting BAT might be the most suitable ligand for their inhibition. The decomposition of the average binding free energy of all the nine complexes into energy terms suggests the binding of the ligands BAT, BUR and ELU in the active site of majority of the proteins is stabilized by Van der Waal's (vdW) energies (Fig. 12). However, another critical energy term average electrostatic (ΔE_ele) seems to be the major contributor towards the stabilization of Top1B-ELU (ΔE_ele = −59.94 ± 11.95 kcal/mol) and VETFS-ELU (ΔE_ele = −42.85 ± 10.55 kcal/mol) complex respectively.

We further employed per residue binding free energy decomposition to determine the molecular interactions of every ligand with the binding pocket residues of the proteins BAT, BUR and ELU (Fig. 13 ). Per residue decomposition analysis suggests the active site residues energetically favour the stable binding of BAT with protein I7L. The free energy decomposition analysis plot suggests the active site residues Arg99, Ile298, Ile316, Leu109, Phe256, Pro110 and Pro264 favour the stable binding of BAT with I7L. The residues Met67, Pro64, Pro110, Pro114, Ser112, Thr111, Tyr46, Tyr238, Val47 favour the binding of BUR to the protein I7L. Similarly, the active site residues Asp258, Leu239, Pro264, Ser112 and Tyr238 favour the binding of ELU with I7L. The residues Pro 110, Ser112 and Pro264 of I7L protein were common to two of the three ligands. This suggests these residues might be important to I7L. The residues Met126 and Phe131 favoured the binding of BAT to the active site of Top1B. The residues Arg223, Gly264, Ile129, Phe127, Phe128, Phe174, Tyr233, Val 262, and Val 263 favoured the stable binding of BUR in the active site of Top1B whereas the pocket residues Arg130, Ile129, Phe131, Val 262 and Val 263 contributed towards the favourable binding of ELU in Top1B's active site. The residues Ile129, Phe 131, Val262 and Val263 were common to two of the three ligands bound to Top1B. Therefore, these residues might be important to Top1B. The residues Lys273, Met 275, Tyr276 and Tyr477 favoured the binding of BAT in the active site of protein VETFS whereas the active site residues Arg262, Met275, Phe271, Tyr276 guided the favourable binding of ELU with VETFS. Residues favouring the binding of BUR with VETFS could not be identified. The residues Met275 and Tyr276 were common to both the ligands. As such, these residues might be critical to VETFS.

All the ligand bound complexes maintained a favourable distance below 10 Å between the centre of mass of the interacting active site residues and centre of mass of their associated ligands with the exception of VETFS-BUR complex in the last 50 ns of simulation (Supplementary Fig. S8). This suggests stable binding modes of these ligands with the active site residues of the three proteins. The distance between the centre of mass of the ligand BUR and the centre of mass of the active site pocket residues of VETFS was far greater than 10 Å. This suggests BUR has been displaced from the active site of VETFS and therefore is no longer bound to it. Therefore, taken together only eight complexes namely, I7L-BAT, I7L-BUR, I7L-ELU, Top1B-BAT, Top1B-BUR, Top1B-ELU, VETFS-BAT and VETFS-ELU have successfully attained a state of dynamic equilibrium. These findings suggest the ligands BAT and ELU have the potential to function as inhibitors of all the three enzymes whereas BUR can potentially inhibit the enzymes I7L and Top1B.