A computational study to identify Sesamol derivatives as NRF2 activator for protection against drug-induced liver injury (DILI)

Drug-induced liver injury can be caused by any drugs, their metabolites, or natural products due to the inefficient functioning of drug-metabolizing enzymes, resulting in reactive oxygen species generation and leading to oxidative stress-induced cell death. For protection against oxidative stress, our cell has various defense mechanisms. One of the mechanisms is NRF2 pathway, when activated, protects the cell against oxidative stress. Natural antioxidants such as Sesamol have reported pharmacological activity (hepatoprotective & cardioprotective) and signaling pathways (NRF2 & CREM) altering potential. A Computational analysis was done using molecular docking, IFD, ADMET, MM-GBSA, and Molecular dynamic simulation of the Schrödinger suite. A total of 63,345 Sesamol derivatives were downloaded for the PubChem database. The protein structure of KEAP1-NRF2 (PDB: 4L7D) was downloaded from the RCSB protein database. The molecular docking technique was used to screen compounds that can form an interaction similar to the co-crystalized ligand (1VX). Based on MM-GBSA, docking score, and interactions, ten compounds were selected for ADMET profiling and IFD. After IFD, five compounds (66867225, 46148111, 12444939, 123892179, & 94817569) were selected for molecular dynamics simulation (MDS). Protein–ligand complex stability was assessed during MDS. The selected compounds (66867225, 46148111, 12444939, 123892179, & 94817569) complex with KEAP1 protein shows good stability and bond retentions. In our study, we observed that the selected compounds show good interaction, PCA, Rg, binding free energy, and ADMET profile. We can conclude that the selected compounds can act as NRF2 activators, which should be validated using proper in-vivo/in-vitro models. Graphical abstract


Introduction
Drug-induced liver injury (DILI) is defined "as a liver injury caused by various medications, herbs, or other xenobiotics, leading to abnormalities in liver tests or liver dysfunction with the reasonable exclusion of other etiologies" [1].Based on the damage site and injury period, it is categorized as hepatitis, cholestatic, or a mixed pattern of injury and chronic or acute [2].According to reports, the DILI incidence rate is 1.3-19.1 per 100,000 people worldwide [3][4][5].The exact incidence rate of DILI is complicated to estimate as most cases are underreported and/or due to the utilization of a wide range of diagnostic criteria [6,7].It can be triggered by a variety of drugs, including antituberculosis, chemotherapy drugs, and antipyretics analgesics.DILI-related problems are one of the major factors in the design of new drugs as well as restriction, withdrawal, or project termination of drug candidates and existing drugs.
The complexity underlying the pathophysiology of DILI has gathered interest from the scientific community in recent years.During the 1950s, scientists uncovered a connection between drug-reactive metabolites and DILI, and it has been an area of interest ever since.Several studies have associated mitochondrial dysfunction, inflammatory responses, and oxidative stress with the occurrence and development of DILI [3][4][5]8].DILI could be linked to direct toxicity from the prescribed drug or its metabolites, or it may be caused by immune-mediated processes.Even though the mechanisms are different, it is interlinked, e.g., drug-induced toxicity leads to hepatocyte damage, which may be exacerbated by successive inflammatory reactions [2].
Most drugs are metabolized inside the liver and then excreted in bile or urine.Drug metabolism involves two reactions, phase-I and phase-II, which have their own enzymes that work in unison for the metabolism and excretion of the drug from our body.The initial step in drug metabolism involves the phase-I enzymes, which react with the drug to form an intermediate product.The intermediate product formed during the phase-I reaction was inactivated in the phase-II reaction and excreted out.Any difference in the working of enzymes during phase-II leads to the accumulation of toxic metabolites.These toxic metabolites may interact with numerous cellular organelles (e.g., mitochondria), resulting in the generation of reactive oxygen species (ROS), which intensifies oxidative stress, which can directly damage enzymes, protein, and DNA inside cells and tissues and also prompt immune-mediated liver injury [9,10].
Nuclear factor erythroid 2-related factor 2 (NRF2) is a member of the cap'n'collar basic leucine zipper family located in the cytoplasm of the cell [11,12].The NRF2 has two motifs, i.e., Aspartic acid-Leucine-Glycine (DLG) and Glutamic acid-Threonine-Glycine-Glutamic acid (ETGE), which bound it with Kelch-like ECH-associated protein 1 (KEAP1) protein.DLG motif of NRF2 binds with SER (363, 508, 555 & 602), ARG (380, 415 & 483), ASN 382, GLN 530, and TYR 525, and ETGE motif binds with SER (363, 508, 555 & 602), ARG (380, 415 & 483), ASN 382, GLN 530, and TYR 525 amino acid residue of KEAP1 protein, respectively [13].Under normal conditions, NRF2 remains bounded with a KEAP1, and this bounded NRF2 is ubiquitination and degraded by the proteasome.Under stress conditions, the bonding between the KEAP1 and the ETGE and DLG motif of NRF2 disrupts.The free NRF2 translocates inside the nucleus, where it gets bound with the antioxidant response element.Initiating the transcription of the genes, including glutathione-S-transferases and heme-oxygenase-1, for the production of phase-II enzymes and antioxidant proteins, decreases the sensitivity toward oxidative stress [14,15].Recently, research has been carried out to identify substances that can destabilize the protein-protein interaction of KEAP1-NRF2, freeing NRF2 for protection against oxidative stress [13,16].The activation of NRF2 has the potential to provide protection against liver damage caused by DILI resulting from oxidative stress, which makes it a potential target for further research.

Protein preparation and grid generation
KEAP1-NRF2 crystal structure (PDB id: 4L7D) with human origin was downloaded from the RCSB protein data bank (https:// www.rcsb.org/ struc ture/ 4L7D).The selected PDB has a Resolution of 2.25 Å, RMSD value of 0.8682 Å, and IC 50 of the co-crystalized Ligand (1VX) was 0.75 ± 0.32 μmol.1VX is an Nrf2 activator and a direct inhibitor of the Nrf2-Keap1 complex formation.It was identified through a high-throughput screen of small molecules and has shown promising binding interactions with the Keap1 protein [16].
Before beginning the molecular simulation studies, the Schrödinger suite module "Protein Preparation Wizard" was used to prepare protein.During the protein preparation process, the protein was imported, refined, reviewed, modified, and minimized.Side chains and missing residue were filled using the prime tool.Catalytically significant residue and active sites were preserved in the protein structure.Beyond 5 Å, water molecules were removed, and stage for heteroatoms was generated.OPLS3e (Optimized potential for liquid stimulation) force field was used for energy minimization and generating low-energy state protein.Molecular modeling was done using this protein.The panel "Receptor grid generation" was used for generating a grid encompassing the ligand at coordinate X (− 22.79), Y (39.11), & Z (-36.57) while keeping all the functional residue inside the grid in default [35][36][37].

Ligand docking
Ligand docking studies were performed using the GLIDE (Grid-based Ligand Docking with Energetics) module.High-throughput Virtual Screening (HTVS), Standard precision (SP), and Extra precision (XP) were used to get different scoring functions.Firstly, HTVS docking was used to dock all the ligands.As HTVS docking computationally lacks explicit water technology and descriptor, the top 5 percent of the ligands were again analyzed using SP and XP modes to avoid false-positive results [35][36][37].

Free ligand binding energy (MM-GBSA) calculation
MM-GBSA method was used to find the absolute Ligand binding affinities of the top 100 ligands.It facilitates in comprehending the binding of protein-ligand which will persist long enough to elicit biological response [38].

Induced fit docking (IFD)
IFD of the top 10 selected compounds were performed where 20 was set as the maximum pose for each ligand, and 0.50 was kept as the van der Waals scaling.After that, prime side chain prediction and energy minimization were conducted, and all residue refinement was done for side chains and ligand poses.The induced fit protein structure underwent extensive XP docking, and the IFD dock score was computed [40][41][42].

Molecular dynamic simulations (MDS)
For investigating the dynamics and function of protein-ligand complexes, MDS is used widely.Ligand docking does not replicate the biological environment where the ligand and protein are dissolved in water.Therefore, MDS was applied to resolve the issue.Depending on the ligand docking score, poses, MM-GBSA, ADMET profile, and IFD data, compounds will be selected for MDS for 100 ns.Before beginning the MDS, the complete system was submerged inside a simple point charge (SPC) solvent model.Throughout the system-building process, the boundary condition's orthorhombic box form was maintained.The Buffer box size calculation method is utilized throughout the building process, and the distance (Å) was a = 10, b = 10 & c = 10, and the angle was α = 90֯ , β = 90֯ , γ = 90֯ , respectively.The System Builder tool was used to prepare the system for MD simulation, and the OPLS3e force field was employed throughout the system preparation process.Minimization tool was used for system minimization.Finally, a simulation time of 100 ns was selected, where every 100 picoseconds 1 frame was captured from the trajectory, and 1000 frames were generated throughout the MDS process.NPT (constant particle number (N), pressure(P) 1.01325 bar, and temperature (T) 300 K) ensemble were used for the production run with OPLS3e as the chosen force field for the MDS process.Lastly, the Simulation interaction diagram tool was used to generate the report of the MDS [40][41][42].

Result and discussion
PDB (4L7D) for KEAP1-NRF2 from human origin was downloaded from the "rscb.org"portal.The selected PDB has an RMSD value of 0.8682.63,345 Ligands were downloaded from the PubChem database.Protein and ligand were prepared at 7.4 ± 0.0 pH.All the procedures are followed as mentioned in the material and method.

Molecular docking and MM-GBSA
Prepared ligands were docked using HTVS docking mode, the ligands having a score more than − 6 kcal/mole were selected and again analyzed in SP docking mode.Again, the top five percent of ligands were subjected to the XP docking mode.XP docking mode provides better accuracy, decreases false-positive results, and provides an improved association between ligand pose and docking score.Finally, 100 compounds with pivotal interaction and a good XP dock score were selected.Ligand structure, MM-GBSA score, and docking score of the top 10 compounds and 1VX (cocrystalized ligand) are shown in Table 1, and 2d interactions are shown in Table 2.
The selected ligands, protein-ligand complex binding free energy, were calculated using the Prime MM-GBSA module.The Binding free energy of the 1VX (co-crystalized ligand) was − 86.39 kcal/mol, and the top 10 compounds ranges from − 102.22 to − 63.35 kcal/mol, shown in Table 1.

ADMET analysis
After the ligand XP docking and MM-GBSA, ADMET properties of the selected ligands and 1VX were calculated using the QikProp tool.In the drug-designing process, drug-like properties and ADME properties are key filtration criteria, and they are presented in Tables 3 and 4. All the top 10 compounds show good percent oral absorption, with compounds 46148111 and 94817569 showing 100% absorption.Other properties like H-bond donor (range: 0.0-6.0),acceptor (range: 2.0-20.0),and QPlogPo/w (range: − 2.0 to 6.5) were all within the range.Also, all the compounds follow the Lipinski's rule of five.
Pharmacokinetic parameters, mainly QPlogKhsa, QPlogBB, QPlogHERG, and QPPCaco were accessed and are shown in Table 4.The QPlogKhsa values of the selected compounds were in the range of − 0.539 to 0.599, within the range or recommended value of − 1.5 to 1.5.The QPlogBB normal range is − 3.0 to 1.2, and all the compounds were in the range of − 2.222 to − 0.193.The QPlogHERG of all the compounds was in the range of − 6.374 to − 4.408.The QPPCaco-2 permeability of all the compounds ranged from 72.3 to 2154.497, with compounds 46148111 and 94817569 showing great permeability, i.e., 757.906 and 2154.497,respectively.Toxicity predicted by "admetSAR" is shown in Table 5.

Induced fit docking (IFD)
Ten compounds were chosen for IFD analysis depending on their interactions, XP docking score, and ADMET.After IFD, five compounds, namely, 66867225, 46148111, 12444939, 123892179, and 94817569, were selected for Molecular dynamic analysis by virtue of their IFD score and bond retention.The interaction image and IFD score of the co-crystallized ligand and selected five compounds are shown in Fig. 1.

Protein-ligand complex root mean square deviation (RMSD) and Lig fit prot
The Lig fit prot and RMSD of protein and protein-ligand complex were computed by aligning the generated frame with the reference frame.Complex-1 was stable in the beginning with a slight drift at 4 ns and became stable till 70-75 ns, where a slight drift is observed, after which it stabilizes toward the end (Fig. 2a).In Complex-2, there is a slight drift observed at 6.5-14 ns, after 14 ns, the complex becomes stable till 95 ns, with a slight drift observed at 95.3 ns, and it becomes stable till the end (Fig. 2b).
In Complex-3, a drift is observed during the initial stage of simulation at 4.8 ns, 10.4 ns, 18.4 ns, & 26.4 ns after that, it becomes stable, and again a slight drift is observed at 73.6 ns, and it became stable toward the end of the simulation (Fig. 2c).In Complex-4, some minor deviation was observed initially, and it became stable throughout the simulation (Fig. 2d).In Complex-5, a deviation was observed at 15 ns & 42 ns, which stabilizes after 43 ns till 100 ns (Fig. 2e).The Complex-6 was stable throughout the entire simulation (Fig. 2f).
respectively.Compound 46148111 of Complex-2 show better contact duration than the rest of the compounds.Contact duration of compound 46148111 was similar to that of 1VX (complex-6).The protein-ligand contact histogram and timeline of selected compounds (complex-1, complex-2, complex-3, complex-4, & complex-5) and co-crystallized ligand (complex-6) are shown in Figs. 5 and 6.

Principal component analysis (PCA)
Principal Component Analysis (PCA) is a highly effective technique for reducing complexity and extracting significant motions in the analysis of Molecular Dynamics (MD) simulations.In this approach, a matrix was constructed based on the trajectories, excluding rotational and translational movements.The essential dynamics protocol was then applied to compute the eigenvectors, eigenvalues, and their projections onto the first two Principal Components (PCs).By diagonalizing the matrix, the eigenvectors and eigenvalues were determined, with the eigenvalues representing the magnitude of the corresponding eigenvectors.The matrix of eigenvectors describes the multidimensional space, providing information about the displacement of atoms in the protein along each direction.Desmond trajectories obtained after MDS was converted using ProDy 2.0 for Visual comparative analysis using VMD's Normal Mode Wizard plugin [43][44][45][46] -6) was in the range of 3.08 to − 0.66 nm & 5.39 to − 2.79 nm, respectively.The range of PC1 and PC2 for all the complexes were found to be less than that compared to that of Co-crystalized ligand complex (Complex-6), which points toward the better stability of selected ligands as compared to that of Co-crystalized Ligand.The scatter plot of PC1 vs PC2 is shown in Fig. 8.

Binding free energy (MM-GBSA dG bind)
Post MDS, binding free energy of the selected complexes was calculated at every 10th frame for the duration of entire simulation, i.e., till the 1000th frame shown in

Total energy
The total energy (Coulomb, van der Waals, bond, angle, and dihedral) of the protein-ligand complex along with metals and ions, water, and hets groups during the molecular dynamic simulation is shown in Fig. 10.The total energy of Complex-1 lies in the range of 333923 to 335016 kcal/mol and the average value was 334518 kcal/ mol.For Complex-2, the total energy lies in range of 335404 to 336255 kcal/mol, with average value of 335879 kcal/mol.The Complex-3 total energy lies in range of 333539 to 334560 kcal/mol, with average value of 334049 kcal/mol.For Complex-4, the total energy lies in range 333598 to 334483 kcal/mol, with average value of 334053 kcal/mol.For Complex-5, the total energy lies in range of 333528 to 334544 kcal/mol, with average value of 334060 kcal/mol.For Complex-6, the total energy lies in range of 350796 to 351776 kcal/mol, with average value of 335879 kcal/mol.

Solvent accessible surface area (SASA)
Solvent accessible surface area (SASA) is defined as protein surface area, which are accessible to solvent.It helps in understanding the behavior of target protein when bounded with a ligand.SASA of the protein-ligand complex during the entire simulation was recorded and shown in Fig.  interaction between the ligand-receptor complexes, similar to that of co-crystalized ligand-protein complex (Complex-6).

Radius of gyration (Rg)
The radius of gyration helps in understanding the compactness of the protein-ligand complex.In all the complexes except Complex-5, minor fluctuations were observed and they were stable during the entire simulation period.In Complex-5, from 14 to 42 ns slight fluctuation was observed but it becomes stable after 43 ns till the end of the simulation (Fig. 12).residue of the KEAP1 protein's active pocket, and they were stable.All the complexes also show good RMSF, PCA, Rg, SASA, RMSD, and total energy.All the selected compounds show good binding free energy, pharmacokinetics, and toxicity profiles and follow Lipinski's rule of Five.
Based on the computational assay, the selected compounds 66867225, 46148111, 12,444,939, 123892179, and 94817569 can serve as an activator of NRF2.However, our result should be validated using a proper in-vivo/in-vitro experimental model of DILI.

Table 1
Ligand structure, Docking score, and MM-GBSA of the co-crystal ligand and top 10 selected compound

Table 4
ADMET properties of selected compoundsQPlogKhsa prediction of binding to human serum albumin, QPlogBB predicted brain/blood partition coefficient, QPlog hERG predicted IC 50 value for blockage of HERG K + channels, QPlogS predicted aqueous solubility, QPPCaco predicted apparent Caco-2 cell permeability in nm/sec, QPlog Po/w predicted octanol/water partition coefficient