GUO Hong-Mei FU Le LI Guang-Ping SHU Mao LIN Zhi-Hua

(School of Pharmacy and Bioengineering, Chongqing University of Technology, Chongqing 400054, China)

ABSTRACT Acetaldehyde dehydrogenase 1A1 is a hopeful therapeutic target to ovarian cancer. In this present work, 3D-QSAR, molecular docking and molecular dynamics (MD) simulations were implemented on a series of quinoline-based ALDH1A1 inhibitors to investigate novel acetaldehyde dehydrogenase 1A1 inhibitors as anticancer adjuvant drugs for ovarian cancer. Two reliable CoMFA (Q2 = 0.583, R2 = 0.967) and CoMSIA (Q2 = 0.640,R2 = 0.977) models of ALDH1A1 inhibitors were established. Novel ALDH1A1 inhibitors were predicted by the 3D-QSAR models. Molecular docking reveals important residues for protein-compound interactions, and the results revealed ALDH1A1 inhibitors had stronger electrostatic interaction and binding affinity with key residues of protein, such as Phe171, Val174 and Cys303. Molecular dynamics simulations further verified the results of molecular docking. The above information provided significant guidance for the design of novel ALDH1A1 inhibitors.

Keywords: 3D-QSAR, molecular docking, molecular dynamics simulation, ALDH1A1 inhibitors;

1 INTRODUCTION

Human aldehyde dehydrogenase (ALDH) gene family encodes more than 10 genes, and so far there have been 19 isozymes[1]. The main functions of the ALDH family is to oxidize endogenous aldehydes produced by various cellular processes of the corresponding carbonic acids[2,3]. Mutations in ALDH genes and subsequent inborn errors in aldehyde metabolism correlate with several diseases, including Sjögren-Larsson syndrome (SLS)[4], type II hyperprolinemia[5],γhydroxybutyric aciduria[6]and pyridoxine-dependent seizures[7]. ALDH enzymes also play important roles in hyperammonemia and alcohol-related diseases, late-onset alzheimer’s disease (AD) and cancer[8-11]. The overexpression of some ALDH (notably ALDH1A) plays an important role in many malignant tumors and cancer stem cells (CSCs)[12]. In ALDH1 substructure family, it mainly has three subtypes,ALDH1A1, ALDH1A2 and ALDH1A3 namely and they are markers for some normal tissue stem cells (SC) and tumor stem cells (CSC)[13,14]. In recent years, some experiments on human cells and mice show that ALDH1A1 is not only a biomarker of cancer, but also related to drug resistance to traditional cancer chemotherapy[2,13].

Ovarian cancer is common malignancy in women and is the leading cause of death from gynecological cancers[15]. Recent studies in ovarian cancer models demonstrated an interesting relationship between ALDH1A1 and cancers cell differentiation[16,17]. Given the significant physiological and pathological roles of ALDH1A1 in ovarian cancer, small molecule inhibitors can effectively act on this enzyme to discover its role in ovarian cancer. However, there are no medically available ALDH1A1 inhibitors. So far, the skeletons of some ALDH1A1 inhibitors have been reported,including indole-2,3-dione-based analogs (Fig. 1) and tricyclic pyrimidinone 2 (CM037, Fig. 1)[18,19]. Most recently,Yang and co-workers found two compounds with similar structures, which both contained a dual rings core of two adjacent arms (NCT-501 and qHTS hit, Fig. 1)[20,21]. In view of the high activity of the two compounds in experiment, they envisaged the possibility of forming a new hybrid series as an example, and described a newly designed series of ALDH1A1 inhibitors based on quinoline-based analogs[22]. After extensive reading of the literature, no literature is found on the corresponding structure and activity of quinoline-based ALDH1A1 inhibitor compounds. Therefore, based on the structure and activity of these reported quinoline-based ALDH1A1 inhibitors, we systematically studied the structureactivity relationship of inhibitors using 3D-QSAR, molecular docking and molecular dynamics simulation and designed several novel ALDH1A1 inhibitors with better inhibitory activity.

Fig. 1. Representative small molecule ALDH1A1 inhibitors, quinoline-based qHTS, and newly discovered quinoline-based hybridization inhibitors

2 MATERIALS AND METHODS

2. 1 Data sources

In this article, the data of ALDH1A1 inhibitors with quinolone-based derivatives were collected from literatures[23],IC50 has been converted to logarithmicpIC50 (−lgIC50). We listed quinolone-based derivatives’ structures and activities(pIC50) in Table 1. The 3Dstructures were obtained from SYBYL-X 2.0 software. To build CoMFA/CoMSIA model,the dataset was divided into a training set (N= 41) and the rest compounds were utilized as a test set marked with “*”.

Table 1. Structures and Respective Experimental pIC50 of ALDH1A1 Inhibitors

10 6-F 7.066 7.620 12 11 6-F 7.569 13 6-F 7.602 14 6-F 7.824 15 6-F 7.886 16 6-F 8.097 17 6-F 7.854 18*6-F 8.222 19 6-F 8.222 20 6-F 7.921 21*6-F 8.155 22 6-F 7.921 23 6-F 8.046 24 6-F 8.155 25 6-F 7.097 26 6-F 8.000 27 6-F 7.959 28 6-F 7.215 29*6-F 7.420 30 6-F 7.959 31 6-F 8.155 32*6-F 8.046 33 6-F 7.796 34 6-F 8.155 35 6-F 7.854 36 6-F 7.854 37 6-F 6-Cl 7.824 To be continued

“*” means test set

2. 2 Methods

2. 2. 1 Molecular optimization and alignment

With Gasteiger-Huckel charge, Tripos force field and Powll energy gradient method, all molecules were optimized in SYBYL-X 2.0 software. All parameters were default except the maximum optimization limit and convergence criterion[24,25].They were set as 10,000 times and 0.005 kcal/mols,respectively. We used the structure obtained by the above method as a subsequent 3D-QSAR analysis.

Molecular alignment is considered to be an important step in the establishment of a 3D-QSAR model. We selected compound 48 (pIC50 = 8.301, Fig. 2) as the template, which has the highest activity. After selecting the common Skeleton(Fig. 2), the superimposed structures of aligned compounds are shown in Fig. 3.

Fig. 2. Structure of No.48

Fig. 3. Alignment of all molecules

2. 2. 2 3D-QSAR model

SYBYL-X 2.0 software was used to establish CoMFA and CoMSIA models. In partial least squares (PLS) analysis, the CoMFA and CoMSIA descriptors are used as independent variables, while thepIC50 value is used as a dependent variable for the development of a 3D-QSAR model. The correlation coefficient (Q2) and the best principal component value (N) of cross validation were determined by leave one method (LOO) for cross validation. We performed non-crossvalidation using previously acquired N values to estimate the general determination factor (R2). In addition, the estimated standard error (SEE) and Fischer statistical values (F) were determined[26,27].

2. 2. 3 External validation

To make the established QSAR model more responsible,external verification is an indispensable step. The method includes Golbraikh-Tropsha’s method and Roy’s method[28-31].The 3D-QSAR model with credible external validation capabilities must satisfy the following criteria:

2. 2. 3. 1 Golbraikh-Tropsha method

2. 2. 3. 2 Roy method

2. 3 Molecular docking

Molecular docking was utilized to investigate interactions between quinoline-based analogs and 5TEI. Surflex-Dock binds ligands to the docking pockets of the receptor protein scored accurately and quickly, thus surveying molecular docking. The crystal structure of protein was downloaded from protein data bank (PDB ID: 5TEI)(http://www.rcsb.org/)[32]. Docking simulations were carried out using a standard Surflex-Dock protocol with default values for adjustable parameters. The binding pose with the top inhibitory activity was selected and the corresponding complexes were output for subsequent MD simulations.

2. 4 ADME/T property prediction and synthetic availability prediction

ADME/T, including absorption, distribution, metabolism,elimination and toxicity, is a critical parameter that is commonly used in clinical trials and the selection of the development of drugs[33,34]. These properties are used to assess the oral bioavailability of five novel radon-based derivatives.

We used the online tool SwissADME(http://www.swissadme.ch/index.php) to evaluate the synthetic availability of new compounds[35]. The synthetic accessibility difficulty scale was 1～10 and the smaller the score, the simpler the synthetic route of the compound.

2. 5 Molecular dynamics (MD) simulation

MD simulation was performed using the AMBER16 package with reference to the official tutorial[36,37]. We obtained the 3Dstructure of 5TEI and ligand complexes output from molecular docking. With the Amber ff99SB force field[38]for receptor and general Amber force field (GAFF)[39]for ligands, the complexes were optimized. All protein inhibitor complex systems were immersed in TIP3PBOX(Buffer ≥ 10.0 Å). The systems were then electrically neutral by adding Na+or Cl-counter ions. The other systems used the same protocol settings.

We optimized discord on the system to eliminate potential space crashes. The system is then gradually heated from 0 to 300 K for the heating phase. Finally, the temperature is kept at 300 K in the next phase. A time step of 2 fs was employed for the entire MD process. Periodic boundary conditions are used to maintain constant temperature and pressure. Using the Langevin dynamics method to adjust the temperature at a collision frequency of 2ps～1s and set the pressure on 1 atm under each anisotropic pressure scale protocol. The particle mesh Ewald (PME) method was employed to deal with longrange electrostatics and the cutoff value range of the realspace interactions was less than 1 nm. The SHAKE method was used to constrain all covalent bonds involving hydrogen atoms. Subsequently, the system performed 100 ns MD simulation and saved the trajectory of the simulated system every 2 ps.

We used MM/GBSA algorithms to process the saved MD simulation trajectories and calculate the binding energy with crystal complexes of different ligands[40,41]. A total of 1000 frames were extracted from the last 90 to 100 ns for calculating the average binding energy. The formulas are as follows:

ΔGBind=Gcomplex –(Gprotein+Gligand)=ΔGsol+ΔGgas

ΔGsol=ΔEGB+ΔESURF;ΔGgas=ΔEele+ΔEvdw

Among them, ΔGBindis a combination of free energy,Gcomplex,GproteinandGligandare related free energy, and ΔGsolrepresents the sum of molecular mechanical energy in the vacuum, and can be further divided into electrostatic contributions (ΔEele) and Van der Wael (ΔEvdw). The term can be calculated using molecular mechanics. ΔESOlis a solventbased energy, including a polar solvent-based energy (ΔEGB)calculated by a generalized natural (GB) approximate model,and a non-polar portion (ΔESURF) overlaps (LCPO) model calculated by fitting the solvent to the surface area (SASA)and two linear combinations. In addition, the energy of each residue is broken down into main and side-chain atoms.Energy decomposition can be analyzed to determine the contribution to key residues to binding.

3 RESULTS AND DISCUSSION

3. 1 3D-QSAR results and analysis

It can be found in Table 2 that each parameter of the COMFA model shows that this model has high reliability. The value ofQ2is 0.583 andR2is 0.967, as well asSEEis 0.059.The ratios of stereo field and electrostatic site are 62.9% and 37.1%, respectively, suggesting that the stereo and electrstatic fields both are indispensable.

Analysis of CoMSIA's 5 fields produced the best model.The predicted and experimental activity values of compounds in the CoMSIA model are listed in Table 2. Among them, our comprehensive analysis selected the best model andQ2=0.640,R2= 0.977, andSEE= 0.052 in the CoMSIA model.

We can see from Table 2 that it has higherQ2than the model and theirQ2is 0.685, 0.731, 0.657, 0.648, 0.676 and 0.708, respectively. Although these models have higherQ2than the template, we believe that they overlooked the influence of hydrogen bonds and hydrophobic fields which were also deemed to be crucial for the activity of the inhibitors. This is why we chose the COMSIA-SEHA model.

In CoMSIA-SEHA model,Q2is not large enough to be standard, but it is an essential condition for the QSAR model to have high predictive ability[27]. These data show a good correlation between experimental and predictive values. The ratios of steric field, electrostatic field, hydrophobic field and hydrogen bond acceptor field are 13.0%, 37.5%, 26.8% and 22.8%, respectively. Those values are listed in Table 2.

3. 2 3D-QSAR results and analysis

It can be found in Table 2 that each parameter of the COMFA model shows high reliability for this model. The value ofQ2is 0.583 andR2is 0.967, as well as SEE is 0.059.The ratios of stereo field and electrostatic site are 62.9% and 37.1%, respectively, suggesting that both the stereo and electrostatic fields are indispensable.

Analysis of CoMSIA's 5 fields produced the best model.The predicted and experimental activity values of compounds in the COMSIA model are given in Table 2. Among them, our comprehensive analysis selected the best model andQ2=0.640,R2= 0.977 andSEE= 0.052 in the CoMSIA model.

We can see from Table 2 that it has higherQ2than the model and theirQ2is 0.685, 0.731, 0.657, 0.648, 0.676 and 0.708, respectively. Although these models have higherQ2than the template, we believe that these models overlooked the influence of hydrogen bonds and hydrophobic fields which were also deemed to be crucial for the activity of the inhibitors. This is why we chose the COMSIA-SEHA model.

In CoMSIA-SEHA model,Q2is not large enough to be standard, but it is an essential condition for the QSAR model to have high predictive ability[27]. These data show a good correlation between experimental and predictive values. The ratios of steric field, electrostatic field, hydrophobic field and hydrogen bond acceptor field are 13.0%, 37.5%, 26.8% and 22.8%, respectively. Those values are summarized in Table 2.

Table 2. CoMFA and CoMSIA’ Statistic Results

3. 3 Analysis of the external validation

The external validation results are shown in Table 3.Comparing the parameters of Golbraikh-Tropsha and Roy method, the 3D-QSAR model is responsible and has good statistical significance.

Table 3. External Validation Method of CoMFA and CoMSIA Models

3. 4 Analysis of the 3D-QSAR model

WhenQ2> 0.5 andR2> 0.9, the model is considered to have reliable predictive power.Q2in CoMFA-SE and CoMSIA-SEHA are 0.583 and 0.640, respectively, andR2in CoMFA-SE and CoMSIA-SEHA are 0.967 and 0.977,respectively. In Fig. 4, the actual and predicted values of all compounds were near the trend line, which indicates the reliability of the model. Combined with the statistical parameters ofPLS, it further shows that the 3D-QSAR model is well predicted and statistically stable.

Fig. 4. pIC50 of predicted versus actual activity of training set and test set. (a) CoMFA (b) CoMSIA

Fig. 5 shows the contour maps of the CoMFA-SE model. It is well known that in stereoscopic fields, the activity of a compound increases in the accumulation of the base group,while the yellow profile means adverse. In Fig. 5a, there is a larger volume of green equipotential region at the position of the small ring of group R1, which indicates that the volume of the group can be increased here, which will be beneficial to increase the activity. In the electrostatic field, the red contour indicates that it is advantageous to increase the negative charge group, while the blue contour means that the increase in the positive charge group is beneficial. In Fig. 5b, there are a certain number of blue and other potential regions in the R1and R2groups, showing that increasing the positive power group will help increase the compound’ activity. The CoMFA model suggests that there should be a bulky group of the R1substituent, and at the same time, the R2substituent should be a positively charged group.

Fig. 5. Three-dimensional equipotential maps of CoMFA with No.48. (a) Stereo and (b) Electrostatic

Fig. 6 depicts the stereo field, electrostatic field, hydrophobic field, and hydrogen bond receptor field of the CoMSIA-SEHA model in turn. The contour maps of CoMSIA are consistent with the CoMFA model and will not be repeated here. In hydrophobic fields, a yellow profile indicates that an increase in hydrophobic groups contributes to an increase in activity, while white means an increase in hydrophilic groups contributes to an increase in activity. In Fig. 6c, the R1and Ra substituents have a certain white isopotential region,indicating that these locations should increase the hydrophilic

group. In Fig. 6d, a purple profile suggests that increasing the hydrogen bond receptor will help increase the activity, while red is the opposite. As can be seen from Fig. 6d, there is a red outline of the position of the R2group, which shows that an atomic or group that can accept the hydrogen bond should be added here. At the same time, the contribution field of the model (Table 2) indicates that the electrostatic field of this model is more important than the stereo one, and the roles of hydrophobic and hydrogen bond acceptor fields cannot be ignored.

Table 4. Errors and Prediction pIC50 in the Training and Test Sets of the 3D-QSAR Model

Fig. 6. Contour maps of CoMSIA with No.48. (a) Stereoscopic (b) Electrostatic (c) Hydrophobic and (d) Hydrogen bond acceptor

3. 5 Results and analysis of molecular docking

Fig. 7 represents the overlap between the co-crystallined molecule CM039 of the protein and the re-docking conformation and the RMSD of the two conformations is 1.32 Å.RMSD value less than 2 Å indicated the reliability of the docking protocol[42].

Fig. 7. 5TEI co-crystal ligand CM039 and re-docking conformation.CM039 shows blue and re-docked conformation as green

Fig. 8a shows 50 active sites where the molecules dock with the 5TEI crystal structure. All compounds are well folded on the active pockets of the receptors, while inhibitor molecules bind together with the receptors. This provides a structural basis of such small molecular compounds to inhibit 5TEI receptors. We selected No.48 as a template for the interaction between the inhibitor and 5TEI according to Table 1.

The docking interactions between No.48 and 5TEI are shown in Fig. 8. As shown in Fig. 8b, No.48 was surrounded by many acid residues, such as Ser121, Asn170, Phe171,Met175, Trp178, Gly294, Tyr297, Cys302, Cys303, Ile304,Tyr457 and Val460, which is in accordance with the former study[32]. Fig. 8c shows the specific interaction between 5TEI and No.48. Therefore, the key residues will be important evaluation criteria in evaluating the interaction between newly designed compounds and target proteins.

Fig. 8. Interaction between compounds and 5TEI Protein. (a) No.1～50 are in the combination pocket, (b) Interaction between compounds with No.48

3. 6 Molecular design and activity prediction of novel derivatives

The model obtained from the QSAR method of docking produces important amino acid residues, provides guidance of the design of new ALDH1A1 inhibitors with more efficient power. According to the information obtained from molecular docking and QSAR model, we selected lowactivity, medium-activity and high-activity compounds as templates to design 5 potential ALDH1A1 inhibitor derivatives. Table 5 lists the structure and predictedpIC50 value of the new ALDH1A1 inhibitor compounds. The results show that the five new compounds designed have betterpIC50 values than their model, of which No.a05 has the highest activity, and the result of molecular docking further improves the reliability of the obtained compound.

Table 5. Structures and Predicted pIC50 of Novel ALDHA1 Inhibitors

Fig. 9 shows No.a05 engages inπ-πinteractions with Tyr297 and van der Waals interactions with Gly125 and Asp122, while the pyridazine ring participates inπ-alkyl interactions with Cys303, Ile304, Phe171, Val460 and Val174. There are two hydrogen bonds: one between the NH base of Cys302 and the vacation oxygen of No.48 (–C=O)with the distance of 2.9 Å, and the other between the NH group of Cys303 and one of the oxygen atoms of the sulfonate group with its distance being 3.1 Å

Fig. 9. Interaction between 5TEI and compound No.a05

3. 7 ADME/T and synthetic availability

We used the DS 3.0 software to predict the novel inhibitors’ ADME/T characteristics, including Human Intestinal Absorption, Aqueous Solubility, Blood Brain Barrier (BBB), Cytochrome P450 2D6 (CYP2D6),Hepatotoxicity and Plasma Protein Binding. The ADME model was developed by using descriptors 2D PSA and AlogP_98 to predict intestinal absorption and blood-brain barrier penetration[43,44]. The descriptors include 95% and 99% confidence ellipses. These ellipses define the areas where compounds that are expected to find good absorption are found. The results showed that the quinoline derivatives a01, a02 and a03 had a 99% confidence in human intestinal absorption and blood-brain barrier (BBB) penetration. The plot of polar surface area and ALog_P for No.48 and its derivatives are represented in Fig. 10.

Fig. 10. Diagram of polar surface area (PSA) and Alog_P98 of 48 and its derivatives

We can see the synthetic availability of new compounds from the online tool SwissADME, and the computational parameters included oral bioavailability (Lipinski’s rule of five) (Tables 6 and 7). All designed compounds are easier to synthesize and conform to Lipinski’s rule of five.

Table 6. No.48 and Its Derivatives Computational Parameters of Pharmacokinetics (ADME)

Table 7. No.48 and Its Derivatives to the Computational Parameters of Oral Bioavailability (Lipinski’s Rule of Five)

3. 8 Analysis of molecular dynamics simulations and free energy

To verify the results of molecular docking and the combination of new inhibitors with 5TEI, MD simulation was carried out, and 100 ns MD simulation was performed on the system of template No.48 by using AMBER software. Fig. 11 shows the RMSD of MD simulation of two complexes, which were No.48 and No.a05. In Fig. 11, the RMSD of 5TEI complex compound 48/a05 fluctuated near ～2.5 and ～3.0 Å after 70 ns, respectively. The RMSD of all complex compounds fluctuated in range less than 2.0 Å after 70 ns.These above statistics shows that the two docking complexes have been effectively combined.

Fig. 11. RMSDs of complex 5TEI with compounds No.48 (a)/a05(b) for 100 ns MD simulations

Fig. 12 shows amino acid residual formed hydrogen bond after the MD simulation. Compared with compound No.48,No.a05 formed two hydrogen bonds with 5TEI and fluorine atom interacts with Ser121 and Asp122 via halogen interaction. We found No.48 interacts Cys303, Phe171, Met175 and Tyr297viaπ-alkyl interaction, while Cyr302 and Trp178 formedπ-sulfur interaction with No.48. Additionally, we observed No.48 engages in hydrogen bond interactions with Cys303, Asn170, Thr129 and Trp178, while No.a05 just engages in Cys302 and Cys303. We agreed that the hydrophobic interaction is the main binding force of the interaction between the compound and the protein, not the hydrogen bond.

Fig. 12. Interaction between 5TEI and No.48(a, b)/No.a05(c, d) after 100ns MD simulation

After MD simulation, MM/GBSA method was used for free energy decomposition to explore the key residues binding to inhibitors in 5TEI. The energy contribution to the residues in active pocket is shown in Fig. 13. Compared with template compound No.48, Novel compound No.a05 had lower binding energy shown in Table 8, indicating that it may have better inhibitory activity, which remained consistent with the results of 3D-QSAR and molecular docking.

Fig. 13. Free energy contribution to 5TEI’ Key residues with No.48(a)/No.a05(b) calculated by MM/GBSA

Table 8. Binding Free Energies and Various Energy Terms (kcal/mol)

4 CONCLUSION

In order to find new and highly effective cancer adjuvant drugs ALDH1A1 inhibitors, based on the CoMFA and CoMSIA methods we established a 3D-QSAR model by using molecular docking and dynamics. The interaction between compounds and proteins was simulated, and the activities of 6 novel ALDH1A1 inhibitors were designed and predicted.From the 3D-QSAR model's contour map analysis, the contribution value of the electrostatic field of the molecule is larger than that of the stereo field, and the effects of hydrophobic and hydrogen bond receptor fields on inhibitor molecular activity cannot be ignored. The introduction to large-volume groups near the R1substituent, or the introduction to positively charged groups near R2replacement bases or the groups capable of accepting hydrogen bonds or the hydrophilic groups near Ra, can help to increase the activity of ALDH1A1 inhibitor compounds. In addition, we also found that the interaction mode of quinoline compounds and proteins is mainly hydrogen bonding interaction. Some amino acids played an important role in combining protein and inhibitor, such as Trp178, Tyr297, Cys302, Val460,etc.The predicted activity values of the newly designed compounds are higher than the respective templates. In addition, these compounds have shown good results of the synthesis feasibility and ADMET evaluation. In summary,through construction, verification and analysis of the 3DQSAR model, combined with the analysis of the mechanism of molecular docking, new ideas and directions for the subsequent development of new cancer-assisted drugs are provided.