CD73, as a surface enzyme anchored on the outside of the cell membrane via glycosylphosphatidylinositol (GPI), can convert the AMP in the tumor cell microenvironment into adenosine to promote the growth of tumor cells. It has been overexpressed in many different types of human tumors, such as gastric cancer, pancreatic cancer, liver cancer and other tumor cells. Therefore, targeted inhibitors of CD73 are considered as potential tumor treatment methods. Due to the low bioavailability of nucleoside CD73 inhibitors, it is necessary to develop new inhibitors. In this study, through molecular docking, three-dimensional quantitative structure-activity relationship (3D-QSAR) and molecular dynamics (MD) simulations, a series of CD73 inhibitors were calculated and studied to reveal their structure-activity relationships. Through molecular docking studies, explore the possible mode of interaction between inhibitors and protein. Subsequently, a 3D-QSAR model was established by comparative molecular field analysis (CoMFA) and comparative molecular similarity indices analysis (CoMSIA). For the best CoMFA model, the Q2 and R2 values are 0.708 and 0.983, respectively, while for the best CoMSIA model, the Q2 and R2 values are 0.809 and 0.992, respectively. In addition, the stability of the complex formed by the two inhibitors and CD73 was evaluated by molecular dynamics simulation, and the results are consistent with the results of molecular docking and 3D-QSAR research. Finally, the binding free energy was calculated by the surface area method (MM-GBSA), and the results are consistent with the activities that Van Der Waals and Coulomb contribute the most during the binding process of the molecule to the CD73 protein. In conclusion, our research provides valuable information for the further development of CD73 inhibitors.
Research Article
Computational Investigation of Adenosine 5′- (α, β-methylene)- diphosphate (AMPCP) Derivatives as Ecto-5′-nucleotidase (CD73) Inhibitors by Using 3D-QSAR, Molecular Docking, and Molecular Dynamics Simulations
https://doi.org/10.21203/rs.3.rs-873328/v1
This work is licensed under a CC BY 4.0 License
Journal Publication
published 08 Jan, 2022
You are reading this latest preprint version
CD73
5′- (α
β-methylene)-diphosphate
3D-QSAR
Molecular docking
Molecular dynamics simulations
As a new method, immune checkpoint inhibitors (ICIs) have caused a revolutionary cancer treatment [1]. However, the drug resistance of tumor cells to ICIs is a problem that cannot be ignored. Therefore, developing a new immune checkpoint inhibitor is very important for immunotherapy [2]. CD39 is a cell surface enzyme, which is highly expressed in tumor cells, exhausted T cells and multi-inhibitory cells [3, 4]. It is a key enzyme for the hydrolysis of ATP and ADP into AMP [5]. CD73 is a surface enzyme anchored outside the cell membrane by glycosylphosphatidylinositol (GPI), it can convert AMP into adenosine and expressed in cancer cells, DC cells, Treg (regulatory T cells), NK cells [6, 7]. Studies have found that CD73 is overexpressed in gastric cancer, pancreatic cancer, rectal cancer, and liver cancer cells [8–11]. The sequential expression of CD39/CD73 pathway can eliminate extracellular ATP in the tumor microenvironment (TME) to generate immunosuppressive adenosine and promote tumor growth [12]. In addition, CD73 also acts as a signal and adhesion molecule, which can regulate the interaction between cells and extracellular matrix components (such as laminin and fibronectin) to mediate the invasion and metastasis characteristics of cancer [13] Therefore, targeted therapy that inhibits CD73 is expected to be a new immunotherapy method.
Structurally, CD73 presents an N-terminal domain and a C-terminal domain connected by a long α-helix. The rotation and bending of these two domains provide the open and closed transition of the conformation to bind to the substrate and released [14]. In recent years, based on the structure of adenosine, a series of competitive nucleotide derivative inhibitors have been designed, the most effective of them is AOPCP, which is also one of ideas to design more effective inhibitors [15]. Although these compounds show potential inhibition to CD73, they have negatively charge under physiological pH because they contain phosphorus atoms so that their bioavailability is very low [16]. In fact, stability studies in rat liver microsomes showed that AOPCP was completely degraded within 10 minutes, while the compound based on it, PSB-12489, degraded 5% within 8 hours [17, 18]. After that, the AB680 was designed by changing the purine ring to pyrazolo-pyridine ring has a good tolerance and has entered the phase І of clinical trials [19]. However, recent studies have found that AB680 competes with AMP in vivo is weaker than AMP, indicating that binding to CD73 is reversible and slow, so that AB680 have a direct implication for dosage [20]. Therefore, it is necessary to develop new CD73 inhibitors to change their bioavailability.
This study used computer-aided drug design methods to study a series of compounds with adenosine 5′- (α, β-methylene)-diphosphate derivatives, and these compounds have been tested for CD73 inhibitory activity [21]. In order to study the relationship between structure and activity, based on these 60 CD73 inhibitors, a three-dimensional quantitative structure-activity relationship (3D-QSAR) model with predictive ability was constructed. Through the analysis of the model, the structure characteristics that affect its activity were found. Secondly, through molecular docking studies, the bioactive conformation of the inhibitor was calculated and the receptor-ligand interaction between the compound and the protein was analyzed. In addition, through molecular dynamics (MD) simulations and binding free energy calculations, more realistic conformational changes were analyzed, and the understanding of the binding energy and mode of interaction of inhibitors and CD73 proteins was deeper. In summary, this study provides a basis for clarifying the structure-activity relationship of such CD73 inhibitors and the mode of interaction with proteins, and for designing more active CD73 inhibitors.
Dataset preparation
In this study, a total of 60 adenosine 5′-(α, β-methylene)-diphosphate derivatives have the same skeleton were selected from the literature, the chemical structures of the compounds with their respective activity values are listed in Table 2[21]. The IC50 active values distributed in a gradient of 0.027–18000 nM were converted into the corresponding pIC50 (-logIC50) values to develop the 3D-QSAR models. The 3D structures in the dataset were constructed and energy optimization was carried out using SYBYL-X 2.1.1 by tripos force field [22] with the gradient convergence criterion and max interations of 0.005 kcal/mol and 10000 respectively and Gasteiger-Hükel charges [23]. The compounds in dataset were divided randomly into a training set of 42 compounds (70%) and a test set of 18 compounds (30%) for the generation of the 3D-QSAR model [24].
Molecular Alignment
Molecular alignment plays a key role in 3D-QSAR study [25]. Herein, two alignment methods including ligand-based alignment and docking-based alignment methods were employed to construct the optimal 3D-QSAR model [26, 27]. For docking-based alignment, the conformation of each molecule in the dataset obtained from molecular docking were selected with the best docking score for alignment. For ligand-based alignment, generally, low energy conformation of the highest bioactive molecule was set as template. Thus, Compound 45 as a template molecule for the rest of molecules to align on, the public skeleton was marked in blue bold part (Fig. 1a). In the process of alignment, the eXtended Electron Distribution (XED) force field was used to calculate the field points while other parameters were set by default [28]. A very accurate and slow calculation method was selected in the process of conformational search to improve the alignment effect.
Construction Of Comfa And Comsia Models
CoMFA and CoMSIA models of CD73 inhibitors were founded by SYBYL 6.9.1 software to describe the relationship between the biological activity and 3D structure of CD73 inhibitor compounds [29]. The fields of CoMFA model were calculated including steric (S), and electrostatic (E)[30], meanwhile the CoMSIA model with the above two fields, also calculated it′s hydrophobic (H), hydrogen bond donor (HBD), and hydrogen bond acceptor (HBA) [31]. Other parameters were set in default.
Partial Least Squares (Pls) Analysis
Evaluation of the relationship between structure and biological activity of compounds by PLS method can quantitate the reliability of the model [32]. The analysis was divided into two parts, the training set database were cross-validated by using the leave-one-out method to calculate the cross-validation correlation coefficient (Q2) and the optimum number of components (ONC) [33]. Then, the non-cross-validation correlation coefficient (R2), Fisher’s test value (F value) and the standard error of estimate (SEE) were calculated through non-cross-validation analysis [34]. Among them, high score (Q2 > 0.5, R2 > 0.9, F > 100, SEE༜0.95) models were considered to be predictive[35, 36].
What‘s more, whether the models were accurate, reliable and capable of prediction, the following formula was used to calculate the R2pred value to evaluate CoMFA and CoMSIA.
R=(Sd-press)/sd
In the above formula, SD represents the sum of the square deviation between the biological activity of the test set and the average activity of the training set, while PRESS represents the sum of the square deviation between the experimental value and the predicted value of the training set [37]. After CoMFA and CoMFA models were developed, contour maps of the field type “StDev*Coeff” represented graphically the results [38].
Molecular Docking Study
Molecular docking is a theoretical method to effectively verify and simulate the binding effect of drugs and target proteins through the geometric relationship between small molecular ligands and proteins, energy matching and recognition [39]. In this study, Glide module of Schrödinger suite software was applied to molecular docking and The X-ray crystal structure (PDB code: 6Z9B) of CD73 was obtained from the Protein Data Bank (http://www.rcsb.org/) with resolution of 2.17 Å. Then the docking target protein must be preprocessed, the crystal structure was prepared by Protein Preparation Wizard in Maestro v9.7 to assign bond orders, add hydrogen atoms and build the missing side chain atoms [40]. In order to get a more convincing molecular docking score, Glide was used to generate a grid based on the co-crystal ligand and the grid size was 10×10×10 Å. Then use the extra precision (XP) algorithm in Glide to complete the molecular docking, and use other parameters as default values to further analyze the docking results [41].
Molecular Dynamics Simulation
Molecular dynamics simulation can track the conformational movement of target proteins and other molecules on the atomic level [42]. Therefore, the Desmond (v3.8) module developed by Schrӧdinger software in 2018 was used for molecular dynamics simulation research to further study the specific interaction between molecules and proteins in the binding process [43]. The two complexes formed by the highly active compound 55 and the lower active compound 5 with CD73 were selected for molecular dynamics simulation.
In the process of molecular dynamics simulation, an orthogonal box filled with simple point charge (SPC) water molecules was first created to embed into these complexes [44]. Afterward, in order to reach the physiological salt concentration, a necessary amount of counter ion must be added to neutralize the system. Then, the energy was minimized under the action of the OPLS-2005 force field, the complex system was relaxed to local energy minimization, the maximum number of iterations was set to 5000, and the convergence threshold was set to 1.0kcal/mol/ Å [45].
After energy minimization, the relaxed system was submitted to a simulation time of 100ns with a time step of 2fs, using a Nose-Hoover thermostat at 300K and using Martyna-Tobias-Klein Barostats at 1.01325 bar to simulate the NPT ensemble [46]. In the simulation process, the system energy was calculated every 1.2ps and the trajectory was recorded every 4.8ps. Finally, in order to determine the stability of the two complexes, a series of data of the two complexes were analyzed, including root mean square deviation (RMSD), root mean square fluctuation (RMSF), protein-ligand interaction and potential energy (U)[47].
Binding Free Energy Calculation
The molecular mechanics generalized Born surface area (MM-GBSA) method combines molecular mechanics energy, non-polar solvation conditions, solvation model, non-polar solvent accessible surface area and van der Waals interaction to calculate protein-ligand complex binding [48].
The compound 55 with higher activity and compound 5 with lower activity were selected to form two complexes which have large difference in biological activity with CD73, and their binding free energy was calculated by Prime using molecular mechanics generalized Born surface area (MM-GBSA), which was embedded into Schrӧdinger software, and three formulas were involved in the calculation [49]:
∆Gbind = Gcomplex – (Gprotein + Gligand)
∆Gbind = ∆EMM - T∆S + ∆Gsol
∆Gbind = ∆EvdW + ∆Eele + ∆GGB + ∆GSA - T∆S
In the above formula, ∆Gbind is the binding energy. The Gcomplex, Gcomplex, and Gligand are the corresponding free energy. Where ∆EMM represents the minimized energy differences between the complex structures of protein–ligand and the sum of the energies of the ligand and protein [50]. TΔS was obtained based on the conformational entropy change caused by ligand binding. Added the electrostatic solvation energy ∆GPB(polar contribution), non-electrostatic solvation energy ∆GSA, and ΔEDISPER together to get the free energy of solvation ΔGsol. ΔEele represents electrostatic energy andΔEvdw represents van der Waals energy.
Molecular docking study
Molecular docking is a computational tool in which receptors and ligands recognize each other to form molecular complexes through energy matching, spatial matching and chemical property matching, and predict the structure of the complex. To evaluate the reliability of the XP algorithm in Glide, the co-crystal ligand was redocked into the crystal structure of CD73 (PDB code: 6Z9B). The root mean square deviation (RMSD)between the native conformation and the redocked conformation was 0.483 Å, which indicates that the XP algorithm in Glide was reliable, and it can be considered that this algorithm can be used for subsequent molecular docking studies. As shown in Fig. 1b, it can be seen that the native conformation and the redocked conformation of the natural ligand overlap each other well. Therefore, all inhibitors in the dataset were docked into the active pocket of CD73, and the docking scores were shown in Table S1.
The compounds in our dataset were derived from four types, including partial compounds obtained by substituting AMPCP as the skeleton and three skeletons obtained by substituting its adenine with different nitrogen heterocycles such as: pyrazolopyrimidine, pyrrolopyrimidine, pyrazolopyridine. By analyzing the structural information of the compounds, it was found that they have similar interactions with protein.
For compounds with AMPCP as the skeleton, take compound 5 as an example to analyze the binding mode of this series of compounds and proteins. As we can see from Fig. 2a, it was found that several hydrogen bonds were formed. It can be seen that most of the hydrogen bonds were formed by two phosphate groups, among which HIS38, ARG395 and ARG354 formed hydrogen bonds with oxygen atoms, and ASN117, ASN245 formed hydrogen bonds with oxygen anions. The existence of these hydrogen bonds ensured the stability of the interaction between the compound and the protein. At the same time, HIS118, ARG395, and ARG354 all formed a salt bridge with the negative oxygen atom on the phosphoric acid group. It can also be seen that the hydroxyl group on the tetrahydrofuran ring formed a hydrogen bond with ASP506, and ASN390 formed a hydrogen bond with the N atom on the pyrimidine ring. In addition, there were hydrophobic effects on the imidazopyrimidine ring, residues PHE500, PHE417 work with the two nitrogen heterocycles to form pi-pi stacking interaction.
For pyrazolopyridine compounds, take compound 55 as an example, as shown in Fig. 2c. ASP36, ARG395, HIP118, ASN117, ARG354, ASN245 formed hydrogen bonds with phosphate groups and ASP506 with hydroxyl groups on tetrahydrofuran. These hydrogen bonds increased the stability of binding to proteins. In addition, two salt bridges are formed with ARG354 and ARG395 through the oxygen atom of the phosphate group. PHE500 and PHE417 formed pi-pi stacking interaction with pyrazole ring and pyridine ring to produce hydrophobic effect.
As for pyrrolopyrimidine and pyrazolopyrimidine compounds, the effective inhibitors (compound 54 and 56) were shown in Fig. 2e and 2g. The interaction between compound 54 and protein was similar to that of compound 55, except that the hydrogen bond formed by compound 54 and ASP36 was changed to combine with HIP118. In addition, the salt bridge formed by oxygen atom and HIP38 was added. On the substituent of R1, LEU formed a new hydrogen bond with the N atom, and the benzene ring formed a pi-pi stacking interaction with PHE500. The compound 56 lacked the hydrogen bond of ASN245 on the basis of 54. In addition, it was found that the nitrogen on the pyrazole ring produces hydrogen bond with ASN390, and the rest of the forces were the same.
From the above interaction analysis, it can be seen that these four kinds of compounds with different frameworks had some interaction in common, and their different biological activities and structural characteristics provided data for the 3D-QSAR research. However, it has been reported that compounds with ribose analogs are more potentially active than compounds with β-fluoroarabinose, and the substitution of a diaza heterocycle for the adenine group may have better activity. Therefore, compounds with pyrazolopyrimidines as the core had better inhibitory activity against CD73, and were selected as objects for further research.
CoMFA And CoMSIA Statistical Results
The accuracy of the model was directly related to the alignment methods, so ligand-based methods and docking-based were used to construct models respectively. Both methods were analyzed by PLS, and the following parameters were used to evaluate the quality of the model, including Q2, ONC, SEE, R2, and F value. The calculated results showed in Table 1 and the alignment results represented in Fig. 3. After analyzing the results, it can be found that the model parameters based on the docking were not within the acceptable range while the parameters of ligand-based model the parameters all meet the study requirements which indicated a high model quality, therefore the model constructed by ligand-based method was applied to the following 3D-QSAR study. The detailed data of CoMFA model were: Q2 of 0.708, ONC of 8, SEE of 0.170, R2 of 0.983 and F value of 241.755, the contributions of steric field and electrostatic field were 78.8% and 21.2%, respectively, which meant that the steric field played a crucial role in the CoMFA model. As for CoMSIA model, the Q2 of 0.809, ONC of 10, SEE of 0.121, R2 of 0.992 and F value 383.332 and the contribution of H-bond acceptor, hydrophobic, steric, and electrostatic fields were19.2%, 36.3%, 16.2% and 28.4% respectively, indicating that the hydrophobic and electrostatic fields were the main contribution in CoMSIA model.
The reliability of the CoMFA and CoMSIA models was verified with 18 molecules in the test set, and the statistical parameters R2pred were 0.987 and 0.991, respectively. Table 2 listed the experimental activities and predicted activities generated by the 3D-QSAR model, and the correlation between them were shown in Fig. 4. It can be observed from the graph that there was a good linear relationship between experimental activities and predicted activities. The above data showed that a reliable and high-quality three-dimensional quantitative structure-activity relationship (3D-QSAR) model has been established.
CoMFA And CoMSIA Contour Maps Analysis
The contour map of the QSAR model helped to visually illustrate the relationship between structural features and biological activity. This relationship used the field type stdev*coeff field to determine the results of CoMFA and CoMSIA, in which the contributions of all favorable and unfavorable regions were set as default parameters. In this study, high-activity compound 55 (IC50 = 0.041nM) and low-activity compound 5 (IC50 = 18000nM) were selected as representative inhibitors to explain the structure-activity relationship of these inhibitors.
The CoMFA steric contour maps of compounds 55 and 5 were shown in Fig. 5a and 5b. The green contour lines indicated that bulky groups in the region were conducive to increased activity, while yellow was the opposite. It can be seen that the green part was mainly distributed near the pyrimidine group in compound 5, indicating that the substitution of large groups in this region will increase the activity of the compound, which explained why the activity of compounds 15 and 49 whose hydrogen atoms on pyrimidines are replaced by cyclohexane was better than that of compound 5, the high activity of compound 55 may be due to the large benzene ring substitution in this part. At the same time, it can also be observed that the chlorine atom on the pyrimidine group of compound 5 was covered by the yellow part, indicating that this position can be replaced by small groups as much as possible to have better activity, which can explain why the activity of compound 40 and 44 with substitution of benzene ring and epoxy ring was not satisfactory.
The CoMFA electrostatic profiles of compounds 5 and 55 were shown in Fig. 6a and 6b. The red outline indicates where the negatively charged group can enhance the biological activity, and the blue outline indicated that the positively charged substituent can increase the biological activity. It can be seen that the blue part was mainly distributed in the adenine part of inhibitors 5 and 55, because this part contains multiple N atoms which were negatively charged, so positively charged groups can be considered for substitution on the R1 substituent.
The assumption of the CoMSIA method was similar to that of the CoMFA method, the non-covalent interaction between the ligand and the receptor was determined by the molecular field. The calculation of the molecular field used five different similarity field calculations: steric field, electrostatic field, hydrogen Bond field (including hydrogen bond donor and hydrogen bond acceptor) and hydrophobic field. The steric and electrostatic contour maps were shown in Fig. S1 and S2 The hydrophobic outlines were shown in Fig. 7a and 7b. The white outlines indicated that the introduction of hydrophilic substituents in these portions can increase activity, while the yellow outlines indicate areas that were favorable for hydrophobicity. The R1 on the pyrimidinyl group of the compound 5 was surrounded by a white part, indicating that if this site is substituted by a hydrophilic group, the activity of the compound will increase. At the same time, there were a lot of white parts around the imidazole ring, indicating that this part should have a hydrophilic group. In addition, R2 was surrounded by a yellow outline, indicating that a hydrophobic substituent was required.
In the HBA region shown in Fig. 8a and 8b, the blue-green contour lines indicated increased activity of the hydrogen bond acceptor, while the orange contour lines were the opposite. The blue-green region was close to the R2 and the imidazolopyrimidine ring, indicating that the introduction of HBA chemical substituents at these positions can help improve biological activity.
As shown in Fig. 7, the structural characteristics of CD73 inhibitors to increase the activity were summarized from the CoMFA and CoMSIA models. It can be seen from the structural activity model that we have established that the imidazopyrimidine ring on the molecular skeleton mainly played a hydrophobic role and it needed hydrogen bond receptor. The R1 part has larger groups, hydrophilic and positively charged groups will have greater activity, R2 if there were small groups, hydrophobic receptor groups will enhance the activity, if R3 were hydrophilic groups existence is crucial, and this should be preserved during structural modification.
Molecular Dynamics Simulation
In order to illustrate the stability of the conjugates of compounds 5 and 55 with CD73 and the contribution of key residues in the active pocket, a 100ns molecular dynamics simulation was used to evaluate the complex. As shown in Fig. S3, 19.1% of the residues are distributed in the α-helix, and 20.2% of the residues were distributed in the β-strand, which ensureed the stability of the protein and the low fluctuation of the RMSD value.
For the complex of CD73 and compound 5, the RMSD value of the protein backbone rose to about 2.5 Å in the first 20 ns, and kept fluctuating. Until the end of the 60ns simulation, the RMSD value stabilized at around 2.2 Å (Fig. 10a). In addition, during the entire simulation process, the RMSD value of the ligand configuration fluctuated frequently in the first 20 ns to 1.5 Å and reached equilibrium at about 80 ns (Fig. S4). These results indicated that the complex reached a steady state at the end of the molecular dynamics simulation. As for the complex with compound 55, the RMSD value of the protein backbone rose to close to 2.0 Å within 30 ns, and stabilized to 1.7 Å at 80 ns Fig. 10a. During the whole process, the RMSD value of compound 55 did not change much, and at the same time, the RMSD value of its ligand kept fluctuating around 1.5 Å, indicating good stability during the simulation process.
The changes in the flexibility of the protein were detected by averaging the atoms of all the residues, and the root mean square fluctuation (RMSF) value of the main chain residues was calculated. The graphs of the two complexes were shown in Fig. 10b. In RMSF, most of the residues fluctuated within the range of 2.5 Å, and very few residues exceeded 3.0 Å. In addition, the overall change trend of the amino acid residues of the two compounds was similar, and the changes of some key residues were even completely the same, indicating the two complexes had similar binding mode.
In the MD simulation process, through the monitored receptor-ligand binding mode, a histogram based on the duration and docking results was obtained. The diagram shows a variety of interactions including hydrogen bonds, hydrophobicity, ions and water bridges.
For compound 5 with lower activity(Fig. 11a), because the phosphate group contains multiple oxygen atoms, it can form hydrogen bonds with multiple residues, including hydrogen bonds with ARG395, HIS38, HIS243, HIS118, and ARG354, respectively, accounting for 113%、65%、96%、99%、97%、93%、88%. In addition, the hydroxyl groups on the tetrahydrofuran ring also formed hydrogen bonds with ASP506, accounting for 36%. It can be seen that hydrogen bonds played a vital role in the interaction between the ligand and the receptor.
For compound 55(Fig. 11b), similar to compound 5, most of the hydrogen bonds were formed with the oxygen atom of the phosphoric acid group, such as HIS243、ASN117、HIS118、HIS38、ARG354、ARG395, and account for 45%、33%、81%、73%、100%、98%、57%、43%、91% respectively. In addition, it was different from compound 5 was ASP506 formed a hydrogen bond with the two hydroxyl groups of tetrahydrofuran in a proportion of 75% and 77%. Meanwhile, ASN390 also formed a hydrogen bond with a hydroxyl groups in tetrahydrofuran in a ratio of 86%, and formed a hydrogen bond with the N atom of the pyridine ring in a ratio of 94%. In addition, due to the hydrophobicity of the pyrazole ring, PHE417 formed Pi-Pi stacking interaction with it and accounting for 71%.
During the entire molecular dynamic simulation process, the rotatable bond in the two molecules was also investigated. A summary of the conformational changes of each rotatable bond (RB) in the ligand in the whole simulation trajectory (0.00 to 100.00 ns) through the ligand torsion diagram. It can be seen that compound 5 had 8 rotatable bonds and compound 55 has 10 rotatable bonds Fig. S5 and S6. This was because compound 55 has a larger group substitution on R1 but compound 5 does not. These bonds increased the possibility of binding between the active site and residues by changing the spatial position, thereby increasing the activity. In contrast, compound 5 had a small number of rotatable bonds, which hindered it′s good interaction with proteins, and also the reason for it′s reduced biological activity.
Binding Free Energy Calculation
The MM-GBSA algorithm was used to calculate the binding free energy of the complex formed by CD73 with the higher activity inhibitor 55 and the lower activity compound 5. The statics and histograms are shown in Table 3 and the Fig. S7. It can be seen from the results that the ΔGBind value of compound 55 (-69.82 kcal/mol) was significantly lower than that of compound 5 (-48.49 kcal/mol), which was consistent with the difference in IC50 value.
The binding free energy was consisted of several types of interactions, among which Coulomb and Van der Waals contribute the most to the total free energy of binding of the two compounds. In addition, the values of these two forces of compound 55 were lower than those of compound 5 because its binding site forms additional interactions with surrounding residues. The ΔGBind_covalent of inhibitors 5 and 55 were 6.89 kcal/mol and 9.25 kcal/mol respectively, both of them were positive, indicating that covalent contact hardly contributed to the binding free energy. The ΔGHbond values of compounds 5 and 55 were − 11.08kcal/mol and − 12.21kcal/mol, respectively, because the number of hydrogen bonds formed with compound 55 and the hydrogen bonds formed with ASN390 were more stable than those of compound 5.
In this paper, a series of CD73 inhibitors (5′-(α, β-methylene)- diphosphate) with the same skeleton were studied by 3D-QSAR, molecular docking, molecular dynamics simulation and binding free energy calculation. Molecular docking studies found that 60 CD73 inhibitors interact with related proteins, revealing several key residues, including ARG395, ARG354, ASN117, ASN245, HIP119, ASP506. The CoMFA and CoMSIA models was calculated based on the ligand-based molecular alignment with high reliability. In the 3D-QSAR study, the Q2 and R2 of the CoMFA model were 0.708 and 0.983, respectively, and the Q2 and R2 of the CoMSIA model were 0.809 and 0.992, respectively. In addition, the contour maps can explain the relationship between the chemical structure of the compound and the biological activity, which provided a certain degree of hints and guidance for the design of new molecules. Besides, molecular dynamics simulations were used to analyze the stability of the two docking results. Through the exploration of the entire simulation process, important residues that contributed to the docking results were found. Subsequently, the binding free energy of the two molecules was calculated using the MM-GBSA method, and the results showed that van der Waals and Coulomb energy made major contributions to these two molecules. Using a variety of calculation methods to clarify the mode of action and chemical properties of CD73 inhibitors, which had guiding significance for the further design of high-efficiency CD73 inhibitors with high anti-cancer activity.
Conflicts of interest
The authors declare that they have no conflicts of interest.
Availability of data and material
All the data and materials are available.
Code availability
Not applicable.
Authors' contributions
The authors of the current manuscript Jiatong Wen, Heng Zhang and Churen Meng all contributed to this work. All authors read and approved the final manuscript.
Ethics approval
All authors agree with ethical responsibilities of authors.
Consent to participate
All authors consent to participate.
Consent for publication
All authors consent for publication.
Acknowledgments
This work was financially supported by National Natural Science Foundation of China (Grant No. U1903122, 81872768), LiaoNing Revitalization Talents Program (XLYC1807118), Shenyang Young Scientific and Technological Innovators Programme (RC200408), Doctoral Scientific Research Foundation of Liaoning Province (2020-BS-129) and Special Fund of Research Institute of Drug Regulatory Science Research Shenyang Pharmaceutical University (2021jgkx010), Educational Committee Foundation of Liaoning Province (2020LJC09), and China Postdoctoral Science Foundation (2020M671384).
[1] Pietzner K, Nasser S, Alavi S, et al. Checkpoint-inhibition in ovarian cancer: rising star or just a dream? J Gynecol Oncol. 2018;29(6): e93.
[2] Li B, Chan HL, Chen P. Immune Checkpoint Inhibitors: Basics and Challenges. Curr Med Chem. 2019;26(17):3009-3025.
[3] Bastid J, Regairaz A, Bonnefoy N, et al. Inhibition of CD39 enzymatic function at the surface of tumor cells alleviates their immunosuppressive activity. Cancer Immunol Res. 2015;3(3):254-65.
[4] Canale FP, Ramello MC, Núñez N, et al. CD39 Expression Defines Cell Exhaustion in Tumor-Infiltrating CD8+ T Cells. Cancer Res. 2018;78(1):115-128.
[5] Antonioli L, Pacher P, Vizi ES, Haskó G. CD39 and CD73 in immunity and inflammation. Trends Mol Med. 2013;19(6):355-367.
[6] Minor M, Alcedo KP, Battaglia RA, et al. Cell type- and tissue-specific functions of ecto-5'-nucleotidase (CD73). Am J Physiol Cell Physiol. 2019;317(6):C1079-C1092.
[7] Clayton A, Al-Taei S, Webber J, et al. Cancer exosomes express CD39 and CD73, which suppress T cells through adenosine production. J Immunol. 2011;187(2):676-683.
[8] Xu Z, Gu C, Yao X, et al. CD73 promotes tumor metastasis by modulating RICS/RhoA signaling and EMT in gastric cancer. Cell Death Dis. 2020;11(3):202.
[9] Chen Q, Pu N, Yin H, et al. CD73 acts as a prognostic biomarker and promotes progression and immune escape in pancreatic cancer. J Cell Mol Med. 2020;24(15):8674-8686.
[10] Miao Yu, Gang Guo, Lei Huang, et al. CD73 on cancer-associated fibroblasts enhanced by the A2B-mediated feedforward circuit enforces an immune checkpoint[J].Nature Communications,2020,11(1).
[11] Ma XL, Hu B, Tang WG, et al. CD73 sustained cancer-stem-cell traits by promoting SOX9 expression and stability in hepatocellular carcinoma. J Hematol Oncol. 2020;13(1):11.
[12] Neo SY, Yang Y, Record J, et al. CD73 immune checkpoint defines regulatory NK cells within the tumor microenvironment. J Clin Invest. 2020;130(3):1185-1198.
[13] Zhu J, Zeng Y, Li W, et al. CD73/NT5E is a target of miR-30a-5p and plays an important role in the pathogenesis of non-small cell lung cancer. Mol Cancer. 2017;16(1):34.
[14] Fini C, Talamo F, Cherri S, et al. Biochemical and mass spectrometric characterization of soluble ecto-5'-nucleotidase from bull seminal plasma. Biochem J. 2003;372(Pt 2):443-451.
[15] Ghoteimi R, Braka A, Rodriguez C, et al. 4-Substituted-1,2,3-triazolo nucleotide analogues as CD73 inhibitors, their synthesis, in vitro screening, kinetic and in silico studies. Bioorg Chem. 2021; 107:104577.
[16] Du X, Moore J, Blank BR, et al. Orally Bioavailable Small-Molecule CD73 Inhibitor (OP-5244) Reverses Immunosuppression through Blockade of Adenosine Production. J Med Chem. 2020;63(18):10433-10459.
[17] S. Bhattarai, J. Pippel, A. Meyer. Structure Guides the Way to Subnanomolar Competitive Ecto-5′-Nucleotidase (CD73) inhibitors for Cancer Immunotherapy, Adv. Ther. 2 (2019) 1900075.
[18] S. Bhattarai, J. Pippel, E. Scaletti, R. 2-Substituted α, βMethylene-ADP Derivatives: Potent Competitive Ecto-5′-nucleotidase (CD73) Inhibitors with Variable Binding Modes, J. Med. Chem. 63 (2020) 2941-2957.
[19] Bhattarai S, Freundlieb M, Pippel J, et al. α, β-Methylene-ADP (AOPCP) Derivatives and Analogues: Development of Potent and Selective ecto-5'-Nucleotidase (CD73) Inhibitors. J Med Chem. 2015;58(15):6248-6263.
[20] Bowman CE, da Silva RG, Pham A, Young SW. An Exceptionally Potent Inhibitor of Human CD73. Biochemistry. 2019;58(31):3331-3334.
[21] Lawson KV, Kalisiak J, Lindsey EA, et al. Discovery of AB680: A Potent and Selective Inhibitor of CD73. J Med Chem. 2020;63(20):11448-11468.
[22] Balasubramanian PK, Balupuri A, Bhujbal SP, Cho SJ. 3D-QSAR-aided design of potent c-Met inhibitors using molecular dynamics simulation and binding free energy calculation. J Biomol Struct Dyn. 2019;37(8):2165-2178.
[23] Purohit D, Saini V, Kumar S, Kumar A, Narasimhan B. Three-dimensional quantitative structure-activity relationship (3DQSAR) and molecular docking study of 2-((pyridin-3-yloxy) methyl) piperazines as alpha 7 nicotinic acetylcholine receptor modulators for the treatment of inflammatory disorders. MiniRev Med Chem (2020) 20:1031-1041.
[24] Gajjar KA, Gajjar AK. CoMFA, CoMSIA and HQSAR Analysis of 3-aryl-3-ethoxypropanoic Acid Derivatives as GPR40 Modulators. Curr Drug Discov Technol. 2020;17(1):100-118.
[25] Sun G, Fan T, Zhang N, Ren T, Zhao L, Zhong R. Identification of the Structural Features of Guanine Derivatives as MGMT Inhibitors Using 3D-QSAR Modeling Combined with Molecular Docking. Molecules. 2016;21(7):823.
[26] Sharma H, Cheng X, Buolamwini JK. Homology model-guided 3D-QSAR studies of HIV-1 integrase inhibitors. J Chem Inf Model. 2012;52(2):515-544.
[27] Abu-Hammad A, Zalloum WA, Zalloum H, Abu-Sheikha G, Taha MO. Homology modeling of MCH1 receptor and validation by docking/scoring and protein-aligned CoMFA. Eur J Med Chem. 2009;44(6):2583-2596.
[28] Gu X, Wang Y, Wang M, Wang J, Li N. Computational investigation of imidazopyridine analogs as protein kinase B (Akt1) allosteric inhibitors by using 3D-QSAR, molecular docking and molecular dynamics simulations. J Biomol Struct Dyn. 2021;39(1):63-78.
[29] Ding L, Wang ZZ, Sun XD, et al. 3D-QSAR (CoMFA, CoMSIA), molecular docking and molecular dynamics simulations study of 6-aryl-5-cyano-pyrimidine derivatives to explore the structure requirements of LSD1 inhibitors. Bioorg Med Chem Lett. 2017; 27(15):3521-3528.
[30] Alzate-Morales J, Caballero J. Computational study of the interactions between guanine derivatives and cyclin-dependent kinase 2 (CDK2) by CoMFA and QM/MM. J Chem Inf Model. 2010;50(1):110-122.
[31] Bhatt HG, Patel PK. Pharmacophore modeling, virtual screening and 3D-QSAR studies of 5-tetrahydroquinolinylidine aminoguanidine derivatives as sodium hydrogen exchanger inhibitors. Bioorg Med Chem Lett. 2012;22(11):3758-3765.
[32] Kouman KC, Keita M, Kre N'Guessan R, et al. Structure-Based Design and in Silico Screening of Virtual Combinatorial Library of Benzamides Inhibiting 2-trans Enoyl-Acyl Carrier Protein Reductase of Mycobacterium tuberculosis with Favorable Predicted Pharmacokinetic Profiles. Int J Mol Sci. 2019;20(19):4730.
[33] Kaur K, Talele TT. 3D QSAR studies of 1,3,4-benzotriazepine derivatives as CCK2 receptor antagonists. J Mol Graph Model. 2008;27(4):409-420.
[34] Gu W, Li Q, Li Y. Environment-friendly PCN derivatives design and environmental behavior simulation based on a multi-activity 3D-QSAR model and molecular dynamics. J Hazard Mater. 2020; 393:122339.
[35] Chen Y, Cai X, Jiang L, Li Y. Prediction of octanol-air partition coefficients for polychlorinated biphenyls (PCBs) using 3D-QSAR models. Ecotoxicol Environ Saf. 2016; 124:202-212.
[36] Yang LZ, Liu M. A Double-Activity (Green Algae Toxicity and Bacterial Genotoxicity) 3D-QSAR Model Based on the Comprehensive Index Method and Its Application in Fluoroquinolones' Modification. Int J Environ Res Public Health. 2020;17(3):942.
[37] Iyer P, Bolla J, Kumar V, Gill MS, Sobhia ME. In silico identification of targets for a novel scaffold, 2-thiazolylimino-5-benzylidin-thiazolidin-4-one. Mol Divers. 2015;19(4):855-870.
[38] Wang M, Wang Y, Kong D, et al. In silico exploration of aryl sulfonamide analogs as voltage-gated sodium channel 1.7 inhibitors by using 3D-QSAR, molecular docking study, and molecular dynamics simulations. Comput Biol Chem.2018; 77:214-225.
[39] Meng XY, Zhang HX, Mezei M, Cui M. Molecular docking: a powerful approach for structure-based drug discovery. Curr Comput Aided Drug Des. 2011;7(2):146-157.
[40] Li M, Wei D, Zhao H, Du Y. Genotoxicity of quinolones:substituents contribution and transformation products QSAR evaluation using 2D and 3D models. Chemosphere.2014;95:220–226.
[41] Pradiba D, Aarthy M, Shunmugapriya V, et al. Structural insights into the binding mode of flavonols with the active site of matrix metalloproteinase-9 through molecular docking and molecular dynamic simulations studies. J Biomol Struct Dyn. 2018;36(14):3718-3739.
[42] Nagamani S, Kesavan C, Muthusamy K. Atom-based and Pharmacophore-based 3D - QSAR Studies on Vitamin D Receptor (VDR). Comb Chem High Throughput Screen. 2018;21(5):329-343.
[43] Rajagopal K, Varakumar P, Aparna B, et al. Identification of some novel oxazine substituted 9-anilinoacridines as SARS-CoV-2 inhibitors for COVID-19 by molecular docking,free energy calculation and molecular dynamics studies. J Biomol Struct Dyn 2020; 12:1-12.
[44] Ru Q, Fadda HM, Li C, et al. Molecular dynamic simulations of ocular tablet dissolution. J Chem Inf Model. 2013;53(11):3000-3008.
[45] Vora J, Athar M, Sinha S, et al. Binding Insight of Anti-HIV Phytocompounds with Prime Targets of HIV: A Molecular Dynamics Simulation Analysis. Curr HIV Res. 2020;18(2):132-141.
[46] Sidler D, Riniker S. Fast Nosé-Hoover thermostat: molecular dynamics in quasi-thermodynamic equilibrium. Phys Chem Chem Phys. 2019;21(11):6059-6070.
[47] Ul Haq F, Abro A, Raza S, et al. Molecular dynamics simulation studies of novel β-lactamase inhibitor. J Mol Graph Model. 2017; 74:143-152.
[48] Gao Y, Wang H, Wang J, Cheng M. In silico studies on p21-activated kinase 4 inhibitors: comprehensive application of 3D-QSAR analysis, molecular docking, molecular dynamics simulations, and MM-GBSA calculation. J Biomol Struct Dyn. 2020;38(14):4119-4133.
[49] Chaudhari HK, Pahelkar A. 3D QSAR, Docking, Molecular Dynamics Simulations and MM-GBSA studies of Extended Side Chain of the Antitubercular Drug (6S) 2-Nitro-6- {[4-(trifluoromethoxy) benzyl] oxy}-6,7-dihydro-5H-imidazo[2,1-b] [1,3] oxazine. Infect Disord Drug Targets. 2019;19(2):145-166.
[50] Jayaraj JM, Krishnasamy G, Lee JK, Muthusamy K. In silico identification and screening of CYP24A1 inhibitors: 3D QSAR pharmacophore mapping and molecular dynamics analysis. J Biomol Struct Dyn. 2019;37(7):1700-1714.
Table 1 Statistic parameters for the ligand-based and dockingbased molecular alignment of CoMFA and CoMSIA models
|
Statistical parameters |
Ligand-based alignment |
Docking-based alignment |
||
|
CoMFA |
CoMSIA |
CoMFA |
CoMSIA |
|
|
Q2a |
0.708 |
0.809 |
0.431 |
0.565 |
|
R2b |
0.983 |
0.992 |
0.974 |
0.991 |
|
SEEc |
0.170 |
0.121 |
0.210 |
0.125 |
|
ONCd |
8 |
10 |
7 |
7 |
|
F valuee |
241.755 |
383.332 |
178.607 |
517.126 |
|
R2predf |
0.987 |
0.991 |
0.995 |
0.996 |
|
Fieldcon tribution(%) |
|
|
|
|
|
Steric |
0.788 |
0.162 |
0.809 |
0.131 |
|
Electrostatic |
0.212 |
0.284 |
0.191 |
0.143 |
|
Hydrophobic |
|
0.363 |
|
0.302 |
|
Donor |
|
|
|
0.291 |
|
Acceptor |
|
0.192 |
|
0.143 |
Table 2 Chemical structures, IC50 values of the actual and predicted by the selected CoMFA and CoMSIA models of CD73 inhibitors
Table 3 Calculated binding free energy (kcal/mol) values of CD73 complexed with compound 5 and compound 55 with MM-GBSA method
|
Contribution |
Compound 5 |
Compound 55 |
|
ΔGBind |
-48.49 |
-69.82 |
|
ΔGBind_coulomb |
-81.00 |
-88.02 |
|
ΔGBind_covalent |
6.89 |
9.25 |
|
ΔGBind_Hbond |
-11.08 |
-12.21 |
|
ΔGBind_lipo |
-6.88 |
-17.07 |
|
ΔGBind_packing |
-7.83 |
-7.75 |
|
ΔGBind_vdW |
-41.90 |
-59.76 |