Structure–activity relationship guided design, synthesis and biological evaluation of novel diaryl urea derivatives as antiproliferative agents

Abstract

Diaryl urea-based compounds have attracted the attention of many researchers due to their potential as anticancer agents. Following our previous study on a series of diaryl urea compounds and implementation of the obtained structure activity relationship (SAR) analysis, a new set of derivatives were designed and synthesized. The synthesized compounds were subjected to evaluation for their in vitro antiproliferative activities against A549 and HT-29 cell lines. Among all, 6a emerged as the most potent antiproliferative agent with an IC₅₀ value of 15.28 and 2.566 µM against HT-29 and A549 cells, respectively. Comparing the activity of the newly designed and synthesized diaryl urea compounds 4a-b and 6a-e with those for the previously reported compounds 8a-b and 9a-f confirmed the importance of the substitution of amide groups instead of ester between the central and distal benzene rings of diaryl urea scaffold. The results of current study revealed that the substitution of proximal and distal benzene rings with chlorine and methyl groups, alongside the linear extension of molecules through the introduction of a methylene spacer group could enhance antiproliferative activity, which is in agreement with previously reported SAR analysis. Molecular docking simulations demonstrated that all designed compounds exhibit binding affinity to VEGFR-2 similar to that observed experimentally for sorafenib. The findings of this study may offer valuable insights for the further development of diaryl urea-based anticancer agents.

Graphical Abstract

Introduction

Metastasis (the capability of spreading to new organs) makes cancer a life-threatening issue. Metastasis of cancer tissue occurs by the growth of the vascular network (blood vessel formation) which is known as angiogenesis [1]. The regulation of angiogenesis relies on vascular endothelial growth factors (VEGFs) [2]. Among these receptors, VEGFR-2 is mainly expressed in endothelial cells [3, 4]. The discovery of VEGFR-2 inhibitors has gained attention as a promising treatment approach for cancer [5]. Sorafenib, a diaryl urea derivative (Fig. 1) as a VEGFR2-TK (VEGFR2- tyrosine kinase receptor) inhibitor was approved for clinical use in the treatment of advanced renal cell carcinoma (RCC) and unresectable hepatocellular carcinoma (HCC) [6, 7]. However, adverse side effects like hypertension and hand-foot skin reaction, as well as unfavorable physicochemical properties limit sorafenib application. Lead optimization represents a critical step in drug discovery before entering a drug candidate into the development phase [8]. The structure–activity relationship (SAR) as an important step in the drug discovery process explains the correlation between the structural features of a series of molecules and their corresponding activities. In the lead optimization process, the elucidated SAR analysis are employed to make structural modifications on lead structures to optimize their physicochemical and activity properties [9]. Encouraged by the results obtained from our previous works (17, 18), we decided to further expand the SAR exploration of diaryl urea derivatives as anticancer agents. To this aim, the introduction of various substituents was performed onto the proximal and distal phenyl rings (phenyl rings next to and distal to the urea group, respectively) of diaryl urea compounds. Figure 1 shows the overall structure of newly designed compounds.

Fig. 1

Chemical structure evolution of diaryl urea derivatives 6a-h based on sorafenib as a reference drug

Nowadays, the utilization of short and efficient synthetic routes is crucial for the synthesis of pharmaceuticals. Such short synthetic pathways enable the preparation of key compounds with minimal energy expenditure and resource utilization, accelerating drug discovery and development activities. In this study, a series of novel diaryl-urea compounds were designed and synthesized using the route shown in Scheme 1. Then, the in vitro antiproliferative potency of synthesized compounds was evaluated against two human cancer cell lines including human non-small cell lung cancer (NSCLC) cell line A549 and human breast cancer cell line HT-29.

Scheme 1.

General synthetic route to target derivatives 6a-g and 9a-f. Reagents and conditions: (a) DCM, Et₃N, r.t., (b) DCM, reflux

Results and discussion

Chemistry

Scheme 1 shows the synthetic route used for the preparation of designed target compounds (6a-f and 4b). Applying a sequential one-pot reaction strategy, the designed compounds were synthesized with high yield and easy workup. Initially, isocyanate derivatives (2a and 2b) and 1,4-diaminophenyl 3 were reacted for 3 h at room temperature in dry CH₂Cl₂ to obtain compounds 4a and 4b. Then the required acyl chloride [5] was added dropwise to afford desired compounds (6a-g) in a single reaction mixture with quantitative yields. The ester derivatives were synthesized reported in our previous study [10, 11].

The suggested route has two advantages including high production yield (> 90%) and simplicity. The established protocol outlines a practical, simple, and rapid one-pot sequential synthesis procedure for preparing new diaryl urea compounds with yields close to 100% [12]. Furthermore, it eliminates the need for tedious work-up and purification processes, making the method advantageous and easily applicable for large-scale production.

In vitro antiproliferative activity

The in vitro antiproliferative activity of the target compounds 4a-b, 6a-g, 8a-b, and 9a-f against HT-29 and A549 cells was evaluated by the MTT assay using sorafenib as the positive control (Table 1). Based on our previously reported SAR analysis for diaryl urea derivatives [10, 12], (i) the presence of an amide group between central benzene ring and distal benzene ring, (ii) presence of benzothiophene ring, (iii) substitution of R₁ in phenyl ring with chlorine and using 4-methylphenyl as R₂, (iv) as well as the elongation of the molecule by introduction of a methylene spacer group between the distal benzene ring and urea group could collectively enhance antiproliferative activity. Accordingly, we designed and synthesized compounds 4a-b and 6a-g, evaluated their antiproliferative activity against HT-29 and A549 cells, and compared the results to those of previously synthesized compounds 8a-b and 9a-f. Among the synthesized compounds, the most potent compound was 6a with benzo[b] thiophene ring (R₂), 4-Cl group on phenyl ring (R₁), and an amide group between central benzene ring and distal benzene ring. The compound 6a showed antiproliferative activity with an IC₅₀ value of 15.28 µM against HT-29 and 2.566 µM against A549 cells comparable to those shown for the positive reference drug sorafenib (IC₅₀ value of 14.01 on HT-29 and 2.913 to A549). The IC₅₀ values for the corresponding ester compound 9a (ester group instead of amide in 6a) were 114.4 and 99.67 against HT-29 and A549 cells, respectively. This result showed that consistent with observed SAR, the substitution of amide group instead of ester leads to increased antiproliferative activity. Interestingly, the antiproliferative activity of all compounds with amide groups between central benzene ring and distal benzene ring was more than corresponding compounds with ester moiety. For example, IC₅₀ values for compounds 4a-b, 6a-e, and 6 g are less than those for 8a-b and 9a-f on both cells. Furthermore, we synthesized compounds 6a-f to investigate the effects of substitutions of proximal and distal rings of diaryl urea derivative as well as the replacement of ester group with amide group, respectively.

Table 1 In vitro antiproliferative activities (IC₅₀) of synthesized compounds (4a-b, 6a-g, 8a-b and 9a-f) against HT-29 and A549 cells.^a

As indicated in Table 1, the antiproliferative activities of some compounds were affected by the substitution of a chloro group (R₁) on the phenyl ring proximal to the urea group. For instance, compounds 6b and 6d (with R₁ = 4-Cl₎, exhibited stronger antiproliferative activity against both cells compared to those derivatives lacking the 4-Cl group, i.e., compounds 6e and 6 g. However, such an enhancing effect of 4-chloro substituted phenyl moiety was not observed for compounds 6c, 4b, and 8b compared to that of 6f, 4a, and 8a (the corresponding compound without 4-Cl group). Also, the comparison of antiproliferative activity of compound 6b with 4-methylphenyl (R₂) to that of corresponding compound 6c with phenyl revealed that substitution of methyl group led to the increased antiproliferative activity. Another strategy that was employed for further structural modification in this study was linear elongation (homologation) in diaryl urea derivatives with an amide group. This modification increased the distance between the proximal and distal phenyl rings due to the addition of a methylene group, resulting in compounds 6d. Comparison between 6c and 6d suggests that benzyl moiety (i.e., elongation) can enhance antiproliferative activity on both cells. Collectively, the results shown in Table 1 indicated that, in both HT29 and A549 cell lines, antiproliferative activity was increased by inserting a chloro group (R₁) on the proximal ring, introducing methyl group on distal ring, and also increasing the distance between the proximal and distal phenyl rings by inserting a methylene group. The obtained SAR on newly synthesized diaryl urea compounds (Fig. 2) is completely consistent with our previous works [10,11,12]. Based on the findings described above, further in silico evaluations were performed on the synthesized compound.

Fig. 2

Graphical representation for SAR of new diarylurea compounds

Interaction of diaryl urea derivatives with VEGFR-2 studied by molecular docking

To explore the interactions between the synthesized derivatives and the kinase domain of VEGFR-2, a molecular docking study was conducted. The binding site for the studied compounds on VEGFR-2 was determined based on the experimentally known binding site for sorafenib (PDB code: 4ASD) [13]. In silico docking calculations were performed using the GOLD program with the ChemPLP scoring function and applying the default settings of the genetic algorithm. The docked poses of the most potent compound 6a and sorafenib are shown in Fig. 3. As shown in figure, 6a established different interactions with the receptor similar to those observed for sorafenib in crystal structure 4ASD [13]. (The PoseView presentations of other synthesized compounds are available in supporting information). The interaction modes of the synthesized compounds with VEGFR-2 obtained from the docking study are detailed and compared to those of sorafenib. The docking results revealed the formation of two hydrogen bonds between 6a and the receptor (Fig. 3). As represented in the figure, one hydrogen bond was formed between a urea nitrogen atom of 6a and the oxygen atom of the carboxyl group in the side chain of Glu⁸⁸⁵ from VEGFR-2. Another hydrogen bond was identified between the urea oxygen of compound 6a and the amide hydrogen of Asp¹⁰⁴⁶ from the receptor. Three aromatic rings of the ligand form different hydrophobic interactions with Phe¹⁰⁴⁷, Phe⁹¹⁸, Asp¹⁰⁴⁶, Leu⁸⁴⁰, Leu⁸⁸⁹, Cyc¹⁰⁴⁵ residues of the receptor. Docking calculations indicated the formation of hydrophobic interactions between 5-nitrobenzo[b] thiophene and central benzene ring of 6a and some residues of the receptor not shown for sorafenib (Fig. 3). These additional interactions may be the reason for the appropriate potency of 6a. Other compounds form different hydrophobic interactions and two H-bonds with the receptor (See supporting information). Taken together, the findings from the docking simulations and MTT antiproliferative assays show a reasonable level of consistency, indicating the reliability of the docking results presented. Nevertheless, experimental studies such as ligand-receptor binding assays and X-ray structure determination are necessary to validate the in silico findings presented.

Fig. 3

PoseView presentation of 6a and sorafenib docked on VEGFR-2

Stability of 6a-VEGFR2 proposed complex assessed by MD simulation

The above described solution for 6a docked into the binding site of VEGFR2 was subjected to MD simulation for 50 ns. The RMSD values and total energy plots for the simulation systems are shown in Fig. 4. The RMSD values for the whole complex and 6a ligand alone fluctuate in narrow range about 1.49 ± 0.33 Å and 1.10 ± 0.08 Å, respectively, which are in close agreement with the same measurements for the complex of VEGFR2 with sorafenib (1.70 ± 0.46 Å) and sorafenib alone (0.71 ± 0.15 Å). The total energy for the complexes of 6a and sorafenib with VEGFR2 are very stable during the MD simulation with an average of -113,865.6127 ± 225.2623 and -113,726.9126 ± 224.5364 kcal/mol, respectively. These results indicated the stability of the proposed docked model for the complex between 6a and VEGFR2. Furthermore, the MD trajectories were used to calculate the binding free energy by applying molecular mechanics (MM) Generalized Born Surface Area (MM/GBSA) and NAB normal mode analysis (NMode) methods implemented in the AMBER package. The obtained ∆G binding (Generalized Born calculation without entropy contribution to the binding) and binding entropy (–TΔS_bind calculated using NMode analysis) values were -53.0846 ± 3.2089 kcal/mol and −19.3750 ± 7.3495 kcal/mol, respectively, leading to the estimated total binding energy of −33.7096 kcal/mol. Applying quasi-harmonic entropy approximation and MMGBSA methods, the calculated ∆G binding was −16.0533 kcal/mol which is due to the difference in the predicted entropy of −37.9748 kcal/mol by the latter method. Applying the same calculations to the results of MD simulation for sorafenib-VEGFR2 complex resulted in binding energy values of −29.5963 kcal/mol (taking into account both GB and NMode analysis) and −18.4358 kcal/mol (using quasi-harmonic entropy approximation and MMGBSA methods). According to these results, the estimated energy values for the binding of 6a and sorafenib to VEGFR2 receptor are very close.

Fig. 4

The results of MD simulations on 6a-VEGFR2 and sorafenib-VEGFR2 complexes. A. RMSD values of 6a-VEGFR2 and sorafenib-VEGFR2 complexes and 6a and sorafenib ligands during 0.5 ns equilibration and 50 ns production MD simulations. B. The total potential energy of 6a-VEGFR2 and sorafenib-VEGFR2 complexes during the MD simulation

Conclusions

A series of novel diaryl urea derivatives were designed based on the results of SAR analysis from our previous study. The designed compounds were synthesized using a facile one-pot sequential reaction method. The synthesized derivatives were then evaluated for their antiproliferative activities against two cancer cell lines, HT-29 and A549. Notably, compound 6a exhibited significantly high activity, with IC₅₀ values of 15.28 µM and 2.566 µM against HT-29 and A549 cells, respectively. The potencies observed for 6a were comparable to the positive reference drug sorafenib, which had IC₅₀ values of 14.01 µM against HT-29 and 2.913 µM against A549 cells. The SAR analysis of diaryl urea derivatives verified the previously described beneficial effects of substituting an amide group instead of ester group. Such modification can be considered as having significant antiproliferative enhancing effects in compounds with diaryl urea scaffold. Furthermore, molecular docking studies revealed that compound 6a binds to VEGFR-2 like sorafenib. The proposed complex with favorable binding free energy was stable during the 50 ns MD simulation. These results provide valuable insights for the design of new sorafenib-related anticancer drugs leveraging the provided SAR analysis.

Experimental

Materials and methods

¹H-NMR and ¹³C-NMR spectra of the compounds were recorded in DMSO-d6 by using Bruker 400 MHz (Bruker Bioscience, Billerica, MA, USA) with tetramethylsilane (TMS) as an internal standard. Agilent 5973 Network Mass Spectrometer system (Agilent Technologies, Inc., USA) was used for obtaining mass spectra. MTT (3-(4, 5-dimethyl-2-thiazolyl)-2,5-diphenyl-2H-tetrazolium bromide) were purchased from Sigma-Aldrich (USA). Without additional purification, all commercially available chemicals and reagents were used.

General procedures for synthesis of target compounds (6a-h)

To prepare target compounds 6a-g and 9a-f, a sequential one-pot reaction strategy was applied which the details of the synthetic pathway was described in our previous work [12]. The structure of the compounds has been characterized by various spectral methods (¹H-NMR, ¹³C-NMR, Mass, IR) for compounds 4b, 6a-6e which explained in detail below and for compound 6f and 6 g reported previously [10].

4b: 1-(4-aminophenyl)-3-(4-chlorophenyl) urea:

¹H NMR (400 MHz, DMSO-d6) δ (ppm): 4.79 (s, 2H, NH2), 6.51 (d, J = 8 Hz, 2H, ArH-6,2), 7.07 (d, J = 8 Hz, 2H, ArH-13,15), 7.29 (d, J = 4 Hz, 2H, ArH-3,5), 7.45 (d, J = 8 Hz, 2H, ArH-16,12), 8.16 (s, 1H, NH-7), 8.62 (s, 1H, NH-9). ¹³C NMR (100 MHz, DMSO-d6) δ (ppm): 114.15 (2C-2,6), 119.43 (2C-16,12), 120.95 (2C-3,5), 124.83 (1C-14), 128.34 (2C-13,15), 128.54 (1C-4), 139.21 (1C-11), 144.22 (1C-1), 152.82 (1C-8).

6a: N-(4-(3-(4-chlorophenyl)ureido)phenyl)-5-nitrobenzo[b]thiophene-2-carboxamide.

IR(KBr): 3302, 1367, 1345, 1308 cm⁻¹. ¹H-NMR (400 MHz, DMSO-d6) δ (ppm): 7.31 (d, J = 4 Hz, 2H, ArH-25,27), 7.48 (t, J = 8 Hz, 4H, ArH-24,28,14,18), 7.68 (d, J = 8 Hz, 2H, ArH-15,17), 8.25 (d, J = 4 Hz, 1H-3), 8.32 (d, J = 8 Hz, 1H-2), 8.49 (s, 1H-9), 8.71 (s, 1H-6), 8.80 (s, 1H, NH-21), 8.92 (s, 1H, NH-19), 10.66 (s, 1H, NH-12). ¹³C NMR (100 MHz, DMSO-d6) δ (ppm): 118.66 (2C-15,17), 119.71 (2C-28,24), 120.09 (1C-2), 120.96 (1C-6), 121.25 (2C-14,18), 124.22 (1C-3), 125.32 (1C-9), 125.99 (1C-26), 128.61 (2C-25,27), 132.68 (1C-13), 135.95 (1C-8), 138.77 (1C-5), 139.05 (1C-16), 144.08 (1C-23), 145.44 (1C-1), 146.15 (1C-4), 152.46 (1C-20), 159.28 (1C-10).). ESI-HRMS (m/z): Calcd. for C₂₂H₁₅ClN₄O₄S [M + H] + : 466.05, found:467.4.

6b: N-(4-(3-(4-chlorophenyl)ureido)phenyl)-4-methylbenzamide.

IR(KBr): 3417, 1634 cm⁻¹. ¹H NMR (400 MHz, DMSO-d6) δ (ppm): 2.38 (s, 3H, CH3-26), 7.32 (d, J = 8 Hz, 4H, ArH-17,19,22,24), 7.43 (d, J = 8 Hz, 2H, ArH-21,25), 7.49 (d, J = 4 Hz, 2H, ArH-6,10), 7.69 (d, J = 8 Hz, 2H, ArH-16,20), 7.88 (d, J = 4 Hz, 2H, ArH-7,9), 8.65 (s, 1H, NH-13), 8.78 (s, 1H, NH-11), 10.08 (s, 1H, NH-3). ¹³C NMR (100 MHz, DMSO-d6) δ (ppm): 21.01 (1C-CH3-26), 118.62 (2C-7,9), 119.70 (2C-21,25), 121.11 (2C-6,10), 125.28 (2C-23), 127.62 (1C-16,20), 128.62 (2C-17,19), 128.88 (2C-22,24), 132.19 (1C-1), 133.74 (1C-5), 135.29 (1C-8), 138.84 (1C-15), 141.38 (1C-18), 152.50 (1C-12), 165.02 (1C-2). ESI-HRMS (m/z): Calcd. for C₂₁H₁₈ClN₃O₂ [M + H] + : 379.11, found:379.11.

6c: N-(4-(3-(4-chlorophenyl) ureido)phenyl)benzamide.

IR(KBr): 3381, 1645 cm⁻¹. ¹H NMR (400 MHz, DMSO-d6) δ (ppm): 7.32 (d, J = 4 Hz, 2H, ArH-22,24), 7.44 (d, J = 8 Hz, 2H, ArH-17,19), 7.50 (d, J = 4 Hz, 2H, ArH-21,25) 7.53 (d, J = 8 Hz, 2H, ArH-16,20), 7.57 (d, J = 8 Hz, 1H, ArH-18), 7.70 (d, J = 4 Hz, 2H, ArH-6,10), 7.96 (d, J = 4 Hz, 2H, ArH-7,9), 8.67 (s, 1H, NH-13), 8.79 (s, 1H, NH-11), 10.18 (s, 1H, NH-3). ¹³C NMR (100 MHz, DMSO-d6) δ (ppm): 118.62 (2C-7,9), 119.71 (2C-21,25), 121.12 (2C-6,10), 125.29 (2C-23), 127.60 (1C-16,20), 128.37 (2C-17,19), 128.64 (2C-22,24), 131.42 (1C-18), 133.65 (1C-1), 135.08 (1C-8), 135.40 (1C-5), 138.83 (1C-15), 152.51 (1C-12), 165.21 (1C-2). Calcd. for C₂₀H₁₆ClN₃O₂ [M + H] + : 365.09, found:365.1.

6d: N-(4-(3-(4-chlorophenyl)ureido)phenyl)-2-phenylacetamide.

IR(KBr): 3415, 1640 cm⁻¹. ¹H NMR (400 MHz, DMSO-d6) δ (ppm): 3.61 (s, 2H, CH2), 7.25 (d, J = 4 Hz, 1H, ArH-18), 7.30–7.33 (m, 6H, ArH-16,17,19,20,22,24), 7.37 (d, J = 8 Hz, 2H, ArH-21,25) 7.47 (d, J = 4 Hz, 2H, ArH-6,10), 7.51 (d, J = 8 Hz, 2H, ArH-7,9), 8.60 (s, 1H, NH-13), 8.75 (s, 1H, NH-11), 10.06 (s, 1H, NH-3). ¹³C NMR (100 MHz, DMSO-d6) δ (ppm): 43.29 (1C-1), 118.78 (2C-7,9), 119.68 (2C-21,25), 119.82 (2C-6,10), 125.25 (1C-23), 126.50 (1C-18), 128.31 (2C-17,19), 128.62 (2C-22,24), 129.10 (2C-16,20), 133.78 (1C-5), 134.97 (1C-8), 136.17 (1C-14), 138.83 (1C-15), 152.48 (1C-12), 168.69 (1C-2). Calcd. for C₂₁H₁₈ClN₃O₂ [M + H] + : 379.11, found:379.2.

6e: 4-methyl-N-(4-(3-phenylureido)phenyl)benzamide.

IR(KBr): 3416, 1642 cm⁻¹. ¹H NMR (400 MHz, DMSO-d6) δ (ppm): 2.38 (s, 3H, CH3-26), 6.96 (t, J = 8 Hz, 1H, ArH-23), 7.28 (t, J = 8 Hz, 2H, ArH-22,24), 7.33 (d, J = 4 Hz, 2H, ArH-17,19), 7.44 (d, J = 8 Hz, 2H, ArH-21,25), 7.47 (d, J = 8 Hz, 2H, ArH-6,10), 7.70 (d, J = 8 Hz, 2H, ArH-7,9), 7.88 (d, J = 4 Hz, 2H, ArH-16,20), 8.62 (s, 1H, NH-13), 8.64 (s, 1H, NH-11), 10.08 (s, 1H, NH-3). ¹³C NMR (100 MHz, DMSO-d6) δ (ppm): 21.02 (1C-26), 118.20 (2C-21,25), 118.48 (2C-7,9), 121.14 (1C-6,10), 121.76 (2C-23), 127.63 (2C-16,20), 128.80 (2C-22,24), 128.89 (2C-17,19), 132.21 (1C-1), 133.60 (1C-8), 135.52 (1C-5), 139.83 (1C-15), 141.38 (1C-18), 152.62 (1C-12), 165.02 (1C-2). Calcd. for C₂₁H₁₉N₃O₂ [M + H] + : 345.15, found: 345.2.

Cell toxicity assay

In vitro MTT assay was conducted to evaluate the cytotoxic activities of the synthesized compounds against HT-29 and A549 cells based on the method described in detail in our previous work [12].

Molecular docking studies

To evaluate the interactions between synthesized compounds and the kinase domain of VEGFR-2, a molecular docking study was performed which is explained extensively in our previous work [12].

Molecular dynamic simulation

The complex between the inhibitor and VEGFR2 (PDB code of 4ASD: the juxtamembrane and kinase domains) obtained by docking calculations was subjected to molecular dynamics (MD) simulation. The simulation was performed using the Assisted Model Building with Energy Refinement (AMBER) suite of programs running on a Linux-based GPU (NVIDIA TESLA K40) work station. The Antechamber package of AmberTools was used to create prmtop (topology and the parameters defining the force field) and inpcrd (containing atom positions) files for the ligand using GAFF force field [14]. Then, producing the input parameter files including frcmode (force field parameters) and lib (library) files, neutralizing the total charge of the system using one Na⁺ ion, and solvating the system with TIP3P water models were performed by using “leap” module of AMBER. Using the “Sander” module, the system was minimized by performing a 50 ps minimization step. Then the system was heated from 0 to 300 K during 50 ps followed by 50 ps of density equilibration with weak restraints on the complex followed by 500 ps of constant pressure equilibration at 300 K. The used box size for MD simulation is 81.6954910, 81.3204940, 87.0509170 Ǻ. All simulations were performed under periodic boundary conditions, and using the SHAKE algorithm on hydrogen atoms, with a 2 fs time step and by applying a Langevin thermostat for temperature control. The equilibrated system was subjected to a 50 ns MD simulation, and the coordinates were written out to obtain the trajectory file of the simulation. Finally, MMGBSA and NAB (Nucleic Acid Builder) normal mode analysis were applied to the trajectory to calculate ligand binding free energy.