Evaluation of the IKKβ Binding of Indicaxanthin by Induced-Fit Docking, Binding Pose Metadynamics, and Molecular Dynamics

Background: Indicaxanthin, a betaxanthin belonging to the betalain class of compounds, has been recently demonstrated to exert significant antiproliferative effects inducing apoptosis of human melanoma cells through the inhibition of NF-κB as the predominant pathway. Specifically, Indicaxanthin inhibited IκBα degradation in A375 cells. In resting cells, NF-κB is arrested in the cytoplasm by binding to its inhibitor protein IκBα. Upon stimulation, IκBα is phosphorylated by the IKK complex, and degraded by the proteasome, liberating free NF-κB into the nucleus to initiate target gene transcription. Inhibition of the IKK complex leads to the arrest of the NF-κB pathway.

Methods: To acquire details at the molecular level of Indicaxanthin’s inhibitory activity against hIKKβ, molecular modeling and simulation techniques including induced-fit docking (IFD), binding pose metadynamics (BPMD), molecular dynamics simulations, and MM-GBSA (molecular mechanics-generalized Born surface area continuum solvation) have been performed.

Results: The computational calculations performed on the active and inactive form, and the allosteric binding site of hIKKβ, revealed that Indicaxanthin inhibits prevalently the active form of the hIKKβ. MM-GBSA computations provide further evidence of Indicaxanthin’s stability inside the active binding pocket with a binding free energy of −22.2 ± 4.3 kcal/mol with respect to the inactive binding pocket with a binding free energy of −20.7 ± 4.7 kcal/mol. BPMD and MD simulation revealed that Indicaxanthin is likely not an allosteric inhibitor of hIKKβ.

Conclusion: As a whole, these in silico pieces of evidence show that Indicaxanthin can inhibit the active form of the hIKKβ adding novel mechanistic insights on its recently discovered ability to impair NF-κB signaling in melanoma A375 cells. Moreover, our results suggest the phytochemical as a new lead compound for novel, more potent IKKβ inhibitors to be employed in the treatment of cancer and inflammation-related conditions.

Introduction

Indicaxanthin, ((2S)-2,3-dihydro-4-[2-[(2S)-2-carboxypyrrolidin-1-yl]ethenyl]pyridine-2,6 dicarboxylic acid), is a betalain pigment belonging to the betalain class compounds (Figure 1). This includes vacuole pigments restricted to fruits and flowers of 10 families of Cariophyllalae plants and a few superior fungi of the genus Amanita of the Basidiomycetes (Slimen et al., 2017). Indicaxanthin has been investigated over the last 20 years by the author’s research group for its chemical and nutraceutical properties. Interestingly the phytochemical has been demonstrated to be highly bioavailable in humans and even permeable to the rat brain-blood-barrier (Tesoriere et al., 2004; Allegra et al., 2015). Moreover, because of its reducing and amphipathic properties, it has been shown to interfere with cellular, redox-dependent signal transduction pathways in several experimental models of inflammatory-related, oxidative stress-dependent pathological conditions. Along these lines, Indicaxanthin has been demonstrated to exert significant reducing, anti-oxidative, anti-inflammatory, spasmolytic, and neuromodulatory effects both in vitro and in vivo (Allegra et al., 2019).

FIGURE 1

FIGURE 1. Indicaxanthin structure.

In addition to its redox-modulating and anti-inflammatory properties, Indicaxanthin has also been shown to exert antiproliferative effects against melanoma cells being also able to significantly impair tumor progression in a mouse model of cutaneous melanoma (Allegra et al., 2018). Mechanistic evaluations have individuated in the inhibition of the nuclear factor-κB (NF-κB) signaling, a key event underlying the antitumoral effects exerted by the pigment.

The NF-κB transcription factor family consists of five different DNA-binding proteins that form a variety of homodimers and heterodimers regulating both the innate and adaptive immune responses (Taniguchi and Karin, 2018). More interestingly, besides its modulating effects on the entire inflammatory response, NF-κB is responsible for accelerating cancer progression, metastasis, angiogenesis, and drug resistance (Zhang et al., 2017). Indeed, in several types of cancer, both in malignant cells and in the tumor microenvironment, NF-κB is constitutively activated, and only rarely is such activation due to NF-κB-related genetic alterations (Staudt, 2010; Ben-Neriah and Karin, 2011).

Melanoma is an aggressive skin cancer, notoriously resistant to current cancer therapies (Van Herck et al., 2021). To effectively prevent chemoresistance development and reduce the risk of unwanted side effects, current pharmacological strategies employ a multi-targeted approach. Given that constitutive or drug-induced upregulation of NF-κB activity is associated with chemoresistance, NF-κB is considered one of the most relevant targets to discover new active compounds (Taniguchi and Karin, 2018). Permanent activation of NF-κB signaling in melanoma has been reported to proceed through the activation of the so-called canonical pathway (Yang and Richmond, 2001; Dhawan and Richmond, 2002; Amiri and Richmond, 2005; Amschler et al., 2010). Through this process, activation of NF-κB depends on the degradation of its specific inhibitors (IκB) consisting of IκBα, IκBβ, and IκBε. Typically, IκBs bind to NF-κB complexes, inhibiting their DNA binding and keeping them in a predominantly, inactive cytoplasmic form. Tumour microenvironmental stimuli can lead to the activation of a large cytoplasmic protein complex: the IκB kinase (IKK). The precise nature of this molecular mechanism remains to be elucidated but it contains IKKα, IKKβ, and IKKγ as the three seminal components. The phosphorylated and thus activated IKK complex is responsible for the phosphorylation of IκB, marking it for degradation through the proteasomal degradation machinery. The free NF-κB dimers (p50-p65) can, then, translocate from the cytoplasm to the nucleus, bind to DNA, and regulate gene transcription. Typical targets within the classical NF-κB signaling include genes encoding pro-inflammatory cytokines, growth factors, chemokines, matrix metalloproteinases, pro-proliferative proteins, anti-apoptotic proteins, pro-inflammatory enzymes, angiogenic factors, and adhesion molecules (Staudt, 2010). As a primary druggable mediator of canonical NF-κB signaling, the IKKβ enzyme inhibition has been the historical focus of drug development pipelines. Thousands of compounds with activity against IKKβ have been characterized, with much demonstrating promising efficacy in pre-clinical models of cancer. However, severe on-target toxicities and other safety concerns associated with systemic IKKβ inhibition have so far prevented the clinal approval of any IKKβ inhibitors (Prescott and Cook, 2018).

Indicaxanthin has been demonstrated to inhibit IκBα degradation in melanoma A375 cells at 100 μM, a concentration at which it impairs NF-κB signaling and inhibits 50% cell proliferation (IC₅₀) (Allegra et al., 2018). In the light of this evidence, we employed an in silico approach to evaluate Indicaxanthin’s inhibitory activity against IKK. To this end, molecular modeling and simulation techniques including induced-fit docking (IFD), binding pose metadynamics (BPMD), molecular dynamics simulations, and MM-GBSA (molecular mechanics-generalized Born surface area continuum solvation) free energy calculation have been performed.

Materials and Methods

Protein and Ligand Preparation

For the purpose study, the 2.83 Å resolution crystal structure of human IKKβ (hIKKβ), which is partially phosphorylated and bound to the staurosporine analog K252a (PDB ID: 4KIK) (Liu et al., 2013) was used. The structure was optimized using the Protein Preparation Wizard in Maestro (Schrödinger, 2017) adding bond orders and hydrogen atoms to the crystal structure using the OPLS3 force field. Prime was used to fix missing residues or atoms in the protein and to remove co-crystallized water molecules. PROPKA was used to check for the protonation state of ionizable protein groups (pH = 7.0). The hydrogen bonds were optimized through the reorientation of hydroxyl bonds, thiol groups, and amide groups. In the end, the system was minimized with the value of convergence of the RMSD of 0.3 Å. Indicaxanthin and staurosporine analog K252a were prepared using LigPrep The force field adopted was OPLS3 (Jorgensen et al., 1996) and Epik 3.9 (Schrödinger, 2017-1) was selected as an ionization tool at pH 7.2 ± 0.2. Tautomers generation was unflagged and the maximum number of conformers generated was set at 32.

Induced-Fit Docking

The induced-fit docking (IFD) is a method for modeling the conformational changes induced by ligand binding developed by Schrödinger (Sherman et al., 2006). This protocol models induced-fit docking of one or more ligands using the steps as also reported in (Tutone et al., 2019). Initial docking of each ligand is performed using a softened potential (van der Waals radii scaling). Then, a side-chain prediction within a given distance of any ligand pose is performed. Subsequently, a minimization of the same set of residues and the ligand for each protein/ligand complex pose is performed. After this stage, any receptor structure in each pose reflects an induced fit to the ligand structure and conformation. Finally, the ligand is rigorously docked, using Glide XP, into the induced-fit receptor structure. The grid boxes for the binding sites of Chain A (inactive form) and B (active form) were built considering the co-crystallized ligand staurosporine analog K252a as a centroid. For the allosteric sites of Chain A, the amino acid residues previously identified by Liu et al. (2018) were considered crucial for centering the docking grid. During the initial docking procedure, the van der Waals scaling factor was set at 0.5 for both receptor and ligand. The Prime refinement step was set on side chains of residues within 5 Å of the ligand. For each ligand docked, a maximum of 20 poses was retained to be then redocked at XP mode.

Binding Pose Metadynamics

Binding pose Metadynamics (BPMD) is an automated, enhanced sampling, metadynamics-based protocol. During the simulations, the ligands are forced to stir in their binding pose. This method showed the ability to reliably discriminate between the correct ligand binding pose and plausible alternative generate with docking or induced-fit docking studies (Clark et al., 2016).

Ten independent metadynamics simulations of 10 ns each were performed using as a collective variable (CV) the measure of the root-mean-square deviation (RMSD) of the ligand heavy atoms concerning their starting position. The alignment before the RMSD calculation was done by selecting protein residues within 3 Å of the ligand. The alpha carbons of these binding site residues were then aligned to those of the first frame of the metadynamics trajectory before calculating the heavy atom RMSD to the ligand conformation in the first frame. The hill height and width were set to 0.05 kcal/mol (about 1/10 of the characteristic thermal energy of the system, kBT) and 0.02 Å, respectively (Clark et al., 2016). Before the metadynamics run, the system was solvated in a box of SPC water molecules followed by several minimizations and restrained MD steps that allow the system to reach slowly the desired temperature of 300 K as well as releasing any bad contacts and/or strain in the initial starting structure. The final snapshot of the short unbiased MD simulation of 0.5 ns is then used as the reference for the following metadynamics production phase. The pivotal idea of BPMD is that ligands that are not stably bound to the target binding pocket will move forward higher RMSD as compared to the stably bound ones if they are exposed to the same biasing force. After the simulation, the stability of the ligand during the course is represented by three scores: PoseScore, PersistenceScore (PersScore), and CompositeScore (CompScore). PoseScore is indicative of the average RMSD from the starting pose. The steepest increase of this value is a symptom that the ligand is not in a well-defined energy minimum and probably it might not have been accurately modeled. PersScore is a measure of the hydrogen bond persistence calculated in the last 2ns of the simulation that have the same number of hydrogen bonds as the input structure, averaged over all the 10 repeated simulations. PersScore covers a range between 0 and 1, where 0 indicates that either the starting ligand pose did not have any interaction with the target or that the interactions have been lost during the simulations, while 1 indicates that the interactions between the staring ligand pose and the last 2 ns of the simulations have been retained. CompositeScore is the linear combination of PoseScore and PersScore, lower values equate to more stable complexes. Each complex, previously obtained, was run on Desmond on a single node with 1 GPU card, taking for a typical system (1 complex = 1 × 10 metadynamics run) 72 h.

MD Simulations

The plain MD simulations were carried out using Desmond 4.9 (Bowers et al., 2007) using the OPLS3 force field (Jorgensen et al., 1996). The complexes were solvated in orthorhombic boxes using the TIP3P water model. Ions were added to neutralize charges. The systems were minimized and equilibrated at a temperature of 303.15 K and a pressure of 1.013 bar. The system was simulated as an NPT ensemble; a Nose–Hover thermostat and Martyna–Tobia–Klein barostat were used. The integration time step was chosen to be 2 fs. To keep the hydrogen–heavy atom bonds rigid, the SHAKE algorithm was used. A 9 Å cutoff radius was set for the short-range Coulomb interactions, and smooth particle mesh Ewald was used for the long-range interactions. For each system, we carried out 100 ns MD, with 1.2 ps detection ranges for energy, and 4.8 ps for the trajectory frames. Visualization and analysis of the MD trajectories were performed using Desmond simulation analysis tools in Maestro.

MM-GBSA Binding Free Energy Calculation

The MM-GBSA approach employs molecular mechanics, the generalized Born model, and the solvent accessibility method to determine free energies from structural information circumventing the computational complexity of free energy simulations wherein the net free energy is treated as a sum of a comprehensive set of individual energy components, each with a physical basis (Genheden and Ryde, 2015). We applied this method to the snapshots extracted from the 100 ns production MD trajectories. Protein-ligand binding free energy using MM-GBSA was calculated as the difference between the energy of the bound complex and the energy of the unbound protein and ligand. In this work MM-GBSA calculations were also achieved in Prime software (Jacobson et al., 2004); the entropy term ‒TΔS was not calculated to reduce computational time. The VSGB solvation model was chosen using OPLS3 Force Field with minimized sampling method.

Results and Discussion

The hIKKβ protomer adopts a trimodular structure that closely resembles that of Xenopus laevis (xIKKβ): an N-terminal kinase domain (KD), a central ubiquitin-like domain (ULD), and a C-terminal scaffold/dimerization domain (SDD). In particular, the selected crystal structure has one protomer in the active conformation with phosphorylated Ser177 and Ser181 (Chain B), and the other protomer is in the inactive conformation with the same residues of serine unphosphorylated in the activation loop (Chain A). Recently, Liu and co. identified a druggable allosteric site between KD and ULD.

To evaluate the binding capability of Indicaxanthin into the hIKKβ, we performed a series of computational studies with increasing accuracy, induced-fit docking (IFD), binding pose metadynamics (BPMD), unbiased molecular dynamics (MD) followed by MM-GBSA (molecular mechanics-generalized Born surface area continuum solvation) free energy calculation. We began the computational studies by selecting the crystal structure of the human IKKβ (PDB ID: 4KIK) (Liu et al., 2013) bound to the staurosporine analog K252a in the KD (Figures 2A,B). Firstly, we optimized the crystal structure by completing and refining the missing loops and residues and optimizing amide groups of asparagine (Asn) and glutamine (Gln), and the imidazole ring in histidine (His); and predicting protonation states of histidine, His, aspartic acid (Asp) and glutamic acid (Glu) and tautomeric states of histidine. Then, the docking studies have been performed centering the docking boxes on the 3D coordinates of K252a, both in the inactive Chain A and in the active Chain B. The RMSD of K252a in Chain B and A was calculated showing values 0.8–0.18 Å respectively (The superimposition of re-docked structures are reported in Supplementary Material). Moreover, another docking box for the allosteric binding pocket was centered on the previously identified residues by Liu and Co. (Liu et al., 2018).

FIGURE 2

FIGURE 2. 3D structure of hIKK forms: (A) the inactive form, Chain A, in purple the KD binding pocket, in red the allosteric binding site; (B) the active form, Chain B, in cyan the KD binding pocket of the inactive Chain A with a volume of 551 Å. This pocket is surrounded by three α helixes (Arg118-Ser127, Leu265-Leu273, and Leu303-His313) of KD and a loop (Thr368-Leu386) between two β sheets of ULD. They suggested that small molecules binding in the pocket between KD and ULD probably interfere with the kinase function by disrupting the interaction between these two domains (Liu et al., 2018).

Induced Fit Docking

Induced-Fit docking (IFD) experiments were first carried out. IFD confers flexibility to the protein side chains, allowing the ligand to adjust and optimize binding interactions within the active site. We performed the IFD of K252a and Indicaxanthin in the two sites previously identified, and the IFD of Indicaxanthin in the allosteric binding site.

In most kinases, the ATP binding site is a narrow hydrophobic pocket located between the N-lobe and C-lobe of the kinase domain (KD) which are linked by a flexible hinge region. This binding site is partly covered by an activation loop comprised of Serine, Threonine, and Tyrosine residues in the unphosphorylated state. The N-terminal side of the activation loop consisting of a highly conserved triplet DLG (Asp166, Leu167, and Gly168) which is involved in the catalytic transfer of the γ-phosphate group (Xu et al., 2011). In IKKβ the residues Glu97, Tyr98, and Cys99 are part of the hinge region. The backbone groups of Glu97 and Cys99 can provide hydrogen bonding interactions with an inhibitor. (Leung et al., 2013).

Docking of K252a in Chain A showed a single H-bond interaction with Cys99 (hinge region), with the carbonyl of lactam ring. Other H-bonds were found between the amine of the K252a lactam group and Glu97 and the hydroxyl group and Glu149. Moreover, an aromatic H-bond was found with Asp166 (DLG triad). In chain B, K252a, as in Chain A, established the same interactions with Glu97, Cys99 (hinge region), and Glu149, with the only difference in another H-bond interaction between the carbonyl group of ether group and Thr23, and an aromatic H-bond with Tyr98. (Figures 3A,B).

FIGURE 3

FIGURE 3. IFD pose of K252a in the inactive Chain A (A); IFD pose of K252a in the active Chain B (B); IFD pose of Indicaxanthin in the inactive Chain A (C); IFD pose of Indicaxanthin in the active Chain B (D). H-bonds interactions are represented in yellow dashes, aromatic H-bonds are represented in light blue dashes, the salt bridges in purple dashes.

We performed the IFD of Indicaxanthin in the three sites identified previously. In the best IFD result, in particular, considering the Chain A (Docking score = −6.166, Table 1), the carboxyl group of the pyrrolidinium ring interacts with Thr23 and Gly24 forming two H-bonds, while the Cys99 simultaneously interacts with the carboxyl groups and the nitrogen of pyridine moiety. The carboxyl group of pyridine in position 11 interacts with Tyr98 (hinge region) with another H-bond and with Lys106 (solvent accessible region) through one salt bridge. The carboxyl groups in the pyridine moiety establish two aromatic H-bonds with Tyr98. Finally, another salt bridge involves the pyrrolidinium nitrogen and the Asp103. In Chain B (Docking score = −7.293), Indicaxanthin showed similar interactions as in Chain A. Two H-bond interactions involve the carboxyl group of the pyrrolidinium ring and Thr23, and Asn28; one of the carboxyl groups of pyridine interacts through H-bond with Tyr98, Lys106. The last interactions were a salt bridge between two carboxyl groups of the pyridine and Lys106 (Figures 3C,D). The IFD study of Indicaxanthin in the allosteric site (Docking score = −6.117) showed a single H-bond interaction was found between one of the carboxyl group of pyridine and His380; another interaction was a salt bridge, between the nitrogen of pyrrolidinium moiety and Asp373 (Figure 4).

TABLE 1

TABLE 1. BPMD and docking scores for the Indicaxanthin and K252a complexes.

FIGURE 4

FIGURE 4. IFD pose of Indicaxanthin in the inactive Chain A allosteric binding site. H-bonds interactions are represented in yellow dashes, the salt bridges in purple dashes.

Moreover, we performed docking simulations with different hIKKβ inhibitors to be considered as reference compounds for the analyses of the binding interactions with key residues from the binding sites and to compare with the predicted binding affinity for Indicaxanthin. We selected three known orthosteric inhibitors: MLN120B (IC₅₀ = 60 nM) (Prescott and Cook, 2018); the imidazo [1,2-a]quinaxoline derivative 6a (IC₅₀ = 324 nM) (Moarbess et al., 2016); LASSBio-1524 (IC₅₀ = 20 µM) (Guedes et al., 2016); and one allosteric inhibitor: BMS345521 (IC₅₀ = 300 nM) (Prescott and Cook, 2018). The docking scores of these known inhibitors are consistent with the experimental IC₅₀, considering that Indicaxanthin showed an IC₅₀ = 100 µM. Additionally, the key residues identified for these known inhibitors are similar to those identified for Indicaxanthin (Cys99, Asp103, Lys106). The allosteric inhibitor BMS345521 showed key interactions with His380, already identified for indicaxanthin, Asn308, Ser127. The 2D and 3D figures of the docked reference inhibitors together with key interactions and docking scores are reported in Supplementary Material.

Binding Pose Metadynamics

Binding pose metadynamics (BPMD) is an automated, enhanced sampling, metadynamics-based protocol, in which the ligand is forced to move around its binding pose. This method showed the ability to reliably discriminate between the correct ligand binding pose and plausible alternative generate with docking studies. In particular, Clark and colleagues introduced this metadynamics plus IFD strategy for accurate and reliable prediction of the structures of protein−ligand complexes at a useful computational cost. Their strategy allows treating the problem in full atomistic details, significantly enhancing the predictive power of IFD methods (Clark et al., 2016).

Firstly, we decided to use BPMD to evaluate the stability of the K252a best poses (in terms of docking scores) obtained from IFD studies into the binding site in Chain A and Chain B to evaluate the reliability of the IFD poses obtained as validation of the docking scoring functions. The results were defined in terms of pose stability based on the PoseScore that is the RMSD of the ligand related to the starting coordinates of the heavy atoms of the ligand. A PoseScore <2 Å is considered stable for the co-crystallized ligand. Moreover, another metric is used to evaluate the results such as the PersScore that is a clue of H-bonding formed between the ligand and the target during the simulations. In general, maintaining 60% of the total H-bonds (PersScore>0.6) is considered a good sign of stable interaction. The linear combination of these two scores provides a third score, the CompScore which is calculated as follows:

CompScore = PoseScore—5 × PersScore.

Lower values of this equate indicate more stable complexes.

The simulation performed on K252a pose in Chain A showed a PoseScore of 1.206. The PersScore showed that the hydrogen bonds were kept for 81% of the simulation time. The interactions by the lactam ring were confirmed to stabilize the molecule, in particular, we observed for 98.2% of the simulation time, interactions between the amine group and Glu97, for 97.3% of the simulation time interactions between the carbonyl oxygen and Cys99. The H-bond interaction between the hydroxyl group and Glu149 was kept for 48.2% of the simulation time. The CompScore of −2.852 confirmed that the starting molecule pose is stable into the active site.

The results of the K252a pose in Chain B showed a PoseScore of 0.909, the PersScore proves that for 63.9% of the simulation time, the H-bonds were maintained. As in Chain A, the interactions by the lactam ring were confirmed, in particular, 100% of the simulation time, the interaction between the carbonyl oxygen and Cys99, and 98.2% of the simulation time interaction between the amine group and Glu97. The same for H-bond interaction between the hydroxyl group and Glu149 was kept for 57.3% of the time. The linear combination of PoseScore and PersScore, CompScore, was −2.284. The result obtained confirmed the accuracy of the identified IFD poses. (Table 1).

Successively the BPDM simulations have been performed at the binding site in Chain A, the best pose of Indicaxanthin reached a steady PoseScore of 3.133, considered stable, while PersScore showed that the hydrogen bonds identified at the start of the metadynamics run were kept for 12% of the averaged time. In particular, the H-bond between carboxyl oxygen and Cys99. The CompScore value was 2.533.

At the binding site in chain B, the averaged RMSD of Indicaxanthin reached a steady PoseScore of 4.631, PersScore of 0.04. In particular, three H-bond were recorded during the ten replicas. H-bond interaction between the carboxyl group of pyrrolidinium moiety and Thr23 was kept for 11.8% of the simulation time, the same group showed an ulterior H-bond with Asn28 for 7.3% of the simulation time. The last interaction was between one of the carboxyl groups of pyrimidine and Lys106 for 2.7% of the simulation time. This interaction was supported by salt bridges between the two inferior carboxyl groups and Lys106. The last score, CompScore was 4.441. The BPMD analysis confirmed that the poses obtained from the IFD are valuable for Indicaxanthin and they can be used to perform a time-dependent study which could give more insights about the binding of Indicaxanthin. The BPDM results of Indicaxanthin into the allosteric binding site give unstable values suggesting that Indicaxanthin is not a stable binder of the allosteric site. However, to gain more insights, we decided to perform also the time-dependent simulation of this latter complex. (Table 1).

Molecular Dynamics Simulations

The running of dynamics simulations of protein-ligand complexes over time could be considered the major accuracy step in computer-assisted drug design. In this study, unbiased molecular dynamics simulations have been performed to explain the stability of Indicaxanthin as an inhibitor against the two forms of hIKKβ, active and inactive, and allosteric inhibitor of the inactive form. Additionally, we tried to understand if the protein target undergoes conformational alteration after interacting with Indicaxanthin. Therefore, starting from the previous IFD poses which BPMD analysis showed to be accurate and reliable except for Indicaxanthin in the allosteric site, five systems have been generated and submitted each for 100 ns in MD simulations (Chain A-Indicaxanthin, Chain B-Indicaxanthin, Chain A allosteric site-Indicaxanthin, Chain A-K252a, Chain B-K252a). The various analysis such as the root mean square deviation (RMSD), the root mean square fluctuation (RMSF), number, and types of protein-ligand contacts have been carried to have a more detailed analysis of Indicaxanthin-target complexes compared to the co-crystallized ligand.

Stability Analysis

The RMSD has been selected as a criterion to evaluate the dynamic stability of ligand-bound systems. The first systems took into consideration are Chain A and Chain B bound to the co-crystallized inhibitor K252a. We decided to perform these analyses to compare the behavior of Indicaxanthin with respect to an inhibitor with experimental evidence in terms of binding interactions. All protein frames are first aligned on the reference frame backbone, and then the RMSD is calculated based on the atom selection, in these cases on the Cα. For these complexes, the RMSD values of protein’s Cα atoms and ligand are reported in Figure 5. The Ligand RMSD indicates how stable the ligand is with respect to the protein and its binding pocket. This plot shows the RMSD of a ligand when the protein-ligand complex is first aligned on the protein backbone of the reference and then the RMSD of the ligand heavy atoms is measured. If the values observed are significantly larger than the RMSD of the protein, then it is likely that the ligand has diffused away from its initial binding site. The Chain A-K252a system reached equilibrium quickly and fluctuated around the average value of 3 Å, the low average value of ligand RMSD = 1.8 Å indicated strong stability of K252a in the binding pocket as expected due to the low number of rotatable bonds and the eight fused rings. The Chain B-K252a behavior is slightly different. The system reached equilibrium after ∼10 ns and the fluctuation of the protein is higher than the previous system analyzed (∼4.3 Å). The same evidence is reported for the average ligand RMSD that is higher (3.72 Å). This first analysis showed that the phosphorylation of the serine residues in the active form (Chain B) could confer more flexibility to the binding pocket in the KD, as confirmed by the higher fluctuations of K252a despite the rigid structure.

FIGURE 5

FIGURE 5. Time evolution of the RMSD (Å) values of backbone atoms (A) and ligand K252a (B) in the active form (Chain B) and the inactive form (Chain A).

For the Chain A-Indicaxanthin and Chain B-Indicaxanthin systems, the RMSD values of protein’s Cα atoms and ligand are reported in Figure 6 showing the behavior of these along the 100 ns of simulations. The Chain A-Indicaxanthin complex, after fluctuating in the first 10 ns, between RMSD values of 2-6Å, reached the equilibrium maintaining the fluctuation around the average value of 3.46Å. But it is worthy to note that ligand RMSD undergoes a biphasic behavior. In fact, in the first 38ns of simulation, the average ligand RMSD value is 6.0 Å, while in the rest of the simulation this value grows about 12 Å, stabilizing around ∼17.5 Å but the ligand continues to occupy the binding pocket acting as a lid. The behavior of the Chain B-Indicaxanthin complex is quite stable if compared to the previously discussed. The system reached rapidly the equilibrium around 7 ns and maintained stable fluctuations along the rest of the simulation with an average RMSD = 5.1 Å. The ligand RMSD, unless a quick cliff around 36–37.5 ns, maintains an average value = 3.2 Å that is similar to the behavior of the co-crystallized K252a. This value is interesting considering the less rigid structure of Indicaxanthin and explaining the good stability of this compound in the KD binding pocket.

FIGURE 6

FIGURE 6. Time evolution of the RMSD (Å) values of backbone atoms (A) and ligand Indicaxanthin (B) in the active form (Chain B) and the inactive form (Chain A).

For the complex of the Indicaxanthin docked into the allosteric site of Chain A, the protein maintains a stable behavior for the whole duration of simulation with an average RMSD value = 3.5 Å. In contrast, the ligand maintains its stability around the allosteric pocket for about 23 ns, but this stability is not maintained over time, and around 40 ns, Indicaxanthin flies out from the allosteric pocket until the end of the simulation. We repeated the simulations two times: in the second simulation the ligand flies out in the first nanoseconds, and in the third, it remains in the allosteric pocket for about 15 ns to fly out in the next nanoseconds (Supplementary Material). These last findings confirmed the BPDM results in which the unstable binding of Indicaxanthin in the allosteric pocket was retrieved.

Residues Mobility Analysis

To examine the structural flexibility effect of K252a and Indicaxanthin upon the Chain A, Chain B and the effect of Indicaxanthin when bound to the allosteric site per residue, the main chain average RMSF of the complexes have been calculated for the entire 100ns of simulation. The residue-wise fluctuation of the complexes was plotted and presented in Figure 7. The plot has been coupled according to the chain bound to K252a and Indicaxanthin, respectively, to compare the differences due to different ligands. As reported, the RMSF plot of Chain A is quite comparable for the residues of the active site with ∆RMSF <1 Å. The same could be observed for Chain B in which the ∆RMSF <0.5 Å. The major fluctuations happen in the allosteric pocket where the presence of Indicaxanthin determines an ∆RMSF>5 Å, in particular in the loop Thr368-Lys386. From this RMSF analysis, it can be concluded that the binding of Indicaxanthin determines similar fluctuations of the Chain A and Chain B backbone atoms with respect to K252a, in particular in the residues involved in the binding pocket. Contrariwise, the RMSF of the allosteric pocket is highly influenced by the presence of the Indicaxanthin.

FIGURE 7

FIGURE 7. The comparative RMSF values for the complexes of backbone atoms of Chain A (A) and Chain B (B) bound to K252a, Indicaxanthin, and of the allosteric pocket.

Protein-Ligands Contact Analysis

Estimation of protein interactions provides a measure of interaction power between the ligands and the target protein. Protein interactions with the ligand can be monitored throughout the simulations. These interactions can be categorized by type and summarized, as shown in the plot represented in Figure 8. Protein-ligand interactions (or “contacts’”) are categorized into four types: Hydrogen Bonds, Hydrophobic, Ionic, and Water Bridges. The stacked bar charts are normalized throughout the trajectory: a value of 1.0 suggests that 100% of the simulation time the specific interaction is maintained. A timeline representation of the interactions and contacts (H-bonds, Hydrophobic, Ionic, Water bridges) summarize the total number of specific contacts the protein makes with the ligand throughout the simulation. The bottom panel of Figure 8 shows which residues interact with the ligand in each trajectory frame. Some residues make more than one specific contact with the ligand, which is represented by a darker shade of orange, according to the scale to the right of the plot.

FIGURE 8

FIGURE 8. Protein-Ligand contacts: K252a in Chain A (A) K252a in Chain B (B).

Analyzing the trajectories of Chain A and Chain B with the co-crystallized ligand K252a, it is interesting to note that the residues involved in the interactions are variable both in terms of the type of contacts and in the time of interactions. As reported in Figure 8, K252a has four key residues involved: Met96, Glu97, Cys99, and Ile165 in both active and inactive forms. The H-bonds with Glu97 and Cys99 stable all 100ns simulations long. Hydrophobic contacts of Ile165 are maintained for the major part of the simulations, while hydrophobic contacts with Met96 are involved in both simulations but just for a little fraction of time. The differences in the two complexes regard the less interaction time of the H-bond with Glu149 in Chain B. In Chain A, an H-Bond with Gln175 is observed during the simulation in a discontinuous fashion, while in Chain B, Hydrophobic contacts are observed with Phe26 for about 50% of the simulation time. Other interactions are involved in both complexes such as with Leu21, Thr23, Val29, Ala42, Tyr98, Val152.

The behavior of Indicaxanthin, when bound to the two forms of IKK, is rather different in terms of interactions, but above in terms of residence in the binding pocket. As previously commented, in Chain A after ∼38 ns Indicaxanthin moves from the deep pocket but remains over it by interacting with H-bonds and ionic interaction with Arg427, Arg575, Arg579, Arg582 acting as a lid. In the first 40ns, the interactions of Indicaxanthin are both ionic and H-bond interactions involving the residues of the pocket such as residues 20–24 and 103–106. Indicaxanthin remains confined into the binding pocket for all the simulation time when Chain B is considered. Even though the residues involved in the protein-ligand interactions are different from those involved in the interaction with K252a. Indicaxanthin interacts for the major of time with Thr23, Asn28, Arg47, and Lys106 through H-bonds, but a key role for the binding stability is due to the ionic interaction with Asp103 and Lys 106 (Figure 9).

FIGURE 9

FIGURE 9. Protein-Ligand contacts: Indicaxantin in Chain A (A) Indicaxanthin in Chain B (B). Legenda is the same as Figure 8.

The stability of the Indicaxanthin’s binding is also mediated by several water bridges that are hydrogen-bonded protein-ligand interactions mediated by a water molecule. The hydrogen-bond geometry is slightly relaxed from the standard H-bond definition. The current geometric criteria for a protein-water or water-ligand H-bond are a distance of 2.8 Å between the donor and acceptor atoms (D—H···A); a donor angle of 110° between the donor-hydrogen-acceptor atoms (D—H···A); and an acceptor angle of 90° between the hydrogen-acceptor-bonded atoms (H···A—X).

Due to the evidence that during the simulation of Indicaxanthin into the allosteric pocket of Chain A, the ligand flies out after about 40 ns of simulation time, we do not comment on the evolution of the interaction in this analysis. Likely, Indicaxanthin does not act as an allosteric inhibitor of IKK.

MM-GBSA and Binding Free Energy Analysis

To understand the biophysical basis of molecular recognition of Indicaxanthin with targets (the inactive form of IKK, and active form of IKK), a molecular mechanics-generalized Born surface area MM-GBSA approach was used. As a comparative analysis, we performed the binding free energy analysis for the K252a complexes. It provides different individual components such as ∆G_vdW, ∆G_Coul, ∆G_Hbond, ∆G_Lipo. The conformational entropy change ‒TΔS can be computed by normal-mode analysis on a set of conformational snapshots taken from MD simulations, but many authors have reported that the lack of the evaluation of the entropy is not critical for calculating the MM-GBSA (or MM-PBSA) free energies for similar systems (Hou and Yu, 2007; Massova and Kollman, 2000; Wang and Kollman, 2001). For these reasons, we did not carry out ‒TΔS calculation in our study. For binding free calculation of Indicaxanthin and K252a/protein systems, 101 frames from 100 ns (each 1 ns) were retrieved to calculate ∆Gbind and individual contributions. A summary of binding components in the binding free energy is reported in Table 2. Furthermore, in Figure 10 it is reported the free energy landscape of Indicaxanthin in Chains A and B. To get a better view of which energy terms have more impact on the inhibitory potency, these individual energy components were compared. As expected from previous analysis of the MD trajectories, Indicaxanthin showed a higher binding affinity for the Chain B concerning Chain A (∆G_bind = −22.2 ± 4.3 kcal/mol, ∆G_bind = −20.7 ± 4.7 kcal/mol, respectively) that it could be justified by the shifting of Indicaxanthin towards the mouth of the binding pocket. It can be seen that ΔG_vdW has the major favorable contribution to the total free energy in both complexes. The same consideration could be stated for ΔG_Hbond and ΔG_Lipo, but there is no great difference among these values for different complexes even though it has minor weight on the total free energy. ΔG_Hbond has a major impact in Chain A, while ΔG_Lipo has a major impact in Chain B. The ΔG_Coul has the same effect in the complexes with a negative contribution to the global free energy which is partially balanced by the ∆G_Solv. As previously stated, Indicaxanthin has been demonstrated to inhibit 50% cell proliferation at 100 µM (Allegra et al., 2018), while K252a showed to be an efficient inhibitor of IKK at IC₅₀ = 0.09–018 µM (Sadler et al., 2004). Comparing the binding free energy values of Indicaxanthin and K252a, they are consistent with the experimental evidence (K252a ∆G_bind = −72.02 ± 4.1 kcal/mol for Chain B, and ∆G_bind = −74.45 ± 4.0 for Chain A) (Table 2).

TABLE 2

TABLE 2. Predicted MM-GBSA free energies (kcal/mol) and individual energy terms of the Indicaxanthin-target complexes and K252a-target complexes.

FIGURE 10

FIGURE 10. The binding energy (∆G expressed in kcal/mol) landscape of Indicacaxanthin bound to Chain A and Chain B.

Conclusion

The stability of Indicaxanthin-hIKKβ complexes compared to K252a, a co-crystallized inhibitor, was assessed by using Induced fit docking, binding pose metadynamics, and molecular dynamics. Finally, MM-GBSA free energy calculations have been performed to establish what form of IKK Indicaxanthin prefers. Induced fit docking results showed that the binding of Indicaxanthin with the active form, the inactive form, and the allosteric site of hIKKβ showed has the strongest stability with the active form. MD trajectories analysis (RMSD, RMSF, and protein-ligand contacts number and along the time) also showed that Indicaxanthin enhanced the stability of the active form at the same level of the known inhibitor K252a. The stability of the inactive form complex with Indicaxanthin is quite similar but it did not reach the quality of the active form. Contrariwise, even though for 40 ns over 100 nsIndicaxanthin can bind the allosteric pocket, it should not be considered as an allosteric inhibitor of hIKKβ.

As a whole, this work shows that Indicaxanthin is able, in silico, to inhibit the active form of the hIKKβ and adds novel mechanistic insights on its recently discovered ability to impair NF-κB signaling in melanoma A375 cells. Along these lines, present results further suggest the molecule as a useful nutraceutical tool in combo-therapy i.e., with other therapeutical agents targeting different checkpoints of melanoma development.

The currently demonstrated ability of Indicaxanthin to inhibit the active form of hIKKβ may, then, suggest the phytochemical as a new lead compound to synthesize novel and more potent IKKβ inhibitors for the treatment of cancer and inflammation-related conditions.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

MA and MT designed experiments. MT and GC carried out the experiments. MT and GC analyzed experimental results. MT and GC analyzed and interpreted the data. MA, AMA, LT, AA, and MT drafted and revised the manuscript. All authors contributed to the article and approved the submitted version.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2021.701568/full#supplementary-material

References

Allegra, M., Carletti, F., Gambino, G., Tutone, M., Attanzio, A., Tesoriere, L., et al. (2015). Indicaxanthin from Opuntia Ficus-Indica Crosses the Blood-Brain Barrier and Modulates Neuronal Bioelectric Activity in Rat hippocampus at Dietary-Consistent Amounts. J. Agric. Food Chem. 63, 7353–7360. doi:10.1021/acs.jafc.5b02612

PubMed Abstract | CrossRef Full Text | Google Scholar

Allegra, M., De Cicco, P., Ercolano, G., Attanzio, A., Busà, R., Cirino, G., et al. (2018). Indicaxanthin from Opuntia Ficus Indica (L. Mill) Impairs Melanoma Cell Proliferation, Invasiveness, and Tumor Progression. Phytomedicine 50, 19–24. doi:10.1016/j.phymed.2018.09.171

PubMed Abstract | CrossRef Full Text | Google Scholar

Allegra, M., Tutone, M., Tesoriere, L., Almerico, A. M., Culletta, G., Livrea, M. A., et al. (2019). Indicaxanthin, a Multi-Target Natural Compound from Opuntia Ficus-Indica Fruit: From its Poly-Pharmacological Effects to Biochemical Mechanisms and Molecular Modelling Studies. Eur. J. Med. Chem. 179, 753–764. doi:10.1016/j.ejmech.2019.07.006

CrossRef Full Text | Google Scholar

Amiri, K. I., and Richmond, A. (2005). Role of Nuclear Factor-Kappa B in Melanoma. Cancer Metastasis Rev. 24, 301–313. doi:10.1007/s10555-005-1579-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Amschler, K., Schön, M. P., Pletz, N., Wallbrecht, K., Erpenbeck, L., and Schön, M. (2010). NF-kappaB Inhibition through Proteasome Inhibition or IKKbeta Blockade Increases the Susceptibility of Melanoma Cells to Cytostatic Treatment through Distinct Pathways. J. Invest. Dermatol. 130, 1073–1086. doi:10.1038/jid.2009.365

CrossRef Full Text | Google Scholar

Belhadj Slimen, I., Najar, T., and Abderrabba, M. (2017). Chemical and Antioxidant Properties of Betalains. J. Agric. Food Chem. 65, 675–689. doi:10.1021/acs.jafc.6b04208

CrossRef Full Text | Google Scholar

Ben-Neriah, Y., and Karin, M. (2011). Inflammation Meets Cancer, with NF-κB as the Matchmaker. Nat. Immunol. 12, 715–723. doi:10.1038/ni.2060

PubMed Abstract | CrossRef Full Text | Google Scholar

Bowers, K. J., Chow, D. E., Xu, H., Dror, R. O., Eastwood, M. P., Gregersen, B. A., et al. (2006). “Scalable Algorithms for Molecular Dynamics Simulations on Commodity Clusters,” in SC '06: Proceedings of the 2006 ACM/IEEE Conference on Supercomputing (Institute of Electrical and Electronics Engineers (IEEE)), 43. doi:10.1109/sc.2006.54

CrossRef Full Text | Google Scholar

Clark, A. J., Tiwary, P., Borrelli, K., Feng, S., Miller, E. B., Abel, R., et al. (2016). Prediction of Protein-Ligand Binding Poses via a Combination of Induced Fit Docking and Metadynamics Simulations. J. Chem. Theor. Comput. 12, 2990–2998. doi:10.1021/acs.jctc.6b00201

PubMed Abstract | CrossRef Full Text | Google Scholar

Dhawan, P., and Richmond, A. (2002). A Novel NF-Kappa B-Inducing Kinase-MAPK Signaling Pathway Up-Regulates NF-Kappa B Activity in Melanoma Cells. J. Biol. Chem. 277, 7920–7928. doi:10.1074/jbc.M112210200

CrossRef Full Text | Google Scholar

Genheden, S., and Ryde, U. (2015). The MM/PBSA and MM/GBSA Methods to Estimate Ligand-Binding Affinities. Expert Opin. Drug Discov. 10, 449–461. doi:10.1517/17460441.2015.1032936

PubMed Abstract | CrossRef Full Text | Google Scholar

Guedes, I. A., Freitas, R. H., Cordeiro, N. M., do Nascimento, T. S., Valerio, T. S., Fernandes, P. D., et al. (2016). LASSBio-1829 Hydrochloride: Development of a New Orally Active N-Acylhydrazone IKK2 Inhibitor with Anti-inflammatory Properties. Chem. Med. Chem. 11, 234–244. doi:10.1002/cmdc.201500266

PubMed Abstract | CrossRef Full Text | Google Scholar

Jacobson, M. P., Pincus, D. L., Rapp, C. S., Day, T. J., Honig, B., Shaw, D. E., et al. (2004). A Hierarchical Approach to All-Atom Protein Loop Prediction. Proteins 55, 351–367. doi:10.1002/prot.10613

PubMed Abstract | CrossRef Full Text | Google Scholar

Jorgensen, W. L., Maxwell, D. S., and Tirado-Rives, J. (1996). Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 118, 11225–11236. doi:10.1021/ja9621760

CrossRef Full Text | Google Scholar

Leung, C. H., Chan, D. S., Li, Y. W., Fong, W. F., and Ma, D. L. (2013). Hit Identification of IKKβNatural Product Inhibitor. BMC Pharmacol. Toxicol. 14, 3. doi:10.1186/2050-6511-14-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, H., Liang, H., Meng, H., Deng, X., Zhang, X., and Lai, L. (2018). A Novel Allosteric Inhibitor that Prevents IKKβActivation. Medchemcomm 9, 239–243. doi:10.1039/c7md00599g

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, S., Misquitta, Y. R., Olland, A., Johnson, M. A., Kelleher, K. S., Kriz, R., et al. (2013). Crystal Structure of a Human IκB KinaseβAsymmetric Dimer. J. Biol. Chem. 288, 22758–22767. doi:10.1074/jbc.M113.482596

CrossRef Full Text | Google Scholar

Moarbess, G., Guichou, J. F., Paniagua-Gayraud, S., Chouchou, A., Marcadet, O., Leroy, F., et al. (2016). New IKK Inhibitors: Synthesis of New Imidazo[1,2-A]quinoxaline Derivatives Using Microwave Assistance and Biological Evaluation as IKK Inhibitors. Eur. J. Med. Chem. 115, 268–274. doi:10.1016/j.ejmech.2016.03.006

CrossRef Full Text | Google Scholar

Prescott, J., and Cook, S. (2018). Targeting IKKβin Cancer: Challenges and Opportunities for the Therapeutic Utilisation of IKKβInhibitors. Cells 7, 115. doi:10.3390/cells7090115

PubMed Abstract | CrossRef Full Text | Google Scholar

Sadler, T. M., Achilleos, M., Ragunathan, S., Pitkin, A., LaRocque, J., Morin, J., et al. (2004). Development and Comparison of Two Nonradioactive Kinase Assays for I Kappa B Kinase. Anal. Biochem. 326, 106–113. doi:10.1016/j.ab.2003.11.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Schrödinger (2017). Schrödinger. New York: Schrödinger Release, 2017–2021.

Sherman, W., Day, T., Jacobson, M. P., Friesner, R. A., and Farid, R. (2006). Novel Procedure for Modeling Ligand/receptor Induced Fit Effects. J. Med. Chem. 49, 534–553. doi:10.1021/jm050540c

CrossRef Full Text | Google Scholar

Staudt, L. M. (2010). Oncogenic Activation of NF-kappaB. Cold Spring Harb. Perspect. Biol. 2, a000109. doi:10.1101/cshperspect.a000109

PubMed Abstract | CrossRef Full Text | Google Scholar

Taniguchi, K., and Karin, M. (2018). NF-κB, Inflammation, Immunity and Cancer: Coming of Age. Nat. Rev. Immunol. 18, 309–324. doi:10.1038/nri.2017.142

PubMed Abstract | CrossRef Full Text | Google Scholar

Tesoriere, L., Allegra, M., Butera, D., and Livrea, M. A. (2004). Absorption, Excretion, and Distribution of Dietary Antioxidant Betalains in LDLs: Potential Health Effects of Betalains in Humans. Am. J. Clin. Nutr. 80, 941–945. doi:10.1093/ajcn/80.4.941

CrossRef Full Text | Google Scholar

Tutone, M., Pibiri, I., Lentini, L., Pace, A., and Almerico, A. M. (2019). Deciphering the Nonsense Readthrough Mechanism of Action of Ataluren: An In Silico Compared Study. ACS Med. Chem. Lett. 10, 522–527. doi:10.1021/acsmedchemlett.8b00558

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Herck, Y., Antoranz, A., Andhari, M. D., Milli, G., Bechter, O., De Smet, F., et al. (2021). Multiplexed Immunohistochemistry and Digital Pathology as the Foundation for Next-Generation Pathology in Melanoma: Methodological Comparison and Future Clinical Applications. Front. Oncol. 11, 636681. doi:10.3389/fonc.2021.636681

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, G., Lo, Y. C., Li, Q., Napolitano, G., Wu, X., Jiang, X., et al. (2011). Crystal Structure of Inhibitor ofκB Kinaseβ. Nature 472, 325–330. doi:10.1038/nature09853

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, J., and Richmond, A. (2001). Constitutive IkappaB Kinase Activity Correlates with Nuclear Factor-kappaB Activation in Human Melanoma Cells. Cancer Res. 61, 4901–4909.

PubMed Abstract | Google Scholar

Zhang, Q., Lenardo, D., and Baltimore, D. (2017). 30 Years of NF-κB: A Blossoming of Relevance to Human Pathobiologyy. Cell 168, 37–57. doi:10.1016/j.cell.2016.12.012

PubMed Abstract | CrossRef Full Text | Google Scholar