Computational Investigation of Structural Basis for Enhanced Binding of Isoflavone Analogues with Mitochondrial Aldehyde Dehydrogenase.

Isoflavone compounds are potent inhibitors against mitochondrial aldehyde dehydrogenase (ALDH2) for the treatment of alcoholism and drug addiction, and an in-depth understanding of the underlying structural basis helps design new inhibitors for enhanced binding. Here, we investigated the binding poses and strengths of eight isoflavone analogues (including CVT-10216 and daidzin) with ALDH2 via computational methods of molecular docking, molecular dynamics (MD) simulation, molecular mechanics Poisson-Boltzmann surface area (MM-PBSA), steered MD, and umbrella sampling. Neither the Vina scoring of docked and MD-sampled complexes nor the nonbonded protein-inhibitor interaction energy from MD simulations is able to reproduce the relative binding strength of the inhibitors compared to experimental IC50 values. Considering the solvation contribution, MM-PBSA and relatively expensive umbrella sampling yield good performance for the relative binding (free) energies. The isoflavone skeleton prefers to form π-π stacking, π-sulfur, and π-alkyl interactions with planar (Phe and Trp) or sulfur-containing (Cys and Met) residues. The enhanced inhibition of CVT-10216 originates from both end groups of the isoflavone skeleton offering strong van der Waals contacts and from the methylsulfonamide group at the 4' position by hydrogen bonding (HB) with neighboring receptor residues. These results indicate that the hydrophobic binding tunnel of ALDH2 is larger than the isoflavone skeleton in length and thus an extended hydrophobic core is likely a premise for potent inhibitors.

Introduction

Aldehydes are a family of reactive organic compounds existing widely in nature.[1] A well-known member is acetaldehyde (AcH), the first metabolite of ethanol, and is responsible for the flushes after consuming alcohol. A high level of aldehydes in the human body may cause DNA damage, cytotoxicity, and carcinogenesis.[2−6] Aldehyde dehydrogenases (ALDHs) are the most important enzymes to metabolize the aldehydes into corresponding carboxylic acid derivatives for the relief of aldehyde stress.[7,8] The activity of ALDHs is essential for protecting the human body from the toxic influences of various exogenous and endogenous aldehydes in vivo.[9,10] Human ALDH superfamily contains 19 NAD(P)+-dependent enzymes displaying similar but not identical functions due to the difference in gene expression and substrate specificity.[11−13] ALDH2, a mitochondrial isozyme of ALDHs, is the most efficient one for the oxidative elimination of AcH via catalyzing it to nontoxic acetate.[14] It also participates in the metabolization of other short-chain aliphatic, aromatic, and polycyclic aldehydes,[15] and relates to a variety of human pathologies such as drug addiction.[16−19]

The number of cancer deaths caused by alcoholism has increased rapidly in recent years, and alcoholism was believed to relate to a number of genetic predisposing factors.[20,21] About 8% of the global population (mostly in East Asia, ca. 40%) carries an ALDH2*2 allele that encodes a nonfunctional ALDH2 enzyme.[22−24] As a result, AcH cannot be eliminated promptly and efficiently, and a high level of blood AcH leads to facial flushing, nausea, headache, and cardiac palpitations.[25] The individuals carrying this allele have a reduced ability of metabolizing AcH, and the resulting discomfort after alcohol intake endows them with a low danger of alcoholism. Moreover, an inhibition of ALDH2 activity was shown to eliminate the reinstatement of alcohol seeking for the rats even when alcohol is absent (that is, no acetaldehyde). This can be attributed to the fact that alcohol-induced dopamine levels in the central nervous system were downregulated by the ALDH2 inhibitor,[26] and similar findings were reported upon exposure to methamphetamine and cocaine.[17,18] These observations indicate that ALDH2 is a promising target for the suppression of heavy drinking and drug addiction.

A few pharmacotherapies are available for the treatment of alcoholism.[27−29] As the first approved medicine, disulfiram targeted both cytoplasmic ALDH1 and mitochondrial ALDH2 and inhibited the former more potently than does the latter.[30] Besides the interaction with ALDHs, however, it was involved in adverse drug reactions on its own leading to toxic hepatitis and severe metabolic disorders.[31−34] Moreover, disulfiram is nonspecific and is able to inhibit other enzymes such as the dopamine β-hydroxylase[35] and phosphoglycerate dehydrogenase (serine synthesis).[36] Different from disulfiram, other approved anti-craving medicines of, for instance, acamprosate, nalmefene, and naltrexone did not target ALDHs; their mechanisms of action refer to a review by Wang et al.[29] Unfortunately, these medicines have poor bioavailability and are often accompanied by a variety of unwanted side effects.

Naturally occurring isoflavones appear as a promising agent against alcohol abuse.[37−39] Daidzin, an active isoflavone, was proved to be a potent and selective inhibitor of ALDH2.[40−42] Keung and co-workers synthesized a number of daidzin analogues and suggested that a potent inhibitor required a free hydroxyl (−OH) group at the 4′ position of the isoflavone skeleton and a straight-chain alkyl substituent with a polar end like carboxyl (−COOH) at the position 7 (Figure a).[37] The compounds with a chain length of 5–10 performed best with an IC50 of ca. 0.05 μM, and daidzin showed a slightly low inhibition (IC50 = 0.08 μM). They then crystallized the ALDH2/daidzin complex, providing a structural basis for the interaction mechanism.[43] Inspired by the co-crystal structure, a highly selective ALDH2 inhibitor, CVT-10216 (also known as GS 455534), was synthesized with an IC50 of 0.029 μM.[26] Design of other selective inhibitors under development such as GS-548351 and its prodrug ANS-6637 for substance (alcohol and other drugs) use disorders were aided by the structure of daidzin in complex with ALDH2 as well,[44] although they do not have an isoflavone skeleton. Besides the suppression of alcohol abuse, preclinical findings in the rats indicated that CVT-10216 played a positive role in a number of neurobiological behaviors such as binge eating[45] and anxiety[46] as well as cocaine[18] and nicotine[47] intake. This is likely due to the downregulated neurotransmitter (like dopamine and serotonin) levels by the inhibition of ALDH2 activity.[17,18,26]

Figure 1

Substituted groups of R1 and R2 at both ends (4′ and 7 positions) of the isoflavone skeleton (a–c) and experimentally determined IC50 values (d) for the investigated inhibitors of ALDH2. MSA-1 stands for the known inhibitor CVT-10216, and OH-3 is for daidzin. A, B, and C indicate the ring groups in the skeleton (a).

Given the excellent profile of isoflavones, an in-depth understanding of the structural basis for the interactions with ALDH2 is helpful for the design of potential inhibitors or activators. Structure-based screening of chemical compounds was attempted by multiple studies,[48−51] while the structural basis for the enhanced binding of isoflavone analogues is yet to be clarified. The enhanced inhibition of CVT-10216 probably originates from the substitution of a methylsulfonamide group for 4′-hydroxyl of the isoflavone skeleton (Figure b). Here, we aim to explore the binding mechanism of isoflavone analogues (including CVT-10216 and daidzin) with ALDH2 at a molecular level and investigate the substitution effects of both ends (4′ and 7 positions) of the isoflavone skeleton on the binding pose and binding strength. A variety of computational approaches were examined such as molecular docking, steered molecular dynamics (SMD), molecular mechanics Poisson–Boltzmann surface area (MM-PBSA), and umbral sampling simulations. We identified the key residues involved in the protein–inhibitor interactions by the binding energy decomposition. An assessment of these approaches is addressed by comparing the relative binding strength with experimental IC50 values. This work provides a molecular basis for the design of new compounds with enhanced inhibition against ALDH2.

Computational Methods

Structural Model

The three-dimensional (3D) structure of human mitochondrial aldehyde dehydrogenase (ALDH2) was extracted from the Protein Data Bank (PDB code: 2VLE) with a resolution of 2.4 Å where ALDH2 binds to the isoflavone daidzin and forms a tetramer.[43] Each subunit contains three domains; two of these domains are for binding with the coenzyme (or cofactor) and substrate, and the third one is the oligomeric domain (residues 140–158 and 486–495, forming a three-stranded antiparallel β-sheet). The cofactor NAD+ is absent in the protein 2VLE,[43] and we obtained an NAD+-bound complex via an alignment of protein backbone atoms with the chain A of protein 1CW3.[52] Molecular structures of the investigated eight inhibitors (MSA- and OH-, –4) are given in Figure . MSA and OH are short for the methylsulfonamide and hydroxyl groups, respectively, and are used to modify the R1 terminal (at the 4′ position) of the isoflavone skeleton (Figure a,b); is for modifying the R2 terminal (at the 7 position, Figure c). MSA-1 (known as CVT-10216), OH-2, OH-3 (daidzin), and OH-4 are synthetic or naturally occurring compounds, and the experimental IC50 values for ALDH2 are given in Figure d.[26,37] Other four inhibitors (OH-1, MSA-2, MSA-3, and MSA-4) are virtual compounds designed for investigating the substitution effect. Via aligning the isoflavone skeleton with OH-3, we obtained the protein/inhibitor/cofactor configurations for the other seven inhibitors and used them as the initial coordination in the following simulations. The alignments for protein backbone and isoflavone skeleton were done using the GROMACS tool “gmx confrms”.[53]

Simulation Protocol

The residues in the ALDH2 dimer interface were reported affecting the stability of catalytic and coenzyme binding domains,[52,54−56] and here we chose a dimer for MD simulations. The protein/inhibitor/cofactor complex was immersed in a cubic box with a length of 9.6 nm, and water molecules and Na+ and Mg2+ ions in the crystal structure were kept in place. Each system is neutral and contains two complexes, two Mg2+ ions, 10 Na+ ions, and approximately 25 000 water molecules. We optimized nine ligands (eight inhibitors in Figure and the substrate acetaldehyde) at HF/6-31G* in the gas phase with Gaussian 09 software[57] and then calculated the restrained electrostatic potential (RESP) charges with the aid of the “antechamber” tool.[58] The Amber99SB-ILDN force field[59] was used to model ALDH2 and ions, and the general Amber force field (GAFF)[60] was chosen for the ligands. Protonation states of titratable protein residues were assigned automatically by the GROMACS utility of “gmx pdb2gmx” at neutral pH,[53] and the assignment of His residues was done via an optimal hydrogen bonding (HB) conformation. Force field parameters of the cofactor NAD+ were taken from the AMBER parameter database collected by Prof. Richard Bryce (http://research.bmh.manchester.ac.uk/bryce/amber).[61,62] The rigid TIP3P model[63] was used to model water molecules using the SETTLE algorithm.[64] All chemical bonds were constrained using the LINCS algorithm,[65] allowing a time step of 0.002 ps. The velocity-rescaling[66] and Parrinello–Rahman algorithms[67,68] were applied to couple the temperature at 298.15 K and the pressure at 1 bar with coupling time constants of 1 and 5 ps, respectively. The particle-mesh Ewald (PME) method was used to calculate the electrostatic interactions,[69,70] and short-range nonbonded interactions were cut off at 1.0 nm.

For all simulated systems, we implemented an energy minimization first to eliminate possible bad contacts, followed by 100 ps NVT and then 400 ps NPT equilibrations. During the equilibration stages, position restraints were exerted on the protein backbone atoms using a harmonic potential with a force constant of 1000 kJ/(mol nm2). Production simulations were extended to 30 ns at NPT without any position restraints. To compute the entropy of the inhibitor in the unbound (free) state, each inhibitor was simulated in a box with a length of 4 nm (containing ca. 2150 water molecules) at NPT for 15 ns, and the last 10 ns was used for entropy calculation via the GROMACS tools of “gmx covar” and “gmx anaeig”. All of the simulations for the dimers and isolated inhibitors were performed with GROMACS 2018 software.[53]

Molecular Docking

We docked the ligands into the binding site of ALDH2 using the Autodock Vina toolkit.[71] The searched space (4 × 4 × 4 nm3) is centered roughly on the binding site of the receptor. In the docking, the ligand was always completely flexible, and two cases of the receptor ALDH2 (rigid or partially flexible) were examined. For a flexible receptor, the side chains of amino acid residues within 0.5 nm of the ligand in the crystal structure of the ALDH2/OH-3 complex were allowed to be moveable. For each ligand, the docking was run 100 times with explicit random seeds, and the best pose with the highest binding affinity for each run was collected for analysis. A preliminary test on the OH-3 (daidzin) system validated that our docking protocol predicted binding poses consistent with the crystal complex (Figure S1 in the Supporting Information, SI). For comparison with the docked poses, we also used the Autodock Vina to score the generated ALDH2/inhibitor complexes from the MD simulations of ALDH2 dimers above.

MM-PBSA Analysis

The molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) is a popular method to estimate the binding energy between the receptor and ligand.[72,73] It is, in principle, more accurate than most of fast scoring approaches in the docking and has a relatively lower computation load than the free energy calculations with an explicit solvent.[74] After 30 ns MD simulations, we stripped water molecules and Na+ and Mg2+ ions and extracted 100 conformations of protein/inhibitor/cofactor complexes from the last 10 ns trajectory with an interval of 100 ps. The “g_mmpbsa” toolkit[73] was used to compute the binding energy (ΔEbind) that includes the contributions from van der Waals (vdW) and electrostatic interactions as well as the polar (ΔGpolar) and nonpolar (ΔGnonpolar) solvation energies. Together with an entropy contribution (−TΔS), one obtained the binding free energy (ΔGbind), as in eq where ΔEMM is the sum of van der Waals and electrostatic contributions. The built-in APBS software[75] was used to compute ΔGpolar, and the nonpolar part was calculated using a solvent-accessible surface area (SASA) model.[73] The MmPbSaStat.py and MmPbSaDecomp.py scripts distributed at http://rashmikumari.github.io/g_mmpbsa/Usage.html were utilized to compute the binding energy and decompose the energy contributions atomically, respectively. Standard deviations for the energies in the MM-PBSA analysis were computed by block averaging via dividing the trajectories into five blocks.

The “g_mmpbsa” toolkit did not include a module for the entropy (ΔS) calculation, which is a challenging and time-consuming issue.[73] Due to the large computational cost, this term is neglected in most cases of practical applications.[74] Here, we attempted to calculate the entropy of inhibitors from the covariance matrices of atomic fluctuations via the quasiharmonic approximation[76] and the Schlitter formula.[77] The entropy change (ΔS) was defined as the difference between the bound and unbound (free) states. The entropy of the receptor and the surroundings were not computed because of the large size of protein and a large number of water molecules. The entropy of inhibitors was used only for quantifying the configurational changes upon complexation and was not added to eq . Instead, the binding energy (ΔEbind) was used for comparing different inhibitors against ALDH2.

Potential of Mean Force (PMF) Calculation

To reduce the computation load, we chose the ALDH2 monomer (chain A) to compute the potential of mean force (PMF, i.e., binding free energy) profiles between the receptor and the ligand. The oligomeric domain (residues 140–158 and 486–495), as well as the αG helix and the loop at residues 463–478 of ALDH2, for the stability of catalytic and coenzyme binding domains[52,54−56] were always position-restrained in the simulations to maintain a crystal-like structure. The protein/inhibitor/cofactor complex was rotated making the isoflavone skeleton parallel to the z-axis (Figure ), and was then placed in a box of 9 × 8 × 10 nm3. Each system contains one complex, one Mg2+ ion, five Na+ ions, and approximately 21 000 water molecules. Steered molecular dynamics (SMD) simulations[78] were carried out at NVT to pull the inhibitor away from the binding site with a rate of 0.005 nm/ps, using a harmonic potential with a force constant of 1000 kJ/(mol nm2). The pulling setup is similar to the ALDH2 dimer simulation, except for the use of a semi-isotropic Parrinello–Rahman barostat.[68] SMD simulations were run for 700 ps and the inhibitor sampled 3.5 nm. The center of the mass distance between the ligand and the ALDH2 oligomeric domain along the z-axis was defined as the reaction coordinate (ξ, Figure ).

Figure 2

Definition of the reaction coordinate (ξ): the center of mass distance along the z-axis between the ALDH2 oligomeric domain (colored in blue) and the inhibitor (stick model colored by elements). The ALDH2 monomer is shown with secondary structures, and the cofactor NAD+ and the inhibitor daidzin (OH-3) are represented with balls as well as the crystals of Na+ and Mg2+ ions.

Along the dissociation process in the SMD simulations, about 70 configurations were chosen with a step of ξ = 0.05 nm and used in the umbrella sampling simulations for PMF calculations. Each configuration (known as umbrella window) was simulated for 1 ns with a similar protocol to the SMD simulations, except that the pull rate was set to zero. Following the same scheme, we carried out SMD and PMF simulations for the eight inhibitors and ALDH2’s substrate acetaldehyde (AcH). After discarding the first 100 ps for equilibration, the weighted histogram analysis method (WHAM) was applied to construct PMFs,[79,80] and the resulting PMFs were set to zero at ξ = 4.5 nm where the protein–ligand interaction vanishes. Statistical errors of PMFs were estimated by the Bayesian bootstrapping of complete histograms.[80] With a cylinder approximation,[81−83] we can compute the binding constant (Ka) by integrating the PMFs (eq ) and then the corresponding thermodynamic parameters (eqs –5)where NA is the Avogadro constant, r(ξ) is the average radius of the cylinder for the sampled volume of the inhibitor movement in the X–Y plane at ξ, πr(ξ)2 is the sampled area, ΔG(ξ) is the PMF profile, R is the ideal gas constant, T is the absolute temperature, and C0 is the standard concentration of 1 mol/L.[82] The integration over ξ in eqs –4 was calculated by a trapezoidal algorithm. PMF simulations were carried out using the GROMACS package (version 4.5.5).[53] For more details of SMD and PMF simulations as well as the calculation of binding thermodynamics, refer to our previous reports.[84−86]

Results and Discussion

ALDH2 Dimer Simulation

MD simulations of the ALDH2 dimer with the eight inhibitors allow equilibration of protein–ligand binding poses upon complexation. The root-mean-square deviations (RMSDs) of the protein backbone from the crystal structure for the ALDH2 dimer and monomers are presented as a function of simulation time (Figure ). RMSDs tend to be stable after 20 ns, and the last 10 ns was selected for data analysis. The ALDH2 dimer displays a large RMSD of 0.25 nm when bound to MSA-1 (Figure a), while a relatively small fluctuation is observed for the binding with OH-3 (RMSD ∼ 0.16 nm, Figure b). The two monomers (chain A and chain B) display a different structural fluctuation, and chain B appears more stable than chain A for all of the inhibitors except OH-4 with a slightly larger RMSD for chain B than chain A by 0.02 nm (Figure and Table S1). A smaller RMSD means a structure much closer to the crystal structure, and we, therefore, choose the monomer with a smaller RMSD for subsequent analysis. The chosen monomer for MSA-1 binding has an RMSD of 0.15 nm, and the RMSDs amount to 0.12 nm, on average, for the other inhibitors (Table S1, SI).

Figure 3

Root-mean-square deviation (RMSD) of the protein backbone from the crystal structure for the ALDH2 dimer (black) and the monomers of chain A (red) and chain B (green) when bound to MSA-1 (a) and OH-3 (b).

The short-range interaction energy between the binding partners contains the contributions from electrostatic and van der Waals (vdW) interactions and is often used as an indicator for the evaluation of relative binding strengths. The vdW part plays a major role in the binding of ALDH2 with the tested inhibitors and contributes 70–90% of the energies (Table ). The substitution of methylsulfonamide (MSA) for 4′-hydroxyl (OH) leads to a decrease (more negative) in the vdW energies (ΔEvdW) by 29, 7, 17, and 28% for the inhibitors with 1, 2, 3, and 4 groups at the 7 position of the isoflavone skeleton, respectively, as revealed by a comparison of OH- with MSA- (–4). Such substitution does not affect the electrostatic contribution (ΔEelec) significantly for the compounds with either glucose (3) or carboxyl (4) groups, whereas more negative values of ΔEelec are observed for the groups of 1 and 2. As a result, the total binding energies (ΔEtotal) become more negative by 13–34% and the substitution leads to a stronger binding (Table ).

Table 1

Short-Range Interaction Energies (kJ/mol) between ALDH2 and the Inhibitor Obtained from MD Simulationsa

compound	ΔEelec	ΔEvdW	ΔEtotal
MSA-1	–47.3 ± 2.9	–218.0 ± 3.6	–265.4 ± 1.7
OH-1	–34.6 ± 2.8	–169.1 ± 1.0	–203.7 ± 3.1
MSA-2	–76.3 ± 6.5	–196.6 ± 2.1	–272.9 ± 4.6
OH-2	–19.3 ± 1.1	–183.9 ± 3.6	–203.2 ± 2.6
MSA-3	–72.5 ± 4.7	–205.0 ± 2.8	–277.5 ± 7.1
OH-3	–71.3 ± 1.4	–175.1 ± 0.5	–246.3 ± 1.1
MSA-4	–65.1 ± 9.3	–192.2 ± 1.4	–257.4 ± 8.9
OH-4	–66.3 ± 10.4	–150.6 ± 1.7	–216.9 ± 11.2

The monomer with a smaller RMSD (Table S1) was used for the calculations (i.e., chain A for OH-4 and chain B for others).

The experimental IC50 values for MSA-1, OH-2, and OH-3 fall in the same range (Figure d),[26,37] implying a binding strength of the order MSA-1 ≥ OH-2 ≥ OH-3 > OH-4. Binding energies for these four inhibitors with ALDH2 are −265.4, −203.2, −246.3, and −216.9 kJ/mol (Table ), respectively. The energies are inconsistent with the IC50 order, and the binding strength of OH-2 appears to be somewhat underestimated, which indicates that the short-range interaction energy between the binding partners is not a good indicator.

Docking

Molecular docking is an efficient and popular method to estimate the binding affinity of protein–ligand complexes for drug screening. We carried out docking calculations with flexible inhibitors and with either rigid or partially flexible receptors. For a flexible receptor, the amino acid residues within 0.5 nm of the ligand can be moveable, which mimics the induced structural fits upon protein–inhibitor binding. Using a flexible receptor results in an increase in the binding affinities by 8% (OH-4) to 30% (MSA-3), and the increment is larger for MSA- than that of OH- (Table ). This indicates that the flexibility of MSA and glucose groups appear to undergo and induce conformational arrangements of the binding partners, giving rise to an enhanced binding. Similar to the observation from short-range interaction energies during MD simulations (Table ), the substitution of MSA for 4′-OH produces an enhanced binding when considering the receptor flexibility (Table ). We note that the substitution effect is not obvious when it comes to a rigid receptor.

Table 2

Calculated Binding Affinities (kJ/mol) of the Inhibitors against ALDH2 with Vina Software

	crystal structurea
compound	rigid	flexible	MDb
MSA-1	–43.5 ± 2.8	–52.9 ± 0.6	–39.8 ± 0.7
OH-1	–44.4 ± 2.9	–50.1 ± 0.7	–38.3 ± 0.5
MSA-2	–37.5 ± 2.7	–44.8 ± 0.8	–29.6 ± 0.5
OH-2	–38.7 ± 1.6	–42.9 ± 0.6	–29.5 ± 0.3
MSA-3	–40.0 ± 1.9	–52.0 ± 0.4	–37.9 ± 0.5
OH-3c	–38.1 ± 2.2	–48.3 ± 2.2	–34.8 ± 0.3
MSA-4	–41.2 ± 3.1	–47.4 ± 0.5	–35.2 ± 0.4
OH-4	–40.3 ± 3.3	–43.6 ± 0.3	–33.9 ± 0.7

Docking with a crystal structure of rigid or partially flexible receptors. The residues within 0.5 nm of OH-3 (daidzin) in the protein 2VLE were set to be moveable in the flexible docking. Each docking was run 100 times with explicit random seeds.

The protein–inhibitor complexes generated from the last 10 ns of MD simulations were scored by Vina.

The complex in the crystal state gives binding affinities of −34.7 and −32.9 kJ/mol for chain A and chain B, respectively.

However, a partially flexible receptor does not necessarily mean good. In the crystal structure, ALDH2 binds to OH-3 with a binding affinity of −33.8 kJ/mol (averaged over two monomers), close to the prediction using a rigid receptor (Table ). A flexible ALDH2 shows a binding affinity of −48.3 kJ/mol with OH-3, a bit too strong compared to the crystal state. Interestingly, equilibrated ALDH2/OH-3 complexes from MD simulations yield a prediction of −34.8 kcal/mol, consistent with that in the crystal complex (Table ). This is likely due to the use of an explicit solvent in the equilibration of protein–ligand binding poses. Docking calculations using either crystal or MD structures predict a weaker binding for OH-2 with ALDH2 than that for OH-3 and OH-4. This is opposite to the inhibition ability (experimental IC50 values in Figure d). These observations, together with the short-range interaction energies (Table ), indicate the necessity of an improved simplification for solvation effects in the binding affinity predictions.

MM-PBSA Calculation

The PBSA method allows a more accurate consideration of solvation contribution to the protein–ligand binding, despite a relatively high computational load. The MM-PBSA calculations were performed using the last 10 ns of MD simulations, and the MM part contains two contributions from bonded and nonbonded interactions for the binding partners.[73] We used a single trajectory to do the analysis, and the bonded contributions therefore amounted to zero. Different from the short-range interaction energies in Table , the nonbonded contributions to the MM part (ΔEMM = ΔEvdW + ΔEelec) were calculated in a vacuum without cut-off.[73] Note that these two cases used identical trajectories for the calculation. The resulting ΔEMM values (Table S2, SI) in the MM-PBSA analysis are more negative than that in Table , and vdW interactions contribute roughly 80% of the MM part. Using ΔEMM as an indicator, the relative binding strength of MSA-1 ≥ OH-2 ≥ OH-3 > OH-4 is well reproduced (Figure a), a better performance than the use of short-range interactions (Table ).

Figure 4

Calculated MM (ΔEMM = ΔEelec + ΔEvdW, (a)) and solvation (ΔGsol = ΔGpolar + ΔGnonpolar, (b)) energies and the total binding energies (ΔEbind = ΔEMM + ΔGsol, (c)) of the inhibitors against ALDH2 by MM-PBSA. The values for the energy decomposition are given in Table S2 (SI).

The solvation part contains polar and nonpolar contributions; the latter favors binding, whereas the former shows the opposite (Table S2). This part amounts to 120–190 kJ/mol, which is on the same order of magnitude as the ΔEMM contribution but with a different sign (Figure b). The total solvation energies are positive, disfavoring the binding, mainly due to the desolvation of binding partners upon complexation. Binding energies (ΔEbind) with ALDH2 from the MM-PBSA calculations are −132.3, −113.0, −105.8, and −88.1 kJ/mol for the inhibitors of MSA-1, OH-2, OH-3, and OH-4, respectively, in excellent agreement with the experimental IC50 order (Table S2 and Figure c). This means that ΔEbind reproduces the relative binding strength accurately and is a good indicator. The enhanced binding for the substitution of MSA for 4′-OH is observed in the MM-PBSA analysis, as revealed by ΔEMM and ΔEbind (Table S2). Although MSA-2 shows a stronger interaction (ΔEMM) with ALDH2 than OH-2, a large polar solvation energy (positive, disfavoring the binding) cancels out most of the ΔEMM and yields a similar binding strength for MSA-2 and OH-2.

For a more in-depth qualitative assessment of how the substituent groups affect binding, we split the inhibitor into three fragments, namely, 4′-group, 7-group, and the isoflavone skeleton, and computed the electrostatic and vdW interactions between the receptor (ALDH2 plus NAD+) and these three fragments (Table S3, SI). For both ends of the inhibitors (4′- and 7-groups), unfavorable electrostatic interactions (ΔEelec) with the receptor are observed (Figure a), whereas the vdW interactions (ΔEvdW) favor the binding for all of the inhibitors (Figure b). For the 7-group, the vdW contributions appear to be larger than or close to the electrostatic interactions (with a different sign in the values, Table S3 and Figure ), thereby favoring the receptor–inhibitor binding (Figure c). However, the vdW contributions from the OH group in the 4′ position are relatively small, and the 4′-OH disfavors binding, as indicated by the positive values of ΔEMM in Figure c. The substitution of MSA for 4′-OH leads to a significant increase in the vdW interaction with the receptor (Figure b), which endows the 4′-group (MSA) with favorable contributions (Figure c). The isoflavone skeleton contributes a lot with favorable interactions for both ΔEelec and ΔEvdW (Figure ).

Figure 5

Contribution from the 4′-group, 7-group, and isoflavone skeleton of the inhibitors to the binding with the receptor (ALDH2 plus NAD+) for the electrostatic (a) and vdW (b) interactions and the sum of both contributions (ΔEMM, (c)). The values for the energy decomposition are given in Table S3 (SI).

Binding Free Energy Profiles

We carried out SMD simulations to generate the dissociation process of the protein–inhibitor complexes and then used umbrella sampling to compute the PMF (i.e., binding free energy) profiles between the protein and inhibitor. When pulling the inhibitor away from the binding site of ALDH2, the external force used did not affect the stability of protein significantly (Figure S2 in the SI). For a complex state, a large force is needed for the dissociation of binding partners, as indicated by a sharp increase in the force in the first 100 ps of SMD simulations for MSA-1 and OH-3 (Figure a). After 500 ps, the ligand escapes from the binding site completely and enters into water, as revealed by the near-zero force (Figure a), the absence of protein–inhibitor hydrogen bonds (Figure b), the increase in the water molecules surrounding the inhibitor (Figure c), and the near-zero interaction energies between the binding partners (Figure d). This highlights the importance of hydrogen bonding and ligand desolvation in the formation of ALDH2/inhibitor complexes. OH-3 tends to form more hydrogen bonds with the receptor ALDH2 than MSA-1 (Figure b), likely due to the glucose unit at the 7 position (Figure ). In the bound state, there are ca. 12 water molecules located within 0.5 nm of the inhibitors (Figure c), indicating a necessity of an explicit consideration of water molecules. In the unbound (free) state, the number of water molecules surrounding MSA-1 and OH-3 is ca. 42 and no significant differences are observed for both inhibitors.

Figure 6

(a) Pulling force, (b) the number of hydrogen bonds between ALDH2 and inhibitor, (c) the number of water molecules within 0.5 nm of the inhibitor, and (d) the short-range interaction energy between the binding partners as a function of the simulation time for the dissociation process of MSA-1 and OH-3 from the binding site of ALDH2 during SMD simulations.

Potential of mean force (PMF) profiles for the formation of ALDH2 complexes with the eight inhibitors as well as its substrate acetaldehyde (AcH) are presented in Figure as a function of the reaction coordinate (ξ). The PMFs level off and amount to zero at ξ = 4.5 nm where the binding partners are completely separated from each other. The well depth of the PMFs correspond to the binding affinity (ΔGwell) and the optimal binding distance (ξwell) between the protein and ligand. ALDH2 has a hydrophobic tunnel for ligand binding (like the isoflavone skeleton), and additional hydrogen bonding (HB) with the ligand at both ends of the tunnel is possible. In the crystal structure, for instance, the glucose unit at the 7 position of daidzin (Figure ) displays HB interactions with Asp457, and the HB network of 4′-hydroxyl with Glu268 is mediated by one water molecule.[43]

Figure 7

Potential of mean force (PMF) profiles of the dissociation process of ligand from the ALDH2 binding site for the eight inhibitors (Figure ) and the substrate acetaldehyde (AcH) of ALDH2.

The optimal binding distances (ξwell) in the PMFs have a large diversity with ξ ranging from 0.76 (OH-4) to 1.54 (OH-2) nm, as shown in Figure and Table . Due to the smallest size of hydroxyl and carboxyl groups, OH-4 is able to penetrate into the binding tunnel deeply. MSA-1 and MSA-2 display a relatively smaller ξ than those of OH-1 and OH-2, respectively, while a larger ξ for MSA-3 and MSA-4 than those of OH-3 and OH-4 is observed, probably due to the relatively hydrophilic groups of 3 and 4 preventing the ligand from a deep penetration into the binding site. This indicates that the hydrophobic tunnel is larger than the isoflavone skeleton in length and additional hydrophobic groups at the 7 position such as the aromatic (1) and alkyl chain (2) can be included in the hydrophobic tunnel. Such finding agrees with experimental observations for evaluating the inhibition of isoflavone derivatives against ALDH2.[37]

Table 3

Calculated Binding Distance (nm) and Binding Energies (kJ/mol) from PMF Profilesa

compound	ξwell	ΔGwell	ΔG0	ΔH0	–TΔS0
MSA-1	0.94	–102(5)	–79	–100	21
OH-1	1.33	–73(8)	–57	–71	14
MSA-2	1.31	–89(6)	–69	–87	18
OH-2	1.54	–87(4)	–67	–86	19
MSA-3	1.40	–65(7)	–45	–64	19
OH-3	1.19	–67(6)	–40	–63	23
MSA-4	1.23	–92(9)	–73	–92	19
OH-4	0.76	–58(6)	–36	–56	20
AcH	1.11	–10(3)	–8	–11	3

ξwell and ΔGwell correspond to the PMF minima. Statistical errors are given in parentheses Binding thermodynamic parameters of ΔG0, ΔH0, and −TΔS0 are computed by eqs –5.

The well depths (ΔGwell) of PMFs yield binding energies of −102, −87, −67, and −58 kJ/mol for MSA-1, OH-2, OH-3, and OH-4, respectively, which are in the same order as the experimental IC50 values (Figure and Table ). The same holds for the binding free energies of ΔG0. An obvious enthalpy gain (ΔH0) is observed favoring binding, while an entropy loss (−TΔS0) cancels out 25% of enthalpy gain and disfavors the complexation (Table ). The substrate acetaldehyde (AcH) displays a weak binding with ALDH2 (ΔG0 = −8 kJ/mol), incapable of competing with the tested inhibitors. In line with the MM-PBSA results (Table S2), there is a tiny difference of 2 kJ/mol in the ΔG0 for MSA-2 and OH-2, although the binding details differ between both inhibitors (Figure and Table ).

Binding Site Identification

Given the fact that the MM-PBSA method was able to reproduce the relative binding strength (Table S2 and Figure ), we further decomposed the calculated binding energy into the contributions per residue to pinpoint the details of protein–ligand interactions and identify key residues in the binding site of ALDH2. The cofactor NAD+ is considered as a residue of the receptor. Energy decompositions for the tested inhibitors are presented in Figures and S3. We identify 17 key residues (16 protein residues plus NAD+), which contribute more than 1 kcal/mol to the binding with at least one of the eight inhibitors tested (Table S4, SI). Val120, Phe170, Phe296, and Phe459 are favorable residues for all of the inhibitors with contributions of −7, −9, −6, and −12 kJ/mol, on average, respectively (Table S4). The solvation contribution of these residues ranges from −1 to 4 kJ/mol, and the ΔEMM interaction therefore plays a major role in the binding. Although the ΔEMM part for Glu268 and Asp457 is negative (favorable), both polar residues show an unfavorable contribution (2–13 kJ/mol) for most of the inhibitors due to the high (positive, unfavorable) solvation energies of 6–32 kJ/mol (Table S4, SI).

Figure 8

Energy contribution per residue to the binding of ALDH2 with the inhibitors of MSA-1 and OH-1 (a) as well as MSA-3 and OH-3 (b). The shown residues have a contribution of ≥1 kcal/mol to the binding with ALDH2 for at least one of the eight inhibitors tested. Dashed lines indicate a value of 1 kcal/mol. The cofactor NAD+ is regarded as a residue of the receptor.

As shown in the two-dimensional (2D) diagrams for the representative binding poses (Figures and S4), the planar residues of Phe170, Phe296, and Phe459 form a π–π stacking with the isoflavone ring structure, and the sulfur-containing amino acids of Met124, Met174, Cys301, Cys302, and Cys303 stabilize the isoflavone skeleton via π–sulfur and π–alkyl interactions. The isoflavone skeleton also interacts with the planar residues of Trp177 and Phe465 via the π–π stacking, as observed for MSA-1, MSA-3, and OH-3 (Figure ) as well as OH-2 and MSA-4 (Figure S4).

Figure 9

Two-dimensional diagrams of receptor–inhibitor interactions for ALDH2 complexes with MSA-1 (a), OH-1 (b), MSA-3 (c), and OH-3 (d). ALDH2 has 500 amino acids, and the residue IDs for the inhibitor and NAD+ are 501 and 502, respectively, in our simulation. The figures were generated with Biovia Discovery studio visualizer software, and the complexes are averaged structures clustered from the last 10 ns simulations.

For the R2 group at the 7 position of the isoflavone skeleton, Keung et al. reported that a hydrophobic alkyl chain with a polar end (like group 2) or a glucose unit (the group 3) lead to potent inhibitors against ALDH2.[37,43] This is due to the vdW contacts of alkyl groups with neighboring protein residues, of which Val120 is of crucial importance for all of the inhibitors, and due to the possible hydrogen bonding of the polar ends with amino acids of, for instance, Asp457 and Glu288 (Figures and S4) with an occupancy of >0.8 (Table S5). Glu288 is located near the protein surface and a long alkyl chain for the R2 group makes the HB interaction possible, as in the binding with MSA-2 (Figure S4). In addition, the alkyl chain at the 7 position of MSA-2 provides favorable vdW contacts with Trp285 and Phe292 residues. The group 1 of the inhibitor CVT-10216 (MSA-1) favors the binding via the π–alkyl interaction between its aromatic ring with Val120 and via vdW contacts with ILE116 and Val458. The polar end (carboxyl group) of MSA-1 forms a hydrogen bond with Phe459, as indicated by the 2D (Figure ) and 3D (Figure ) interaction networks for the ALDH2/MSA-1 complex.

Figure 10

Representative binding pose of the receptor (ALDH2 plus NAD+) with the inhibitor MSA-1. Protein residues of Glu268 and Glu288 are capable of hydrogen bonding with the 4′-group of OH-4 and the 7-group of MSA-2, respectively (Figure S4, SI).

For the R1 group at the 4′ position of the isoflavone skeleton, the water-bridging HB network between 4′-OH and Glu268 in the ALDH2/OH-3 crystal structure is not detected for the inhibitor OH-, while a direct HB with Glu268 or the cofactor NAD+ is formed for OH-1 (Figure ) and OH-4 (Figure S4), respectively. The solvation contributions of Glu268 appear sensitive to the binding poses and it prefers to interact with the MSA group over the OH, as indicated by more negative ΔEMM values (Table S4). The substitution of MSA for 4′-OH appears to offer more vdW contacts, additional π–sulfur interactions, and one or two HBs with neighboring residues such as Asn169, Lys178, Gly270, Cys301, Cys302, and NAD+. Cys-involved HBs are transient and not stable, as indicated by a small occupancy of <0.01 (Table S5). The HB occupancy for other three residues amounts to 0.3–0.7, implying a relatively strong HB interaction. Although Lys178 hydrogen bonds with the 4′-group in MSA-4 (Figure S4) and has a favorable ΔEMM, it still disfavors the binding in the energy decomposition (Figure S3) with a contribution of 5 kJ/mol due to a high (positive, unfavorable) solvation energy (Table S4).

The cofactor-binding domain is adjacent to the domain of ligand binding. The energy decomposition reveals that NAD+ disfavors the inhibitor binding with a contribution of 1–5 kJ/mol because of a high (positive) solvation energy, although the ΔEMM part of NAD+ is negative (favorable, Table S4) and it may participate in the HB interaction with the R1 group at the 4′ position (Figures and S4). NAD+ tends to offer favorable interactions with the isoflavone skeleton, while it is unfavorable (positive ΔEMM) when interacting with both 4′- and 7-groups of the inhibitor (Table S3).

Note that the binding details differ from case to case. The substituent groups at the 4′ and 7 positions of the inhibitor appear not to affect the global orientation of the inhibitor in the bound states significantly but have an influence on the depth of ligand penetration into the ALDH2 binding site (Table and Figure S5). For instance, Phe459 contacts with the B-ring of the isoflavone skeleton in OH-1 via π–π interactions (Figures and 9), and the substitution of MSA for 4′-OH allows a deeper penetration of MSA-1 than OH-1. As a result, Phe459 changes to interact with the A- and C-rings of the skeleton and forms a HB with the 7-group of MSA-1 (Figures and 10). With a small group (carboxyl), OH-4 is able to enter the binding tunnel of ALDH2 deeply, thereby forming a HB with Glu268 (Figure S4).

Method Assessment

We examined eight indicators (i.e., binding energies from different methods) to assess the method performance in reproducing the relative binding strength of inhibitors against ALDH2. Four inhibitors of MSA-1, OH-2, OH-3, and OH-4 with experimental IC50 data (Figure d) were used as a reference. The Spearman rank correlation coefficient (SRCC) is equal to 0.4 for the Vina docking using either rigid or partially flexible receptors (Table ), indicating a poor performance in reproducing relative binding affinities. The receptor flexibility was argued affecting the performance of virtual screening. Equilibrium MD simulations with explicit water may generate a reasonable complex allowing structural arrangements of binding partners. For our inhibitors, however, it still fails (SRCC = 0.4) when scoring the generated complexes using Autodock Vina software.[71] The short-range interaction energy (ΔEvdW + ΔEelec) obtained from equilibrium simulations with an explicit solvent shows a poor performance as well. A full consideration of the nonbonded interactions (i.e., without cutoff; ΔEMM) between protein and ligand leads to improvement (SRCC = 0.8). Considering that the IC50 values for OH-2 and OH-3 are in the same range, ΔEMM without cutoff reproduces the relative binding strength reasonably. These five indicators can be obtained with a low computation load (Table ), which is required for the structural-based virtual screening. The poor performance is mainly due to the neglect or inappropriate treatment of solvent contribution such as the desolvation of binding partners and the resulting changes in enthalpy and entropy of pure solvent.

Table 4

Assessment of Different Methods for Reproducing the Relative Binding Strength

method	indicator	SRCCa	note	computational loadb
equilibrium MDc	ΔEvdW + ΔEelec	0.4	short-range interactions	**
Vina dockingd	ΔEbind	0.4	rigid receptor	*
	ΔEbind	0.4	partially flexible receptor	*
Vina scoringd	ΔEbind	0.4	using MD-generated complexes	**
MM-PBSA	ΔEMM	0.8	ΔEvdW + ΔEelec (without cutoff)	**
	ΔEbind	1.0	ΔEMM + ΔGpolar + ΔGnonpolar	***
nonequilibrium MDe	ΔGwell	1.0	well depth of PMF profiles	*****
	ΔG0	1.0	binding free energy from PMF	*****

Four inhibitors with available IC50 (Figure d) are used for the Spearman rank correlation coefficient (SRCC) calculation.

Computational cost for different methods is rated on a scale of one (cheap) to five (expensive) asterisks (*).

Simulation of protein–ligand complexes.

The scoring function is identical whereas the scored complexes differ.

Simulation with an external force such as SMD and PMF calculations for the formation/dissociation process of protein–ligand complexes.

Consideration of the solvation effect results in a good prediction of the relative binding strength (SRCC = 1.0, Table ) in the MM-PBSA calculation. The MM-PBSA analysis for MSA-1 and OH-4 was replicated three times using independent MD simulations with different initial velocities. The total binding energies (ΔEbind) show consistency, validating that MSA-1 has a stronger binding than OH-4 (Table S6). Similar findings are observed for the identified key residues from the three replicas, while the energy contribution per residue might be different when using a different trajectory (Table S6). This indicates that a more reliable identification of key residues for the ligand binding may necessitate a consideration of different inhibitors with similar structures (as done in this work) or multiple trajectories.

The expensive method of PMF calculation has a good performance as well (Table ). PMF profiles reflect a combination of enthalpy and entropy contributions to the protein–ligand complexation from the binding partners and the surrounding solvent molecules, and generate absolute binding free energies that can be compared, in principle, with experimental observations. The PMF results indicate a total entropy loss disfavoring the binding with a contribution of ∼20 kJ/mol. Our calculated configurational entropy changes (−TΔS, Table S7) of the inhibitors upon complexation are unfavorable and range from approximately 30 (OH-4) to 150 (MSA-2) kJ/mol, on the same order of magnitude as the binding energy of ΔEbind (Table S2). We tested two methods (quasiharmonic approximation and Schlitter formula) for the entropy calculations; both methods give almost identical ΔS, although there is a large discrepancy in the absolute entropy values (Table S7). The entropy loss of ligand is likely accompanied by an entropy loss of the receptor, and most of the entropy loss would be canceled out via an entropy gain of solvent water molecules via, for instance, desolvation of the binding partners.[84−87] Entropy calculation, in particular for solvent molecules, is challenging and is a difficulty faced by virtual screening.[88] The total entropy change is expected to be small for the ligands with similar structures,[73] and one often chooses to use the relative binding energies instead, as done for ΔEbind in MM-PBSA.

Conclusions

A variety of computational methods was utilized to investigate the binding poses and binding energies of eight isoflavone analogues with human mitochondrial aldehyde dehydrogenase (ALDH2). We focused on two potent inhibitors of CVT-10216 (MSA-1) and daidzin (OH-3) and aimed to explore the substitution effects of both ends of the isoflavone skeleton on the binding with ALDH2. The method with a low computational load like molecular docking failed to reproduce the relative binding strength of isoflavone inhibitors against ALDH2. Equilibrium MD simulations with classical force fields in an explicit solvent were able to generate reasonable protein–ligand complexes. Together with MM-PBSA analysis, a good performance was obtained in reproducing the relative binding energies. The PMF calculation also performed well but with a high computational load. Considering the solvation contribution via either an implicit (PBSA) or explicit solvent was therefore required for investigating the inhibition of isoflavone analogues against ALDH2. Key residues were identified via the energy decomposition, which provides a structural basis for further design of inhibitors with enhanced binding.