Computational Studies of P-Glycoprotein Polymorphisms in Antiepileptic Drug Resistance Mechanisms

The treatment of epilepsy using antiepileptogenic drugs is frequently complicated by drug resistance, leading to drug failure in more than one-third of cases. Human P-glycoprotein (hPGP), coded by MDR1 and belonging to the ABC superfamily, is a membrane efflux transporter that has been identified as an epileptogenic mediator. The ability of hPGP to bind a broad spectrum of substrates could be implicated in the emergence of drug resistance in epilepsy treatment. Single-nucleotide polymorphisms (SNPs) in MDR1 could compound the aberrant changes in hPGP activity causing and enhancing drug resistance. Bioinformatics approaches were used to assess the functional impact of 20 missense MDR1 polymorphisms and of these, five MDR1 polymorphisms were prioritised for further study. The structures of the wildtype and five mutant hPGP were modelled using the mouse PGP structure as the template. Docking studies of the wildtype and mutant PGP with four standard FDA-aprroved anti-epileptic drugs were carried out. Our results revealed that the drug binding site with respect to the wildtype protein was constant. However the hPGP mutant proteins displayed a repertoire of binding sites with stronger binding affinities towards the drug. Our studies indicate that specific polymorphisms in MDR1 could drive conformational changes of PGP structure, facilitating novel contacts with drug-substrates and eventually transporting the drug out of the cell, leading to pharmacoresistance.


Introduction
Epilepsy is a chronic neurological condition affecting more than 50 million people worldwide and 1-2 % of the population. 1 The recurring limitation in the treatment protocol of epilepsy is the failure of drug-response in more than one-third of cases. This has been observed with all the FDA-approved drugs for epilepsy, beginning with phenobarbital, the first drug approved for treating epilepsy in 1912. Of the > 30 FDA-approved drugs for epilepsy, all of them are known to be limited by the problem of drug resistance. Membrane efflux transporters could act to decrease the effective intracellular drug concentration by binding the drug and ejecting it out of the cell. This would indirectly enable pharmacoresistance. PGP (HGNC nomenclature: ABCB1) is an ATP-coupled efflux pump protein well-known as an epileptogenic mediator. 2 It is known to be highly expressed in the blood-brain barrier, which is pharmacologically crucial for the bioavailability of drugs which act on the central nervous system. 3 Thus PGP overexpression might be linked with the low intracellular drug concentration in cortical cells observed in epilepsy treatment. Alternatively, PGP might mediate epileptogenesis via its broad substrate specificity. Undesirable PGP activity might be enhanced by mutations that modulate PGP's polyspecific binding ability. The role of mutations in PGP gain-of-function phenotypes has not been exhaustively studied in the literature. 4,5 The structure of the mouse PGP contains 12 transmembrane helices and two nucleotide-binding domains. 6 The drug-binding pocket is large and hydrophobic, spanning multiple trans-membrane domains. 7 Despite the fact that little information is available on the pharmacological role of the missense polymorphisms of MDR1, clearly the key to epilepsy treatment would involve control of the epileptogenic mediator proteins, including hPGP. Structural studies of hPGP might enhance our current understanding of the role of PGP in drug resistance mechanisms, and provide any evidence of the relationship between specific MDR1 haplotypes and altered drug pharmacokinetics. There is an increasing recognition that there is a genetic component to most acquired epilepsies, in addition to simple Mendelian epilepsy. Here, we have attempted to investigate in silico the effect of polymorphisms on the altered activity of PGP leading to variant-dependent sensitivity to drugs. By conducting a systematic study in this direction, we hope that our results would be useful in throwing light on the PGP-facilitated mechanisms of drug resistance in cancer cells as well. 8

Polymorphism analysis
The hPGP sequence was retrieved from UniProt (acc. no. P08183). A PSI-BLAST search was performed using hPGP as query and target database as Vertebrates, with a E-value of 0.001 until convergence. 9 Multiple alignment of all the hits was performed using ClustalX and manually edited. 10 The dbSNP was to used to identify hPGP SNPs with the search term: "human [orgn] AND missense AND PGP". The hits were assessed for the functional impact of polymorphisms using the curated multiple alignment obtained above. Three different tools were used: SIFT, 11 , PolyPhen2, 12 and PhD-SNP 13 . Consensus of these predictions was used to evaluate the functionally important SNPs.

Homology modeling
The template structures were retrieved using a Blast search of hPGP against the PDB database. 14 ClustalX was used to align the template and hPGP (i.e, target). Modellerv9 was used for modelling and energy-minimisation. 15 For each target, five separate models were generated and the model with the least DOPE score was chosen as the best model. The structure of a mutant protein could be obtained by modeling in the mutation on the wild-type structure, however this would not model any global effects due to the mutation. In order to fully account for the effects of the mutation, we modelled the mutant proteins independently of the wildtype protein.
Molprobity was used to validate the models obtained.

Protein and ligand preparation
Autodock4.2 suite of tools was used for carrying out the docking simulations of hPGP variants and anti-epileptic drugs. 17,18 Hydrogen atoms missing in the protein were added. This was followed by the addition of partial charges to the atoms. The protein was then converted to PDBQT format. The SMILES notation of the drugs of interest were retrieved from Pubchem. 19 The PDB co-ordinates of the drugs of interest were generated from their SMILES representation using OpenBabel. 20 To generate the conformers of each drug, we used MGLtools by calculating the number of bond torsions in the 3D structure. The ligand was then converted to the PDBQT format as well using AutoDock Tools. Target affinity maps for each atom type in the ligand were generated by autorgrid by defining a uniform grid box centered in the hPGP internal cavity. This procedure was repeated for each target-ligand pair, for a total of 6 × 4 = 24 times.

Docking
We employed the Lamarckian genetic algorithm with default parameters for docking search, with 2,500,000 cycles per run, and 10 runs per receptor-ligand pair. The binding mode with the least binding energy was defined as the best pose. The ten poses obtained for each receptor-ligand pair were clustered at 2.0Å r.m.s. to validate the convergence to the best pose. The docked complex was then loaded in PDBQT format, converted to PDB coordinates using OpenBabel, and finally visualized using Rasmol2.7. 21 The differential affinity of the mutant for a given ligand relative to the wildtype was estimated as the difference in the binding energies, i.e. ∆∆G mut = ∆G bind,mut − ∆G bind,wt .

Results and Discussion
Nearly 500 hPGP SNPs were retrieved, however most of these were unannotated, and we obtained a set of 20 hPGP SNPs for further study, none of whose functional effects were known in the literature ( Table 1). The results of our assessment of functional impact by various approaches are summarised in Table 2. Most of the SNPs were determined to be neutral, not disease-causing or deleterious. Five SNPs were predicted to be functionally important by at least one of the tools, as shown in Table 2. Table 3 provides the representative structures of P-glycoprotein in the PDB. Of these homologous hits, the mouse structures cover the full length of the hPGP. Some mouse structures co-crystallised with a ligand might not be representative of the native PGP conformation. Even though 4Q9H is of slightly better resolution, we superimposed 4Q9H with the 3G5U structure and observed that right-half of the 'inverted-V' of 4Q9H was displaced relative to the right-half of 3G5U (shown in Fig. 3), which rendered 4Q9H unsuitable for modeling the full hPGP structure. The alignment between the hPGP and 3G5U is very good, showing > 87.5 % sequence identity between the two sequences and good sequence coverage (shown in Fig. 2). 3G5U was used as the template structure for homology modeling. The target structures of the hPGP wildtype and the five variants were independently modeled and energy-minimised five times each, and the best model was used for further studies (in Table 4). Phenobarbital was the first antiseizure drug used in 1912, followed by phenytoin in 1938. Today > 30 drugs are FDA-approved in the treatment of epilepsy,  yet all of them face pharmacoresistance and more than one-third of epilepsy cases remain untreatable. In addition to phenobarbital and phenytoin, two other common antiepileptic medications, namely valproate and carbamazepine, were the drugs studied here as hPGP ligands (Table 5).
Docking between each hPGP protein (wildtype + 5 mutants) and ligand was carried out. Ten docking runs were performed per receptor-ligand pair. Each run provides one low-energy docked conformation of the respective receptor-ligand pair.
The site corresponding to the lesat-energy binding mode was taken as the binding site of the ligand with the receptor. To ascertain convergence to the lowest-energy binding mode, the ten runs of each receptor-ligand binding conformations were clustered at 2.0Å r.m.s. The lowest-energy binding modes showed good convergence. A sample visualization of the binding site is shown in Fig. 4 and a sample energyhistogram of docked receptor-ligand conformations is shown in Fig. 5. Of the 24 receptor-ligand pairs, 21 had energy-histograms showing the least binding energy (±0.2 kcal/mol) as the single most probable conformation. It was clear that the docking procedure resulted in convergence to the optimum receptor-ligand conformation in most of the cases. The location of the binding site within the internal cavity of PGP was variable across the receptor-ligand pairs. Cutaways of the individual binding pockets representing the best poses are shown in Figs. 6,7,8,9,10. The hPGP residues binding the ligand in each hPGP-drug pair represent the drug-specific binding pocket. These residues were defined at a contact distance of < 4.5Å from the drug in the bound conformation. These residues contributed to stabilizing the docked complex by forming hydrogen bonds and Van der Waals interactions with the substrates. The groups of contacting residues specific to each docked complex are shown in Table 6. The following observations were made: (1) The binding site of wildtype hPGP with all drugs except phenobarbital is invariable. This is an eight-residue binding pocket in the internal cavity lined   (4) Three out of the five hPGP variants exhibited a different binding pocket for each drug, i.e, no two drugs bound in the same pocket for the given hPGP variant. For the other two hPGP variants, three distinct binding pockets were seen. (5) The R593C hPGP mutant bound phenytoin very close to the mutation site (Thr594). This suggested that the polymorphis could perhaps play a direct role in altering phenytoin pharmacokinetics, but we found that the differential ligand affinity was small (∆∆G mut = −0.04 kcal/mol). Further clarity on these observations could be obtained on an examination of the estimated binding free energies of the wildtype and variant hPGPs with the different   Table 7 shows these binding energies along with the predicted differential ligand affinity which is estimated by ∆∆G mut = ∆G bind,mut −∆G bind,wt . The following observations were made: (1) Of the 20 estimated differential ligand affinities, all but three were negative. This implied that in each case the variant hPGP was binding the corresponding drug with a stronger affinity than the wildtype hPGP. (2) The binding energies of phenobarbital with the variant hPGPs were all lower than with the wildtype hPGP. (3) The maximum range of differential response was observed with valproate (−1.25 kcal/mol < ∆∆G mut < 0.58kcal/mol). (4) The least range of differential response was observed with phenytoin. Carbamazepine showed an intermediate range of differential response. The following conclusions pertaining to MDR1 polymorphisms could be drawn: (1) Wildtype hPGP tends to be neutral with respect to binding anti-epileptic drugs. This was shown by: (i) Exposure of a constant binding pocket (ii) Relatively higher binding energy with the drug (2) On the other hand, mutant hPGPs tend to assist in the development of drug resistance. This was shown by: (i) Mutant hPGPs bound each drug in a different location in the internal cavity. Variability in location affords a better search of the optimal binding modes (ii) Consistently lower binding energies imply stable complexes with the drugs which could then be pumped out of the cell coupled with ATP expenditure.
The following conclusions pertaining to the drugs were drawn: (1) The threshold of development of resistance to phenobarbital might be achieved with all the MDR1 SNP alleles studied (2) The patient to patient variability in treatment response and the development of pharmacoresistance was predicted to be maximum in the case of valproate administration (3) The patient to patient variability in treatment response and the development of pharmacoresistance was predicted to be minimum in the case of phenytoin administration (4) Drug-to-drug variability was seen, with preferred binding site for valproate and the absence of a preferred binding site in the case of phenytoin.
The following recommendations for personalized medicine could be drawn: (1) The R593C polymorphism of hPGP showed a more favourable binding profile (negative ∆∆G mut ) with respect to all the four drugs studied, especially valproate. This would result in a gain-of-function phenotype. (2) Similarly, A999T showed a gain-of-function phenotype with respect to all the drugs studied, especially phenobarbital and valproate. (3) E566K showed a gain of function with respect to all drugs carbamazepine. Our studies indicated that carbamazepine might be more effective than other drugs in treating patients with E566K MDR1 polymorphism and suspected PGPmediated drug resistance. (4) Similarly G185V showed a gain-of-function phenotype with respect to pheno- Table 7. Free energy of binding (∆G bind ) of each docked receptor-drug pair. The predicted differential ligand affinity is given by ∆∆Gmut = ∆G bind,mut − ∆G bind,wt . All values in kcal/mol. barbital. (5) Similarly Q1107P showed a gain-of-function phenotype with respect to valproate. Our studies indicated that phenytoin might be more effective than other drugs in treating patients with Q1107P MDR1 polymorphism and suspected PGP-mediated drug resistance.

Conclusions
It was clear that the wild-type hPGP had a preferred binding pocket for most of the anti-epileptic drugs, but this behavior was disrupted by polymorphisms. The polymorphisms played a role in relocating the optimal drug-binding cavity for a higher affinity. We concluded that MDR1 polymorphisms played a key modulatory role in hPGP activity potentially leading to the development of pharmacoresistance in epilepsy treatment. MDR1 polymorphisms were likely to mediate the epileptogenic hPGP phenotype. Gain-of-function process in hPGP is a candidate target in treating drug-resistant epilepsy. The course of treatment with respect to MDR1 polymorphisms would include prescribing drugs with minimal differential ligand affinity. Our study has analysed in silico the impact of MDR1 polymorphisms. Experimental validation of our work would ensure that our findings are applicable to achieving pharmacosensitive response to epilepsy treatment. Furthermore our methodology is extendable to studies investigating the effect of genetic polymorphisms on disease phenotypes.