Pharmacoinformatic and molecular docking studies reveal potential novel antidepressants against neurodegenerative disorders by targeting HSPB8

Abstract

Charcot–Marie–Tooth (CMT) disease is an inherited peripheral neuromuscular disorder characterized by length-dependent and progressive degeneration of peripheral nerves, leading to muscular weakness. Research has shown that mutated HSPB8 may be responsible for depression, neurodegenerative disorders, and improper functioning of peripheral nerves, resulting in neuromuscular disorders like CMT. In the current work, a hybrid approach of virtual screening and molecular docking studies was followed by homology modeling and pharmacophore identification. Detailed screening analyses were carried out by 2-D similarity search against prescribed antidepressant drugs with physicochemical properties. LigandScout was employed to ascertain novel molecules and pharmacophore properties. In this study, we report three novel compounds that showed maximum binding affinity with HSPB8. Docking analysis elucidated that Met37, Ser57, Ser58, Trp60, Thr63, Thr114, Lys115, Asp116, Gly117, Val152, Val154, Leu186, Asp189, Ser190, Gln191, and Glu192 are critical residues for ligand–receptor interactions. Our analyses suggested paroxetine as a potent compound for targeting HSPB8. Selected compounds have more effective energy scores than the selected drug analogs. Additionally, site-directed mutagenesis could be significant for further analysis of the binding pocket. The novel findings based on an in silico approach may be momentous for potent drug design against depression and CMT.

Keywords:

• bioinformatics
• computer-aided drug design
• HSPB8
• heat shock protein (HSP)
• molecular docking
• modeling
• pharmacoinformatics
• virtual screening
• neurodegenerative disorders
• antidepressants

Introduction

The superfamily of ten human small HSPs (sHSPs) has been characterized on the basis of a conserved C-terminal part called the α-crystallin domain.^1,2 This domain is involved in hetero- and homodimer formation, which is observed in many HSPs. Large oligomer formations may also occur, based on the interactions of sHSP dimer amino-terminal regions.^3,4 The formation of oligomers by sHSP has been linked with their chaperoning functional properties.⁵ Some sHSPs are restricted to specific cell-type expression, eg, HSPB9 and HSPB10 are present in male germ cells, while HSPB4 (αA-crystallin) is associated with the eye lens. HSPB5 (αB-crystallin), HSPB1 (HSP27), and to a lesser extent HSPB3, HSPB7 (cvHSP), HSPB2 (MKBP) and HSPB8 (HSP22) are expressed in a variety of cell tissues and types. The latter seven sHSPs are abundant in muscles.^6–10

sHSPs, a large protein family with low molecular weight (12–43 kDa), are ubiquitously expressed in all phyla.¹¹ This protein family contains a 90- to 100-residue-long α-crystallin domain, which is present at the C-terminal part of the molecule.^12,13 sHSPs tend to change in composition and size of oligomers, with flexible structure.^14,15 This significant feature obscures sHSP-structure investigation for the long, stable oligomer structure of Hsp16.5 of Archaea Methanococcus¹⁶ and Hsp16.9¹⁷ of wheat. To date, structures of certain mammalian members of sHSPs have been revealed.^18–20 Ten sHSPs are encoded in the human genome.^1,2 It is supposed that human sHSPs may play a role in the accumulation of harmful protein aggregates by interacting with partially unfolded proteins.^12,13 HSPs maintain the soluble state of unfolded proteins and prevent their precipitation.²¹ This family is present in animals,²² plants,²³ fungi, eukaryotes,²⁴ bacteria,²⁵ and Archaea.¹⁶ Multiple homologues of sHSPs are present in the organisms: up to 15 in plants,²⁶ 13 in zebrafish,²⁷ ten in humans,²⁷ and only a few in yeast and bacteria.²⁴ In higher eukaryotes, sHSPs are tissue specific and distributed subcellularly.^6,28 Mutations in sHSPs are associated with some severe pathologies, such as neurodegenerative diseases, cataract, Charcot–Marie–Tooth (CMT) disease, distal hereditary motor neuropathy (HMN), and desmin-related myopathy.²⁹ Prominent sHSP-specific activity is mediated by binding of other proteins and protecting them from aggregation. This property of sHSPs may inhibit the aggregation of many pathogenic proteins, responsible for Huntington’s disease, dementia with Lewy bodies, Alzheimer’s disease, and Parkinson’s disease. Moreover, filament formation in the mutated and wild-type Parkinson’s forms of α-synuclein can be stopped in the presence of HSPB8.³⁰

Stressful life events are considered a major predisposing risk factor for the development of depression. Depressive individuals usually show profound neuroendocrine alterations (hyper- or hypocortisolism) arising from dysregulation of the hypothalamic–pituitary–adrenocortical axis. The hippocampus has a significant role in stress-related mood disorders. Acute versus chronic stress effects on the integrity of hippocampal circuitry in influencing the vulnerability to or resilience against neuronal injury are poorly understood. Filipović et al concluded that chronic isolation stress induces profound alterations in the functionality of the hypothalamic–pituitary–adrenocortical axis. In contrast to acute stress, chronic isolation strongly induces a reduction in parvalbumin immunoreactivity in γ-aminobutyric acidergic interneurons in all regions of the hippocampus.³¹ The efficacy of antidepressant drugs used in short-term major depression treatment considered as under the control of genetic variations. Numerous pharmacogenetic predictors of antidepressant efficacy have been studied in the last two decades, but reliability is still contradictory. An association between HSP70-coding genes and antidepressant response in mood disorders has been reported, such that genetic variants within the genes coding for HSP70-family proteins may affect the action of antidepressants and thus their therapeutic efficacy.³²

HSPB8 is linked with breast cancer cells, glioblastoma, protection from apoptosis in melanoma, and cell proliferation.^33–35 However, HSPB8 antiapoptotic activity may differ in cancer types. The overexpression of HSPB8 in melanocytes enhances the transformation of tumor during the G₁ cell-cycle stage. Alternatively, HSPB8 opposes Sam68 in glioblastoma cells and facilitates cell-cycle progression from G₁ to S by enhanced expression of cell-cycle regulatory A (PCNA) and cyclin E proteins.³⁴ HSPB8 activates the Akt–PI3K pathway in cardiac cells and supports their survival and growth.³⁶ The level of HSPB8 in some tumor tissues is lower than in normal cells, while restoring normal concentration by cell transfection and DNA demethylation leads to apoptosis induction in p38 MAPK- and caspase-dependent pathways.³⁷ Therefore, an HSPB8-specific role cannot be considered simply anti- or proapoptotic.

The early expression of sHSP orthologues (HSPB1, HSPB2, HSPB5, HSPB6, and HSPB8) during heart development has been reported in a wide range of animal species, including human,^38,39 pig,⁷ mouse,⁴⁰ zebrafish,⁴¹ sea squirt,⁴² and Drosophila melanogaster.⁴³ The expression of HSPB8 is upregulated by ischemic stress.^44,45 It has been inferred in experimental transgenic mice that HSPB8 is involved in cardioprotection via a nitric acid-dependent mechanism in an ischemia context.⁴⁶ HSPB8-specific protective activity is similar to ischemic preconditioning, which involves metabolic adaptation, inhibition of apoptosis, upregulation of HSPs, and activation of cell-cycling pathways.³⁶ Overall, HSPB8 upregulation is correlated with increased expression of the inducible isoform of nitric oxide synthase.⁴⁴ HSPB8 is found in the hypertrophic heart, and its upregulation in the transgenic mouse heart results in cardiac hypertrophy. Knockout HSPB8 mice under stress conditions and pressure overload show various abnormalities, comprising increased mortality, faster transition into heart failure, accumulation of interstitial collagen, increased cardiomyocyte length, impaired contractile function, ventricular dilation, and hypertrophy.⁴⁷

HSPB8 is expressed in the brain, where its expression is not limited but concentrated into the neurogenic niche of the hippocampus (www.brain-map.org; Hspb8-sagittal-b04-0153). The expression of HSPB8 is upregulated in response to hypoxia⁴⁸ and also expressed in cultured hippocampal neurons.⁴⁹ The expression of HSPB8 is in the hippocampus of adult rats and mice, but HSPB8 messenger RNA expression depends on the developmental stage.^49,50 Its expression in adults increases considerably, whereas at postnatal and embryonic stages it is modestly or not at all expressed in the hippocampus of rats.⁴⁹ It is also expressed in blood vessels and neurons.⁵¹ It is also associated with neuronal survival by interacting with Bag3 to induce macroautophagic removal of misfolded proteins.^52,53 HSPB8 mutations are involved in the hereditary peripheral neuropathy of CMT neuropathy type 2.^54,55 The missense mutations occur in HSPB8 crystalline-domain neurodegeneration in cultured motor neurons without inducing cell death.⁵⁵ Moreover, alterations in the crystalline domain lead to CMT neuropathy type 2 and distal HMN. However, the effect of missense mutation is specific only to motor neurons.⁵⁵

CMT is an inherited peripheral nerve disorder that is characterized by length-dependent and progressive degeneration of peripheral nerves, leading to weakness of muscles and wasting in distal hands, feet, and limbs. Neuropathological and neurophysiological defects in sensory and motor nerves create sensory deficits, wheelchair dependence, walking disabilities, and foot deformities. Depending on the severity of sensory and motor deficiency, CMT is grouped with hereditary sensory and autonomic neuropathies and distal hereditary motor neuropathies. Mutations in HSPB8 have been reported to induce CMT2L (OMIM 60867354), neuronopathy, and distal HMN2A (OMIM 15859056).

The present work demonstrates pharmacophore-based virtual screening to reveal novel inhibitors. Pharmacophore-based molecular libraries were screened by 2-D similarity search against recommended US Food and Drug Administration-approved antidepressants (CMT and HMN2A drugs). The inclusive in silico analyses may provide evidence for a reliable framework that could assist medicinal chemists in the design and development of lead molecules and targeting potential antidepressant drugs.

Materials and methods

In the present work, structure prediction, molecular dynamic (MD) simulation, sequence comparison, library screening, pharmacoinformatics, and docking analyses were performed on an HP workstation. The amino acid sequence of HSPB8 (196 residues) was retrieved for homology modeling, as the gene is a suspected candidate of depression, neurodegenerative disorders, HMN2A,⁵⁶ and CMT2L.⁵⁴ The sequence was retrieved in FASTA format from UniProt Knowledgebase with the accession number Q9UJY1. The retrieved sequence of HSPB8 was subjected to a protein–protein Basic Local Alignment Search Tool search against the Protein Data Bank⁵⁷ for the identification of a suitable template. A 3-D structure of the α-crystallin domain of human HSPB1 and HSPB6 (ID 3Q9P) at a resolution of 2.0 Å was selected as a suitable template, with 55% identity, 37% query coverage, and E-value of 2e-22. The automated protein-modeling program Modeller 9v10⁵⁸ was utilized to predict the 3-D structure of HSPB8 by satisfying spatial restraints. The threading approach was also utilized to generate an effective 3-D structure. Geometry optimization followed by energy minimization was performed by utilizing UCSC Chimera 1.6⁵⁹ for a total of 750 steps (step size 0.02 Å), employing a conjugate-gradient method followed by protonation of wild-type histidines using the Amber ff98 method.⁶⁰ The structure was heated to 303 K and simulated for 30 nanoseconds using Gromacs 4.6.5. For input parameters to Gromacs, HSPB8 was provided and MD simulation performed with total steps of 15,000,000 with 0.002 time steps at an average temperature of 303.3 K. The trajectory file was written for every 1,000 MD steps. The ensemble used for MD simulation was NTP (constant number of atoms, temperature, and pressure), with an average grid for the x-axis, y-axis, and z-axis to 8.965 Å. For analysis of MD trajectories, the Visual Molecular Dynamics tool was utilized. Root mean square deviation (RMSD), RMS fluctuation (RMSF), and average RMSD per frame was calculated by loading HSPB8 with a trajectory file. Graphs for RMSD, RMSF, and average RMSD per frame were plotted for better induction of results and analysis of MD simulation. The evaluation tools Anolea,⁶¹ Errat,⁶² Rampage,⁶³ and ProCheck⁶⁴ were applied to assess the quality of the predicted HSPB8 model. The predicted structure was further evaluated by the MolProbity⁶⁵ server. Finally, Ramachandran outliers and poor rotamers were corrected by the WinCoot⁶⁶ tool.

Numerous online servers and tools (ChemDraw,⁶⁷ Mcule,⁶⁸ Vega ZZ,⁶⁹ Osiris Property Explorer,⁷⁰ Molinspiration,⁷¹ Cresset, PyMol, Discovery Studio, AutoDockTools,⁷² and UCSF Chimera 1.6) were utilized to design novel compounds that might potentially inhibit HSPB8 by interacting with a 3-D model. Docking analyses were performed by AutoDock 4.2⁷³ and AutoDock Vina. The hydrogen polar atoms were added to the receptor molecule. The total docking runs were set to 100 for each docking experiment. The grid size was set at 60×60×60 Å in the x-, y, and z-axis, respectively, with 0.575 Å grid spacing to cover the whole receptor. The genetic algorithm implemented in AutoDock 4.2 was employed as the main search protocol, while other parameters were set to default values. The number of H-bond acceptors, rotatable bonds, and H-bond donors were obtained using Molinspiration, Mcule, and PubChem.⁷⁴ The online tool Osiris Property Explorer was utilized to estimate the possible mutagenic risks, reproductive or tumorigenic behavior, and to calculate the drug-like properties of novel designed molecules. Lipinski’s rule of five was analyzed by Mcule and Cresset servers. Drug-likelihood values calculated by Osiris software are positive when fragments of designed molecules are frequently present in the approved drugs. The Osiris program was used to estimate the mutagenesis of novel molecules, and no mutagenic risks were detected. The properties of already reported drugs (fluoxetine, paroxetine, fluvoxamine, and ethacrynic acid) of depressant, neurodegenerative disorders HMN2A and CMT2L were used for library screening and designing novel molecules. The energy minimization and geometry optimization of 3-D structures of designed molecules were performed by ChemDraw Ultra, Vega ZZ, and LigandScout.⁷⁵ The absorption, distribution, metabolism, excretion, and toxicity (ADMET) properties were analyzed by the AdmetSAR server.⁷⁶

Pharmacophoric screening of compounds was performed by employing LigandScout. The three compound libraries (ZINC, Drug, and DrugLike) were screened against HSPB8 structure. The procedure of screening was performed using the LigandScout’s default parameters. The compounds with high pharmacophore score (>70) were extracted, and docking studies were performed on the top hits. The interacting residues were elucidated by employing AutoDock 4.2. The grid box was used to define the screening site. Three pharmacophore-based compounds with optimum binding energies were selected and visualized by LigPlot,⁷⁶ Discovery Studio 4.0, and UCSF Chimera 1.6. Results were analyzed by using Discovery Studio, LigPlot, UCSF Chimera 1.5.2, and AutoDock 4.2.

Results and discussion

The objective of the study under consideration was based on HSPB8’s relation to depression and the neurodegenerative disorders HMN2A and CMT2L, and its computational analysis for designing, identifying, and evaluating novel inhibitors. The top five optimally aligned templates with query coverage, maximum identity, E-values, and total scores are given in Table 1. The 3Q9P template showed better evaluation results. The 3-D model of HSPB8 was generated by homology modeling. Query coverage and identity were not satisfactory for the generation of 3-D structure by comparative modeling. The evaluation tools also validated that the predicted structures had errors and were not satisfactory for further analyses.

Table 1 Five BLAST aligned templates of HSPB8 with identity, query coverage, and E-values

PDB ID	Query coverage	E-value	Maximum identity
3Q9P	37%	2e-22	55%
2KLR	84%	2e-18	34%
2WJ7	38%	1e-17	44%
3L1G	38%	1e-17	45%
2Y1Y	38%	7e-17	44%

Abbreviations: BLAST, Basic Local Alignment Search Tool; PDB, Protein Data Bank.

The threading approach and multiple-template comparative modeling approach were utilized for 3-D structure prediction of HSPB8. To build the 3-D structures, various web servers (3D-Jigsaw, Swiss Model, M4T, and I-TASSER) and Modeller 9.10 were employed. As a result, various evaluation tools were utilized for the comparison of predicted structures and suitable optimal predicted model of HSPB8 was selected. The HSPB8 3-D structure was comparatively modeled using the threading approach.

The predicted model of HSPB8 was subjected to MD simulation by employing appropriate solvent molecules, ensemble, and temperature. Energy minimization was performed on the scrutinized structure to remove steric constraints by relaxing the system. The gradual velocities of the system were gradually increased through heating to the 303.3 K temperature. In 30 nanoseconds’ simulation of protein, it was observed that the protein changed its structure dramatically. The changes in structure were analyzed through RMSD, RMSF, and average RMSD per frame of protein. According to the RMSD graph, the protein seemed to attain higher RMSD at the start of the simulation.

RMSD increased to 14 Å in the first 8 nanoseconds, while stabilization of the protein was observed afterward till the end of the simulation. The line of linearity fell at 9.6476 Å, which depicted high RMSD attainment during simulation (Figure 1). The overall RMSD change according to the Pearson product-moment coefficient was 0.5242 Å, demonstrating protein stabilization after abrupt change in structure.

Figure 1 Root mean square deviation (RMSD) versus time.

Notes: HSPB8 showed an increasing trend in RMSD in the first 2 ns. From 2 to 10 nanoseconds, RMSD value continuously increased, while it stabilized for the next 20 nanoseconds.

The RMSF graph demonstrated per-residue fluctuations throughout simulation. The observed RMSF values showed that residue fluctuations were less while fluctuation observed in residues from 1–17 and 167–196 that can also seen in protein structure at last nanoseconds (Figure 2). Overall protein residues did not show high fluctuation, fluctuations lasted only up to 8 nanoseconds, and afterwards the protein was stabilized with minimum fluctuation, as can be observed in the RMSD graph.

Figure 2 Root mean square fluctuation (RMSF), showing residue fluctuations.

Note: Starting and ending residues showed higher fluctuations during molecular dynamic simulation.

Average RMSD per frame demonstrated per-residue fluctuation, including side chains, to give an overall RMSD of a frame. Average RMSD is different from RMSD in the sense that RMSD provides information for fluctuations of only backbone atoms from reference coordinates of protein, while average RMSD calculates changes of coordinates of protein from the reference structure, including complete residues with side chains. While average RMSD per frame values were being compared with RMSD and RMSF values, it could be analyzed that the protein showed higher fluctuation only at the start of the simulation, while fluctuation decreased toward the end of the simulation. Based on average RMSD, starting frames of simulation showed much higher RMSD values, while afterward stabilization of the protein was observed (Figure 3). The simulated minimized 3-D structure of HSPB8 was visualized by UCSF Chimera 1.5.2 (Figure 4).

Figure 3 Average root mean square deviation (RMSD) values per frame.

Note: Average RMSD seemed to be higher at the start, then fell, depicting stabilization of the structure.

Figure 4 Predicted HSPB8 structure simulated at 30 nanoseconds.

The utilized evaluation tools indicated the reliability and efficacy of predicted HSPB8 structure. The Ramachandran plot showed the presence of 93.3% residues in favored region, and only 15 of 196 residues were observed in the outlier region. Errat showed an overall quality factor of 77.444 for structure.

Subsequently, poor rotamers and outliers were corrected to refine the HSPB8 predicted model. Experimental analysis elucidated that selected drugs used in this study (Figure 5) exhibited significant abilities (Table 2) to cure depression, neurodegenerative disorders, and CMT. However, our docking studies of selected drugs exposed variations in their binding energies. Initially, docking studies were performed with 100 runs and 15 poses saved, of which the best pose with the lowest binding energy was selected for each compound. Our results indicated that the four selected compounds (fluoxetine, paroxetine, fluvoxamine, and ethacrynic acid) efficiently bound to HSPB8 (Table 2).

Table 2 Compounds investigated in this study

Ligand properties	Fluoxetine	Paroxetine	Fluvoxamine	Ethacrynic acid
IUPAC name	N-methyl-3-phenyl-3-(4-(trifluoromethyl)phenoxy)propan-1-amine	(3S,4R)-3-[(2H-1,3-benzodioxol-5-yloxy)methyl]-4-(4-fluorophenyl)piperidine	2-{[(E)-{5-methoxy-1-[4-(trifluoromethyl)phenyl]pentylidene}amino]oxy}ethanamine	[2,3-dichloro-4-(2-methylenebutanoyl)phenoxy]acetic acid
Estimated free energy of binding, kcal/mol (AutoDock Vina)	−7.1	−7.8	−6.4	−7.2
Estimated free energy of binding, kcal/mol (AutoDock 4)	−7.37	−7.5	−6.37	−6.16
Estimated inhibition constant, Ki (µM)	3.97	3.2	11.6	30.34
Final intermolecular energy, kcal/mol	−8.88	−8.11	−8.73	−7.11
Ligand efficiency	−0.34	−0.31	−0.31	−0.32
Torsional free energy, kcal/mol	1.92	1.10	3.29	1.92
Molecular weight	309.32	329.36	318.33	303.13
Hydrogen bond acceptor	5	5	7	4
Hydrogen bond donor	1	1	1	1
Rotatable bonds	7	4	10	6
Rule of five (violation)	0	0	0	0
ClogP	3.96	3.82	3.01	3.22
Solubility	−3.56	−4.49	−3.7	−4.17
Blood–brain barrier (probability)	BBB^+ (0.9830)	BBB^+ (0.9869)	BBB^+ (0.9775)	BBB^− (0.6921)
Human intestinal absorption (probability)	HIA^+ (1)	HIA^+ (1)	HIA^+ (1)	HIA^+ (0.9591)
Caco2 permeability (probability)	Caco2^+ (0.8867)	Caco2^+ (0.5195)	Caco2^− (0.5149)	Caco2^+ (0.5597)
CYP450 2D6 inhibitor (probability)	Inhibitor (0.8932)	Noninhibitor (0.9231)	Inhibitor (0.8932)	Noninhibitor (0.9231)
Carcinogens (probability)	Noncarcinogen (0.8089)	Noncarcinogen (0.9046)	Noncarcinogen (0.6370)	Noncarcinogen (0.8344)
Acute oral toxicity (probability)	Category III (0.7805)	Category III (0.5546)	Category III (0.5789)	Category III (0.7998)
Aqueous solubility (logS)	−3.9179	−3.3848	−3.0217	−3.3381
Rat acute toxicity (LD50, mol/kg)	2.6048	2.8239	2.6997	2.4505
Polar surface area	26.14	–	56.84	63.60
Solvent accessibility surface area (Å^2)	464.067	496.845	490.579	472.769
Drug likeness	0.05	2.16	−11.4	−2.84
Drug score	57%	64%	24%	19%
Binding residues (AutoDock 4)	Ser57, Ser58, Trp60, Pro61, Gly62, Thr63, Thr114, Lys115, Val152, Thr153, Phe155, Leu186, Asp189, Ser190, Gln191, Glu192	Met37, Ser57, Ser58, Trp60, Thr63, Thr114, Lys115, Asp116, Gly117, Val152, Val154, Leu186, Asp189, Ser190, Gln191, Glu192	Ser57, Ser58, Trp60, Pro61, Gly62, Thr63, Thr114, Lys115, Val152, Leu186, Ser190, Gln191, Glu192	Trp60, Pro61, Gly62, Thr114, Lys115, Asp116, Gly117, Val152, Thr153, Val154, Phe155, Leu186, Ser190, Gln191

Abbreviations: IUPAC, International Union of Pure and Applied Chemistry; LD₅₀, median lethal dose.

Figure 5 2-D structures of investigated drugs.

Notes: (A) Fluoxetine; (B) fluvoxamine; (C) paroxetine; (D) ethacrynic acid.

Subsequently, four selected drugs were analyzed on the basis of binding energy values, and ADMET and drug properties (Table 2) indicated that the selected drugs exhibited significant biological properties. Therefore, these drugs may be considered potential agents against depression, neurodegenerative disorders, and CMT. The lowest binding energy value (−7.8 kcal/mol by AutoDock 4 with Vina) was observed for paroxetine. These four drugs were used to generate the pharmacophore.

Libraries were screened by LigandScout, and the 21 compounds with the lowest pharmacophore scores (>70) were selected for further analyses. AutoDock 4 and AutoDock Vina were employed for docking analyses, and the three complexes with the lowest binding energy were chosen (Figure 6). The current work gave reliable results and conclusions by docked analyses of selected compounds against HSPB8. In an effort to investigate the novel hits, three scrutinized molecules from a combination of selected libraries were revealed (Table 3). It was found that the majority of compounds were bound at the binding residues between Ser57 and Ser190. The docking analyses revealed that known (Figure 7) and novel molecules bound at the same types of residues, sharing a quite similar binding pocket. Another observation was that most molecules occupied the space between Ser57 and Gly117.

Table 3 Bioinformative details and ligand properties of top four screened compounds

Ligand properties	PB-765894052	PB-411001374	PB-411001436
Estimated free energy of binding, kcal/mol (AutoDock Vina)	−8.5	−7.1	−7.1
Estimated free energy of binding, kcal/mol (AutoDock 4)	−4.05	−6.3	−6.3
Estimated inhibition constant, Ki (µM)	1.07	24.25	38.08
Final intermolecular energy, kcal/mol	−6.01	−8.09	−7.14
Ligand efficiency	−0.12	−0.25	−0.24
Torsional free energy, kcal/mol	2.2	1.65	1.65
Molecular weight	505.71	373.56	372
LogP	2.49	1.87	0.25
Hydrogen bond acceptor	9	8	10
Hydrogen bond donor	5	6	7
Rotatable bonds	9	8	9
Rule of five (violation)	0	0	0
ClogP	−5.2	−0.95	−2.57
Solubility	−2.53	−1.87	−1.6
Blood–brain barrier (probability)	BBB^− (0.5624)	BBB^+ (0.6250)	BBB^+ (0.6707)
Human intestinal absorption (probability)	HIA^+ (0.6418)	HIA^+ (0.8951)	HIA^+ (0.9045)
Caco2 permeability (probability)	Caco2^− (0.6904)	Caco2^− (0.5959)	Caco2^− (0.6007)
CYP450 2D6 inhibitor (probability)	Noninhibitor (0.87)	Noninhibitor (0.83)	Noninhibitor (0.8532)
Carcinogens (probability)	Noncarcinogen (0.7482)	Noncarcinogen (0.7687)	Noncarcinogen (0.6637)
Acute oral toxicity (probability)	Category III (0.5645)	Category III (0.6150)	Category III (0.6188)
Aqueous solubility (logS)	−3.3358	−2.6424	−2.6113
Rat acute toxicity (LD50, mol/kg)	2.5859	2.5595	2.5258
Polar surface area	134.19	112.16	92.85
Solvent accessibility surface area (Å^2)	718.053	574.367	546.72
Drug likeness	1.44	5.78	3.55
Drug score	65%	89%	88%
Binding residues (AutoDock 4)	Phe21, Arg29, Leu30, Asp33, Arg55, Leu56, Ser57, Ala59, Pro70, Arg71, Ala75, Thr76, Phe79, Lys113, Glu192	Ser57, Ser58, Gly62, Thr63, Lys113, Thr114, Lys115, Asp116, Glu120, Leu186, Asp189, Ser190, Gln191, Glu192	Phe8, Arg15, Leu16, Trp60, Phe79, Val81, Phe155, Glu167, Pro169

Abbreviation: LD₅₀, median lethal dose.

Figure 6 2-D structures of novel molecules.

Notes: (A) PB-765894052; (B) PB-411001374; (C) PB-411001436.

Figure 7 The potential binding interactions of investigated drugs through LigPlot.

Notes: (A) Ethacrynic acid; (B) fluoxetine; (C) fluvoxamine; (D) paroxetine.

The highest-ranked hits with the lowest binding energy values from all screened libraries were elucidated, namely PB-765894052. It was observed that Met37, Ser57, Ser58, Trp60, Thr63, Thr114, Lys115, Asp116, Gly117, Val152, Val154, Leu186, Asp189, Ser190, Gln191, and Glu192 of HSPB8 exhibited high binding affinity with all the docked ligands. In an attempt to evaluate the occurrence of detailed interactions among HSPB8-active sites and analyzed ligands, plots of ligand and residue interactions were generated by utilizing UCSF Chimera 1.5.2 (Figure 8).

Figure 8 Binding pocket and interacting residues of novel molecules.

Notes: (A) PB-765894052, (B) PB-411001374, (C) PB-411001436.

The HSP family members share significant features: a highly dynamic structure that plays a chaperone-like role, increased expression under stressful conditions, into large oligomers consisting in many cases of subunits, monomeric molecular mass range (12–43 kDa), and formation of the conserved α-crystallin domain structure with about 80–100 residues in the C-terminus.⁷⁷ HSPB8 is highly demanding, due to its significant role in pathological and physiological processes. The mutations link with the development of CMT and distal motor neuropathy type II.^55,78 HSPB8 has also been detected in hereditary cerebral hemorrhage with amyloidosis of the Dutch type and in cerebral amyloid angiopathy and senile plaques in patients with Alzheimer’s disease.^50,79,80 HSPB8 has a conserved α-crystallin domain.

The α-crystallin domain extends from residues 86 to 176, while the N-terminus extends from residues 1 to 85. The HSPB8 has 21.6 kDa mass.^9,77,81 Fluvoxamine is a widely used drug for depression treatment. The analgesic effects of fluvoxamine have been documented for diabetic neuropathic pain in rats, and are limited due to its administration as a single dose.⁸² Fluvoxamine has effectively used for ameliorates cardiac hypertrophy, σ₁-receptor and dysfunction via σ₁-receptor stimulation. Patients in old age usually receive statin drugs for the secondary and primary prevention of cerebrovascular and cardiovascular events.

Patients commonly receive selective serotonin-reuptake inhibitors (SSRIs) as antidepressant drugs for the treatment of anxiety, depression, and other conditions. SSRIs are linked with pharmacokinetics through drug interactions linked to cytochrome CYP450 metabolic pathway inhibition. SSRIs’ specific effect and statin metabolism on CYP450 enzymes suggest that paroxetine, citalopram, and escitalopram are safe with statins. Fluoxetine is a commonly used medication for the treatment of depression. Fluvoxamine used with the combination of lovastatin, simvastatin, or atorvastatin, and the risk can be minimized by monitoring patients and utilizing lower doses of statin.⁸² Recent research has suggested a possible effect of fluoxetine on bone. This drug may improve the quantity of bone regeneration in calvarial defects of rats.⁸³ Gene-expression changes of β₁- and β₂-adrenoceptors were studied in the left and right ventricles and atria after treatment with fluoxetine in stress-induced depression in adult male rats. Fluoxetine treatment led to increased expression of β₁-adrenoceptor messenger RNA, which suggests that therapy with fluoxetine plays a significant role in animals with heightened sympathetic nervous activity.⁸⁴ It was investigated whether paroxetine contained adverse effects in comparison with other antidepressant agents, and no clear evidence has yet been reported that paroxetine is worse or better in comparison to antidepressant agents give response to increase the treatment.⁸⁵ Patients with depression, neurodegenerative disorders, or CMT disease lose quality of life and feel mentally depressed. Mentally relaxed patients recover from CMT earlier compared to depressed patients. The drugs analyzed are used as antidepressant agents, and the identified novel hits may be considered better antidepressants on the basis of overall analysis by targeting HSPB8.

Bioinformatics utilizes in silico analyses of biological problems by using mathematical and statistical approaches,⁸⁶ and structural bioinformatics predicts and analyzes 3-D structures of macromolecules. Structural bioinformatics has solved numerous biological problems and led to novel computer-aided designed molecules against neurological disorders^87,88 and cancer.^89–91 Computer-aided drug design helps scientists to design effective active molecules against drugs. Drug design and development is not only a costly procedure but also time-consuming. Therefore, we applied computational approaches and methodologies for pharmacophore generation in a drug-discovery procedure. This emerging tendency has colossal significance in reducing the time phase, by amplifying the design of pharmacophore molecules with minimal side effects for specific targets of disease and better biological activity.

Chemical structures are assessed for better oral bioavailability and to be efficacious drug-like compounds subjected to Lipinski’s rule of five.⁹² According to Lipinski’s rule of five, molecules should have not violated more than one of the following conditions: logP no higher than 5, molecular weight no more than 500, no more than ten hydrogen-bond acceptors, and no more than five hydrogen-bond donors.⁹³ ADMET properties of molecules within an organism were predicted for selected compounds by utilizing the AdmetSAR online server. Eight mathematical models (blood–brain barrier penetration, aqueous solubility [logS], carcinogens, acute oral toxicity, Caco2 permeability, cytochrome P450 2D6 inhibition, rat acute toxicity, and human intestinal absorption) were utilized to predict the ADMET characteristics of molecular chemical structures. A variety of toxicities were predicted that are often used in drug design (Table 3). These toxicities help in evaluating pollutants, metabolites, and intermediates, along with adjusting the dose range for animal assays.

Polar surface area (PSA) has an inverse relationship with human intestinal absorption. Therefore, cell-wall permeability through a relationship between permeability and PSA has been revealed. PSA values for massive molecular transport through the membrane showed lower PSA at a standard limit of 140 Ų. The prediction of aqueous solubility (defined water at 25°C) showed that selected compounds are soluble in water. LogP is the lipophilicity measurement, which is the solubility ratio of a compound in octanol compared to its solubility in water. It was observed that novel compounds showed lower logP values compared to known drugs, and followed Lipinski’s rule of five, showing better oral bioavailability. The process of excretion that purges the compounds from the human body also depends on logP.⁹² Orally administered drugs and novel compounds must be absorbed by the human intestine. The analyzed human intestinal absorption results showed that all the selected drugs and novel compounds can be easily absorbed by the intestine. From the selected four compounds, paroxetine was found to be a noninhibitor of cytochrome P450 2D6, which indicated that it may be well metabolized in phase I metabolism. It was also observed that all the identified novel compounds showed noninhibitor properties of cytochrome P450 2D6 and could be metabolized in phase I metabolism. Cytochrome P450 2D6 is one of the significant enzymes involved in drug metabolism.⁹² Toxicity-risk assessment and carcinogenicity was performed for the selected drugs and novel compounds, and showed that all the compounds were noncarcinogenic.

In the current study, in silico analyses were carried out, and the 3-D structure of HSPB8 was modeled and simulated at 30 nanoseconds. The simulated model had a good degree of accuracy, especially at the active site of the protein. Comparative molecular docking studies were analyzed by AutoDock 4 and AutoDock Vina, which revealed the ligand–receptor interactions of the most representative drugs. The selected drugs showed binding in the binding domain and critical binding residues (Met37, Ser57, Ser58, Trp60, Thr63, Thr114, Lys115, Asp116, Gly117, Val152, Val154, Leu186, Asp189, Ser190, Gln191, Glu192), observed in molecular docking studies of selected drugs with HSPB8. It was found that these residues also interacted with novel compounds (Figure 9), revealed by the binding region. The novel compounds also showed effective drug likeness and drug score compared with selected drugs. The molecular docking studies suggested that effective ligand molecule must have the least binding energy and also it has to satisfy the drug-likeness parameters. Given this, and the satisfaction of the parameters of drug likeness, ADMET properties, drug score, and lowest docking energy, it is suggested that paroxetine and novel molecules are potential drug molecules. It stands to reason that the discovered compounds have the propensity to be good candidate molecules.

Figure 9 All the analyzed ligands bound at the same binding pocket in HSPB8.

In conclusion, the present study suggested that novel compounds and selected drugs are effective in depression, neurodegenerative disorder, and CMT treatment. Though various divergences exist among trail-methodology studies, baseline population- and bioinformatics-based pharma-coanalysis seem to be sufficient to conclude that paroxetine may be the better option for treatment. Further studies and synthesis of novel compounds considering these findings may expect similar response rates.

Acknowledgments

The work was sponsored by a CAS-TWAS President’s Fellowship for International PhD Students. The authors are thankful to Mirza A Hammad, and Rana Adnan Tahir, Bioinformatics Research Laboratory, Department of Bioscience, COMSATS Institute of Information Technology, Sahiwal, Pakistan, for their kind guidance in interpreting and suggestions for improving the results. The project was funded by a grant from the National Natural Sciences Foundation of China (NSFC31370028 to Lei Du).

Disclosure

The authors report no conflicts of interest in this work.