In Silico Repositioning of Dopamine Modulators with Possible Application to Schizophrenia: Pharmacophore Mapping, Molecular Docking and Molecular Dynamics Analysis.

We have performed theoretical calculations with 70 drugs that have been considered in 231 clinical trials as possible candidates to repurpose drugs for schizophrenia based on their interactions with the dopaminergic system. A hypothesis of shared pharmacophore features was formulated to support our calculations. To do so, we have used the crystal structure of the D2-like dopamine receptor in complex with risperidone, eticlopride, and nemonapride. Linagliptin, citalopram, flunarizine, sildenafil, minocycline, and duloxetine were the drugs that best fit with our model. Molecular docking calculations, molecular dynamics outcomes, blood-brain barrier penetration, and human intestinal absorption were studied and compared with the results. From the six drugs selected in the shared pharmacophore features input, flunarizine showed the best docking score with D2, D3, and D4 dopamine receptors and had high stability during molecular dynamics simulations. Flunarizine is a frequently used medication to treat migraines and vertigo. However, its antipsychotic properties have been previously hypothesized, particularly because of its possible ability to block the D2 dopamine receptors.

Introduction

Drug repurposing (also known as drug repositioning, rediscovering, rescuing, or redirecting) investigates whether it is possible to develop new therapeutic applications for drugs currently in use.[1] It reduces time, production costs, and stages of development of new drugs, all of which usually take over 10–15 years and represent an investment of around several hundred million dollars.[2] Interestingly, psychopharmacology is one of the medical areas where drug repositioning has more potential to be applied. New therapeutic applications of some drugs and formulations currently in use were observed during their clinical administration to treat other diseases. This situation has occurred relatively frequently with medicaments used in neuropsychiatry, strengthening the approach to drug reuse in this field of medicine.[3]

Over the last years, schizophrenia patients have been markedly characterized,[4] and the analyses have been centered on one of the meaningful approaches, a sharp increase in presynaptic dopamine function. Therefore, targeting dopaminergic neurotransmission is continuously used for treating and approaching schizophrenia, aside from other mental disorders and psychiatric illnesses such as Tourette’s syndrome, Parkinsonian disorders, Huntington’s disease, attention deficit hyperactivity disorder (ADHD), and bipolar disorder.[5] The physiological cues of dopamine are intervened by five intimately related yet functionally distinct G protein-coupled receptors (GPCRs) that are sorted based on their sequence and pharmacological similitudes in two main subfamilies: the D1-class receptor (D1 and D5) and the D2-class receptor (D2, D3, and D4) receptors.[6] When dopamine (DA) is liberated from presynaptic axonal terminals, it might interact with D1-like receptors (D1Rs), which are positively linked to the production of cAMP and the stimulation of adenylyl cyclases, or with D2-like receptors (D2Rs), whose activation inhibits adenylyl cyclases and interrupt the generation of cAMP. Additionally, Ca2+ levels and several intracellular signaling processes are notably modulated by DRs. Neuronal activity and synaptic plasticity, for instance, are outstandingly affected by DA and in large part by its dopamine receptor (DR) binding.

The dopamine receptor dysfunction is, ostensibly, associated with human disorders in a vast number of multiple genetic studies. Variants in the D2, D3, and D4 dopamine receptor (DR) genes have been linked to schizophrenia or involved in the reaction to antipsychotic medications. These divergent variants have undeniably led to a more complex understanding of psychopharmacology. Haloperidol, for example, is relatively selective for dopamine D2, D3, and D4 receptors, whereas the pharmacology of olanzapine looks intricate.[6−8]

Schizophrenia is a neurodevelopmental disorder with pronounced disorders of thought, emotional disbalances, hallucinations, and mood changes. The annual prevalence is around 1.4 to 4.6 per 1000 inhabitants with an incidence between 0.16 and 0.42 cases per 1000 inhabitants.[9] The most prominent antipsychotic drugs used over the long term have been clozapine, aripiprazole, or olanzapine even though some of these drugs have little to null specificity with multiple side effects, and some could even contribute to a dramatic decline in the production of white blood cells. In addition, orthostatic hypotension, weight gain, and other extrapyramidal symptoms are commonly diagnosed in patients prescribed such medication.[10] Therefore, to address the gaps in the treatment, it is necessary to implement new drugs to minimize the known side effects and also to try to promote drug adherence. Evidence from systematic analysis has stated that patients with schizophrenia, generally, are over-sensitive to dopamine psychostimulants due to the considerable proportion of active D2-like dopamine receptors. As a result of this circumstance, psychotic reactions have been studied as an integrated state of D2 (high) dopamine supersensitivity.[11,12]

In some former reports, certain advancements have emerged in computational studies applied to drug repositioning in schizophrenia. Given the great value and the vast spectrum of the in silico tools used to understand drug design, it is not a coincidence that many new genome-wide associations,[13] machine learning approaches,[14] and interactome relations[15] are becoming so widespread in the later years. The central objective of the present study does open the door to far analysis of information, which is already conserved in one of the largest publicly available clinical trials databases (ClinicalTrials.gov, 18 January 2018 data freeze). The aim is to show how the clinical trials database information could be used and studied with cheminformatics techniques. In a recent revision of the abovementioned database, emerging data was collected and reviewed for about 88 repurposing candidates, stated across 231 clinical trials for their therapeutic actions in schizophrenia.[1] These compounds were tested and reported in terms of their pharmacological effects. They were included with at least one primary intervention using a regulatory-approved drug that is not already designated for schizophrenia but assessed in a target of schizophrenia. They were designed in categories depending on their therapeutic use in schizophrenia (nootropic agents 25%, hormonal agents 17%, monoamine neurotransmitters agents 15%, compounds that suppress inflammation 14%, compounds with anti-convulsant/mood disorders/sedative action 12%, metabolic/cardiovascular 11%, infectious disease 3%, and others 3%). The present study aims to perform computational calculations with 70 of these drugs (chosen for not being proteins or antibodies) interacting with D2-like receptors regarding the hypothesis of the dopaminergic system in schizophrenia.

The repurposing candidates were retrieved and then analyzed with different robust computational methodologies. To test and predict the conformations of the ligands when interacting with receptors, docking programs can mimic a similar structure to the one observed in crystallographic protein/ligand complexes. It was assumed that some limitations would arise from the attempt to rank these conformations as the top docking solutions in conjunction with the fact that it would promote an incomplete understanding of the pharmacology of dopamine receptors and human diseases. To tackle these issues, virtual screenings based on the best pharmacophore matching and molecular dynamics simulation were the leading centers of the present study.[16−19] A set of crystallographic ligands already bonded to DRs were inspected, and those certain frequent 3D shared pharmacophore features were grouped and then tested on the molecules identified as repurposing dopamine modulators.[20]

Selectively targeting DRs has been widely researched by computational methods such as molecular dynamics, pharmacophore model, molecular docking, and machine learning.[21−23] Molecular docking is guided by the target, where each compound is coupled and assessed against the characterized target structure (e.g., protein, DNA, RNA), whereas a pharmacophore model could be assembled whether or not in the presence of the biomolecular structure and by selecting the key points of the interaction, such as chemical, steric, and electronic features required on the compound. Each of these techniques has some attributes and scopes that make them profitable and correlative resources to search for molecules that inhibit or modulate different molecular targets in a disease.[24] One of the most highlighted merits of molecular dynamics is the possibility of examination of thermodynamic properties and other molecules that are determining on the free protein binding sites, as well as it would explain the evolution through the time of the system subjected to forces between atoms and molecules for a fixed period, to characterize at the atomic level, the ligand–target conformational and energetic landscapes.

Pharmacophore models coupled to molecular docking and molecular dynamics have been reported in many drug repositioning and targeting studies.[25−27] Kagami et al.[28] identified an inhibitor of a promising target of the human infection malaria, the purine nucleoside phosphorylase of Plasmodium falciparum, using these techniques. Two inhibitors were proposed and supported using molecular docking, shared pharmacophore features, and molecular dynamics simulation. In another study by Kumar et al.,[29] the protein–ligand interactions of different complexes were applied to Alzheimer’s disease, and in silico repositioning of antipsychotic drugs was pointed out. Dilly et al.[30] conducted in silico studies of naproxen and its analogs to demonstrate that they could have antiviral activity against the influenza A virus. Teli and Rajanikant[31] performed in silico repositioning using drugs approved by the FDA as inhibitors of PHDs, which are compounds of the family of dioxygenases called HIF prolyl-4-hydroxylases. Other examples include computational methods, which have been increasing as the first strategies on the drug design process in pharmaceutical companies, making them almost mandatory during the procedure.[32−35]

Blood–brain barrier (BBB) penetration and gastrointestinal (GI) absorption are prime factors in neurological drugs, and these two properties are related to lipophilicity (WLOGP) and polarity (tPSA) molecular structure computation.[36] Since the in vitro and in vivo BBB tests are difficult, time-consuming, and low cost, both physicochemical features (WLOGP and tPSA) were calculated as predictive models to mimic the behavior of the drug within the human body, especially to determine a prediction of BBB penetration and possible side effects with the metabolism of the compounds tested. The following physicochemical-related structural properties are canonical and must be met by any new drug: (i) lipophilia; XLOGP3 [−0.7, +5.0], (ii) molecular weight [150-500] (g/mol), (iii) topological polar surface area (TPSA) [20, 130] Å2, (iv) solubility (log S) < 6, (v) carbons in sp3 hybridization >0.25, and (vi) flexibility <9 rotatable bonds.[104] The compounds that fall within the previously mentioned ranges will probably have good oral bioavailability, i.e., because of their structural configuration it will be easier for them to reach their therapeutic target (channels, transporters, receptors, or other protein macromolecules) in their active form.

Results and Discussion

The increase in new cases of schizophrenia and its strong complications have shown why it is classified among one of the leading causes of long-term disability.[37] Typically it emerges in late adolescence and late adulthood, which is why, although huge efforts have been overcome to prolong the life expectancy of people with schizophrenia, currently there is no effective way to prevent the cognitive impairments, and the causes are still unknown. Amid advancements in applied schizophrenia research, yet it is not well understood which should be the main focus, for example, environmental influences (e.g., obstetric complications) or genetic factors along with social factors (e.g., poverty) are strong connections and contributing factors.[38,39]

The pathogenic mechanisms underlying schizophrenia are unknown as well. One of the key factors that have been exhibited during the past years is that schizophrenia subjects have a notably higher amount of D2-like receptors compared to the control people.[40] Dopaminergic synapse and neuroactive ligand–receptor are some of the most representative and admitted pathways characterized in the illness, even though other pathways are being studied.[41,42]Figure shows a brief overview of dopaminergic synapse,[43] in red are highlighted D2-like dopamine receptors. Dopamine receptors, also called G protein-coupled receptors or rhodopsin-like receptors, are distinctively characterized by having seven-transmembrane regions interconnected through three extracellular and intracellular loops. All DRs traverse the membrane seven times, and they are recognized as serpentine receptors because of how they wind back and forth across the membrane.[44] Due to the fact that DRs are considerably intervened by the activation of G proteins, most of their effects are acknowledged as G protein-coupled. In the context of drug therapy, drugs that block D2-like receptors are substantially recommended for the treatment of schizophrenia and other psychotic disorders; meanwhile, drugs that stimulate D1-like or D2-like receptors help with the motor symptoms that result from the decay of dopamine-related neurons of some diseases like Parkinson.[45]

Figure 1

Brief overview of the dopaminergic synapse pathway [KEEG copyright permission obtained: Material: Dopaminergic synapse, pathway map04728. (2021/ KEGG: Kyoto Encyclopedia of Genes and Genomes database resource/Bioinformatics Center, Institute for Chemical Research, Kyoto University and Human Genome Center, Institute of Medical Science, University of Tokyo).[105−107]

To analyze different binding modes of drugs reported as modulators of D2-like receptors and co-crystallized in with D2, D3, and D4 dopamine receptors, a 3D pharmacophore hypothesis was used to analyze the shared features in the binding site of the crystallographic ligands risperidone, eticlopride, and nemonapride bounded to D2, D3, and D4 protein receptors, respectively.

One widely used antipsychotic is risperidone. It is considered a second-generation antipsychotic and is one of the main drugs administrated during the treatment of schizophrenia. Risperidone has been recognized as one of the leader drugs due to its properties of reducing the overactivity of D2R receptors in the brain.[46,47] Eticlopride is sorted as a dopamine agent and neurotransmitter agent, and it is a substituted benzamide with a powerful affinity for D2-like but not D1-like receptors.[48] Nemonapride is another antipsychotic drug approved in Japan for schizophrenia, and it belongs to the benzamide class as well. It has a potent affinity to D2, D3, and D4 receptors.[49]

Figure shows the pharmacophores established between the drugs risperidone, eticlopride, and nemonapride with D2, D3, and D4 receptors, respectively. Risperidone and the D2 dopamine receptor interact with three hydrophobic points: two owing to the 6-fluorobenzo[d]isoxazole region where the fluor is interacting with ALA89A, OLA6, PHE370A, and PHE165A and the phenyl group interacts with PHE378A, ALA89A, PHE165A, and TRP37A, while the third hydrophobic interaction occurs with the methyl group of the 4H-pyrido[1,2-a]pyrimidin-4-one. Additionally, risperidone presented one positive ionizable area located in the piperidine. Eticlopride interacts with the D3 receptor with four hydrophobic areas and one positive ionizable region. The hydrophobic interactions were in the 4-chloro-2-ethyl-5-methoxyphenol and the ethyl group bounded to pyrrolidine. The involved aminoacids in the hydrophobic interactions of eticlopride were VAL111A, PHE346A, ILE183A, VAL350A, VAL189A, PHE345A, TRP342A, TYR373A, and THR369A. The pyrrolidine ring of eticlopride interacts with ASP110A through one positive area due to the formation of a salt bridge between them. This is a result of the high likelihood of the pyrrolidine ring being charged at physiological pH.[50] Nemonapride interacts with the D4 dopamine receptor with three hydrophobic regions. The hydrophobic zones were located in the 1-benzyl, the 2-methyl bounded to pyrrolidine, and the 2-chloro-5-methoxy-N-methylaniline. The involved aminoacids in the hydrophobic interactions of nemonapride were LEU111A, PHE91A, VAL87A, TYR390A, THR386A, PHE362A, VAL116A, VAL193A, and LEU187A. Nemonapride presented two types of interaction with Asp115A,: a positive ionizable area with the pyrrolidine ring and a H-bond interaction with the formamide group. It is worth noting that all the presented drugs had at least three hydrophobic interactions and one positive ionizable feature into the established pharmacophore.

Figure 2

Established pharmacophore between (A) risperidone and D2 dopamine receptor (PDB: 6CM4-8NU), (B) eticlopride and D3 dopamine receptor (PDB: 3PBL-ETQ), and (C) nemonapride and D4 dopamine receptor (PDB: 5WIU-AQD).

The pharmacophore feature model was considered the first strategy to search into the database of 70 possible repurposing drugs (Table ). The discovery and deeper understanding of such interactions presented a crucial starting point for this study.

Table 1

Information of Prospective Repurposing Drugs Selected for the Present Study

drug information
drug	mechanism of action	disease area	indication
allopurinol	xanthine dehydrogenase inhibitor	rheumatology	gout
citalopram	serotonin transporter inhibitor	neurology/psychiatry	depression
clonidine	adrenergic receptor alpha-2 agonist	cardiology	hypertension
cycloserine	alanine racemase inhibitor	infectious disease	tuberculosis
cysteamine	cystine hydrolytic enzyme	metabolism	cystinosis
dextroamphetamine	dopamine transporter releasing agent	neurology/psychiatry	ADHD
dextromethorphan	glutamate (NMDA) receptor subunit epsilon 1 antagonist	pulmonary	cough suppressant
dipyridamole	3′,5′-cyclic phosphodiesterase inhibitor	cardiology	coronary artery disease (CAD)
donepezil	acetylcholinesterase inhibitor	neurology/psychiatry	Alzheimer’s disease
duloxetine	serotonin transporter inhibitor	neurology/psychiatry	depression
esomeprazole	potassium-transporting atpase inhibitor	gastroenterology	gastroesophageal reflux disease
estradiol	estrogen receptor alpha agonist	endocrinology	contraceptive
eszopiclone	GABA-A receptor; anion channel positive allosteric modulator	neurology/psychiatry	insomnia
famotidine	histamine H2 receptor antagonist	gastroenterology	heartburn
fingolimod	sphingosine 1-phosphate receptor agonist	neurology/psychiatry	multiple sclerosis
flunarizine	voltage-gated T-type calcium channel blocker	neurology/psychiatry	migraine headache
fluvoxamine	serotonin transporter inhibitor	neurology/psychiatry	obsessive–compulsive disorder (OCD)
gabapentin	voltage-gated calcium channel modulator	neurology/psychiatry	epilepsy
galantamine	acetylcholinesterase inhibitor	neurology/psychiatry	senile dementia
guanfacine	adrenergic receptor alpha-2 agonist	cardiology	hypertension
isradipine	voltage-gated L-type calcium channel blocker	cardiology	hypertension
lamotrigine	sodium channel alpha subunit blocker	neurology/psychiatry	epilepsy
levetiracetam	synaptic vesicle glycoprotein 2A modulator	neurology/psychiatry	epilepsy
levodopa	dopamine D3 receptor agonist	neurology/psychiatry	Parkinson’s disease
linagliptin	dipeptidyl peptidase IV inhibitor	endocrinology	diabetes mellitus
lorazepam	GABA-A receptor; anion channel positive allosteric modulator	neurology/psychiatry	epilepsy
losartan	type-1 angiotensin ii receptor antagonist	cardiology	hypertension
mecamylamine	Nach receptor (a3/b4) negative allosteric modulator	cardiology	hypertension
meclofenamic acid	arachidonate 5-lipoxygenase inhibitor	rheumatology	rheumatoid arthritis
memantine	glutamate (NMDA) receptor negative allosteric modulator	neurology/psychiatry	Alzheimer’s disease
metformin	5′-AMP-activated protein kinase activator	endocrinology	diabetes mellitus
methotrexate	dihydrofolate reductase inhibitor	oncology	acute lymphoblastic leukemia (ALL)
minocycline	bacterial 70S ribosome inhibitor	infectious disease	respiratory tract infections
mirtazapine	serotonin 2a (5-HT2a) receptor antagonist	neurology/psychiatry	depression
modafinil	dopamine transporter inhibitor	neurology/psychiatry	shift work disorder (SWD)
nicotine	Nach receptor (a4/b2) agonist	neurology/psychiatry	smoking cessation
nitroglycerin	soluble guanylate cyclase activator	cardiology	angina pectoris
nitroprusside	soluble guanylate cyclase activator	cardiology	hypertension
ondansetron	serotonin 3a (5-HT3a) receptor antagonist	gastroenterology	nausea
oxcarbazepine	sodium channel alpha subunit blocker	neurology/psychiatry	epilepsy
oxybate	GABA-B receptor agonist	neurology/psychiatry	narcolepsy
papaverine	phosphodiesterase 4 inhibitor	cardiology	myocardial infarction
pentosan polysulfate	fibroblast growth factor 2 binder	urology	interstitial cystitis (IC)
pergolide	dopamine receptor agonist	neurology/psychiatry	Parkinson’s disease
pramipexole	D2-like dopamine receptor agonist	neurology/psychiatry	Parkinson’s disease
pravastatin	HMG-coa reductase inhibitor	endocrinology	hypercholesterolemia
prednisolone	glucocorticoid receptor agonist	ophthalmology	conjunctivitis
pregabalin	voltage-gated calcium channel modulator	neurology/psychiatry	peripheral neuropathy
pregnenolone	nuclear receptor subfamily 1 group I member 2 agonist	rheumatology	rheumatoid arthritis
progesterone	progesterone receptor agonist	obstetrics/gynecology	infertility
pyrimethamine	dihydrofolate reductase inhibitor	infectious disease	malaria
raloxifene	estrogen receptor beta modulator	orthopedics	osteoporosis
ramelteon	melatonin receptor agonist	neurology/psychiatry	insomnia
rasagiline	monoamine oxidase B inhibitor	neurology/psychiatry	Parkinson’s disease
reboxetine	norepinephrine transporter inhibitor	neurology/psychiatry	depression
riluzole	sodium channel alpha subunit blocker	neurology/psychiatry	amyotrophic lateral sclerosis (ALS)
roflumilast	phosphodiesterase 4 inhibitor	pulmonary	COPD
selegiline	monoamine oxidase B inhibitor	neurology/psychiatry	Parkinson’s disease
sertraline	serotonin transporter inhibitor	neurology/psychiatry	depression
sildenafil	phosphodiesterase 5A inhibitor	cardiology	hypertension
simvastatin	HMG-coa reductase inhibitor	endocrinology	hypercholesterolemia
tiagabine	GABA transporter 1 inhibitor	neurology/psychiatry	epilepsy
tolcapone	catechol O-methyltransferase inhibitor	neurology/psychiatry	Parkinson’s disease
topiramate	GABA-A receptor; anion channel positive modulator	neurology/psychiatry	epilepsy
tropisetron	serotonin 3a (5-HT3a) receptor antagonist	gastroenterology	nausea
valacyclovir	human herpesvirus 1 DNA polymerase inhibitor	infectious disease	virus herpes simplex (HSV)
valproic acid	4-aminobutyrate aminotransferase inhibitor	neurology/psychiatry	epilepsy
varenicline	Nach receptor (a4/b2) agonist	neurology/psychiatry	smoking cessation
vorinostat	histone deacetylase 1 inhibitor	oncology	cutaneous T-cell lymphoma (CTCL)
vortioxetine	serotonin 1a (5-HT1a) receptor agonist	neurology/psychiatry	depression

The creation of the shared feature pharmacophore was extracted by an alignment and interpolation of the merging pharmacophores of the native ligands with DRs. The model, which is presented in this research (Figure ), only contains four features commonly aligned in the pharmacophore. The four-point pharmacophore model was selected as the input for the virtual screening against the 70 compounds of this research.

Figure 3

Construction of the pharmacophore-shared feature model. (A) Three hydrophobic interactions are represented by the yellow spheres and one positive ionizable area, which is represented as a blue astral center. (B) Four-point pharmacophore model based on shared features of the superposed structures 8NU-risperidone (red), ETQ-eticlopride (green), and AQD-nemonapride (magenta), crystallographic ligands of proteins 6CM4, 3PBL, and 5WIU.

The algorithm of screening, pharmacophore fit, retrieved a total of six structures as the best match of the pharmacophore model. All the selected structures had a minimum of three features located in the same disposition of the model. The selected compounds of the pharmacophore input were A - linagliptin, B - citalopram, C - flunarizine, D - sildenafil, E - minocycline, and F - duloxetine (Figure ). Each of the six drugs has one positive ionizable area that is susceptible to protonation and two hydrophobic points positioning in the same conformation and 3D space the input pharmacophoric model. The atoms or groups classified as positive ionizable are susceptible to be protonated, such as the 3-amine-piperidine region in linagliptin, the 3-(dimethylamino)propyl group in citalopram, the piperazine group in flunarizine, the 1-methylpiperazine group in sildenafil, the 5-dimethylamino group in minocycline, the N-methylamine group in duloxetine. The hydrophobic areas were in the aromatic, saturated, unsaturated, and fluoridated regions. The physicochemical properties of the selected six drugs are shown in Table .

Figure 4

Pharmacophore screening results of the four-point shared pharmacophore. The selected compounds as a perfect match were (A) linagliptin, (B) citalopram, (C) flunarizine, (D) sildenafil, (E) minocycline, and (F) duloxetine. The yellow spheres represent the hydrophobic interactions, and the blue areas represent the positive ionizable areas.

Table 2

Structure and Physicochemical Properties of the Drugs Selected as the Best Match of the Pharmacophore Screeninga

XLOGP3: Octanol–Water Partition Coefficient Atomistic method, optimal range between −0.7 and +5.0. TPSA: Topological polar surface area, optimal range between 20 and 130 Å2. Ali LogS: Water solubility, implemented by the optimal range between.[104] Fraction Csp3: Saturation, the fraction of carbons in the sp3 hybridization not less than 0.25. Flexibility: Rotatable bonds, no more than 9 rotatable bonds.

To confront the results of the pharmacophore screening, molecular docking calculations with the D2, D3, and D4 dopamine receptors were conducted and they are shown in Table . The exactness and precision of the docking process during the simulation were validated by redocking analysis of the crystallographic ligands (Table ). To avoid inexact results, the root-mean-square deviation (RMSD) was rooted in the threshold value of 2.0 Å: this process is termed “self-docking”. By self-docking it was shown that the binding sites had the optimal size and fidelity of the crystallographic complexes. These results show that the docking program, AutockVina, adopted an accurate binding pose during the whole simulation, and it has the potential to compute with confidence the computational ligand-binding affinities of 70 compounds to D2, D3, and D4 dopamine receptors.

Table 3

Molecular Docking Scores for Drugs with D2, D3, and D4 Receptorsb

Docking score (kcal/mol) calculated using AutoDock Vina.

GI absorption: High, PSA < 142 Å2 and log P between −2.3 and + 6.8. Low: outside the high GI range. BBB Permeant: Yes: PSA < 79 Å2 and log P between +0.4 to +6.0. No: outside the high BBB permeant range. The drugs selected in the pharmacophore screening are presented in bold (see the complete table in the Supporting Information).

Table 4

Self-Docking Results and RMSD of Crystallographic Ligands of the Proteins 6CM4, 3PBL, and 5WIU

dopamine receptor	crystallographic ligand code	RMSD (Å)	docking score (kcal/mola)
D2 (6CM4)	8NU	1.244	–8.9
D3 (3PBL)	ETQ	1.232	–4.9
D4 (5WIU)	AQD	1.442	–6.9

Docking score (kcal/mol) calculated using AutoDock Vina.

Consequently, the key element of the drugs used to treat schizophrenia is the crossing of the blood–brain barrier. The study hereafter is guided into the selection of the repurposing molecules with the best docking scores and in conjunction with a high probability of being absorbed by the human intestine and a high possibility of penetration of the BBB. Predictive GI as well as BBB penetration were added at the bottom of Table . Interestingly, flunarizine showed better docking scores in comparison to crystallographic native ligands and the six drugs selected in the pharmacophore-based model. From the six pharmacophore matchings, both flunarizine and citalopram were simultaneously predicted with high GI absorption and BBB permeation. Flunarizine was coupled with −9.9, −9.2, and −11.0 kcal/mol, respectively, with D2, D3, and D4 dopamine receptors, and citalopram was coupled with −7.8, −7.2, and −8.2 kcal/mol to D2Rs. Although native ligands must have a high affinity to D2R targets, the docking scores exposed in Table were not better than those accomplished by flunarizine.

Complexes of the docked structures of flunarizine with D2, D3, and D4 receptors are compared in Figure . Overall, the key components of flunarizine were the phenyl, 4-fluorobenzyl groups, and tertiary amine in the piperazine ring. In particular, this piperazine ring is plausible to be charged at physiological pH and it could form a salt bridge to the D2R. In accord with the shared pharmacophore input, the aforementioned groups established hydrophobic regions, polar contacts, and salt bridges with the three DRs. The main binding sites of flunarizine with the D2-like receptors are described in Figure , most of the contact residues in flunarizine were identical to the native ligands, and even some distances in those conserved amino acids are smaller in flunarizine than in risperidone, eticlopride, and nemonapride.[46,51] In general, hydrogen bonds are divided into a strong hydrogen bond when the involucrated atoms are N–H···O, N–H···N, and O–H···O, and when the hydrogen bond is built between, C–H···O, it is considered as a weak hydrogen bond. Although the type of involved atoms is a paramount factor of the hydrogen bonds, other aspects must be included as discrimination elements. The distance between acceptors and donor atoms is another aspect that is recommended to contemplate. For instance, weak hydrogen bonds are usually among 3.0–3.7 Å, and strong hydrogen bonds should measure over 2.0–3.9 Å. Besides those facts, hydrophobic interactions could be separated into five different classes, and they could be between carbon–carbon, carbon–halogen, or aromatic carbon–sulfur with distances slightly longer, usually ≤4.0 Å. Salt bridges are established with two regions of opposite charge, for instance, between a positively charged nitrogen and a negatively charged oxygen, and the average distance between is approximately 3.0–4.0 Å.[52]

Figure 5

Docked structures of the flunarizine and the native ligands with D2-like receptors. (A) Amino acids that were involucrated in the pharmacophore between flunarizine (deep purple) and D2R (light blue), amino acids that were involucrated in the pharmacophore between the native ligand, risperidone (green), and D2R (lightblue). (B) Amino acids that were involucrated in the pharmacophore between flunarizine (deep purple) and D3R (bright orange), amino acids that were involucrated in the pharmacophore between the native ligand, eticlopride (blue), and D3R (bright orange). (C) Amino acids that were involucrated in the pharmacophore between flunarizine (deep purple) and D4R (light cyan), amino acids that were involucrated in the pharmacophore between the native ligand, nemonapride (brown), and D4R (light cyan).

Five key residues in the D2 dopamine receptor were selected to measure their distances to the ligands (Figure A): PHE-377, TRP-67, TYR 396, and PHE 377, which are aromatic amino acids, and ASP-81, which is negatively charged. Aspartic acid amino acid (ASP) is reported as one of the most frequent acceptors in weak hydrogen bond interactions, C–H···O interactions,[52] which had a reduction from 4.0 Å with risperidone to 3.8 Å with flunarizine. In terms of D3, six amino acids were chosen to analyze the ligand binding, and nonpolar, aromatic, and polar amino acids were compared (Figure B). VAL-86, a nonpolar amino acid, was markedly reduced to 3.2 Å, establishing a stable interaction of carbon-halogen type.[48] Representative binding amino acids of the binding pocket between the D4-flunarizine and D4-nemonapride[53] complexes were selected (Figure C). Some interactions slightly increased in the D4-flunarizine complex, for instance, ASP-115 increased from 3.0 to 3.4 Å and SER-196 increased from 3.4 to 4.1 Å. Yet, other interactions, especially those related to aromatic amino acids, were reduced (e.g., PHE-91 and PHE-363). So, hydrophobic interactions were the main driving force in flunarizine–D2Rs complexes.

To study the stability of the flunarizine and native ligands with the DRs, the protein–ligand complexes were analyzed by molecular dynamics simulations, and RMSD, RMSF, and radius of gyration (Rg) were examined (see Figure ). The ligand topologies were generated using the CHARMM General Force Field due to the fact that it groups a huge scope of chemical groups presented in drug-like molecules, as well as many heterocyclic scaffolds.[54] The CHARMM36 force field was used to study the protein.[55−60] Another field, OPLS3e, has been studied by other authors to generate the protein field of DRs simulations, however, significant differences have not been reported between the two fields, and CHARMM36 has been primarily recommended.[61,62] Aiming to generate the first report of the molecular dynamic simulation of DRs with flunarizine and its comparison with the risperidone, eticlopride, and nemonapride, 30 ns was selected as the running time of the present simulations. Each parameter was adjusted to the final purpose of this article. To avoid the molecule becoming unstable and any error during the integration, it is recommended the use a short time step. Thus, selecting the time, time step, and the number of steps was a critical step of the present molecular dynamic study. Simulation time and the time step are linear in molecular dynamics simulation; for the specific objective of our study, the simulation time was 30,000 ps (30 ns) with 15,000,000 steps and a time step of 2 fs.[63,64]

Figure 6

Molecular dynamics simulations of DRs unbound and the protein complex with crystallographic ligands and flunarizine: (A) RMSD, (B) RMSF, and (C) radius of gyration plots.

To measure the difference between the backbone of a protein from an initial structure conformation to the final position, the RMSD analysis can be used. When a protein is unstable, its fluctuations in the RMSD protein are high. Regarding the atomic positions in the trajectory, root-mean-squared fluctuations (RMSF) may be used to calculate their mean square fluctuations, and the radius of gyration calculates the closeness of the protein structure from the start to the end of simulations.[65,66]

The RMSD analysis and tracing during the whole simulation with the selected compound, flunarizine, show that the RMSD value starts with an increasing trend with the three DRs (Figure A). The complex of flunarizine with the D2 (D2_FLUN) dopamine receptor had a low RMSD course in comparison to the risperidone complex during all the running times, showing that the flunarizine molecule binds more stably to the targeting site of the D2 receptor. The complex flunarizine–D3 dopamine receptor (D3_FLUN), at the initial stage of the molecular dynamics simulation, had a gradual decline until 1.25 ns and then a sightly increase up to 2.5 ns, which are not statistically significant, and the trend of the complex shows a lower RMSD in comparison to the complex with eticlopride (D3_EQT). Similarly, in the simulation of the D4 dopamine receptor, comparable trajectories were observed for the D4 dopamine receptor, in its unbounded manner, and complexes with flunarizine (D4_FLUN) and nemonapride (D4_ADQ), despite the fact that there was a peak increase in the RMSD between 15 and 20 ns for the complex D4_ADQ. In all cases, the binding of the compound flunarizine showed a better stabilization to the complexes than the native co-crystalized ligands.

The dynamic behavior of individual amino acid residues is shown in Figure B. The RMSF plots showed a matching pattern regarding the flexibility influence of the crystallographic compounds with the flunarizine at different levels. The regions of the DRs that are most susceptible to a 300 K temperature are the residues 50–75, 100–130, 140–160, 200–210, 215–250, 300–340, and 380–400 in the D2 dopamine receptor; residues 200–250, 275–300, and 325–375 in the D3 dopamine receptor; and residues 25–40, 50–60, 90–110, 140–160, 175–180, 220–260, 270–290, 300–325, and 360–380 in the D4 dopamine receptor. The radius of gyration showed a similar pattern with the complex during the simulation. The complex of flunarizine with the D2 dopamine receptor had a variation of radius of gyration between 2.8 and 3.0 nm during all the simulations. The complex of flunarizine with the D3 dopamine receptor had a variation of 2.4–3.0 nm, and the complex of flunarizine with the D4 dopamine receptor had a variation of 2.7–3.2 nm. These analyses evidenced that flunarizine had a similar pattern in the protein structure compactness with the three DRs. The RMSD analyses indicate a more stable complex between flunarizine with DRs in comparison to the crystallographic compounds, risperidone, eticlopride, or nemonapride.

To analyze a set of diverse conformations that resulted from the molecular dynamics simulations, snapshots from the 30,000 ps simulation (ignoring ions and water) at each 5000 ps intervals (Figures S1–S3) were created from the GROMACS simulation results. Figure S1 shows the structures of flunarizine–DRs complexes at different stages, beginning from 5000 to 30,000 ps, and each stage was generated with the gmx trajcon command of GROMACS.[67−70]

The main contacts of flunarizine were the phenyl, 4-fluorobenzyl, and tertiary amine in the piperazine ring as well as in the docked structures. Flunarizine, during the whole simulation, mainly interacts with nonpolar and aromatic amino acids and, in some frames, with established polar contacts with THR and CYS amino acids. Due to the fact that flunarizine has six rotatable bonds and four H-bond acceptors, it could adopt multiple shapes within the pocket.

Considering all the results, it is clear that flunarizine not only showed the greatest affinity energy and stability with the D2, D3, and D4 receptors in comparison to the native ligands risperidone, eticlopride, and nemonapride, but also the key features and distances were significantly ameliorated in flunarizine outcomes. The shown results of this study are supported by previous research where the clinical antipsychotic efficacy of flunarizine has been questioned, especially by its D2 blocking properties.

When it comes to how flunarizine has been involucrated into the schizophrenia spectrum, there is just one report of flunarizine with schizophrenia, and very promising results were retrieved although it was the first approach of the evaluation of its antipsychotic properties. Bisol et al.[71] emphasized that flunarizine could show good efficacy and tolerability for the treatment of schizophrenia with an uncommon profile. It was worth noting that the study scored the patients by PANSS (Positive and Negative Syndrome Scale), along with the Clinical Global Impressions-Improvement (CGI-I) scale and the Extrapyramidal Symptom Rating Scale (ESRS). The key points of the study were some facts such as the important reduction in the PANSS with the flunarizine group (21%) and compared to those using haloperidol (19%). Other motivators such as the long half-life (2–7 weeks), low cost, and low inference of extrapyramidal symptoms justified advanced investigation, particularly in psychiatric patients with low adherence to treatment. The present study is the first computational support of flunarizine as a possible treatment in schizophrenia by its high affinity and stability with DRs compared with antipsychotic drugs risperidone, eticlopride, and nemonapride, and the results presented here should be contrasted with more clinical trials and experimental data.[72−79] Further studies should be carried out to establish whether the possible antipsychotic effect of flunarizine is due to its interaction with D2 receptors or as a calcium channel antagonist or both.

Conclusions

An in silico drug repositioning study was conducted. Considering that D2-like receptors (D2, D3, and D4 receptors) are one of the main pharmacological targets for treating schizophrenia, a hypothesis of shared pharmacophore features was formulated for this study. The common features of risperidone, eticlopride, and nemonapride with DRs were studied and screened with a total of 70 promising compounds that are not currently indicated for schizophrenia. Linagliptin, citalopram, flunarizine, sildenafil, minocycline, and duloxetine were the drugs that perfectly matched the pharmacophore input. Blood–brain penetration and human intestinal absorption prediction results were analyzed and compared with the pharmacophore results. From the six drugs selected in the pharmacophore-share input, flunarizine showed the best docking score with the D2, D3, and D4 dopamine receptors (D2-like receptors) and high stability in comparison to all the tested ligands. Molecular docking and molecular dynamics simulations were used to corroborate the binding score and the complex stability of flunarizine–D2-like receptors. Flunarizine has been hypothesized in some previous studies by its D2-like dopamine blocking properties. These results are the first computational reports where it is supported that flunarizine should be studied as a D2-like dopamine blocker. Flunarizine has many advantages in comparison to other antipsychotics. Further studies must be performed with it in the field of psychopharmacology. Low induction of extrapyramidal symptoms, a long half-life (2–7 weeks), and low cost are some facts that justified further investigation, particularly in psychiatric patients with low adherence to other treatments.

Materials and Methods

Selection and Optimization of the Compounds

The suggested compounds in this study were selected using the following inclusion criteria: (1) having been tested in clinical studies (ClinicalTrials.gov)[1] with a pharmacological target approved for schizophrenia by the FDA but that had not yet been approved to treat schizophrenia and (2) have not been reported side effects, alone or in combination with other drugs. Thus, 70 (out of a group of 88 possible drugs reported in 231 clinical studies) were selected, which were the only nonprotein and non-antibody compounds (Table ). The compounds’ 3D structures were downloaded from the PubChem database in “.sdf” format, and then they were submitted to quantum chemistry models, such as ab initio calculations with the DFT method with the B3LYP multiparameter functional and the 6-31G Gaussian basis set, using the program Gaussian 09.[80] The resulting optimized geometry was converted to “.mol2” format, and subsequently polar hydrogens, Gasteiger charges, and free bond rotation were incorporated to the molecules in “.pdbqt” format file using Open Babel 2.4.1.[81]

Selection and Optimization of the Targets

The selected pharmacological targets for schizophrenia were the dopamine receptors D2, D3, and D4.[53,82] The criteria of selection to crystallographic proteins were as follows: (1) organism, human as a source organism; (2) resolution at a minimum of 3.0 Å; (3) structures with identified and known binding sites with natural ligands of dopamine crystallographic receptors; and (4) wild type protein and without variants or mutations.[83]

The crystallographic structures of the used proteins in the present study were obtained from the Protein Data Bank, whose codes are PDB ID: 6CM4 (D2 receptor), 3PBL (D3 receptor), and 5WIU (D4 receptor). The 3D crystallographic structures of PDB proteins co-crystallized with the ligands 6CM4-8NU (risperidone), 3PBL-ETQ (eticlopride), 5WIU-AQD (nemonapride) were curated using the Swiss model web server to select reliable and complete three-dimensional protein structures without missing residues.[105−107] The 3D structures were optimized with the program UCSF Chimera,[84] which means the water molecules and other nonprotein ligands were removed and assigned AMBER force field charges as it is reported from the literature that it used a gradient convergence criteria of 0.005 kcal/mol. The structures were placed in their minimum energy conformation and were stored in “.pdbqt” format.

Pharmacophore Modeling and Screening

Crystal structures of complexes (D2 dopamine receptor with risperidone, D3 dopamine receptor with eticlopride, and D4 dopamine receptor with nemonapride) were superposed to create a shared feature pharmacophore using LiganScout 4.4. A four-point pharmacophore model was selected by combining all the original features and screened against 70 repurposing drug candidates, the algorithm of screening was the pharmacophore-fit scoring function, the conditions were a maximum number of omitted 1 feature, and the minimum required feature was 3. The screening mode was the matching of all query features and Multi-Thread execution mode. The interactions considered were hydrophobic interaction (represented as a yellow single sphere, distance constraint 1.0–5.9 Å), hydrogen bond donors (represented as a green dotted arrow, distance constraint 2.2–3.8 Å), a hydrogen bond acceptor (represented as a red dotted arrow, distance constraint 2.2–3.8 Å), positive ionizable interaction (represented as a blue astral center, distance constraints 1.5–5.5 Å with negative ionizable and 1.5–10.0 Å with aromatic ring), among others.[85]

Gastrointestinal Absorption and Blood–Brain Penetration Parameters

The SwissADME server and admetSAR online database were used to check ADME parameters of the drugs, absorption, distribution, metabolism, and excretion, along with some physicochemical properties. The blood–brain penetration and human intestinal absorption estimation were computed by the lipophilicity (WLOGP) and polarity (tPSA) of the molecules using the SwissADME web server.[36,86−88] This model is highly reliable because it is supported by a big library of molecules that were curated by the literature, patents, and database cross-checks. The guidelines for good absorption[89] are commonly accepted by a rectangular limit of log P between −2.3 and +6.8 with the PSA lower than 142 Å2. The molecules with a high possibility of ingress in the CNS had moderately polar (PSA < 79 Å2) and relatively lipophilic (log P from +0.4 to +6.0).[87,90−92]

The Combination of these predictive models had shown an accurate predictive model of molecules being absorbed by the human intestinal absorption (HIA) and blood–brain barrier (BBB), and it is of great support for lead optimization.[88,93−98]

Molecular Docking

With the ligands and proteins in “.pdbqt” format, molecular docking computation was performed in triplicate using AutoDock Vina,[24] the flexible side chains protocol, and an iterated local search global optimizer algorithm. Thus, flexible docking of the 70 compounds with each of the molecular targets (D2, D3, and D4) was automatically performed, through a bash script, running on a Precision 7920 workstation using a Linux Ubuntu 18.04. The docking site was selected upon the location of the co-crystalized ligands found in each PDB file: for the D2 receptor, the XYZ coordinates chosen to be the center of the grid box were 9.9, 5.8, and −9.5 Å, respectively; for the D3 receptor, the XYZ coordinates were 0.08, −14.8, and 10.4 Å, respectively; and for the D4 receptor, the XYZ coordinates were −17.4, 15.2, and −16.7 Å, respectively. For each docking, a grid box (that is, a grid or defined space large enough to allow the docked ligand to rotate freely) of x 25, y 25, and z 25 Å was constructed with a number of bonding modes equal to 20 and exhaustivity of 25. The pose with the most negative affinity (kcal/mol) was stored. Redocking estimates were studies as an accuracy process of the docking scores in the binding poses.

Molecular Dynamics Simulations

GROMACS 2018.1 was used for the molecular dynamics simulation of DR structures (PDB: 6CM4, 3PBL, 5WIU) and the ligands (8NU, ETQ, AQD, and flunarizine).[99−101] CHARMM36 force field was used for protein and ions, meanwhile the TIP3P model for water.[57,102] CHARMM General Force Field (CGenFF)[54] was used to retrieve the ligand parameters (legacy version 1.0.0). A dodecahedron box[103] and other periodic boundary conditions were optimized to perform the simulations. To solvate and neutralize the system, ions were added. A total of 5000 steps were performed in terms of energy minimization using the steep descent method of having a stable conformation. Canonical ensembles (NVT) and isobar isothermal ensembles (NPT) were performed, respectively, with a constant temperature of 300 K for 100 ps for NVT followed by a constant temperature of 300 K and a constant pressure of 1 atm per 100 ps for NPT. Molecular dynamics simulation was performed for 30,000 ps. The first step was minimized and then equilibrated with restraints on the ligand, its heavy atoms, and protein backbone atoms, followed by production runs. The Verlet cutoff scheme was used, a gid box with a nstlist of 10 was used, and rcoulomb and rvdw were fixed to 1.0 nm. The energies were saved every 10.0 ps with an update log file every 10.0 ps. The simulation covered 15 million steps with 30,000 ps (30 ns). Based on the average of these results, longer (100 ns) MD simulations should be employed in the following study. The RMSD, RMSF, and radius of gyration (Rg) were calculated by g-rmsd, g-rmsf, and g-Rg, respectively. Different snapshots were generated by the g-trjconv command, using 5000 ps intervals. The graphics for RMSD and RMSF were designed using the QtGrace v0.2.6 program.