Computer aided identification of sodium channel blockers in the clinical treatment of epilepsy using molecular docking tools.

UNLABELLED: Phenytoin (PHT) and Carbamazepine (CBZ) are excellent sodium channel blockers administered in clinical treatment of epileptic seizures. However, the narrow therapeutic range and limited pharmacokinetics of these drugs have raised serious concerns in the proper management of epilepsy. To overcome this, the present study attempts to identify a candidate molecule with superior pharmacological profile than PHT and CBZ through In silico approaches. PHT and CBZ served as query small molecules for Tanimoto based similarity search with a threshold of 95% against PubChem database. Aided by MolDock algorithm, high affinity similar compound against each query was retrieved. PHT and CBZ and their respective similar were further tested for toxicity profiles, LC 50 values and biological activity. Compounds, NSC403438 and AGN-PC-0BPCBP respectively similar to PHT and CBZ demonstrated higher affinity to sodium channel protein than their respective leads. Of particular relevance, NSC403438 demonstrated highest binding affinity bestowed with least toxicity, better LC 50 values and optimal bioactivity. NSC403438 was further mapped for its structure based pharmacophoric features. In the study, we report NSC403438 as potential sodium channel blocker as a better candidate than PHT and CBZ which can be put forth for pharmacodynamic and pharmacokinetic studies. ABBREVIATIONS: AEDs - Antiepileptic drugs, BLAST - Basic Local Alignment Search Tool, CBZ - Carbamazepine, GEFS+ - Generalized Epilepsy with Febrile Seizures Plus, GPCR - G Protein Coupled Receptor, Nav - Sodium channel with specific voltage conduction, PDB - Protein Data Bank, PHT - Phenytoin, PIR - Protein Information resources, SAVES - Structural Analysis and Verification Server, VGSC - Voltage-gated Sodium channels.

Background

International League Against Epilepsy [ILAE] defines epilepsy as a disease condition characterized by an enduring predisposition to generate seizures by neurobiological, cognitive, psychological and social consequences [1]. The worldwide prevalence of active epilepsy ranges from 4 to 10 per thousand populations [2] and in India almost 5.59 per thousand suffers epileptic attack once in their life time [3]. The electrogenic property of an individual neuron forms an important marker for hyper excitability of neuronal circuits which are dependent on the functional properties of ion channels like Na+, K+, and Ca2+ in the membranes. VoltageGated Sodium Channels [VGSCs], in particular, are the mediators of intrinsic neuronal excitability and are central to most important determinants of pathophysiology of epileptic seizures and execute initiation of action potentials, synaptic transmission and neurotransmitter release.

The major role of sodium channels in epileptic pathogenesis is reflected from the studies wherein mutations on chromosome 21q were coupled to progressive myoclonus epilepsy of the Unverricht–Lundborg type. The affected gene was later identified to code for a sodium channel protein, which evidently indicated sodium channels may play a crucial role in this syndrome [4]. Additionally, involvement of sodium channels in epilepsy came from yet another compelling evidence wherein gain-of-function mutations on the β1 sodium channel subunit gene -SCN1b on chromosome 19q that modified gating capability of channels was linked to generalized epilepsy with febrile seizures plus (GEFS+) type [5]. The absolute necessity for intact sodium-channel function in normal neurophysiology was underscored by the finding that knockout mice with specific ablations of sodium channels (Nav1.2, Nav1.5 or Nav1.6) develop seizures [6, 7, 8]. Given such accumulating evidences from past 15 to 20 years paved way in rational designing of sodium channel blockers like Phenytoin (PHT), Carbamazepine (CBZ) which are still regarded as best- prescribed medication amongst currently available antiepileptic drugs (AEDs) [9, 10]. Sodium channel blockers act by preventing the repetitive firing of the axons by stabilizing the inactive form of channel. In addition, presynaptic and postsynaptic blockade of sodium channels of the axons causes stabilization of the neuronal membranes, prevents post tetanic potentiation, and finally limits seizure activity [11].

Despite the target hypothesis of AED treatment, around 30% of the patients still continue to have uncontrolled seizures (drug resistance) and drug toxicity which may be due to structural or functional change at the site of drug action or alteration in the drug pharmacodynamics [12, 13, 14]. Hence, development of novel sodium channel blockers with an optimal efficacy that can enhance improved clinical outcome is undoubtedly an important medical demand. In the view of above, the present study was designed to identify novel compounds through structure and ligand based strategies which we anticipate would have superior channel blocking potential endowed with optimal ADMET features than customary drugs - Phenytoin and Carbamazepine.

Methodology

Selection of multiclass sodium channel blockers:

Table 1 (see supplementary material) shows two established sodium channel blockers selected as leads for virtual screening.

Molecular Modeling of SCN1A protein:

Protein sequence of voltage-gated sodium channel alpha subunit SCN1A was retrieved from NCBI with accession number: AAK00217.1. Crystallographic and solution NMR structures of the protein were available in PDB [4JPZ, 2KBI] [15, 16] however structures were interrupted with missing residues and loops. Therefore, complete structure was predicted with computational homology modeling. Similarity searching was performed against PDB database for finding an appropriate template for homology modeling using BLOSUM- 62 matrix enabled BLAST program. Top 10 templates were used for the alignment against SCN1A. Threading alignment with a normalized Z-score >1 was considered optimal. The entire top 10 template alignment file (.ali) was used for building loops using MODELLER program. The FASTA was converted to PIR using EMBL׳s Readseq algorithm. Structure similarity was performed by using the profile.build (), an inbuilt command in MODELLER. The build_profile.py was then used for the local dynamic algorithm to identify homologous sequences against SCN1A target sequence. The final model thus obtained (Figure 1a) from MODELLER was further used for structure validation using by Procheck using UCLA׳s SAVES server.

Figure 1

A) Secondary structure view of the modeled SCN1A protein. The top view of the protein depicts an ion pore (green) surrounded by four domains; B) The ion-conducting aqueous pore (solid green cylinder) calculated by V3 (Volume Calculation and Extraction Procedures) web server by probing measures.

Prediction of channel in the protein:

3V web server was employed to predict the channel through comprehensive analyses of the internal volumes considering difference as large as possible probe radius and the solvent radius (typically 1.5 Å for water) [17]. The volumetric representations of the channel are provided in Table 2 (see supplementary material).

Similarity search, preparation of protein and compounds :

Similarity search was supervised by Binary Finger Print Based Tanimoto similarity equation to retrieve compounds similar to PHT and CBZ with similarity threshold of 95 % against NCBI׳s Pubchem compound database. All the structures were optimized through OPLS 2005 force field algorithm [18] embedded in the LigPrep module of Schrödinger suite, 2013 (Schrodinger. LLC, New York, NY) [19]. The three-dimensional structure of the modeled structure was processed by removing all bound crystal water molecules and adding hydrogen bonds. Explicit hydrogen, bond orders, disulphide bonds, hybridizations and charges were assigned wherever missing. The resulting structure was energy minimized using OPLS- 2005 force field by protein preparation wizard of Schrödinger suite 2013.

Toxicity screening and bioactivity prediction of compounds:

All the similar compounds retrieved were screened for its drug ability by lipinksi filters. The toxicity screening was achieved by using LAZAR toxicity prediction server [20]. Biological activity of the ligands was predicted using Molinspiration webserver (© Molinspiration Cheminformatics 2014). LC 50 was predicted using T.E.S.T. Version 4.1 (2012, U.S. Environmental Protection Agency) software. The complete ADMET properties was calculated using admetSAR [21].

Optimization of virtual screening parameters:

Molecular docking program- Molegro Virtual Docker (MVD) which incorporates highly efficient PLP (Piece wise Linear Potential) and MolDock scoring function provided a flexible docking platform [22– 24]. The leads (PHT and CBZ) and similar chemical structures were docked in predicted channel of protein. Docking parameters were set to 0.20Å as grid resolution, maximum iteration of 1500 and maximum population size of 50. Energy minimization and hydrogen bonds were optimized after the docking. Simplex evolution was set at maximum steps of 300 with neighborhood distance factor of 1. Binding affinity and interactions of similar compounds with protein were evaluated on the basis of the internal ES (Internal electrostatic Interaction), internal hydrogen bond interactions and sp2-sp2 torsions. Post dock energy of the ligand-receptor complex was minimized using Nelder Mead Simplex Minimization (using non-grid force field and H bond directionality) [25]. On the basis of MolDock − rerank score best interacting compound was selected respective to each parent compound.

Molecular docking program- Molegro Virtual Docker (MVD) which incorporates highly efficient PLP (Piece wise Linear Potential) and MolDock scoring function provided a flexible docking platform [22– 24]. The leads (PHT and CBZ) and similar chemical structures were docked in predicted channel of protein. Docking parameters were set to 0.20Å as grid resolution, maximum iteration of 1500 and maximum population size of 50. Energy minimization and hydrogen bonds were optimized after the docking. Simplex evolution was set at maximum steps of 300 with neighborhood distance factor of 1. Binding affinity and interactions of similar compounds with protein were evaluated on the basis of the internal ES (Internal electrostatic Interaction), internal hydrogen bond interactions and sp2-sp2 torsions. Post dock energy of the ligand-receptor complex was minimized using Nelder Mead Simplex Minimization (using non-grid force field and H bond directionality) [25]. On the basis of MolDock - rerank score best interacting compound was selected respective to each parent compound.

Softwares, Suites and Webservers used :

For virtual screening Pubchem database was used to search and prepare library of similar chemical compounds. All the chemical structures were drawn in MarvinSketch 5.6.0.2, (1998- 2011, Copyright © ChemAxon Ltd). Ligands were optimized with LigPrep module of Schrodinger suite 2013. Protein was processed and refined with protein preparation wizard of Schrodinger suite 2013 (Schrodinger. LLC, 2009, New York, NY). Flexible molecular docking of the compounds with target was completed using Molegro Virtual Docker 2010.4.0.0. Accelrys Discovery Studio® Visualizer 3.5.0.12158 (Copyright© 2005-12, Accelrys Software Inc.) was used for molecular visualizations. LAZAR online server was used to predict In silico toxicity. T.E.S.T software (2012, U.S. Environmental Protection Agency) and Molinspiration web server (© Molinspiration Cheminformatics 2014) were respectively used for predicting LC50 and bioactivity of the compound. ADMET profiles were calculated using admeSAR (Laboratory of Molecular Modeling and Design. Copyright @ 2012, East China University of Science and Technology, Shanghai Key Laboratory for New Drug Design,)

Results

Number of similar compounds screened with ≥95 similarity corresponding to each PHT and CBZ has been listed in Table 3 (see supplementary material). The procheck results revealed the modified status of the modeled structure of SCN1A protein (Figure 1A). In final model 97.6 percentages of overall amino acids were in allowed region of Ramachandran plot, validating the model in close proximate to experimental quality. Evident from rerank score, compound NSC403438 (CID: 345700) (Figure 2a) showed 1.47 folds better interaction than PHT whereas; compound AGN-PC-0BPCBP (CID: 57199333) (Figure 2b) was 1.29 folds better interacting compound than CBZ (Table 3). In further investigation we observed that compound NSC403438 was marginally (1.4 folds) better interacting drug than AGN-PC-0BPCBP Table 3 (see supplementary material). NSC403438 not only showed better interaction against the channel than its parent compound PHT, but also showed superior binding affinity than CBZ and rest of the virtually screened molecules (Table 3). The overall interaction profile of PHT, CBZ and their respective similar NSC403438 and AGN-PC-0BPCBP is shown in Table 4 (see supplementary material). The superior rerank score of NSC403438 can be probably attributed to its optimal electrostatic and hydrogen bond interactions Table 4 (see supplementary material). The LC 50 values at 96 hour interval were predicted to be 1.6 folds superior for NSC403438 than its parent compound PHT; similarly AGN-PC-0BPCBP had 2.4 folds better LC 50 values than CBZ Table 5 (see supplementary material). In addition the similar compounds identified against their parents showed enhanced bioactivity providing a clue for target specificity Table 5 (see supplementary material). Except for AGN-PC-0BPCBP, all the compounds were safe and predicted to be non-carcinogen and non-mutagenic Table 6 (see supplementary material). Further, results from ADMET prediction revealed that NSC403438 to be better drug like compound compared to AGN-PC-0BPCBP Table 7 (see supplementary material). In the present study, we were able to identify similar compounds to have better pharmacological profile than their parents, however, CBZ similar - AGN-PC-0BPCBP failed to pass carcinogenic filter in cell lines of DBS Hamster. Taking this fact into consideration, only PHT similar - NSC403438 was mapped for its structure based pharmacophoric features. Comprehensively shown in figure 3a, the molecule demonstrated van der Waals interactions with His 614, Pro 611, Arg 500, Leu 716, Arg 501 and Asp 606 and electrostatic interactions with Gln 554, His 553 Ser 607, 550 and 603, Val 610, Asn 499 and Arg 498. Further π- π interactions were observed with Arg 501 and 498. The ligand binding pattern of NSC403438 in the channel site is shown in Figure 3b. Electrostatic and hydrophobic interactions of NSC403438 in the channel is shown in Figure 4a & Figure 4brespectively.

Figure 2

Structure of A) NSC403438 (CID: 345700) - compound similar to PHT; B) AGN-PC-0BPCBP (CID: 57199333) - compound similar to CBZ

Figure 3

A) Interactions of NSC403438 in the channel of the SCN1A protein. Residues circled in green participate in van der Waals interaction with the ligand while residues in pink forms electrostatic interactions. Hydrogen bonds are shown as green arrows between ligand and residues Arg 498 and 501; B) Binding pattern of NSC403438 in the channel. The pink lines represent various interactions like electrostatic, van der Waals, stearic, hydrogen bonding and hydrophobic interactions that enable energetically favourable binding of the ligand in the channel.

Figure 4

A) NSC403438 deeply embedded in the channel surrounded by highly electropositive residues; B) The channel harboring NSC403438 is shown with hydrophobic intensities. The hydrophobic intensities of the binding site ranges from - 3.00 (least hydrophobic area - blue shade) to 3.00 (highly hydrophobic area –brown shade).

Discussion

Since, serendipitous discovery of potassium bromide in midnineteenth century, AEDs have emerged as the most effective treatment for hyperexcitable neuronal network. These anticonvulsants were the mainstays of seizure treatment until the 1990s, when newer AEDs with good efficacy, fewer toxic effects, better tolerability and without the need for blood level monitoring were developed [26].

Regardless of progress achieved in understanding the neuropharmacology in treatment of epilepsy, the mortality and morbidity associated with this disease have not appreciably declined [27] and it becomes quite obvious to explore for innovative avenues to epileptic therapy. Most AEDs have a narrow therapeutic window and many patients tolerate and need serum concentrations above the usual therapeutic range and some even experience adverse effects [25]. Secondly, narrow therapeutic range of AEDs is apparent from the fact that serious idiosyncratic effects like skin rashes which later advances to Stevens-Johnson syndrome have been observed within several weeks or months of initiating the drug dose [28]. The third type of adverse drug effect documented was the cumulative toxicity which usually occurs over years of treatment [25].

Published literature has shown that although the effect of PHT is optimally voltage dependent with IC50 ranging from 600 µM to <10 µM, in a holding-potential range of –90 to –50 mV [29, 30, 31, 32, 33]. Binding of PHT to channels was much slower and required longer depolarizations (up to seconds) [34, 35] to produce blocking effects which otherwise hinders emergency treatment. In contrast CBZ blocks the channels ~5 times faster, but with lower affinity, suggesting that antiepileptic effects might be expected to be more pronounced with short-duration discharges [36]. Considering the narrow therapeutic range of PHT and CBZ, we in possible attempt endeavored to overcome these setbacks by identifying compounds bestowed with high affinity, better kinetics, least toxicity and optimal bioactivity.

Two compounds viz., NSC403438 and AGN-PC-0BPCBP similar to PHT and CBZ were identified to have higher affinity than their respective parents. NSC403438 and AGN-PC- 0BPCBP showed 1.47 and 1.29 folds higher affinity than their respective parent compounds– PHT and CBZ. Further, NSC403438 proved to be still better than AGN-PC-0BPCBP with augmented affinity order of 1.14 folds. The superior affinity of NSC403438 can be attributed to its excellent interaction profile especially in terms of electrostatic and Hbonding interactions. Further, NSC403438 was the only compound in the study to demonstrate 2 hydrogen bonds and 8 electrostatic interacting residues participating in the interaction.

The ADMET profiles of these compounds revealed that NSC403438 was better compound and most likely drug gable compared to AGN-PC-0BPCBP. The predicted bioactivity as well as the LC 50 values of NSC403438 and AGN-PC-0BPCBP was quite appreciable. The LC 50 value at 96 hour interval was predicted to be 1.6 folds superior for NSC403438 than its parent compound PHT; similarly AGN-PC-0BPCBP had 2.4 folds better LC 50 values than its parent compound CBZ. In addition the similar compounds identified against their parents show enhanced bioactivity providing a clue for target specificity. The toxicity profiles of all the four compounds were although appreciable, however, AGN-PC-0BPCBP demonstrated to be carcinogenic in DBS Hamster cell line therefore it was not further analyzed for pharmacophoric mappings.

From the present study we predict compound NSC403438 can overcome the narrow therapeutic window of PHT and CBZ and can be to put forth for pharmacodynamic and pharmacokinetic experiments for better clinical outcomes in the successful treatment of epilepsy.

Conclusion

The narrow therapeutic range of PHT and CBZ necessitates an urgent need in the treatment of epilepsy. Therefore in the given view, through In silico strategies we identified NSC403438 to have better therapeutic properties endowed with least toxicity than hitherto administered AED׳s -PHT and CBZ.