Network metrics, structural dynamics and density functional theory calculations identified a novel Ursodeoxycholic Acid derivative against therapeutic target Parkin for Parkinson's disease.

Parkinson's disease (PD) has been designated as one of the priority neurodegenerative disorders worldwide. Although diagnostic biomarkers have been identified, early onset detection and targeted therapy are still limited. An integrated systems and structural biology approach were adopted to identify therapeutic targets for PD. From a set of 49 PD associated genes, a densely connected interactome was constructed. Based on centrality indices, degree of interaction and functional enrichments, LRRK2, PARK2, PARK7, PINK1 and SNCA were identified as the hub-genes. PARK2 (Parkin) was finalized as a potent theranostic candidate marker due to its strong association (score > 0.99) with α-synuclein (SNCA), which directly regulates PD progression. Besides, modeling and validation of Parkin structure, an extensive virtual-screening revealed small (commercially available) inhibitors against Parkin. Molecule-258 (ZINC5022267) was selected as a potent candidate based on pharmacokinetic profiles, Density Functional Theory (DFT) energy calculations (ΔE = 6.93 eV) and high binding affinity (Binding energy = -6.57 ± 0.1 kcal/mol; Inhibition constant = 15.35 µM) against Parkin. Molecular dynamics simulation of protein-inhibitor complexes further strengthened the therapeutic propositions with stable trajectories (low structural fluctuations), hydrogen bonding patterns and interactive energies (>0kJ/mol). Our study encourages experimental validations of the novel drug candidate to prevent the auto-inhibition of Parkin mediated ubiquitination in PD.

Introduction

Parkinson's disease (PD) ranks second amongst neurodegenerative disorders after Alzheimer's, affecting millions with a higher prevalence in males of mean age ∼ 55 years [1], [2]. Clinically, PD is characterized by several motor (bradykinesia, tremor, muscle rigidity) and non-motor (weight loss, anxiety, depression, sleep dysfunction) disorders [3]. The pathological hallmarks of PD are the loss of dopaminergic neurons in the Substantia Nigra pars compacta (SNpc) and aggregation of α-synuclein (SNCA) that develops into insoluble fibrillar aggregates called Lewy Bodies (LB) [4]. Reduction in dopamine levels results in dopaminergic cell loss and intraneuronal inclusions leading to motor symptoms [5], [6]. Early diagnosis and prophylaxis of PD remains a primary challenge for clinicians as the long latency between first dopaminergic cell damage to the onset of clinical symptoms (∼70–80 % damaged dopaminergic neurons) poses great challenges in its therapeutic interventions. Hence, it is crucial to find molecular biomarkers reliable enough to distinguish PD from other disease conditions that can be exploited for theranostic purposes [7]. The conventional treatment regimen involves the usage of dopaminergic drugs (Levodopa) that replace the action of dopamine in the body. However, these drugs stabilize the motor symptoms and not their underlying causes [8]. Available data sources and compound database screening highlighted several repurposed drugs as alternate therapeutics against PD. However, several side-effects due to neurodegeneration in brain cells upon nilotinib administration, unclear action mechanism of simvastatin and risks of hepatocellular carcinoma from Ursodeoxycholic Acid (UDCA) over-dosage widened the gap in designing effective treatment strategies [9].

Modern-day pharmaceutics has tremendously advanced with the advent of artificial intelligence (AI) through machine learning (ML) and deep learning (DL) algorithms [10], [11], [12]. DL and ML have been a boon to physicians in designing high precision treatment strategies in surgery, neurodegenerative disorders and can also predict upcoming diseases thereby contributing to increased overall survival of patients [13], [14], [15]. Our research team has been instrumental in predicting alternate drug-targets/biomarkers [16], [17], [18], [19], [20] and novel inhibitor molecules [21], [22], [23], [24] against nosocomial pathogens through genomics, systems biology and structural biology approaches. Direct influencers of SNCA-mediated insoluble fibrillar aggregate formation can be efficient drug-targets for PD. The extent of involvement of key genes/proteins in insoluble fibrillar aggregate formation (α-synuclein aggregation) and the structural suitability of the most suitable therapeutic target was evaluated. The specific goals were to establish the association of the most suitable therapeutic target in the modulation of α-synuclein aggregation through gene interaction network (GIN), functional-association studies (enrichment analyses) and network-centrality values (hubness metrics). Subsequently, the conformational stability and druggability of the potential target were assessed through molecular dynamics and energy parameters. An amalgamation of GIN, virtual-screening, structural dynamics and density functional theory (DFT) calculations was found to be efficient while addressing the hypothesis. Evaluation of therapeutic targets and candidates in PD is limited due to experimental tediousness using animal models, time, cost and ethical clearance norms. The proposed integrated computational approach is hence advantageous in all aforesaid aspects while proposing experimentally translatable leads against similar human disorders.

Materials and methods

Data curation, GIN construction and functional enrichment analysis

GIN analysis was performed to identify the key functional influencers (genes/proteins) of SNCA-mediated insoluble fibrillar aggregate formation which is the hallmark of PD. The network associated annotation data was also used for determining suitable drug-target in the interactome. PD specific human genes identified from several literary sources were subjected for GIN analysis from the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING v11.5) database, which is a comprehensive data repository for both experimentally validated and computationally predicted interactions. STRING retrieves interaction data from several low and high throughput experiments, co-expression analyses, curated pathway databases, computational predictions (gene fusion, phylogenetic co-occurrence, chromosomal neighborhood) and text mining [25]. The interactions from these seven channels are scored that ranges from zero to one based on the probability of the association being true. STRING integrates these scores into a final confidence score (CS) imparting a meaningful biological association. The CS is further categorized (based on the level of experimentation) into four classes of highest confidence (1.0 ≥ CS ≥ 0.9), high confidence (0.9 > CS ≥ 0.7), medium confidence (0.7 > CS ≥ 0.4) and low confidence (0.4 > CS ≥ 0) scores [26], [27], [28]. STRING designates each interacting entities (genes/proteins/metabolites) as 'node(s)' while the interactions are coined as 'edge(s)'. In the present study, interactions from experimentally evident (highest CS) to computationally predicted (low CS) were considered with statistical significance of P-value < 0.05.

To study the functional association of the interacting genes in the network, STRING possesses an inbuilt plug-in for performing functional enrichment analysis (FEA). To decipher the role(s) of each interacting entities imparting specific biochemical and physiological functions in the respective pathway(s), FEA is further segregated into gene ontology (GO) terms, domain descriptors (InterPro, Pfam, SMART and UniProt keywords) and pathway enrichments (KEGG and Reactome). GO terms are further classified into three independent ontologies namely Biological Process (BP), Molecular Functions (MF) and Cellular Components (CC). BP highlights the physico-chemical transformation of gene or gene products (expressed protein) while the intervening biochemical activities are defined by MF and CC refers to the sub-cellular compartments where the gene/protein is active.

Clustering and topological analysis

The interactome retrieved from STRING was further visualized in an user-friendly open-source visualization software Cytoscape v3.9.1 [29] which is used extensively for data integration, biological network modeling and construction. Cytoscape offers diverse inbuilt plug-ins that aids in better understanding of the network characteristics and analysis [30], [31].

Clustering analysis was performed using molecular complex detection (MCODE) algorithm which detects highly efficient and densely interconnected nodes in the GIN. The clustering is performed based on the three stages of vertex weighting, complex prediction and optionally post-processing. The MCODE algorithm calculates local neighbourhood density for each node that traces the densely interconnected gene/protein (seed protein) from the network. The tool thereby extrapolates from the seed in all direction to encapsulate the dense region and finally represent it graphically in the form of clusters. Genes/proteins presumably function as a group (operon or regulon) rather than individually. Therefore, each nodes in a cluster is likely to impart similar functionality and plays active roles in disease progression and virulence. Each cluster is designated with a score based on the density, size and connectivity of the nodes in the respective cluster [32].

A wide array of topological network metrics are computed using NetworkAnalyzer an inbuilt plug-in of Cytoscape which computes various simple to complex topological parameters aiding in better understanding of a biological network. The topological parameters are calculated for both directed and undirected networks which include the number of edges, nodes, degree of distribution, betweenness and closeness centralities, average shortest path length, density, radius and diameter of the network [33]. Identifying central/hub-genes in an interactome can be best established by betweenness centrality as this topological parameter measures the influence of each node on other nodes in the interactome. The parameter also computes the flow of information from one node to another, detecting the node that acts as a bridge between other nodes. This clearly highlights the importance of the node as well as loss of functionality in its absence. The parameter hence identifies hubs that connect functionally significant genes in the interactome. Clustering coefficient signifies the average degree of node neighborhoods. The nodes thus identified tend to form tightly interconnected clusters characterized by high degree of edge counts [34], [35], [36], [37], [38]. Hence, these topological parameters are prioritized for identification of the hub-genes playing crucial roles in influencing the entire interactome.

Protein structure modeling, secondary and tertiary structure validations

The 3D coordinates of potential therapeutic target were either retrieved from RCSB-PDB or else modeled by extensive dual-step modeling approach integrating both knowledge based (homology, threading and iterative) and molecular dynamics algorithms as previously described [19]. In the present study, uncharacterized protein stretches in template were cured using automated template based modeling platform SWISS-MODEL and python script based standalone software MODELLER 10.0 [39], [40]. MODELLER executes comparative protein structure modeling upon satisfying the spatial restraints (stereochemical, homology-derived, manually curated and statistical preferences) using combination of conjugate gradients and molecular dynamics with simulated annealing algorithm. This model building algorithm of MODELLER generates 3D conformers analogous to NMR spectrometry derived structures [40].

The Ramachandran outliers, poor rotamers, steric clashes and torsions of the target protein were refined by furnishing and repacking the side chains through GalaxyRefine server [41]. The overall structure relaxation was performed by short molecular dynamics simulation for 0.6 ps under CHARM22 force-field with Coulomb potential energy, FACTS solvation free energy and Lennard-Jones interaction energy. Later, the overall energy of the protein target was minimized in-vacou by carrying out 2000 steps each for steepest-descents and conjugate-gradients algorithms with GROMOS96-43B1 force-field mechanics using standalone Swiss-PDB Viewer v4.1.0 software [42].

Secondary structural assessments were performed using PSIPRED protein structure prediction server [43] while the tertiary structure was characterized using ProSA-web [44], HARMONY [45] and ProTSAV [46] servers. Protein Structure Analysis (ProSA)-web is a dynamic server that predicts erroneous 3D protein structures by assessing the protein backbone (Cα-atoms) potentials in terms of Z-score and energy plots. Z-score signify the overall model quality upon calculating the deviations amongst the most stable conformer and random conformers of the protein. The energy plot reveals local model quality by considering each amino acid residues with energy functions, where the positive values signify erroneous conformations and negative values indicate proper folding with minimal errors [44]. For further assessment of protein structure quality, HARMONY was deployed which scrutinizes the sequence-structure compatibility based on their solvent accessibility, backbone, and hydrogen bonds at the residue-level and predicts misfolds through propensity-calibration curve while local errors are assessed through substitution curve [45]. Protein Tertiary Structure Analysis and Validation (ProTSAV) is a robust server that integrates ten online tools (Procheck, ProSA-web, ERRAT, Verify3D, dDFire, Naccess, MolProbity, D2N, ProQ, PSN-QA) to assess various Ramachandran plot parameters, side chain conformations, planar peptide bonds, polarity of residues and several atom–atom outliers. ProTSAV combines the results from all these servers and provides a statistically significant unified quality assessment Z-score ranging from 0 to 1. Structures with lower Z-scores are suitable for post analysis while higher Z-scores signifies misfolded/erroneous structures; therefore misfits for further analyses [46].

Delving deeper, the functional domains retrieved from FEA were further validated using the InterPro [47], [48] and Pfam [49], [50] databases. The native disorders in protein domains and their functional dynamicity were predicted through the DISOPRED3 server, which screens the PSI-BLAST profile to predict the probability of disorderness of each residue [51]. On tertiary level, the global folding free energy change (ΔG) of protein domains were assessed and estimated using a structure based approach to predict stability weakness and strengths through the SWOTein server. The structural descriptors were expressed in terms of inter-residue spatial distance (ΔGdist), inter-residue interaction with solvent (ΔGacc) and backbone torsion angles (ΔGtor). The negative values indicate stability strengths and positive values conferred stability weakness of the protein domains/residues. The distance descriptor indicates the free energy change from the tertiary interactions established by a residue with its neighbouring residues in spatial conformations while the accessibility and torsion components indicates the energy required to interact with the solvent molecules and the energy required to maintain the backbone torsion angles to lie in the favoured regions respectively [52].

Residue level propensity and stability analysis through coarse dynamics and MDS

Residue-level backbone dynamics in the form of N—H S bond-order restraints and coarse-grained dynamics refinements were performed through the DynaMine [53] and CABSflex2.0 [54] servers respectively as described previously [19], [20], [55]. The DynaMine web-server predicts the atomic bond-order restraints according to the molecular reference frame derived from experimentally derived NMR chemical shifts. The N—H S values ranges from 0 to 1, where scores above 0.8 infers structural rigidity while scores below 0.8 represents flexibility [53]. The coarse dynamics represents residue-level fluctuation concerning the most stable protein conformation. The protein stability was estimated from the root means square fluctuation (RMSF) trajectory upon merging coarse dynamics simulation and consensus protein fluctuations in an aqueous environment using default parameters of all-atom force-field with explicit water model, for a timeframe of 10 ns. The simulation was carried out under default settings of conformational distances (minimum: 3.8 Å; maximum: 8.0 Å) and gap = 3 which represents the minimum distance between previous and next amino acid residue to be restrained [54].

The optimized protein was finally subjected to MDS analysis for estimating protein stability upon mimicking an aqueous cellular environment using GROMACS 2020.2 package. Protein topology was built using CHARMM36-Feb2021 force-field and simple-point charge (TIP3P) water model. Protein solvation was performed upon centering the protein inside a cubic box of uniform dimension (1.0 nm) followed by system neutralization with the addition of requisite counter ions (Na+/Cl-). Energy minimization was carried out using steepest descent algorithm for 50,000 steps and convergence tolerance force of 1000 kJ/mol nm−1. Further, equilibration of the system was achieved under standard NVT (constant number of particles, volume and temperature) and NPT (constant number of particles, pressure and temperature) ensembles for 100 ps. Later, MD production was carried out for a time-scale of 75,000 ps using particle-mesh Ewald electrostatics summation and Parrinello-Rahman extended coupling ensemble as previously described [19]. The trajectory data were further visualized using the Grace software.

Virtual screening, pharmacokinetic filtering and ligand optimization upon DFT simulations

Based on text mining, therapeutic target database and DrugBank screening, conventional inhibitor against the target protein was considered as the reference molecule for screening commercial analogues from the ZINC database using SwissSimilarity [56] server. The server aids in rapid ligand based virtual screening from enormous small molecule/ligand libraries. Based on the structural homology descriptors, the similarity is quantified by Tanimoto coefficient which is the ratio between the common positive bits to the total positive bits between the reference drug and predicted molecular fingerprint. The lead molecules were funneled down based on Tanimoto coefficients and pharmacokinetic parameters from SwissADME web-tool that provides a rational outlook to the diverse physicochemical properties, drug-likeliness and medicinal chemistry parameters besides accessing the permissible ranges of absorption, distribution, metabolism and excretion (ADME) [57]. In-silico toxicity predictions including oral toxicity assessments (LD50, toxicity class) and toxicity endpoints (hepatotoxicity, acute toxicity, mutagenicity, carcinogenicity, immunotoxicity, adverse pathways and toxicity targets) were assessed from the ProTox-II server using pharmacophore modelling, molecular homology and ML algorithms [58].

The atomic valences of the shortlisted 3D leads were satisfied upon the addition of polar H-atoms, while spatial restraints were enhanced by universal force-field mechanics using the steepest-descent algorithm. These optimisations were carried out in a powerful standalone Avogadro v1.2.0 platform used extensively for molecule editing, visualization and analysis [59]. Further, the chemical reactivity and stability of the leads were ascertained through DFT simulations using Gaussian-09 software for quantum chemical calculations and the results were subsequently visualised in GaussView v6.0.16 software [60], [61]. Molecular geometry optimization, electron density mapping and frontier orbital (HOMO-LUMO) calculations were achieved using Becke's three parameter exchange–correlation function (B3) conjugated with Lee-Yang-Parr (LYP) and 6–311++G(d,p) basis level for estimating its minimal energy configurations [62], [63].

Molecular docking and stability analysis upon MDS

Site-specific molecular docking was performed concerning the active-site residues which were predicted using the CASTp server that scans the surface topography of the protein. The server simultaneously predicts the superficial, deep-seated pockets and cross channels which might serve as potential active-site for drug binding [64]. Therefore, subsequent druggability validation of these grooves is predicted through the DoGSiteScorer server uses supervised ML algorithm for druggability estimation. The druggability score ranges from 0 to 1 where higher scores represent better druggability profile of the concerned groove [65].

The validated drug-target devoid of unwanted hetero-atoms, crystallographic water molecules was subjected to AutoDock 4.2 and its embedded tools that uses semi-empirical free energy scoring function to predict the binding energies of docked poses with high reproducibility [66]. The free-end residues of protein moiety were stabilized by adding polar H-atoms, requisite Kollman charges and merging non-polar H-atoms. The ligand torsions of optimized reference drug and shortlisted leads were similarly fixed upon merging non-polar H-atoms and addition of Gasteiger charges. The active-site of protein was centered inside affinity grid-box of uniform dimension (60x60x60Å3) with grid-point spacing of 0.375 Å. Thereafter, autogrid4 and autodock4 programs were executed that respectively scanned and generated the docked protein–ligand complex. Lamarckian and genetic algorithms were opted for selecting complex with the least binding energy (BE) and represented as Mean ± SD [66]. The 2D conformer of the protein–ligand docked complexes were visualized in UCSF-Chimera v1.9, an efficient molecular visualization and editing program [67]. The 3D docked poses were viewed in Discovery Studio Visualizer v20.1.0.19295 which provides an user friendly protein modeling, pharmacophore analysis and drug designing platform [68]. The classical drug and shortlisted leads were docked (in triplicates) separately with the proposed drug-target [23], [24], [55].

Based on the binding affinities and inhibition constants (IC), screened protein-inhibitor complexes were subjected to MDS analyses. GROMACS (GROningen MAchine for Chemical Simulation) suite is a highly efficient, open-source and flexible MDS suite integrating tutorials for simulating simple protein in water to highly complex simulations incorporating membrane protein, biphasic system, protein–ligand complex and construction of virtual sites. In the present study, stability of therapeutic target-ligand/lead complexes was performed in an aqueous environment for a timeframe of 75,000 ps using GROMACS 2020.2 package. Protein topology was built as described previously (section 2.4) while ligand topology was built using the CGenFF server. The complex was centered inside a dodecahedron box of uniform edge-distance (1.0 nm) following solvation with simple point-charge water model and system neutralization with requisite counter ions (Na+/Cl-). Steepest-descents algorithm (50,000 steps) and convergence-tolerance force (1000 kJ/mol nm−1) were considered for energy minimization of the system. Thereafter, two cycles of equilibrations were opted with constant volume (NVT) ensemble for 100 ps using leap-frog integrator for attaining desired temperature (300 K). Secondly, constant pressure (NPT) ensemble for 100 ps using Parrinello-Rahman barostat was applied for attaining desired pressure (1 bar) upon applying motion-equations to the box vectors. Long-range electrostatic interactions were treated using particle-mesh Ewald algorithm with a cubic interpolation of order 4.0 and Fourier spacing of 0.16 nm. Finally, the system was subjected to MD production for 75,000 ps timescale with integration timescale of 2 fs and sampling of simulated trajectories were carried out every 10 ps [21], [69], [70].

Results

GIN Analysis, hub-protein identification and validation

The retrieved interactome profile encompassed 49 genes specific to PD machinery having 525 interactions among them. Both experimentally validated and computationally predicted interactions were considered with a confidence score (CS) ranging from 0.40 to 1.00 and statistical significance of P-value < 0.05. All the 49 genes interacted with each other making 86, 47, 136 and 256 interactions with the highest, high, medium and low CS respectively [Supplementary Table S1]. FEA resulted in 658 GO terms (525 BP; 65 MF; 66 CC); 26 pathway enrichments (05 KEGG; 21 Reactome) and 77 domain descriptors (29 InterPro; 10 Pfam; 05 SMART; 33 UniProt keywords) [Supplementary Table S2].

Clustering analysis of 49 genes resulted in four densely interconnected clusters (C1-C4). The clusters C1 (20 nodes; 182 edges), C2 (06 nodes; 12 edges), C3 (14 nodes; 29 edges) and C4 (04 nodes; 06 edges) obtained MCODE scores of 19.16, 4.80, 4.46 and 4.00 respectively as illustrated in Fig. 1(a). The functionalities associated with the clustered proteins were highly similar as evident from nine C1; three C3 and four C4 proteins being involved in the regulation of apoptosis [GO:0042981], a major biochemical pathway in PD [Fig. 1(b)]. The cluster C3 maintained the oxidative stress, kinase activity, ion homeostasis [GO:0043549; GO:0045859; GO:0071900] besides strongly interacting with the cluster C2 [Fig. 1(c)]. The cluster C1 was mainly involved in the regulation of mitochondrial oxidative stress [GO:0006979; GO:1903202], maintaining appropriate levels of α-synuclein, protein ubiquitination [GO:0031397; GO:0031398], regulation of phosphorylation and neurotransmitter transport [GO:0001505; GO:0006836; GO:0046928; GO:0051580] [Fig. 1(d)]. Cluster C4 regulated protein folding mechanisms [GO:0006457; GO:0042026; GO:0061077], while C2 denoted regulation of translational activity, cellular and biochemical metabolism [GO:0006417; GO:0006413; GO:0006446] [Supplementary Table S2].

Fig. 1

Gene Interaction Network Analysis: (a) MCODE clustering analysis of 49 PD associated genes (b) Hub genes and associated genes involved in apoptosis [GO:0042981; GO:2001233] (c) Hub genes and associated genes involved in oxidative stress [GO:2000377; GO:0006979; GO:0034599] (d) Hub genes and associated genes involved in mitophagy [GO:0010821; GO:1903599; GO:0051881; GO:1903146].

Topological analysis highlighted LRRK2 exhibiting the highest number of direct interactions (39) followed by PARK7, SNCA, PARK2 and PINK1 with 38, 36, 33 and 30 interactions respectively. The tendency of a gene to be considered as the hub molecule was majorly sorted based on betweenness centrality, clustering coefficient and degree of interaction. From the comparative analysis of topological metrics, CSs of direct interactors and functionality, LRRK2, PARK2, PARK7, PINK1 and SNCA favoured multiple criteria and thus can be considered as hub-genes. Based on the above analyses and interaction confidence (reliability), a scoring metrics was generated for the top five hub-genes, all belonging to the C1 cluster, as tabulated in Table 1. With an objective to interfere with the α-synuclein aggregation, it was observed that PARK2 (Parkin) made direct association with SNCA with CS > 0.990, the highest amongst all other predicted hub-genes. Thus, Parkin was prioritized as an important candidate for targeted therapy and was subjected for structural evaluations.

Table 1

Topological parameters, scoring metrics (based on CS) and sub-cellular locations of PD hub-genes from GIN analysis.

HUB GENES	TOPOLOGICAL PARAMETERS			CONFIDENCE SCORES (CS)					SUBCELLULAR LOCATION
	Betweenness Centrality	Clustering Coefficient	Degree	LRRK2	PARK2	PARK7	PINK1	SNCA
LRRK2	0.82	0.24	39						Membrane
PARK2	0.36	0.29	33	0.987					Cytoplasm
PARK7	0.86	0.74	38	0.958	0.979				Mitochondria
PINK1	0.34	0.32	30	0.966	0.999	0.988			Mitochondria
SNCA	0.55	0.27	36	0.971	0.992	0.990	0.922		Membrane

In-silico drug-target modeling and refinement

The unavailability of a well-defined 3D Parkin structure prompted us to cure the stretches of missing residues using dual step modeling method. BLASTp search of Parkin (UniProt ID: 060260) against RCSB-PDB highlighted 5C1Z with query coverage and sequence similarity of 100 % and 79.78 % respectively. Thereafter, Parkin was modeled using dual-step modeling approach using SWISS-MODEL and MODELLER using 5C1Z as a template. The modeled structure was energy minimized and refined resulting in high Ramachandran favoured regions (98.5 %) and low poor rotamers (0.8 %). The refined 3D conformer of Parkin was finally submitted to Protein Model Database (PMDB) bearing PMDB-ID:PM0084200 [Fig. 2(a)]. Gross structural characteristics explained in terms of global model quality exhibited Z-score of −8.87 [Fig. 2(b)].

Fig. 2

Structural Analysis of Parkin: (a) Optimized modeled structure (b) Global model quality (c) Local model quality (d) Atomic-level fluctuations (normalized B-factor), solvent accessibility and secondary structural analysis (e) Propensity plot (f) ProTSAV heat-map (g) Extent of disorderness.

Domain configurations of Parkin

After the identification of parkin as a targetable key regulatory protein, the conformational dynamics (structural stability/ dynamics) was evaluated to validate the quality of the structure used. Subsequently, drug-pocket was characterized prior to screening novel therapeutic agents against it. It was observed that the local energetics curve of functional domains of Parkin was positioned well below the threshold cut-off (0.00) [Fig. 2(c)]. Inherent thermal motions expressed as normalized B-factors revealed average atomic-level fluctuations of −0.07 Å2. Secondary structure prediction highlighted 63.65 % solvent accessible and 36.35 % buried regions, while the structure is majorly composed of coils (62.80 %), β-sheets (27.10 %) and α-helix (10.10 %) [Fig. 2(d)]. Parkin laid upon the HARMONY propensity-calibration curve [Fig. 2(e)] while substitution curve revealed minimum intersections between the forward and reverse sequences (graph not shown). The overall protein structure quality assessed from ProtSAV heatmap revealed the structure within RMSD 2–5 Å [Fig. 2(f)], whereas the average domain disorderness was estimated to be ∼ 0.20 positioning below the threshold cut-off (0.50) [Fig. 2(g)].

The structural descriptors based on several kinetic and thermodynamic constraints deduced RING0 domain having high energy requirement for inter-residue distances (ΔGdist = -0.257 kcal/mol) and backbone torsion angles (ΔGtor=+0.683 kcal/mol) when compared with RING1 (ΔGdist = -0.611 kcal/mol; ΔGtor=+0.108 kcal/mol) and RING2 (ΔGdist = −0.349 kcal/mol; ΔGtor=+0.149 kcal/mol) domains. In terms of solvent accessibility potential, the interaction of residues in RING0 (ΔGacc = −0.031 kcal/mol) and RING1 (ΔGacc = −0.305 kcal/mol) domains were energetically favoured. [Fig. 3(a)]. The N—H S2 bond-order restraints based on NMR chemical shifts displayed rigid backbone dynamics of Parkin protein (average S2=∼0.80), while RING0, RING1 and RING2 domains exhibited average backbone dynamics of S2 = 0.82, 0.87 and 0.80 respectively [Fig. 3(b)]. The coarse-grained simulation revealed the average residual-level root means square fluctuations of RING1, RING0 and RING2 domains to be 0.68 Å, 0.97 Å and 1.07 Å respectively, whereas the RMSF trajectory for the entire protein was computed to be 1.22 Å [Fig. 3(c)].

Fig. 3

Stability Analysis of Parkin: (a) Global folding free energy of thermodynamic and kinetic constraints of Parkin (b) Backbone stability profile of Parkin (c) Residue-level fluctuation profile of Parkin.

Stability analysis of Parkin through MDS

The MDS analysis profile of unbound parkin was taken as a standard for evaluating the relative deviation in the structures of inhibitor-bound complexes. MDS trajectories revealed a stable root means square deviation (RMSD) curve concerning protein backbone attaining equilibrium at 0.68 ± 0.09 nm that was maintained throughout the studied time-frame [Fig. 4(a)]. The RMSF values of RING0 (0.22 ± 0.08 nm), RING1 (0.18 ± 0.06 nm) and RING2 (0.21 ± 0.06 nm) as well as for the overall protein (0.27 ± 0.09 nm) were noted [Fig. 4(b)] while a stable radius of gyration curve maintained at 2.73 ± 0.04 nm [Fig. 4(c)]. The conformational coherence between proximal backbone residues was calculated to be 4.89 ± 0.65 nm (>2.0 nm) [Fig. 4(d)]. Persistent intermolecular (protein-solvent) hydrogen bonds (∼1215) [Fig. 4(e)] and steady trajectory of solvent accessible surface area (SASA; 282.44 ± 6.01 nm2) were observed for the Parkin protein [Fig. 4(f)]. The system energetics represented in terms of potential energy (−2.3e + 06 kJ/mol) [Fig. 4(g)] and total energy (-1.8e + 06 kJ/mol) were computed [Fig. 4(h)].

Fig. 4

MDS analysis of unbound Parkin: (a) RMSD curve (b) Residue-level RMSF plot (c) Rg trajectory (d) Minimum distance amongst proximal backbone residues (e) Number of intermolecular (protein-solvent) hydrogen bonds (f) SASA trajectory (g) Potential energy curve (h) Total energy curve.

Lead sorting, ADMET screening and chemical reactivity prediction analysis

Ligand-based virtual screening resulted in 400 commercial analogues of the classical drug UDCA (PubChem ID: 31401) due to its reported activity against Parkin protein [Supplementary Table S3(a)]. The analogues were initially sorted upon filtering the isomers resulting in 58 unique lead molecules [Supplementary Table S3(b)]. ADMET screening was thereby performed narrowing down the data based on positive blood–brain barrier (BBB) permeability, as the leads were destined to penetrate the human brain cells [Supplementary Table S3(c)]. Further, refinements were based on high bioavailability score, high LD50 value, high GI absorption, total polar surface area (20–130 Å2), low synthetic accessibility and minimal Lipinski violations [Supplementary Table S3(d-e)]. The screening resulted in four lead molecules namely Molecule-258 (ZINC5022267), Molecule-297 (ZINC3846933), Molecule-309 (ZINC3846931) and Molecule-371 (ZINC5289889) to be showing favourable scores [Table 2]. We observed the LD50 (2000 mg/kg) and bioavailability (0.56) profiles of UDCA ranged lower compared to the four shortlisted drugs with ∼ 63 % and ∼ 52 % higher LD50 and bioavailability profiles respectively. Amongst all, Molecule-258 and Molecule-271 displayed no reported toxicity and therefore their efficacy was critically scrutinized for declaring potential therapeutic candidates.

Table 2

Pharmacokinetic profiles of UDCA and four shortlisted lead molecules.

Molecule	TPSA	GI Absorption	BBB Permeability	Lipinski Violations	Bioavailability Score	Synthetic Accessibility	LD50	Remarks
UDCA	77.76	High	No	No	0.56	4.93	2000	Hepatotoxicity, MMP
Molecule-258	57.53	High	Yes	No	0.85	4.14	3265	–
Molecule-371	54.37	High	Yes	No	0.85	2.86	3265	–
Molecule-297	57.53	High	Yes	No	0.85	3.91	3265	Hepatotoxicity, MMP
Molecule-309	57.53	High	Yes	No	0.85	3.91	3265	Hepatotoxicity, MMP

 = Total Polar Surface Area; = Gastro-Intestinal; = Blood Brain Barrier; = Lethal Dose;Mitochondrial Membrane Potential.

The structures of four shortlisted leads and reference compound (UDCA) were subjected to stereochemical enhancements with improved torsions for deciphering their chemical reactivity parameters through DFT structure optimization simulations. The energy difference (ΔE = ELUMO-EHOMO) between the frontier orbitals of LUMO and HOMO in Molecule-258, Molecule-297, Molecule-309, Molecule-371 and UDCA were computed to be ΔE = 6.93 eV, ΔE = 6.75 eV, ΔE = 6.71 eV, ΔE = 5.85 eV and ΔE = 7.36 eV respectively [Fig. 5]. The optimized structures of UDCA along with four shortlisted leads with atom numbering scheme and DFT profile of Molecule-309 are provided in Supplementary Figures S4(a-e) and S5(a) respectively.

Fig. 5

DFT simulations highlighting frontier molecular orbital (HOMO-LUMO) and electron density map of: (a) UDCA (b) Molecule-258 (c) Molecule-297 (d) Molecule-371.

Druggability prediction of Parkin and molecular docking analysis

Surface topology scan of Parkin revealed a groove (Area: 1176.92 Å2; Volume: 1024.15 Å3) encompassing the active-site residues from RING0 (residues: 239–241) and RING1 (residues: 259, 263, 266–267, 270–271, 274, 276, 288–292) domains. The pocket revealed high druggability potential (0.81) and therefore was considered as the active-centre for performing molecular docking analysis. The surface view of Parkin highlighting its RING0, RING1 and RING2 domains are illustrated in Fig. 6(a). Intermolecular interactions of Parkin_UDCA complex (BE = -7.28 ± 0.05 kcal/mol) displayed the highest binding affinity with minimum IC of 4.62 µM. Amongst the four shortlisted leads, Molecule-258 (BE = -6.57 ± 0.1 kcal/mol; IC = 15.35 µM) had better binding profiles when compared to Molecule-297 (BE = -6.12 ± 0.25 kcal/mol; IC = 34.66 µM), Molecule-371 (BE = -5.44 ± 0.2 kcal/mol; IC = 34.66 µM) and Molecule-309 (BE = -5.32 ± 0.3 kcal/mol; IC = 125.65 µM) [Fig. 6(b)].

Fig. 6

Molecular Docking Profiles: (a) 3D conformer of Parkin highlighting its RING domains (b) Binding energies and Inhibition Constants of docked complexes (c) Intermolecular interaction profile of UDCA (d) Intermolecular interaction profile of Molecule-258 (e) Intermolecular interaction profile of Molecule-297 (f) Intermolecular interaction profile of Molecule-371.

Docking pose of Parkin_UDCA complex revealed two H-bonds with the electronegative O1-I239 and O2-R396 residues. The complex was further stabilized by van der Waals (vdW) interactions with solvent-accessible residues K48-E49, T240, T242 and several non-canonical alkyl-π interactions between electropositive atoms of cyclohexane-cyclopentane (C5-C8, C10, C12-C13, C15-C16, C19) ring and A397, A402, C241 residues [Fig. 6(c)]. The intermolecular interactions of Parkin_Molecule-258 complex were highly similar to that of UDCA as K48-E49 residues from Ubl substrate recognition domain formed alkyl-π interactions and H-bond respectively besides forming a second H-bond with O2-D394 residue. The RING0 domain residues T240 and T242 formed vdW interactions while C241 formed alkyl-π interactions with the cyclohexane rings consisting C4-C7, C9-C10 and C12-C15 atoms similar to the parkin-UDCA complex [Fig. 6(d)]. Parkin_Molecule-297 displayed three H-bonds between O1-T242, O3-K48 and O2-D394 residues, alkyl-π interactions (C4-C6, C8, C10-C11 atoms) with C241 and vdW interactions with T242 residues of RING0 domain [Fig. 6(e)]. Though safe toxicity profile was obtained for Molecule-371, no correlation could be drawn from the intermolecular interacting residues of Parkin_Molecule-371 complex with UDCA or other studied complexes. A sharp conformational change (dihedral angle −81.49⁰) was observed in C12-C15 ligand backbone atoms that deviated from the parent molecule (dihedral angle 179.71⁰) in the docked pose [Fig. 6(f)]. In Parkin_Molecule-309 complex, although involvement of K48, T240-C241 residues in forming H-bonds and T242 in vdW interactions were displayed, the lowest binding affinity (∼27 % lower than UDCA) was recorded [Supplementary Figure S5(b)]. Hence, Molecule-309 and Molecule-371 were excluded from the current study whereas Molecule-258 and Molecule-297 complexed with Parkin protein were prioritized for stability analysis through MDS analysis.

Stability analysis of Parkin-inhibitor complexes through MDS

Based on the binding profiles and intermolecular interaction patterns of Molecule-258 and Molecule 297 with Parkin, stability of these complexes was further ascertained upon MDS and compared with Parkin_UDCA complex. The RMSD curves indicated that Parkin complexed with UDCA (0.61 ± 0.1 nm), Molecule-258 (0.63 ± 0.12 nm) and Molecule-297 (0.55 ± 0.09 nm) attained equilibrium which was maintained throughout the studied timeframe [Fig. 7(a)]. The positional fluctuations of all residues in simulated system were measured in terms of RMSF concerning the protein backbone atoms. Parkin_Molecule-297 complex exhibited relatively lower RMSF profiles (0.24 ± 0.10 nm) when compared with UDCA (0.28 ± 0.16 nm), Molecule-258 (0.30 ± 0.17 nm) complexes. However, when the residual fluctuations were considered for the interacting RING0 domain, least fluctuations were recorded in Parkin_Molecule-258 (0.2 ± 0.05 nm) followed by Parkin_Molecule-297 (0.23 ± 0.08 nm) while maximum residual disturbances were noted in Parkin_UDCA (0.27 ± 0.12 nm) complexes. Further, the integration of UDCA in active groove of RING0 resulted in relatively high RMSF curves (0.27 ± 0.12 nm) when compared to the fluctuation profiles of Parkin_Molecule-258 (0.20 ± 0.05 nm) and Parkin_Molecule-297 (0.23 ± 0.08 nm). Higher fluctuations were also noted in the RING1 and RING2 domain while the protein was bound with UDCA (11.11 %; 47.62 %), Molecule-258 (33.33 %; 23.81 %) and Molecule-297 (11.11 %; 14.29 %) respectively as compared with the unbound protein [Fig. 7(b)]. The Rg trajectory concerning Cα atoms denoted an average value of 2.56 ± 0.02 nm for Parkin_Molecule-297 while Parkin_UDCA (2.62 ± 0.03 nm) and Parkin_Molecule-258 (2.63 ± 0.04 nm) maintained similar Rg trajectory [Fig. 7(c)]. The intermolecular H-bonding patterns from MDS further validated the molecular docking analysis revealing that UDCA formed three H-bonds, while Molecule-258 and Molecule-397 each formed two stable H-bonds with the Parkin protein [Fig. 7(d)]. The SASA plot of Parkin_UDCA (∼277.25 nm2) and Parkin_Molecule-258 (∼268.49 nm2) revealed similar convergence profile upto 60,000 ps while Parkin_Molecule-297 displayed SASA trajectory with reduced hydrophilicity (∼262.8 nm2) that was constantly maintained throughout the simulation [Fig. 7(e)]. The free energy of solvation revealed that Parkin_Molecule-258 (-51.82 kJ/mol) complex was the most energetically favourable compared to Parkin_Molecule-297 (-45.93 kJ/mol) and Parkin_UDCA (-45.36 kJ/mol) [Fig. 7(f)]. The interaction energy revealed energetically stable (<0kJ/mol) profiles for Parkin_UDCA (-91.55 kJ/mol), Parkin_Molecule-297 (-53.12 kJ/mol) and Parkin_Molecule-258 (-48.55 kJ/mol) complexes [Fig. 7(g)] while the total energy of all the three simulated systems converged at ∼ -1.288e + 06 kJ/mol [Fig. 7(h)].

Fig. 7

MDS analysis of Parkin-Inhibitor Complexes: (a) RMSD curve (b) Residue-level RMSF plot (c) Rg trajectory (d) Number of intermolecular (protein-inhibitor) hydrogen bonds (e) SASA trajectory (f) Free energy of solvation (g) Interaction energy profile (h) Total energy curve.

Discussion

The present study was dedicated to propose potential therapeutic interventions against PD by means of integrated network and structural biology approaches. The constructed GIN highlighted distinct biomolecular clusters whereas FEA and topological analyses resulted in the better understanding of the network and eventual identification of the hub-genes. The hubs were screened based on two crucial topological parameters namely betweenness centrality and clustering coefficient. Betweenness centrality identifies the central point in the interactome while the clustering coefficient indicates the tendency of a node to cluster. Each gene cluster (C1-C4) was enriched with similar GO terms denoting specific biochemical pathways significant to PD. The E3-ubiquitin ligase PARK2 (33 interactions) leads to ubiquitination and proteasomal degradation while maintaining mitochondrial homeostasis. PINK1 and PARK7 are the key players in mitochondrial dysfunction machinery being a target of reactive oxygen species generated from electron transport system [GO:0010821; GO:1902958; GO:0098779]. Regulation of mitochondrial organization and autophagy [GO:0010821; GO:1903599] is of utmost necessity in preventing PD, as increased oxidative stress [GO:2000377; GO:0034599] results in mitochondrial toxicity and dysfunction.

Simultaneously, excessive α-synuclein in PD patients reduces dopamine uptake [GO:0051583; GO:0014059]. PARK2 has substantial functional interaction with PINK1 (0.999), which prevents mitochondrial oxidative stress in association with PINK1 facilitating the mitochondrial translocation of PARK2 for Parkin-mediated mitophagy [R-HSA-5205685] and Ubiquitin-mediated proteolysis [GO:0032434; GO:0032436; GO:0031398]. From the comparative analysis of the CS, topological metrics, functionality and direct interactors, LRRK2, PARK2, PARK7, PINK1 and SNCA favoured multiple criteria and thus can be considered as hub-genes. These five genes were moreover involved in regulatory responses to drugs [GO:0042493; GO:2001023], hence implying the genes to be considered as candidate therapeutic biomarkers. From the scoring metrics, PARK2 was singled out due to its highly experimentally validated functional association (CS > 0.990) with SNCA involved in α-synuclein aggregation (as evident from Table 1). Moreover amongst the five hub-genes, PARK2 (expressed product being Parkin protein) interacted with other hub-genes with the highest average CS of 0.989 followed by PARK7 (0.979), LRRK2 (0.971), PINK1 (0.969) and SNCA (0.946). Therefore, Parkin protein was noted to hold the central point and served as the bridge among other hub-genes in the interactome. Therefore, Parkin was adjudged as a better candidate for targeted therapy.

Text mining revealed that Parkin has already been patented to be a potent diagnostic biomarker in PD progression [71]. On the other hand, a recent study summarizing therapeutic targets at various stages of PD did not suggest the role of Parkin as therapeutic biomarker [72]. Hence, the present study focused on exploring the potency of Parkin as a therapeutic target. With an aim to design novel therapeutic leads against Parkin, extensive database and literature mining revealed the inhibitory activity of UDCA, Nilotinib and Simvastatin against Parkin (PARK2). Parkin inhibition as supported by literary evidence could be easily targeted being present in the cytoplasm or Golgi-complex of human brain [73], hence Parkin structure modeling and refinement were executed. The N-terminal of Parkin comprises of Ubl domain (amino-acid residues: 1–76), followed by RING0 (amino-acid residues: 141–255), RING1 (amino-acid residues: 256–327), IBR (amino-acid residues: 328–378) and RING2 (amino-acid residues: 328–378) domains at C-terminal [74]. Enzyme E1 activates ubiquitin via energy-dependent process that is further transferred to the cysteine residue of E2 conjugating enzyme. Activated ubiquitin is subsequently attached by E3-ubiquitin ligase (Parkin) to the lysine residue of the target protein (α-synuclein), thereby degrading the LB [75]. However, literary evidence suggests auto-inhibition of the catalytic RING2 domain by RING0 domain [74], [76] while several in-vitro studies revealed the high ubiquitin potency of Parkin upon deletion of RING0 domain [77], [78]. Therefore, subsequent analyses were carried out to understand the stability and dynamics of the Parkin-RING domains.

The global model quality revealed a well-poised structure of Parkin amongst all experimentally characterized proteins revealing minimal erroneous structural folds [Fig. 2(b)]. Negative ProSA-web energetics curve and HARMONY propensity-calibration plot indicated nominal local errors and structural misfolds revealing high protein stability [Fig. 2(c, e)]. Normalized B-factor revealed minimum atomic-level fluctuations due to thermal mobility and positional disorders [79]. The RING0 domain unique to Parkin is primarily composed of coils (64.35 %) and appeared to be the most solvent accessible (∼67 %) amongst other RING domains being actively involved in auto-inhibition of the catalytic RING2 domain [Fig. 2(d)] [74]. Protein structural quality as interpreted from ProtSAV heatmap established Parkin to be ideal for protein–ligand interaction analysis [46] [Fig. 2(f)]. Although low domain disorderness was noted in the overall structure, slight disorderness was observed in the residues of RING0 (0.142) in comparison to RING1 (0.076) and RING2 (0.038) domains could possibly arise due to the intramolecular interactions of the RING0 residues with other domains [Fig. 2(g)]. Kinetic and thermodynamics-based structural parameters revealed high energy contribution of the RING0 domain which in turn signifies functionally active residues involved in the interaction with other domains, macromolecular and micromolecular species [52]. The residual-level fluctuations of Parkin domains gave an elaborate perspective about their stability and functionality. Rigid backbone dynamics of overall Parkin protein implied restricted movements of the atomic bond vectors that were evident from the evolutionary conserved domains of E3 ubiquitin ligases. Due to macromolecular interaction in executing auto-inhibition mechanism of RING2 domain by RING0 domain, the latter acts as a scaffold resulting in high rigidity (S2 = 0.82) of the domain [76]. The E2-binding domain (RING1; S2 = 0.87) and catalytic domain (RING2; S2 = 0.80) had reduced flexibility being highly evolutionary conserved among other E3-ubiquitin ligases [Fig. 3(b)]. The coarse-grained dynamics simulation revealed the average residual-level RMSF of RING1, RING0 and RING2 domains to be lower than the entire RMSF trajectory by 44.60 %, 20.33 % and 12.54 % respectively [Fig. 3(c)]. These fluctuations justify the overall stability and rigidity of the Parkin domains, further complemented by the reduced flexibility profiles of the N—H S2 bond-order restraints.

MDS revealed a stable RMSD trajectory for Parkin as well as minimal positional fluctuations further complementing the N—H S2 backbone dynamics and coarse-grained dynamics curves [Fig. 4(a-b)]. The compactness in folding patterns of Parkin was ascertained by a stable Rg curve further validating the domain stability strengths represented by the negative folding free energy curves of ΔGdist from SWOTEIN server [Fig. 4(c)]. The conformational coherence between proximal backbone residues indicated stereochemically compact tertiary structure conformation with negligible steric hindrance further strengthening the findings from ProSA-web, HARMONY and ProTSAV servers [Fig. 4(c-d)]. Parkin formed a uniform layer of solvation upon interaction with the exposed residues as evident from the formation of persistent intermolecular hydrogen bonds and steady trajectory of SASA profile [Fig. 4(e-f)]. The system energetics signified environmentally favoured trajectories for the entire timeframe establishing the stable conformation of the Parkin protein in aqueous environment [Fig. 4(g-h)] [80], [81]. The MDS analyses thus validated the overall stability profile of the therapeutic target.

Ligand-based virtual screening of anti-Parkin drug UDCA analogues resulted in four lead molecules namely Molecule-258, Molecule-297, Molecule-309 and Molecule-371 which were sorted based on BBB permeability, high bioavailability, less toxicity, passive absorption by GI-tract, optimal total polar surface area, easily synthesizable and no Lipinski’s violations. The reference drug UDCA responded negatively to BBB permeability which could be a probable reason for increased risk of hepatocellular carcinoma and Tox21 stress response mitochondrial membrane potential (MMP) toxicity due to high dosage administration for effective drug binding [9]. Out of the four leads, Molecule-258 and Molecule-271 exhibited safe toxicity profiles while Molecule-297 and Molecule-309 reported similar toxicity profiles of MMP and hepatotoxicity similar to that of UDCA. Frontier molecular orbital (HOMO and LUMO) calculated from DFT structure optimization simulations gave functional insights about the chemical stability and reactivity of the leads. Molecule-258 was found to be chemically more reactive (by ∼ 6 %) than UDCA and possessed better stability profiles when compared to the ΔE profiles of Molecule-297, Molecule-309 and Molecule-371. Electron density map of UDCA highlighted four electronegative oxygen atoms (O1-O4) as potent centres for electrophilic attack (red zones) while other leads possessed three electrophilic centres over O1-O3 atoms. The C—C ligand backbone corresponded to high positive potential apt for nucleophilic attack (blue zones) [Fig. 5]. The electrochemical nature of these leads contributing to various intermolecular interactions with Parkin was further interpreted in molecular docking analysis.

The residue T240, being a key factor for E2 binding from the highly druggable groove encompassing RING0 and RING1 domains, was considered as the active-centre for performing molecular docking analysis [82]. As illustrated in Fig. 6(a), RING0 lies in close vicinity to RING2 domain, and therefore, molecular docking was carried out in order to inhibit the RING0 domain, thus restoring the catalytic function of RING2 domain. From the present study, it was clearly evident that Molecule-258 and Molecule-297 showed similar interaction profiles as compared to the classical drug [Fig. 6(c-e)], where two residues (K48-E49) from Ubl domain and three (T240-C241-T242) residues from RING0 domain played predominant roles in the interaction profiles. Though closely placed, no residues from RING2 domains contributed to intermolecular interactions with the lead molecules thereby leaving the catalytic RING2 domain undisturbed. Though safe toxicity profile was obtained for Molecule-371, no correlation could be drawn from the intermolecular interacting residues of Parkin_Molecule-371 complex with UDCA or other studied complexes [Fig. 6(f)]. Amongst the four UDCA analogues, Molecule-258 and Molecule-297 exhibited superior binding profiles with Parkin, while Molecule-309 displayed the least binding affinity [Fig. 6(b)]. Distinct conformational change of Molecule-371 in the docked pose (by 261.2⁰) was detected which might be a probable reason contributing to the disparity in the binding profile. Due to these following reasons, Molecule-309 and Molecule-371 were not considered as anti-PD lead molecules in the current study. Therefore, we can conclude that Molecule-258 and Molecule-297 could be proposed as novel therapeutic leads in targeting specifically the Parkin-RING0 domain thereby restoring the functionality of RING2 domain.

MDS revealed the stability and compactness of the protein-inhibitor complexes as predicted from the RMSD, RMSF and Rg trajectories. Parkin_Molecule-297 complex was slightly rigid (relatively lower fluctuations) while minimal backbone deviations were observed in Parkin_UDCA and Parkin_Molecule-258 complexes. Overall the structural convergence of all three RMSD trajectories was observed after 52,500 ps implying stability of the complexes. Higher fluctuation profiles in the RING0 domain of Parkin_UDCA, Parkin_Molecule-258 and Parkin_Molecule-297 complexes by 22.72 %, 9.09 % and 4.54 % respectively were observed when compared to the unbound Parkin protein. Therefore, integration of the proposed UDCA analogues had minimal residue-level fluctuation than UDCA confirming their stability. Moreover, the higher fluctuations in the catalytic RING2 domain of Parkin-inhibitor complexes further elucidate its flexibility and dynamicity, which might be due to the restoration of the RING2 domain. Through high structural compactness was displayed by the Parkin_Molecule-297 complex, structural convergence in both Parkin_Molecule-258 (2.55 ± 0.01 nm) and Parkin_Molecule-297 (2.55 ± 0.02 nm) was evident after 65,000 ns and was maintained throughout the timeframe [Fig. 7(a-c)].

Molecular docking analysis and MDS were well correlated giving a clear insight about the protein-inhibitor interaction and catalytic functions in the studied timeframe. From the intermolecular H-bonding analysis, it was clearly evident that three H-bonds formed by UDCA and two H-bonds formed each by Molecule-258 and Molecule-397 with Parkin corresponded to ∼ 96 % and ∼ 99 % of total stable H-bonds formed within the complexes [Fig. 7(d)]. The knowledge of exposed and hydrophobic buried core interpreted from secondary structural analysis was further analyzed upon studying the SASA profiles signifying the compactness of the hydrophobic core [Fig. 7(e)]. Conformational geometry interpreted from the Rg, SASA trajectories and energy curves revealed uniform concurrence between Parkin_UDCA and Parkin_Molecule-258 complexes with minimal alterations and energetically favoured global conformation of Parkin [Fig. 7(f-h)]. Finally based on all these parameters, we can conclude that the inhibitor Molecule-258 displayed stable interaction profiles comparable to that of UDCA and therefore can be concluded as novel therapeutic lead molecule inhibiting the RING0 domain of Parkin. Though equipotent interaction and simulation profiles were displayed for Molecule-397, less binding affinity (by ∼ 16 %) and notable toxicity profiles (Hepatotoxicity, MMP) were recorded and therefore not recommended for therapeutic purposes in the present study.

Conclusion

Constraints in early detection and targeted therapy evoked us to investigate therapeutic solutions against PD. Parkin (PARK2) was identified and prioritized as therapeutic biomarker amongst LRRK2, PARK7, PINK1 and SNCA from high network associated parameter scores. Structural elucidation and validation of Parkin revealed an active groove (centering T240 residue) with high druggability in the RING0 domain, which plays active role in auto-inhibition of catalytic RING2 domain. Therefore, commercial analogues of conventional UDCA were virtually screened revealing four leads with safe toxicity and favourable pharmacokinetic profiles. The shortlisted leads were further optimized through DFT simulations prior to site-specific molecular docking analysis. Intermolecular interaction profiles and molecular dynamics simulation revealed Molecule-258 (ZINC5022267) as a potent candidate which can be subjected to further experimental validations and can be potentially used against PD.

Funding

The research was funded by the Indian Council of Medical Research, Govt. of India through research grant [IRIS-ID: 2020–0690].

CRediT authorship contribution statement

Aniket Naha: Formal analysis, Investigation, Writing – original draft. Sanjukta Banerjee: Investigation, Data curation, Writing – original draft. Reetika Debroy: Investigation, Data curation. Soumya Basu: Investigation, Data curation, Visualization. Gayathri Ashok: Formal analysis, Data curation, Visualization. P. Priyamvada: Data curation, Visualization. Hithesh Kumar: Validation, Data curation. A.R. Preethi: Validation, Data curation. Harpreet Singh: Funding acquisition. Anand Anbarasu: Supervision, Funding acquisition. Sudha Ramaiah: Supervision, Funding acquisition, Project administration.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.