Molecular design of anticancer drugs from marine fungi derivatives.

Heat shock protein 90 (Hsp90) is one of the most potential targets in cancer therapy. We have demonstrated using a combination of molecular docking and fast pulling of ligand (FPL) simulations that marine fungi derivatives can be possible inhibitors, preventing the biological activity of Hsp90. The computational approaches were validated and compared with previous experiments. Based on the benchmark of available inhibitors of Hsp90, the GOLD docking package using the ChemPLP scoring function was found to be superior over both Autodock Vina and Autodock4 in the preliminary estimation of the ligand-binding affinity and binding pose with the Pearson correlation, R = -0.62. Moreover, FPL calculations were also indicated as a suitable approach to refine docking simulations with a correlation coefficient with the experimental data of R = -0.81. Therefore, the binding affinity of marine fungi derivatives to Hsp90 was evaluated. Docking and FPL calculations suggest that five compounds including 23, 40, 46, 48, and 52 are highly potent inhibitors for Hsp90. The obtained results enhance cancer therapy research. This journal is © The Royal Society of Chemistry.

Introduction

The heat shock protein family, including Hsp27, Hsp70, Hsp90, etc., are crucial elements in cell homeostasis.[1-3] In particular, Hsp90 is a dedicated chaperone, and makes up approximately (ca.) 2% of the total protein content and reacts with ca. 200 various proteins. These proteins are a vital part of cell signalling and responses to stress. Hsp90 is thus associated with six stages of cancer,[3,4] and consequently the enzyme is associated with cancer cell survival.[4] Therefore, Hsp90 is distinguished as a target for cancer therapy.[5] Thus, several studies have been executed to discover highly potential inhibitors to prevent biological activity of Hsp90.[6-14]

Nature has always been considered as the main source of therapies; providing folk remedies and therapeutics agents for the treatment of a wide spectrum of diseases. It is estimated that among all the medicines provided on the international market, more than 60% of them are derived from or inspired by natural products, mainly originating from terrestrial plants.[15] On the other hand, in the last 50 years, with the advances in new technologies and engineering in marine science, scientists are increasingly shifting their focus toward marine organisms. Numerous publications in this field have proven that marine organisms provide many natural compounds with biological activities ranging from antiviral to anticancer for the pharmaceutical industries.[16,17] In addition, the rare and unique marine chemical structures could be a scaffold for developing new drugs with greater efficacy and specificity for therapeutics.[18] During the last decades, a reduced number of novel compounds isolated from macro-organisms were observed, but microorganisms such as marine fungi and bacteria are a promising source with several useful drug discoveries already.[19,20] In this light, marine fungi are important not only from the perspective of new drugs but also having the advantage of sustainable production of large quantities of compounds with reasonable costs and large-scale cultivation. Based on the information above, in this study, a database of compounds originating from marine fungi will be investigated to look for Hsp90 potential inhibition activity.

Recently, computer-aided drug design (CADD) has been widely used to decrease both cost and time for developing a novel therapy.[21,22] In CADD, the thermodynamics metric ΔG, corresponds to the binding free energy between protein and ligand, and is able to compute over atomistic simulations. This metric is associated with the experimental data via formula ΔGbind = RT ln(ki), where R is the gas constant, T is the absolute temperature, and ki is the inhibition constant. Normally, the experimental data IC50 is assumed to be equivalent to ki in order to characterize the experimental binding free energy ΔGEXP.[23-25] Because the metric discloses the inhibitor efficiency,[21] accurate and speedy determination of ΔG is very important in the discovery of potential inhibitors for a protein target.[26-29] In this work, we have first benchmarked the performance of the free energy approaches including molecular docking and FPL calculations. The obtained results in the estimation the ΔG of marine fungi derivatives is thus more reliable. Compounds 23, 40, 46, 48, and 52 were suggested as highly possible inhibitors for preventing the biological activity of Hsp90. Our observations can enhance cancer treatment research.

Materials and methods

Structure of ligands and Hsp90

The three-dimensional structure of Hsp90 was downloaded from the Protein Data Bank (PDB) with the identity (ID) of 3K99.[10] The resolution of Hsp90 is high at 1.79 Å. The available inhibitors for Hsp90 were also obtained from the PDB database with IDs: 2QF6,[9]2QG0,[9]2QG2,[9]3K97,[10]3K99,[10] 3 R4M,[11]3R4N,[11]3R4O,[11]3RLR,[12]4NH7,[13]4NH8,[13]4O05,[14]4O07,[14]4O09,[14] and 4O0B.[14] The 3D structures of 63 ligands (E1–E63) were obtained via experimental investigations on marine fungi samples (Aspergillus sp. and Penicillium sp.).

Molecular docking simulations

The molecular docking simulations were employed to evaluate both binding affinity and position of ligands to Hsp90. Here we employed Autodock Vina (Vina) 1.1,[30] Autodock4.2 (AD4),[31] and GOLD packages[32] to finish the task force. The parameters of the docking simulations via AD4 and Vina were referred to in previous studies.[24,33] Whereas the GOLD parameters were selected as default options. The best docking structure is selected as the docking configuration forming the largest binding affinity.

Vina docking[30]

The docking simulation was performed using parameters referred to in previous studies.[33-35] In particular, the exhaustiveness is chosen as 8 and the maximum of different energy between docking configurations is selected as 7 kcal mol−1. The grid size of the Vina docking is selected as 20 × 20 × 20 Å with the center of the grid chosen as the center of mass of the experimental ligand (cf.Fig. 1).

Fig. 1

Computational modelling of molecular docking and FPL simulations.

AD4 docking[31]

AD4 was executed with docking parameters referred to in previous studies.[33] Particularly, grid size was selected as 60 × 60 × 60 with the spacing of 0.333 Å (cf.Fig. 1). The Autogrid4 was executed to generate the docking grid. Ligands were docked to the Hsp90 with the genetic algorithm (GA)/Lamarckian GA (LA)/Monte Carlo Simulated Annealing (SA) number of evaluations was 250.000, equivalent to the short option. Moreover, the GA run is 10, the population size is 150, and the number of generations is 27.000.

GOLD docking[32]

The molecular docking simulation was executed using the GOLD docking package referring to previous work.[36] The grid center was selected as well as the AD4/Vina grid center. All atoms within the range of 12 Å would be considered during docking simulations. The GOLD docking package was performed using two scoring functions as well as ChemScore and ChemPLP.

Fast pulling of ligand (FPL) simulations

GROMACS version 5.1.5 (ref. 37) was employed to mimic the unbinding process of a ligand, which was pulled via an external force out of the binding cavity of Hsp90. The Hsp90 and counteracted ions were presented using the Amber99SB-ILDN force field (FF).[38] The TIP3P was engaged in treating water molecules.[39] The general Amber force field (GAFF)[40] was utilized to parameterize the ligand via AmberTools18 with the charged and geometric parameters obtained via QM calculations using the B3LYP functional and 6-31G(d,p) basis set.[41] In particular, the atomic charge was evaluated using the restrained electrostatic potential (RESP) method.[42] The ACPYPE[43] protocol was then utilized to convert the ligand topology to GROMACS formats. The combined FF was selected since they are high coincidence for a free energy assessment.[44,45]

The Hsp90 + ligand complex was initially inserted into a rectangular periodic boundary condition box with a size of 5.8 × 5.8 × 8.5 nm, which is large enough for unbound ligands out of the enzymic binding cavity (cf.Fig. 1). The soluble complex thus comprises of ca. 27 500 atoms involving the Hsp90, ligands, water molecules, and ca. 7 Na+ ions. The steepest descent approach was initially utilized to minimize the energy of the system. Then, NVT and NPT imitations were followed to relax the system with a length of 0.1 ns for each simulation. It should be noted that the Cα atoms of Hsp90 were positionally restrained using a weak harmonic potential. Finally, the inhibitor was stressed to mobilize out of the HSP binding cavity by utilizing an external harmonic force with a pulling speed of k = 0.005 nm ps−1 and a spring constant of ν = 600 kJ mol−1 nm−2. The metrics k and ν were chosen as referred to in previous work.[27] During FPL calculations, the displacement of the ligand and the value of pulling force, F, were recorded every 0.1 ps and were used to estimate the ligand-binding affinity.[27] The FPL calculations were carried out 8 times for each system.

Analyzed tools

The ligand protonation state was predicted using the Chemicalize tools (http://www.chemicalize.com), a website application of ChemAxon. The binding position of the ligand to Hsp90 was analyzed using the free version of Maestro 2020.[46]

Results and discussion

Molecular docking simulation

In CADD, the binding pose and binding affinity of ligands to enzymic targets are generally obtained via molecular docking simulations, a fast and low computing resources technique. In this work, we have tried Vina,[30] an open-source docking protocol, widely used for this purpose to complete this task force. Unfortunately, Vina calculations adopted uncorrelated results in comparison with the respective experimental data with correlation coefficient Rvina = 0.40 (Table 1). AD4 was thus employed to perform the docking simulation since they have a different scoring function to Vina:.[30,31] Vina uses a full empirical scoring function, while AD4 uses a semi-physical/empirical based scoring function.[30,31] We thus expected that AD4 probably is suitable for completing the task force and the obtained docking result would be appropriate. Unfortunately, the obtained correlation coefficient is poor with a value of RGAAD4 = 0.33 (Table 1). It is consistent with previous work.[33] Although we have changed search algorithms from GA to Lamarckian GA (LA) and Monte Carlo Simulated Annealing (SA), the poorly correlated results were observed with a value of RLAAD4 = 0.35 and RSAAD4 = 0.25 (Table 1).

The obtained values of the docking simulations

No.	PDB ID	ΔGvina	ΔGGAAD4	ΔGLAAD4	ΔGSAAD4	ΔGChemScoreGOLD	ChemPLP	ΔGEXPa
1	2QF6	−10.2	−8.8	−8.8	−8.8	−8.8	74.7	−8.91
2	2QG0	−9.6	−8.7	−9.2	−10.4	−7.0	66.4	−7.85
3	2QG2	−10.2	−9.9	−9.4	−10.0	−7.9	75.2	−7.41
4	3K97	−8.3	−7.7	−7.9	−8.1	−7.3	68.3	−10.98
5	3K99	−9.2	−8.2	−8.2	−8.2	−7.9	61.4	−9.91
6	3 R4M	−7.8	−6.6	−6.6	−7.5	−5.1	49.7	−8.13
7	3R4N	−8.3	−7.1	−7.4	−7.6	−5.4	61.2	−9.45
8	3R4O	−10.1	−8.8	−8.7	−9.2	−6.0	73.2	−11.19
9	3RLR	−9.6	−8.2	−8.2	−8.2	−6.4	70.5	−10.33
10	4NH7	−11.3	−10.0	−9.9	−10.6	−11.8	98.9	−11.53
11	4NH8	−10.7	−9.7	−9.8	−10.0	−11.8	98.3	−11.70
12	4O05	−12.5	−11.5	−11.4	−11.5	−11.4	96.4	−10.20
13	4O07	−12.5	−11.9	−11.9	−12.0	−11.5	95.5	−10.13
14	4O09	−11.6	−10.7	−10.7	−10.8	−11.0	93.8	−10.57
15	4O0B	−12.3	−11.5	−11.7	−11.7	−11.1	94.2	−11.39

The experimental binding free energy was obtained through the inhibition constant. The unit of energy and force are in kcal mol−1 and pN, respectively.

The GOLD docking package,[32] a commercial docking program, was then utilized to evaluate the binding affinity and position of ligands to the Hsp90 enzyme. Fortunately, the docking energy provided by the GOLD docking package forms appropriate correlations with the respective experiment with a value of RChemScoreGOLD = 0.51 and RChemPLPGOLD = −0.60 (cf.Table 1 and Fig. 2). The ChemPLP scoring function is much better than ChemScore, and the ligand having a larger ChemPLP score will form a larger binding affinity to Hsp90. Therefore, the GOLD docking package with the ChemPLP scoring function was executed to predict the binding affinity of marine derivatives to Hsp90.

Fig. 2

Association between the ChemPLP score and the experimental binding free energy ΔGEXP.

FPL simulations

Although the docking simulation adopted appropriate results compared with experimental data,[47-49] the docking results are required to refine a more accurate approach.[35,50,51] Because the molecular docking approach normally uses many constraints to enhance the computing speed. The popular constraints of molecular docking simulations are known as lacking receptor dynamics, limiting the number of trial position of ligands, using an implicit water model, etc. Normally, molecular dynamics and/or steered-molecular dynamics simulations are thus utilized to refine the docking outcomes.[28,52,53] FPL simulations were employed to validate the docking results as they form corresponding results with the respective experimental results, and are a low cost computing resource.[54]

Furthermore, it should be noted that FPL simulations are a very efficient free energy approach, which normally forms appropriate results compared to the respective experiments with low computing resources. FPL has thus been successfully applied to various targets before.[27,35,54] In this context, we first benchmarked the performance of the FPL scheme over 15 Hsp90 systems. The obtained results are shown in Table 2 and Table S2† of the ESI file. In particular, the mean of pulling work (W) diffuses in the range 58.3 to 208.5 kcal mol−1, whereas the mean of rupture force (FMax) is in the range 580.7 to 1566.1 pN. The graph of W and FMax in the time revolution is described in the ESI file.† The computed metrics are in good agreement with the respective experimental data,[9-14] because the Pearson correlation coefficient is of RW = −0.78 and RF = −0.81 (Table 2 and Fig. 3). Although the difference between RW and RF is quite small, the rupture force was selected to rank the ligand-binding affinity to Hsp90. Therefore, using linear regression, the predicted binding free energy ΔGPreFPL can be calculated via the rupture force as follows

The obtained values of the FPL calculations in comparison with the respective experiments

No.	Name	F Max a	W b	ΔGEXPc
1	2QF6	580.7 ± 90.9	67.7 ± 12.3	−8.91
2	2QG0	764.3 ± 119.3	76.5 ± 12.9	−7.85
3	2QG2	676.2 ± 78.5	71.1 ± 8.7	−7.41
4	3K97	1065.5 ± 66.3	113.9 ± 10	−10.98
5	3K99	831.4 ± 52.5	82.7 ± 4.5	−9.91
6	3 R4M	605.7 ± 81.7	58.3 ± 9.4	−8.13
7	3R4N	651.0 ± 37.6	69.8 ± 4.4	−9.45
8	3R4O	1140.8 ± 128.1	126.3 ± 17.2	−11.19
9	3RLR	910.8 ± 88.1	98.8 ± 8.6	−10.33
10	4NH7	1357.1 ± 146.7	167.6 ± 28	−11.53
11	4NH8	1566.1 ± 124.2	208.5 ± 26.1	−11.70
12	4O05	1313.2 ± 88.2	164.4 ± 20.3	−10.20
13	4O07	1363.3 ± 131.0	175.7 ± 25.9	−10.13
14	4O09	1348.6 ± 102.6	163.9 ± 18.6	−10.57
15	4O0B	1318.9 ± 112.6	157.5 ± 20.8	−11.39

The obtained value of the mean rupture force FMax.

The recorded metric of the pulling work W.

The experimental binding free energy ΔGEXP was estimated via the reported ki. The calculated error was the standard error of the average. The unit of force and energy is in pN and kcal mol−1, respectively.

Fig. 3

Association between the average of the FMax and the experimental binding free energy ΔGEXP.

The metric ΔGPreFPL was investigated by the RMSE with linear regression, forming RMSEF = 0.77 kcal mol−1, which implies that FPL calculations can categorize ligands having similar binding affinity to Hsp90. Overall, FPL simulations are an efficient approach to predict the ligand-binding free energy of Hsp90.

Marine fungi derivatives inhibit Hsp90

As mentioned above, marine fungi derivatives are vital compounds that are promising candidates for preventing the biological activity of Hsp90. The GOLD docking package using the ChemPLP scoring function was utilized to estimate the binding affinity and binding pose of marine fungi derivatives preliminarily calculated using this application. It should be noted that trial compounds were extracted from marine fungi using experiments. The obtained results are reported in Table 3 and S1† in the ESI file. In particular, the ChemPLP scores of these compounds vary in the range from 31.99 to 77.8. The compounds having ChemPLP scores larger than 70.0 are proposed as potential inhibitors for preventing Hsp90 (cf.Table 3).

The obtained values of the molecular docking and FPL calculations

No.	Name	ChemPLP	F Max a	W b	ΔGPreFPLc	k Pre i
1	48	86.3	627 ± 127.1	70.6 ± 9.2	−8.62	0.29
2	25	77.8	555.3 ± 60.5	61.7 ± 8.2	−8.38	1.25
3	40	76.9	711.7 ± 77.9	75.7 ± 13	−8.90	0.54
4	23	72.6	827.6 ± 195	95.9 ± 22.3	−9.29	0.67
5	46	71.7	669.9 ± 97.2	78.1 ± 9.8	−8.76	0.85
6	52	70.6	753.3 ± 110	77.5 ± 10.5	−9.04	0.43

The obtained value of the mean rupture force FMax.

The recorded metric of the pulling work W.

The predicted binding free energy ΔGPreFPL was estimated using formula (1). The calculated error was the standard error of the average. The unit of force, energy, and inhibition constant are pN, kcal mol−1 and μM, respectively.

The binding positions of the six best compounds docking to the Hsp90 enzyme were analyzed using the PyMOL package.[55] Interestingly, only ligand 23 forms a hydrogen bond to Hsp90, whereas the other compounds just adopt hydrophobic contacts. This observation implies that the hydrophobic interaction dominates over electrostatic interactions in the binding process of ligands to Hsp90. This is consistent with the correlation between experimental binding free energy with the metric C log P, R = −0.64, which means that more hydrophobic inhibitors form stronger binding ligands (Fig. 4).

Fig. 4

The 6 binding positions of the marine fungi compounds docking the best to Hsp90, which were estimated using the GOLD docking program with ChemPLP scoring function.

FPL calculations were then applied to refine the docking results of 63 marine fungi derivatives to Hsp90. The obtained results are shown in Table 3 and S3† of the ESI file. The average of FMax falls in the range 555.3 to 827.6 pN, whereas the mean W falls in the range 57 to 94 kcal mol−1. Using linear regression via formula (1), the predicted binding free energy of 6 marine compounds was obtained as listed in Table 3. The predicted metrics ΔGPreFPL, are in the range −8.38 to −9.29 kcal mol−1, which is equivalent to the predicted inhibition constant kPrei ranging from 0.29 to 1.25 μM. Interestingly, five ligands adopt a sub-micromolar inhibition constant implying that they probably are promising inhibitors for blocking the biological activity of Hsp90. In addition, it should be noted that the docking conformation differs from the native-binding conformation, resulting in a decrease of FPL accuracy.[56] The MD simulations with a length of 20 ns for each trajectory, were carried out to let ligands to reach a stable-binding conformation. The complex reached equilibrium states after 5 ns of MD simulations (cf. Table S4† of the ESI). The last snapshot of MD simulations was then used as the initial conformation for FPL calculations. The obtained results were reported in Table S5 and S6† of the ESI file. Interestingly, the FPL results reported in Table 3 and S5† are in good agreement with each other. Additional evaluation using in vitro and/or in vivo studies should be executed to validate the obtained results.

CPU time consumption

Each Hsp90 + inhibitor system was mimicked over 8 independent FPL calculations as well as the performed calculations of the different targets.[35,51] In particular, different FPL calculations have the same starting structure and one FPL simulation involves 0.1 ns of NVT, 0.1 ns of NPT, and 0.5 ns of SMD simulations. 5.6 ns of MD simulations were thus carried out to predict the binding free energy of a ligand to the Hsp90 target. It should be noted that a Dual Xeon E5-2683 V3 server is able to perform ca. 90 ns of MD simulation each day for the Hsp90 + inhibitor complex. Therefore, the FPL approach can be used to estimate the binding free energy of a ligand to Hsp90 within 1.5 hours. It is quite a low computing resource requirement, allowing us to be able to accurately and rapidly evaluate the binding free energy of several ligands to Hsp90.

Conclusions

We have established that an amalgamation of GOLD docking and FPL simulations can accurately and rapidly calculate the binding free energy of a ligand to Hsp90. The ChemPLP scoring function is better than ChemScore one in docking a molecule to Hsp90, since RChemScoreGOLD = 0.51 is significantly smaller than RChemPLPGOLD = −0.60. Besides that, we also tried AD4 and Vina docking to complete the task force, but a poor correlation to the respective experiments was observed. Moreover, FPL calculations are an efficient protocol to refine the docking result as they form a good correlation coefficient RF = −0.81 and RMSEF = 0.77 kcal mol−1 with a low computing resource requirement. Furthermore, based on the efficiency of the approaches, a shortlist comprising 5 compounds from marine fungi derivatives, including 23, 40, 46, 48, and 52 were indicated as promising inhibitors for preventing the biological activity of Hsp90 because they adopt a sub-micromolar inhibition constant. Further evaluation using in vitro and/or in vivo studies should be executed to validate the obtained results.

Conflicts of interest

The authors declare no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.