(-)-Kusunokinin as a Potential Aldose Reductase Inhibitor: Equivalency Observed via AKR1B1 Dynamics Simulation.

(-)-Kusunokinin performed its anticancer potency through CFS1R and AKT pathways. Its ambiguous binding target has, however, hindered the next development phase. Our study thus applied molecular docking and molecular dynamics simulation to predict the protein target from the pathways. Among various candidates, aldo-keto reductase family 1 member B1 (AKR1B1) was finally identified as a (-)-kusunokinin receptor. The predicted binding affinity of (-)-kusunokinin was better than the selected aldose reductase inhibitors (ARIs) and substrates. The compound also had no significant effect on AKR1B1 conformation. An intriguing AKR1B1 efficacy, with respect to the known inhibitors (epalrestat, zenarestat, and minalrestat) and substrates (UVI2008 and prostaglandin H2), as well as a similar interactive insight of the enzyme pocket, pinpointed an ARI equivalence of (-)-kusunokinin. An aromatic ring and a γ-butyrolactone ring shared a role with structural counterparts in known inhibitors. The modeling explained that the aromatic constituent contributed to π-π attraction with Trp111. In addition, the γ-butyrolactone ring bound the catalytic His110 using hydrogen bonds, which could lead to enzymatic inhibition as a consequence of substrate competitiveness. Our computer-based findings suggested that the potential of (-)-kusunokinin could be furthered by in vitro and/or in vivo experiments to consolidate (-)-kusunokinin as a new AKR1B1 antagonist in the future.

Introduction

Chemotherapy is commonly used in the treatment of cancer. Chemotherapeutic agent development is ongoing, so as to increase treatment efficacy and reduce adverse effects.[1] The major source of synthetic chemotherapeutic agents is active plant-derived chemicals. Various types of natural products have been successfully developed into practical chemotherapeutic drugs.[2]

Lignans, a dimeric phenylpropanoid from natural products, provide a wide range of biological functions. In particular, flaxseed lignans have been revealed to have anticancer effects by multiple targets of the hallmarks of cancer.[3] Well-known cases include lignan-based anticancer drugs, namely, etoposide and other podophyllotoxin semisynthetic derivatives.[4,5] Another reported lignan, arctigenin and its glycoside, arctiin, hindered the cancer cell growth.[6,7] Arctigenin prevented the osteoarthritis progression via targeting of the PI3K/AKT pathway.[8] Arctiin possessed antiangiogenic effects on human myeloma cells via inhibition of STAT3 phosphorylation.[9] In addition, arctigenin and arctiin also exhibited anti-inflammatory effects through the inhibition of reactive oxygen species (ROS) production and COX-2.[10]

(−)-Kusunokinin, one of the reported lignan compounds, exhibited anticancer activity on various human cancer cell lines, including breast cancer, colon cancer, and cholangiocarcinoma.[11] This compound inhibited tumor growth with no side effects and exhibited the normal complete blood count parameters and clinical chemistry of the renal and liver function test in the breast cancer rat model.[12] The anticancer mechanism of (−)-kusunokinin could occur through the partial binding with CSF1R, resulting in the suppression of AKT and its downstream proteins, including c-myc, cyclinB1, and cyclinD1.[12,13]

Computational modeling of (−)-kusunokinin was carried out to investigate its protein target based on a similar case of bursehernin,[13,14] a structural isomer of (−)-kusunokinin. Colony stimulating factor 1 receptor (CSF1R) was reported as a possible target of these compounds, in which the stabilized π–π interaction with aromatic amino acids played a crucial role. This observation suggested that the selectivity of (−)-kusunokinin depended on the π–π stacking with aromatic amino acids at the binding pocket, such as tyrosine, phenylalanine, or tryptophan. In addition, the (−)-kusunokinin structure could perform interactions via various hydrogen-bond acceptor counterparts, namely, an oxygen atom(s) of the methoxy group (−OCH3) or/and an oxygen atom at the ketone group on the γ-butyrolactone ring.[13,14]

As we had already noticed some clues in the interaction and specificity of (−)-kusunokinin, we thus speculated that another protein containing an aromatic-rich domain could behave with (−)-kusunokinin in the same manner. Moreover, due to the partial binding with CSF1R, the downstream molecules of the CSF1R and AKT signaling pathways were selected for screening for the possibility of (−)-kusunokinin target proteins in this study. Herein, we used the molecular docking technique to predict an additional target protein of (−)-kusunokinin. The time-dependent conformational behavior of the (−)-kusunokinin–target complex was elucidated using molecular dynamics simulation. Then, the binding manner of (−)-kusunokinin was also evaluated and compared with the reported inhibitors. Finally, our present work provided more information regarding the target protein of (−)-kusunokinin with aldo-keto reductase family 1 member B1 (AKR1B1) and the downstream molecules of CSF1R and AKT. (−)-Kusunokinin bound this protein at the anion binding pocket. These findings could become constructive for furthering the (−)-kusunokinin development phase as a specifically targeted anticancer inhibitor.

Methods

Protein Structure Preparation

Structures of 41 proteins associated with the CSF1R and AKT pathways were retrieved from the Protein Data Bank (https://www.rcsb.org/) in a PDB format file. All atoms were selected based on the highest occupancy. Crystallographic waters were removed and polar hydrogen atoms were added using AutoDockTools 1.5.6 (ADT).[15] The protein structures were written into a PDBQT format file.

Ligand Preparation

The structure files of (−)-kusunokinin and the other ligands used were obtained from the PubChem database. (−)-Kusunokinin was taken from a structure with CID 384876. The crystal-bound (native) ligands were directly separated from the crystallography structure using the Visual Molecular Dynamics (VMD) package.[16] All ligands used in this study are summarized in Table S1. The structure was then converted into a PDB file using the Biovia Discovery Studio 2019 Client.[17] Hydrogen atoms were added to all ligands. Gasteiger charges and default atom parameters were automatically assigned from ADT. Finally, the structure was written into the PDBQT file format using ADT.

Molecular Docking Parameters

Molecular docking studies were performed using the AutoDock4 program. During the process, the protein structure was set as a rigid molecule with a flexible ligand. The docking site on the protein target was defined by a grid box of 120 × 120 × 120 cubic angstrom (Å3), with the grid spacing at the center of the protein structure. Other parameters followed the default values in ADT. Fifty genetic algorithm (GA) runs with a population size of 200 were executed for conformational sampling. The conformation with the lowest binding energy was considered as the best pose. The interaction of the best docked pose was compared with the native ligands in the crystallographic structure.

Molecular Docking Study and Binding Energy Calculation

The reproducibility of the AutoDock4 program was validated via a self-redocking experiment of the native reference ligand. The all-atom root-mean-square deviation (RMSD) cutoff of less than 2 Å between docked and native crystal structures was set as the criterion for a successful procedure. The AutoDock4 program score was recorded from the binding free energy between the ligand and the protein.[18,19]

Molecular Dynamics Simulation

The AKR1B1 structure selected was the PDB structure code 1US0 (a resolution of 0.66 Å). All water and other solvents were removed. Only AKR1B1 with the cofactor dihydro-nicotinamide-adenine-dinucleotide phosphate (NADPH) was retained. All polar hydrogen in the best ligand docked pose was eliminated. AMBER16 force field was exploited to model the drug–AKR1B1 complex, similar to previous reports of AKR1B1 molecular dynamics study.[20−23] The hydrogen atom was instead replaced via the leap module implemented in the AMBER16 package.[24] An ionizable side chain of the amino acid was considered at pH 7. The side chain in lysine or arginine contained +1e charge, while the glutamate and aspartate side chains contained −1e charge. The imidazole ring in histidine was neutral. No disulfide was observed in the AKR1B1 structure. All ligand molecules, namely, (−)-kusunokinin, AKR1B1 inhibitors, and a substrate (PGH2), were parameterized using the Gaussian16 package,[25] based on calculations of restrained electrostatic atomic potential (RESP). NADPH parameters were adopted from a previous report.[26]

The structure of ligand-free AKR1B1 or AKR1B1 with the best docked conformer was solvated by TIP3P water with a distance of 14 Å from the surface, leading to approximately 12000 water molecules. Sodium ion (Na+) was used to neutralize the system, and 18 NaCl pairs were included to yield a 0.1 M NaCl solution. The aforementioned AKR1B1–ligand complex was used as an initial coordinate for the molecular dynamics process. The AKR1B1 solution was first equilibrated in a canonical (NVT) ensemble using Langevin dynamics at 310 K (37 °C). The harmonic restraint potential was applied to the positions of AKR1B1, NADPH, and the ligand, with a force constant of 200, 100, 50, 20, and 10 kcal/mol-Å2, respectively. Each restrained NVT ensemble lasted for 200 ps using a time step of 1 fs. The system was switched into an NPT ensemble via the Berendsen algorithm at 310 K and 1.013 bar.

The molecular dynamics (MD) simulation was run for 150 ns with a time step of 2 fs. An MD trajectory with 2500 equidistant snapshots was obtained. The last 60-ns-MD conformation (1500 snapshots) was taken for conformation and binding analysis. Root-mean-square deviation (RMSD) calculation and structure visualization were conducted using the VMD package. Root-mean-square fluctuation (RMSF) was computed, and the pattern of distance geometry was acquired using the cpptraj module.

Binding Free Energy Evaluation

Average binding free energy (ΔG) was calculated using the molecular mechanics/generalized Born-surface area (MM/GBSA)[27−30] method from MD trajectories. Average molecular mechanics energy was obtained from the AMBER functional formula.[18,28,31]

Results and Discussion

(−)-Kusunokinin Target Screening

From the selected 41 CSF1R and AKT associated proteins, the top 10 for the binding energy of (−)-kusunokinin are shown in Figure . All binding energies of (−)-kusunokinin and the native ligands with 41 selected proteins are also summarized in Table S2. (−)-Kusunokinin exhibited the best predicted binding potential toward aldo-keto reductase family 1 member B1 (AKR1B1), with a binding energy of −11.11 kcal/mol. The other candidate was the dual specificity mitogen-activated protein kinase kinase2 (MEK2), with a binding energy of −10.14 kcal/mol.

Figure 1

Binding energy of (−)-kusunokinin with protein candidates in the CSF1R and AKR1B1 pathways. A number denotes the binding energy in kcal/mol.

In spite of the closeness in the predicted binding energy between AKR1B1 and MEK2, AKR1B1 piqued our interest not only for its lower binding energy but also due to some supporting data regarding its lignan-based inhibitor. AKR1B1, an NADPH-dependent enzyme, is a target for diabetic complication treatment through the regulation of the polyol metabolic glucose pathway.[32] Recently, AKR1B1 was found in various human cancers,[33] and it could be the novel target for cancer aggressiveness reduction.[34−36] Arctigenin and arctiin exhibited potency in the aldose reductase (AKR1B1) inhibition.[37−39] Since the lignans share a feature for the aromatic ring structure, the π–π stacking with the aromatic binding site of AKR1B1 would then be expected to be like CSF1R.[13,14]

Preliminary Atomistic View of (−)-Kusunokinin on AKR1B1 Binding

(−)-Kusunokinin was suitably inserted in the AKR1B1 active pocket cleft, which was surrounded by several hydrophobic and aromatic residues (Figure ). The site was categorized into three apparent pockets: an anion binding pocket (Tyr48, His110, and Trp111), a small hydrophobic pocket (Trp20, Val47, and Tyr48), and a large hydrophobic specific pocket (Trp79, Trp111, Phe115, Phe122, Trp219, Cys298, Leu300, Cys303, and Tyr309). These active pockets resembled previously reported AKR1B1 binding with β-aminophenylpropanoic acid derivatives.[23]

Figure 2

(−)-Kusunokinin binding site on AKR1B1. The site is composed of anion binding pocket (orange), small hydrophobic pocket (green), and large hydrophobic pocket (yellow).

Later, a molecular docking study of (−)-kusunokinin against AKR1B1 was performed and compared with the other available AKR1B1 substrates,[40,41] commercial carboxylic acid aldose reductase inhibitors (ARIs), commercial hydantoin ARIs, and some suggested potential ARIs, respectively.[42−44] The result is summarized in Table . Interestingly, the binding energy of (−)-kusunokinin (−11.11 kcal/mol) was in a similar range to the reported ARIs, from −6.92 to −11.54 kcal/mol. Moreover, (−)-kusunokinin provided a better predicted energy than substrates (PGH2 and retinol), substrate analogs (UVI2008 and PGA1), and some reported potential ARIs. (−)-Kusunokinin was likewise inserted in the active pocket with a similar position and orientation to the well-known ARIs (Figure ).

Table 1

Predicted Binding Energy of Docked Ligands with AKR1B1

ligands	ΔG (kcal/mol)	ligands	ΔG (kcal/mol)
(−)-kusunokinin	–11.11	suggested potential ARIs
carboxylic ARIs		oleanolic acid	–10.90
zenarestat	–11.54	AD5467	–10.84
zopolrestat	–11.21	γ-mangostin	–9.63
sulindac	–11.19	arctiin	–9.85
lidorestat	–10.92	arctigenin	–9.55
IDD1219	–10.86	10C	–9.32
epalrestat	–10.57	caffeic acid phenethyl ester	–9.08
IDD740	–10.52	substrates
IDD594	–9.83	Retinol	–9.91
IDD552	–9.53	PGH2	–9.71
tolrestat	–9.38	substrate analogs
tolmetin	–9.22	UVI2008	–10.71
alrestatin	–8.92	PGA1	–10.01
hydantoin ARIs
minalrestat	–10.88
ranirestat	–10.41
fidarestat	–7.50
sorbinil	–6.92

Figure 3

(−)-Kusunokinin and ARI interaction with AKR1B1. Interaction of ARIs include epalrestat (orange), zenarestat (blue), minalrestat (red), and (−)-kusunokinin (lime green). Tyr48 and His110 bound ligands using hydrogen bonds, while Trp20 and Trp111 were mainly responsible for π–π interactions. The π–π interactions are shown as black dashed lines, and the hydrogen bonds are shown as red dashed lines.

Trp20 and Trp111 sustained ARI conformation by forming π–π stacking with the aromatic rings of ARIs. Furthermore, Tyr48 and His110 also formed hydrogen bonds with the AKR1B1-bound ligand. These interactions, found in most ARIs, are concordant with previous reports that Trp20, Tyr48, His110, and Trp111 are present in the anion binding pocket.[32] The results gave a hint that the aromatic ring-containing compounds could be the potential ARIs due to π–π stacking formation at the active pocket.

Since molecular docking was performed on the basis of rigid proteins with a flexible ligand, the study of molecular dynamics simulation was later carried out to investigate the time-dependent characteristics of AKR1B1 and its bound compounds.

Molecular Dynamics Trajectory Analysis

Molecular dynamics simulation was performed to reflect the more realistic dynamic behaviors of the drug-free AKR1B1 or drug-AKR1B1 complex. All of the MD trajectories of all systems became stable after 90 ns (Figure a). The flexibility pattern of ligand–AKR1B1 complexes was similar mostly throughout the structure. A distinctive exception was observed in the AKR1B1 residue 210 to 230 in PGH2-AKR1B1 and UVI2008-AKR1B1, which was distinctively more flexible (Figure b).

Figure 4

MSD and RMSF of MD trajectory from the AKR1B1 simulation. The interesting simulations were the AKR1B1 and AKR1B1 structures with six docked compounds. RMSD was plotted from 150 ns simulations (a). RMSF was computed from a 90 to 150 MD trajectory (b).

The next point was to investigate whether the bound ligand can affect the AKR1B1 conformation. The distance between a center of geometry in each amino acid and an origin point was plotted to give the pattern of distance. The reference structure was ligand-free AKR1B1 MD simulation. The same distance pattern indicated a similar conformation of the assigned amino acid. However, a different distance pattern pinpointed a different conformation from the AKR1B1. Herein, similar patterns among the seven AKR1B1 simulations were observed (Figure ). The conclusive finding was that (−)-kusunokinin and other ARIs did not affect the AKR1B1 conformation.

Figure 5

Distance pattern among ligand-free AKR1B1 and drug-AKR1B1 MD simulations. The distance pattern of all simulations was similar to the ligand-free AKR1B1, giving a hint that a bound ligand showed no significant effect on AKR1B1 conformation.

AKR1B1 binding affinity of (−)-kusunokinin was considered from an average MM/GBSA binding free energy with respect to other experimentally known ARIs, which were zenarestat, minalrestat, and epalrestat (Table ). MM/GBSA calculations revealed that the binding free energy of (−)-kusunokinin (−42.23 kcal/mol) was comparable to other selected ARIs, ranging from −29.63 to −38.96 kcal/mol, and suggested that (−)-kusunokinin would bind with AKR1B1 like other inhibitors.

Table 2

MM/GBSA Binding Energy of (−)-Kusunokinin and Selective ARIs toward AKR1B1

ligand	MM/GBSA binding energy (kcal/mol)
(−)-kusunokinin	–42.23 ± 0.09
zenarestat	–29.63 ± 0.10
minalrestat	–36.55 ± 0.07
epalrestat	–38.96 ± 0.08

An oxygen atom at the ketone group on the γ-butyrolactone ring of (−)-kusunokinin interacted with His110 via hydrogen bonding. Meanwhile, the aromatic ring on 1,2-dimethoxybenzene was stabilized by π–π stacking with Trp111 throughout the MD progression (Figure a,b). The γ-butyrolactone ring acted equivalently to the carboxylic acid group of epalrestat and zenarestat to access into a key catalytic His110 in an AKR1B1 anion binding site. The hydantoin ring of minalrestat, however, formed temporary hydrogen bonds with His110. This observed finding is concordant with the reported role of His110 with carboxylic acid or hydantoin ARIs.[45]

Figure 6

Binding interaction of (−)-kusunokinin and other inhibitors/substrates with AKR1B1. Six compounds were chosen for MD simulation: (−)-kusunokinin (a), zenarestat (b), minalrestat (c), epalrestat (d), UVI2008 (e), and prostaglandin H2 (PGH2) (f). The black, red, and green dotted lines represent hydrogen bonds, π–π interaction, and π–alkyl interaction, respectively.

Furthermore, Trp111 reportedly acted as a key residue in ARI specificity.[32,46] Our study showed that an aromatic ring of (−)-kusunokinin was stabilized by π–π stacking with Trp111 (Figure a), similar to selected ARIs (Figure b–d) and UVI2008 (an AKR1B1 substrate analog) (Figure e). Although no aromatic ring was found in the PGH2 structure, Trp111 was also responsible in PGH2 binding, using π–alkyl interaction with the PGH2 hydrocarbon chain (Figure f). Besides, Tyr48 could temporarily form hydrogen bonds with (−)-kusunokinin and ARIs.

Comparative Behavior of (−)-Kusunokinin with ARIs and AKR1B1 Substrate

AKR1B1 requires NADPH cofactor to catalyze the various substrates, including prostaglandin H2 (PGH2). PGH2 is converted into PGF2α and relays the signals in the NF-κB pathway to regulate the epithelial–mesenchymal transition (EMT) process on various cancers. The catalytic imidazole ring in His110 acts as a proton donor for the dioxabicycloheptane ring of PGH2 to produce PGF2α,[40][40] as shown in Figure a.

Figure 7

His110 catalytic role in PGH2 conversion. PGH2 conversion to PGF2α using proton transfer from His110 (a). Position alignment between (−)-kusunokinin (red) and dioxabicycloheptane ring of PGH2 (blue) at the ARK1B1 catalytic site (His110) (b).

The molecular dynamics simulations showed that the PGH2 dioxabicycloheptane ring made contact with a transferable proton of NADPH and His110, as expected. The stable hydrogen bond formation with a catalytic His110 due to (−)-kusunokinin, as well as epalrestat, directly prevented PGH2 access into the catalytic site (Figure b). However, no stable hydrogen bond with His110 was noticed in the other ARIs. Instead, aromatic residue Trp111 was found to be the principal stabilizer of (−)-kusunokinin, carboxylic and hydantoin ARIs. Mutual interaction with Trp111 would be a rather major contributor to AKR1B1 inhibition because the ligand hindered PGH2 access into the catalytic site, no matter whether or not the compound could form a hydrogen bond with His110 (Figure ).

Figure 8

Proposed binding model of ARI in AKR1B1 anion binding pocket. (−)-Kusunokinin (a), carboxylate ARIs (b), hydantoin ARIs (c); red dashed lines: hydrogen bond; orange arrow: π–π interaction.

The binding behaviors of (−)-kusunokinin were closer to carboxylic ARI than hydantoin ARIs. In Figure , (−)-kusunokinin could bind AKR1B1 via its carbonyl oxygen in the γ-butyrolactone ring and aromatic ring (Figure a). The γ-butyrolactone ring resembled a carboxylic group carboxylic ARI. Carboxylic ARIs have exerted carboxylate groups into anion binding pockets,[47] including His110, similar to the carbonyl oxygen of (−)-kusunokinin (Figure a,b). This structural information could also be supported by ARI-reported lignans (arctigenin and arctiin)[37−39] because their core structures contain the γ-butyrolactone ring and aromatic entity (Figure ).

Figure 9

(−)-Kusunokinin, arctigenin, and arctiin structures.

Although many ARIs were synthesized, several ARIs were withdrawn in consequence of the severe adverse effects of hepatotoxicity[48−50] and hypersensitivity.[51,52] Only commercial epalrestat is used in Asian countries for diabetic neuropathy treatment[53] and metastatic triple-negative breast cancer in clinical trials.[54] However, its poor efficiency and adverse effects remain a problem in other countries.[55] Hence, the search for new, efficient ARIs with low toxicity is important.

Conclusions

Related to (−)-kusunokinin action in CSF1R and AKT pathways, AKR1B1 was found as a tentative (−)-kusunokinin target. Predicted binding affinity from (−)-kusunokinin was investigated and compared with other enzyme inhibitors and substrates. Interestingly, (−)-kusunokinin exhibited distinctive binding affinity toward AKR1B1 with respect to the reported AKR1B1 inhibitors (epalrestat, zenarestat, and minalrestat) and substrates (UVI2008 and prostaglandin H2). No conformational alteration in AKR1B1 was found due to compound binding. Interchangeable structural features in (−)-kusunokinin reflected an equivalency to AKR1B1 inhibitors. An aromatic ring acted as an AKR1B1 binding counterpart through π–π stacking with Trp11. An obstruction in AKR1B1 function was also expected because of the hydrogen bond formation of the γ-butyrolactone ring with a crucial catalytic His110. Due to persuasive in silico information, further experiments are encouraged to determine the AKR1B1 inhibitory effect of (−)-kusunokinin.