Computational analysis and predictive modeling of polymorph descriptors

Background A computation approach based on integrating high throughput binding affinity comparison and binding descriptor classifications was utilized to establish the correlation among substrate properties and their affinity to Breast Cancer Resistant Protein (BCRP). The uptake rates of Mitoxantrone in the presence of various substrates were evaluated as an in vitro screening index for comparison of their binding affinity to BCRP. The effects of chemical properties of various chemotherapeutics, such as antiviral, antibiotic, calcium channel blockers, anticancer and antifungal agents, on their affinity to BCRP, were evaluated using HEK (human embryonic kidney) cells in which 3 polymorphs, namely 482R (wild type) and two mutants (482G and 482T) of BCRP, have been identified. The quantitative structure activity relationship (QSAR) model was developed using the sequential approaches of Austin Model 1 (AM1), CODESSA program, heuristic method (HM) and multiple linear regression (MLR) to establish the relationship between structural specificity of BCRP substrates and their uptake rates by BCRP polymorphs. Results The BCRP mutations may induce conformational changes as manifested by the altered uptake rates of Mitoxantrone by BCRP in the presence of other competitive binding substrates that have a varying degree of affinities toward BCRP efflux. This study also revealed that the binding affinity of test substrates to each polymorph was affected by varying descriptors, such as constitutional, topological, geometrical, electrostatic, thermodynamic, and quantum chemical descriptors. Conclusion Descriptors involved with the net surface charge and energy level of substrates seem to be the common integral factors for defining binding specificity of selected substrates to BCRP polymorph. The reproducible outcomes and validation process further supported the accuracy of the computational model in assessing the correlation among descriptors involved with substrate affinity to BCRP polymorph. A quantitative computation approach will provide important structural insight into optimal designing of new chemotherapeutic agents with improved pharmacological efficacies.


Background
The computational tools intended for quantitative assessment of protein-ligand interactions are based on several factors including protein-ligand docking, molecular dynamic simulation and free energy calculations [1]. To better define the role of binding affinity in forming a protein-ligand complex, a structural characterization for putative human off-targets was recently performed on Nelfinavir, a potent HIV-protease inhibitor with pleiotropic effects in cancer cells [2]. In this experiment, they have adapted numerous computational models that integrated molecular dynamic simulation, free energy calculations with ligand binding site comparison and biological network analysis.
There are two integral screening approaches that could help identify and characterize the substrates and inhibitors of the efflux proteins and/or transporter system; the measurement of binding affinity and toxicity analysis of substrate compounds [3]. There was a report that drug resident time and uptake amount are better correlated with drug efficacy than the binding affinity [4][5][6], suggesting that lead optimization could be efficiently accomplished with analyzing the drug uptake profiles. Although numerous methodologies have been proposed for drug-target screening strategies based on binding affinity [7,8], there are no efficient computational tools available for the accurate estimation of the drug uptake profiles from the point of the molecular structures. In this study, the uptake rates of Mitoxantrone in the presence of various substrate compounds were examined as an in vitro screening index that could help to characterize the binding properties of chemotherapeutic drugs to tumor cells or efflux proteins.
Breast cancer resistant protein (BCRP) also known as ABCP or MXR or ABCG2 is a member of transporter super family ATP binding cassette (ABC) proteins. BCRP is known to affect the therapeutically available concentrations of various clinical agents [9][10][11]. Since the BCRP effluxes a wide range of structurally diverse xenobiotic compounds from cells [12], the broad distribution of BCRP not only renders less complete distribution of drugs but also causes a poor response of cells to chemotherapeutics [13][14][15]. BCRP in conjunction with P-gp expression at target sites affected the pharmacokinetic profiles of substrates and inhibitors [16]. Subsequently, the therapeutically available concentrations of certain agents increased in BCRP knock-out animal models that were highly prone to Mitoxantrone induced toxicity [17].
The in vitro studies on the BCRP efflux system have demonstrated that some cell lines displayed erratic efflux profiles of doxorubicin and rhodamine 123, and these observations were attributable to 482 nd position in the amino acid sequence consisting of arginine, glycine or threonine residues which are susceptible to numerous posttranslational modifications [18][19][20]. Three polymorphs, namely 482R (wild type) and two mutants (482G and 482T) of BCRP, have been identified, and alterations in their expressions and functions were reported [21]. Wild-type BCRP and its variants were markedly expressed in human embryonic kidney (HEK) cells [22].
The present study was intended to establish the relationships between chemical properties involved with the uptake rates of structurally diverse substrates and BCRP polymorphs. To achieve this goal, we have designed the computational model consisting of numerous molecular descriptors. The uptake rates of Mitoxantrone by BCRP were examined in the presence of various pharmacological classes of ABC transporter inhibitors, such as antiviral (i.e. Erythromycin, Foscarnet), antibiotic (i.e. Ciprofloxacin, Febendazole, Novobiocin, Quercitin), calcium channel blockers (i.e. Verapamil, Diltiazem, Nifedipine, Qunidine), anticancer (i.e. Mitroxantrone, Acyclovir, FTC, Phenethyl ITC, Raloxifene, Rodamin 123, Saquinavir, Tamoxifene), antifungal agents (i.e. Ketoconazole), hormones (i.e., Estradiol) and immunosuppressant (Cyclosporin) [16,23]. It was hypothesized that any changes in uptake rates of Mitoxantrone are due to competitive binding of these substrates to BCRP.
In the development of a computational model for prediction of structural specificity of substrate compounds to BCRP, three dimensional structures of the substrates were built using AMPAC with Graphical User Interface (Semichem, Shawnee Mission, KS). AMPAC used Austin Model 1 (AM1) for the quantum mechanical semiempirical calculations of interactive energy. CODESSA can generate the numerical values for molecular descriptors, whereas the heuristic method (HM) preselects appropriate molecular descriptors. The multiple linear regression (MLR) is capable of deriving the linear QSAR based on them. The final outcomes were labeled as a characterization of compounds using derived properties X from AM1 calculations and regression with MLR, using measurements Y as response. The knowledge on such descriptors that determine substrate specificity to binding receptors is critical to delineate the drug interaction with BCRP polymorphs and the mechanisms behind their action. The outcomes of this study ultimately lead us to discover efficient new drugs with enhanced chemotherapeutic efficacies.

Results and discussion
The uptake rates of Mitoxantrone by HEK cells The effects of various substrates on the uptake rates of Mitoxantrone by HEK cells were evaluated to determine their binding capacity to BCRP polymorphs. The uptake rates of Mitoxantrone (expressed per mg of protein) in the presence of various substrates were converted to the percentage uptake rate of Mitoxantrone in the absence of the substrates ( Figure 1). The transcellular permeation profiles of the substrate compounds showed a similar trend to those of the uptake profiles, but statistical significance of the latter is much greater than the former. There are several important findings from this study. It was suggested that changes in the uptake amount of Mitoxantrone in the presence of various substrates are due to their influence toward BCRP efflux polymorphs. Although an indirect approach associated with the uptake rate may not be the best option to predict the relationship with the chemical structures, the changes in the uptake rate of Mitoxantrone in the presence of substrate compounds could serve as a valid indicator for the drug affinity to BCRP.
Relationships between the uptake rates and chemical properties I. 482R Polymorphs The relationships between the uptake rates and the chemical structures of substrates were analyzed and quantitatively expressed as the concentration of a substrate required to exert a biological response. As shown in Table 1, four descriptors, HOMO-1 energy, Max atomic orbital electronic population, ESP-Max net atomic charge and average electrophilic reaction index for a O atom, mainly contribute to the linear relationship profiles in a QSAR model for 482R polymorph with the experimental coefficient value R 2 of 0.9740 (F = 56.17 s 2 = 0.024 Q 2 = 0.8561: Table 2), which are indicative of a close correlation among them. As shown in Figure 2, the prediction power of the test substrates was also within 5% of the experimental value, which further corroborated the proposed model.
It was noted that the dynamics of the aromatic cores and the alkyl tails can affect the electronic properties.
The highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) are just two of molecular energy levels, and efficient energy conversion is an integral step in practical applications of substrates for charge transport and transfer processes through efflux proteins or transporters.
Another important descriptor seems to be overall net charge and surface charge of substrates. As the atomic charge of the molecule increases, the modulation of the uptake rate of Mitoxantrone decreases (X-axis is expressed as a negative absolute value), indicating that the charge in the molecules were inversely correlated with the uptake rates by BCRP 482R polymorphs. This finding is in a good agreement with the previous report that ABC transporters have a greater affinity to positively charged molecules which serve as an electron accepting functional group [27].
Phenethyl Isothiocyanate (PEITC), which showed a higher affinity to 482R, has an electropositive property originated from thiocyanate functional groups which may contribute to its strong affinity to 482R. On the other hand, Substrates, such as caffeine, estradiol and verapamil didn't significantly affect the uptake rate of Mitoxantrone by 482R polymorphs.

II. 482G Polymorphs
The results of the study with 482G polymorph were most properly expressed in a linear QSAR equation consist of 4 descriptors, which are ESP-Max net atomic charge, Max SIGMA-PI bond order, Max 1-electron reaction index for a C atom, and ESP-FPSA-1 Fractional  Tables 1 and 2. The predicted value for the substrate molecules was also within 5% of the experimental value as shown in Figure 3. Similar to the results of 482R polymorph, it was also demonstrated that affinity of the 482G polymorph to BCRP decreases, as the charge on the molecule increases. The correlation coefficient between two variables (i.e., charge and binding affinity) for 482G polymorph is less significant than that of 482R. Among substrates, Phenethyl Isothiocyanate (PEITC) has the highest binding affinity to 482G. The binding affinity of Rhodamine123 to 482G polymorph is lower than that to 482R. Febendazole and Riboflavin showed a low binding affinity to 482G, mainly due to the presence of multiple double bonds and aromatic rings. Substrates, such as Rhodamine123 (123%), Estradiol (114%) and Foscarnet (111%), did not significantly affect the Mitoxantrone uptake rate by 482R polymorph.

III. 482T Polymorph
As shown in Table 1 and Min e-n attraction for a C-C bond. The correlation value for the modulated uptake rate was 0.8268 (F = 71.6, s2 = 0.027, Q 2 = 0.5617) and the predicted power of the test molecules was within 5% of the experimental outcomes, indicating that there is a good correlation among selected parameters. As shown in Figure 3, differing from the results of 482R and 482G polymorph, there is a positive relationship between surface charge and binding affinity; as a charge on the molecule increases, its binding affinity to 482 T polymorph increases. The results of this study suggest that the presence of the charged residue significantly affects the affinity of test substrates even though it does not substantially contribute to the specificity to each BCRP polymorph. The steric factors are likely to play a vital role in the relationship between substrate property and binding affinity, even though their contribution to the binding affinity of substrates to BCRP is much less than the charged residue.
482T polymorph has a lower affinity towards substrates, such as Caffeine, Diltiazem, Raloxifene, Quinidine and Verapamil. Fumitremorgin C (FTC) showed the highest impact on the uptake rate of Mitoxantrone by 482T, which is probably due to the presence of the charge species on FTC. Since the side effects of FTC (i.e., neurotoxicity) arises from stereo chemical constraints on the conformation of the diketopiperazine D ring, the replacement of the proline moiety (E ring) by an acyclic substituent might allow the adjacent diketopiperazine ring to assume a new conformation with the less charge species that renders the diastereoisomeric mixtures of FTC analogues less neurotoxic than native FTC [28,29].
Quercitin, which is known as the most active reactive oxygen species (i.e., peroxynitrite) scavenger among the structural analogues of flavonoids, had a significant impact on the uptake rate of Mitoxantrone in all three polymorphs, probably due to the fact that Quercitin has high resonance and donates electrons on the oxygen atom even though it lacks nitrogen atoms on the aromatic ring. It was reported that molecule oxygen atom on C-4 of C ring of Quercetin carries the largest excess charge, whereas charge accumulation on the hydroxy groups at the same ring is not considerably large [30,31]. On the other hand, estradiol did not affect the uptake rate of Mitoxantrone in all three polymorphs, indicating that estradiol is not a major substrate for BCRP.
The results of this study underline importance of complement regulatory proteins in the biologic systems that outline the binding capacity of exogenous compounds and subsequent their uptake rates.

IV. The model validation process
The Cross Validation process was carried out to confirm the predicting power of the QSAR model.   The results of this study suggested that the binding affinity of substrates to specific receptors or efflux proteins is not entirely dependent on a particular variable or an individual group of variables in the relationship, but it is rather affected by various combinations of variables. This finding further supported the robustness of the QSAR approach in predicting the outcome from the medical database by avoiding a spurious association within a set of variables. It was also suggested that descriptors involved with the drug uptake profiles will give new insights on chemical modifications that can lead to designing new chemotherapeutic agent with improved pharmacological properties.
To design a pharmacologically active drug for sitespecific activity is a challenging task that begins with rationally identifying the targets to which that drug binds. There are a number of computational approaches in designing efficient therapeutics based on target identification and lead optimization. As the specific binding to active targets may have a profound impact on the overall pharmacological activity [32], the effects of a variety pattern of protein binding reflected in the uptake profiles by cell membrane on therapeutic efficacy of drugs could be adapted as a primary screening means. The importance of protein binding has already been validated by less specific protein kinase inhibitors which attack tumors through multiple mechanisms [33]. This strategy has been effective to more than one type of cancer therapies.
A thorough understanding of drug binding interactions and their relationship with the biological activity requires high-throughput computational biology approaches. Computational techniques that identify competitive binding substrates and their inhibition range in cellular networks have been intensively developed, but their scales are very limited to initial assessment process [6,34]. Moreover, the qualitative description of the chemical entity currently available showed a limited predictive power due to the high dynamic nature of molecular structures and complicated responses from biological systems including complex efflux pathways. The mathematical modeling approaches, such as ordinary differential equations and pi-calculus, have limitations in that they require a large number of kinetics parameters to simulate the dynamic behavior of the biological system to which chemical entities bind [35,36]. Therefore, a functional dynamic model based on the qualitative descriptors defined from the competitive uptake profiles is integral for parameter optimization and dose regimen specification of new drug entities.
In recent years, major efforts have been placed on identifying and characterizing ABC transporters. They are expressed at the major barriers within the body (e.g., intestine, blood-brain barrier, placenta, kidney, and liver), where they lowered the uptake rate or enhanced the clearance of drugs [37]. Breast Cancer Resistant Protein (BCRP) is one of the most recently discovered members of ABC transporters. BCRP is a homo-dimer and consists of 655 amino acids containing a nuclear binding domain and a membrane spanning domain. BCRP shares broad similarities with bacterial, yeast, insect and other mammalian ABC transporter proteins [38]. In normal human tissues, BCRP was detected at higher levels in the placenta and at lower levels in the brain, prostate, ovary, colon, testis, liver, small intestine, kidney and heart [22,39]. Among normal tissues, BCRP is expressed in sycytiotrophoblasts of placenta, epithelium of small intestine, colon, liver, ducts and lobules of breast, and haemopoietic stem cells [22,40].
ABC transporters comprise various efflux proteins, some of which exert abnormal responses to exogenous compounds due to the presence of polymorphs. Nine polymorphs have been identified for MDR-1 [41]. In a clinical trial, ABCG2 polymorphism has a vital role in delineating the effective dose in chemotherapy [42]. Genetic polymorphism in ABC transporters influences numerous diseases including hypertension [43], lung cancer [44] and colon cancer [45]. The mutation of polymorphs discovered from human BCRP serves an integral criterion for the differential affinity of substrates to BCRP [46]. It was found that the affinity of BCRP to substrates can be modulated by altering the substrate specificity of multi-drug transporters. Several BCRP variants from direct DNA sequencing of the BCRP gene have been reported [47]. It was also demonstrated that single nucleotide polymorphisms (SNP) of BCRP produced individual variations in the pharmacokinetics and toxicity profiles of BCRP substrates. BCRP G34A (Val12-Met) and C421A (Gln141Lys) polymorphisms occurred at high frequency in most ethnic populations and have been associated with the expression and activity of BCRP protein [48]. It has distinctive features including racial differences; for instance, BCRP V12M, Q141K, P269S and Q126Stop were detected in Korean at frequencies of 23, 28, 0.2 and 1.9%, respectively [49].
This study was undertaken to define various physicodynamic and chemical properties of substrates to BCRP polymorphs and elucidate the rationales behind their efficacies. The structurally diverse compounds were evaluated for elucidation of the BCRP polymorph mediated uptakes and establishment of the relationships with their molecular structures. The results of this study suggested that the chemical properties of exogenous compounds significantly influence BCRP polymorphs mediated uptake rate. All the compounds tested in this study are either substrates or inhibitors of at least one of BCRP polymorphs. The analysis on the chemical properties of the substrates based on the combining AMPAC/ CPDESSA approach could help us to identify integral descriptors that should be mirrored by interactions with receptor proteins [50].
In this study, the consistent appearance of surface charge, electrophilic reactivity indices and molecular orbital energy descriptors obtained from the sets of chemotherapeutic substrate compounds supported the proven concept that charged electrophiles with the high energy level affected the affinity to BCRP polymorphs. The quantum chemical descriptors of the substrates, such as atomic orbital electronic population and bond order, also significantly contribute to its affinity to BCRP polymorph. It is known that the docking analysis of descriptors provides a qualitative representation of ligand and protein interactions in the QSAR model, even though the selection of docked conformations is often complicated due to its sensitivity to the scoring function. The results of this study demonstrate that substrate compounds containing net charged radicals can activate efflux proteins or peptides in the complement system. There is also a close correlation between descriptors and molecular weight, especially for bulky groups. The steric contour analysis indicates that the addition of bulky groups in the active region reduces the binding affinity.
It is possible that subjects with these polymorphisms may have different levels of single nucleotide polymorphisms (SNP) expression level and cellular localization and, consequently, varying degrees of efflux capability to model compounds [51]. Further studies are needed to determine which level and sites of SNP mainly contribute to the specificity of BCRP bindings. The findings in this study provide

Conclusions
In summary, the chemotherapeutic effects of the known substrates were classified based on their binding affinity to BCRP. The computational approach with the sequential approaches of Austin Model 1 (AM1), CODESSA program, heuristic method (HM) and multiple linear regression (MLR) was performed to derive QSAR model and its predictive power was validated. The BCRP mutations may induce conformational changes as manifested by the altered uptake rates of Mitoxantrone by BCRP in the presence of other competitive binding substrates that have a varying degree of affinities toward BCRP efflux. At the practical level, the use of a computational The data are expressed as mean +/− SD, p < 0.05, each experiment performed in quadruplicate.
structural approach will help the scientists identify the best compound and its linear dose range with improved pharmacological efficacy, eliminating the need to perform multiple assays over a wide range of concentrations in defining the binding affinity and uptake rates.

Methods
HEK (human embryonic kidney) cell lines transfected with each polymorph (i.e., 482R, 482G and 482T) were kindly donated by NIC (NIH, Bethesda, MD) [53]. The minimum essential medium was purchased from  Cell culture preparation HEK cell lines transfected with 482R, 482G and 482T plasmids were grown in the minimum essential medium supplemented with 10% FBS, 50 IU/ml penicillin, 50 μg/ml streptomycin, 4 mM L-glutamine and 100 nM Geneticin. Cells were incubated in 75 mm 2 plastic culture flasks at 37°C supplemented with 5% CO 2 /95% air. BCRP expression in HEK cells was confirmed by RT-PCR using the method previously described [54].
A drug uptake study HEK cells were trypsinized and loaded in 24 well plates at a seeding density of 2 × 10 5 cells per well. Cell viability was maintained by providing a fresh medium every other day until they reached confluence. Cells were exposed with radioactive Mitoxantrone (~100 μM) in the presence or absence of various inhibitors at specific concentrations or 10 μM Fumitremorgin C (FTC) (as a positive control) at 37°C for 5 min. The drug uptake process was stopped by washing the cells with 1 ml of ice-cold DPBS for 3 times, followed by lysis with the Triton X/0.1 M NaOH solution. A cell digest (100 μl) was taken and diluted to 5 ml with 30% Scintisafe ™ (Fisher Scientific, NJ). The cumulative amount of Mitoxantrone in diluted samples was determined using Beckman Coulter Counter and expressed as the percentage amount of the control. Mitoxantrone accumulation was normalized for cellular protein and presented as the percentage of the control, where the control represents cells treated with Mitoxantrone in the absence of any inhibitors. Data were expressed as mean +/− SD, p < 0.05 and by one way ANOVA.
A drug permeability study HEK cell lines were prepared as described previously [54]. After confluence was achieved, the TEER value of cells was measured to verify the presence of tight junction. The growth medium was replaced and washed with PBS. The radioactive Mitoxantrone (~100 μM) in the presence or absence of inhibitors was added on the basolateral side of the transwells. The samples were collected from the apical side at predetermined time intervals for up to 120 min. Apparent drug permeability (Papp) value was calculated using the formula Papp = (dQ/dt)/(A x D 0 ), where (dQ/dt) is the linear appearance rate of drug in the apical side, A is the cross-sectional area of the Transwell insert and D 0 is the initial concentration of the compound in the baso-lateral compartment [55]. The experiment was repeated for 4 times for each inhibitor (N = 4).

The selection of datasets
As shown in Table 4, the substrates selected according to the different pharmacological categories were evaluated for their effects on the flux rates of Mitoxantrone. The experimental sets (25) were divided into 18 training sets and 7 test sets (i.e., for validation purpose). The structures of the substrates in the training set are shown in Figure 5. The uptake rates for Mitoxantrone or those in the presence of the substrates obtained from the experimental transport study were considered to be directly proportional to the affinity and/or permeation rate of drugs with BCRP.

A QSAR approach
This approach was intended to select the most suitable descriptors among the descriptor sets for defining binding specificity of each polymorph. Three dimensional structures of the substrates were built using AMPAC with Graphical User Interface (Semichem, Shawnee Mission, KS) [56]. For structural classification of descriptors, AMPAC used Austin Model 1 (AM1), which is the Hamiltonian widely utilized in the quantum mechanical semi-empirical calculations of interactive energy. AM1 achieved energy minimization of the gradient norm (0.05 kcal/mol) using 20 simplex iterations followed by 1000 steps of Powell minimization [57]. CODESSA (Comprehensive Descriptors for Structural and Statistical Analysis; Semichem, Shawnee Mission, KS) is an advanced, full featured Quantitative Structure/Activity Relationship (QSAR) program that connects information from AMPAC to experimental data [58]. CODESSA preselects each subset of structure descriptors, which include constitutional, topographical, geometrical, electrostatic, thermodynamic and quantum chemical properties, as shown in Table 5. CODESSA can generate the numerical values for up to 600 molecular descriptors which can be used for the regression analysis [59]. The heuristic method (HM) preselected appropriate molecular descriptors and derived the linear QSAR model based on them. Those descriptors that exist for all molecules in the training set were included, whereas those descriptors whose values did not vary throughout the training set were excluded from the regression. The number of descriptors in the final QSAR models was usually less than one third of the number of molecules in the data set [60]. HM allows us to obtain the best QSAR based on F and t test values, which were set in such way that descriptors having more than 0.99 and less than 0.8 correlations were excluded as they may generate an over-optimistic regression.
Molecular descriptors selected by the heuristic method (HM) in CODESSA were used as inputs to perform multiple linear regression (MLR), which is the simplest method that builds a single regression equation for a given data set. For each regression analysis, the goodness of fit was evaluated by examining the number of molecules (N), coefficient of determination (r 2 ), cross-validated standard error (Q 2 ) and value of the F-statistic (F). The Q 2 value obtained through the leave-one-out algorithm (LOO) cross-validation procedure reflects the stability of the model through perturbation of the regression coefficients with the acceptability criterion of 0.5 in most CoMFA studies [61].

The model validation process
The linear QSAR model was cross validated using the error values acquired in the prediction process [62]. To optimize the validation outcomes, a relationship between the experimentally obtained values and computed uptake rates from the test set was established using the heuristic method. The expected values for each compound were calculated and plotted to elucidate their correlations with the experimental values.
Each descriptor was assigned a number through the web based random number generator (www.random.org) for cross validation. Individual descriptors were sequentially incorporated into the regression process to monitor the error value. The following two-fold cross-validation scheme was implemented.
1. The experimental data (25) were divided into 2 subsets; 18 training sets and 7 test sets for the model validation. 2. The data from the selected subsets were categorized into the descriptors listed in Table 5. 3. The odd ratio for each descriptor was calculated to find the most potential contributors to binding property with defined risk weights. 4. The significance of the final regression is determined by comparing prediction Absolute Relative Errors (AREs) which are obtained using the absolute value of [(Actual Output -Predicted Output)/Actual Output] to the test subset. An estimated p value less than 0.05 was considered to be significant for the consistency of the model. The model with the lowest prediction error generated through the cross validation process was chosen to represent the best outcome for each BCRP polymorph.