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

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.

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 Ca 2+ in the membranes. Voltage-Gated 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.  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. 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). . 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.  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 MolDockrerank score best interacting compound was selected respective to each parent compound.  The hydrophobic intensities of the binding site ranges from -3.00 (least hydrophobic area -blue shade) to 3.00 (highly hydrophobic area -brown shade).

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  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 4b respectively.

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.