Computational study on the mechanism of small molecules inhibiting NLRP3 with ensemble docking and molecular dynamic simulations

Abstract

NLRP3 (Nucleotide-binding oligomerization domain, LRR and pyrin domain-containing protein 3) is a pivotal regulator of inflammation, with strong implications in gout, neurodegenerative diseases, and various inflammatory conditions. Consequently, the exploration of NLRP3 inhibitors is of great significance for the treatment of diseases. MCC950, NP3-146, compound (3), and YQ128 are four highly bioactive NLRP3 inhibitors that show great potential; however, their mechanism of action is currently limited to targeting the ATP binding region (NACHT site) of the NLRP3 protein. To gain deeper insights into the defining features of NLRP3 inhibitors and to develop more potent inhibitors, it is imperative to elucidate the interaction mechanism between NLRP3 and these inhibitors. In this study, we employ a comprehensive computational approach to investigate the binding mechanism between NLRP3 and representative inhibitors. Utilizing the molecular mechanics/generalized Born surface area (MM/GBSA) method, we calculate the binding free energy and pinpoint the key residues involved in the binding of the four inhibitors to NLRP3. The decomposition of binding free energy by the MM/GBSA method reveals that the residues Val71, Arg195, Ile255, Phe419, Arg422, and Met505, situated around the binding pocket, play a crucial role in conferring the high bioactivity of NLRP3 inhibitors. Furthermore, pharmacophore analysis of the four NLRP3 complexes indicates that the primary interaction between the inhibitors and NLRP3 was mainly hydrophobic interaction. Our study provides a profound understanding of the interaction mechanism between NLRP3 and its inhibitors, identifies the key residues involved, and provides theoretical guidance for the design of more efficient NLRP3 inhibitors.

Introduction

Inflammation is a common physiological and pathological process in the body, and the inflammatory corpuscles play an crucial regulatory role in inflammatory response [1]. To date, five types of inflammasomes have been identified: NLRP1, NLRP3, NLR family, CARD domain containing 4(NLRC4), Interstitial Pneumonia with Autoimmune Features(IPAF) and Absent In Melanoma 2(AIM2) [2]. Among these, NLRP3 is a key regulatory protein that participates in the formation of the NLRP3 inflammasome complex [3]. It has been observed that the activation of NLRP3 and its associated molecular regulatory signal pathways are closely associated with the onset and progression of various diseases, including nonalcoholic steatohepatitis (NASH), gout, porphyria-associated periodic syndrome, inflammatory bowel disease, and neurodegenerative diseases. These findings have garnered widespread attention and represent a frontier and hot area of clinical drug research and development [4, 5]. The NLRP3 sensor protein comprises three domains: The N-terminal pyrin domain (PYD), the nucleotide binding oligomerization domain (NACHT) with ATP catalytic function, and the C-terminal leucine-rich repeat (LRR) [6]. Under normal conditions, NLRP3 NLRP3 expression is low, and it is maintained in a self-inhibition state by the LRR region (Fig. 1). Upon recognition of pathogen-associated molecular patterns (PAMPs) or danger-associated molecular patterns (DAMPs) by the host, nuclear factor-κB (NF-κB) can be activated, initiating the transcription of NLRP3. The activated PYD of NLRP3 binds to the PYD domain of ASC (apoptosis-associated speck-like protein containing a CARD), leading to the formation of ASC specks through oligomerization. Subsequently, the CARD domain of ASC interacts with the CARD domain of pro-caspase-1. This oligomeric form of pro-caspase-1 undergoes self-cleavage to produce caspase-1, which then catalytically processes pro-IL-1β and pro-IL-18 into their active forms, ultimately triggering inflammation [7].

Fig. 1

A The structure of NLRP3, the inhibitor binds to the region of NACHT, purple cartoon represents NLRP3, blue sticks represent ligands MCC950 and ADP; (B) The functional area division of NLRP3; (C) The inhibitors studied in this work

At present, numerous NLRP3 small molecule inhibitors have been reported, with MCC950 being a notable example [8], MCC950 efficiently and selectively inhibits NLRP3, with an IC₅₀ of inhibiting IL-1β at the cellular level of 7.5nM [9]. It has also been shown to reduce mortality in newborn mice in a cryopyrin-associated periodic syndrome (CAPS) mouse model, demonstrating the compound’s efficacy [10]. Other NLRP3 inhibitors with good activity include NP3-146 [11], compound(3) [12] and YQ128 [13] are all NLRP3 inhibitors with well activity. They also exhibit significant differences in biological activity. Dekker et al. [11] reported the crystal structure of NP3-146 and NLRP3 complex for the first time, providing valuable insights into the molecular-level interaction mechanism between NLRP3 and its inhibitors. Prior to this structural insight, researchers had shown that these inhibitors primarily target the ATP-binding region of the NLRP3 protein (i.e., the NACHT domain), thereby inhibiting its ATP hydrolysis activity and subsequently preventing the oligomerization and assembly of the inflammasome complex. Therefore, studying the interaction mechanism between NLRP3 and its representative inhibitors is of great importance for understanding the key features of NLRP3 inhibitors and for the development of more effective NLRP3 inhibitors.

Molecular dynamics (MD) simulation and binding free energy calculation tools are widely utilized in drug design [14,15,16,17,18,19], providing valuable insights into the dynamic structure of protein-ligand interactions [20,21,22,23,24]. Based on MD trajectories, numerous methods have been developed to predict binding free energy, including free energy perturbation (FEP) [25], thermodynamic integration (TI) [26] and molecular mechanics/generalized Born surface area (MM/GBSA) method [27]. Compared to FEP or TI, MM/GBSA method offers a favorable balance between computational speed and accuracy, making it a popular choice for various protein-ligand systems, such as kinase systems, viral proteins, and other functional proteins [28,29,30,31].

Therefore, in this study, a comprehensive computational strategy was employed to elucidate the binding mode of NLRP3 and its representative inhibitors. Initially, the structures of the complexes formedby the four NLRP3 inhibitors— MCC950, NP3-146, compound (3) and YQ128 (Fig. 1) — with NLRP3 were obtained through molecular docking methods of varying precision. Subsequently, these four systems were subjected to conventional molecular dynamics simulations for 1 microsecond. Finally, by analyzing the trajectories from these simulations in conjunction with free energy calculations and energy decomposition, key residues that play crucial roles in the binding of inhibitors to NLRP3 were identified. Currently, there is limited research on NLRP3 computational studies; thus, gaining direct insights into the key interactions between NLRP3 and its inhibitors at the protein level will significantly advance the development of NLRP3 inhibitors.

Materials and methods

Ensemble docking

The docking of NLRP3 and ligands was performed using Schrodinger 2015. The crystal structure of the NLRP3 complex of the drug NP3-146, which is the focus of our study, has been determined. Therefore, we retrieved the initial structure of the NLRP3/NP3-146 complex from the RCSB Protein Data Bank (PDB ID code:7ALV [11]). Missing residues were added and aligned using Schrodinger 2015. Subsequently, the structure of the NLRP3/NP3-146 complex was prepared using the Protein Preparation Wizard, which involved adding side chains to residues, assigning protonation states, hydrogen atoms, and relaxing the side chains of the proteins.

To obtain a more accurate initial pose, induced Fit Docking (IFD) was performed in Schrodinger 2015 [32]. The protein molecule was minimized with an RMSD cutoff of 0.20 Å, and the centroid of the residues was automatically generated. The initial docking for each ligand was carried out using Glide. Residues within 5.0Å of the ligand position were carefully selected for the side chain optimization, ensuring that the structure and conformation adapted to each pose of the protein. The ligand was rigorously docked into the protein structure suitable for induction, and as a result, an IFD score was generated for each output pose.

Systems setup and molecular dynamics (MD) simulations

The partial charges of the four ligands were calculated at the HF/6-31G(d) level of theory and then fixed using the RESP methodology [33,34,35]. Each receptor-ligand complex was subsequently parameterized using the AMBER14SB [36] and GAFF force fields. The complexes were solvated with the TIP3P water model [37] in a 10Å cubic box using Leap, and Cl- ions were added to neutralize the net + 2 charge of the system(the salt concentration is about $\:2.5\times\:{10}^{-3}$M).

All of the MD simulations were performed using the AMBER20 package. Initially, a steepest-descent minimization scheme was applied to the systems for 40,000 steps. The systems were then gradually heated from 0 to 310 K over 100 ps in the NVT ensemble, with weak harmonic restraints of a constant force of 10 kcal/mol·Å² applied to the C and N atoms of the protein backbone. Subsequently, the restraints were gradually decreased over 0.9 ns from 10 to 0.01 kcal/mol·Å². Finally, 1 µs MD simulations were conducted at a temperature of 310 K and a pressure of 1 atm without any restraints. Throughout the simulation process, short-range non-bonded interactions were computed using a cutoff distance of 10Å, while long-range electrostatic interactions were treated using the Particle Mesh Ewald (PME) algorithm [38]. The SHAKE algorithm [39] was used to constrain all bonds involving hydrogen atoms, and the time step was set to 2 fs.

Thermodynamic calculation

The Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) approach [40, 41], commonly used to elucidate receptor-ligand interaction mechanisms, was employed to estimate the binding free energies for the protein-ligand complex [23, 42,43,44,45,46,47,48]. Within the MM/GBSA framework, the binding free energy is decomposed into several components. For the free energy calculations, 500 snapshots were extracted from the final 200ns of the MD trajectory for free energy calculations. These calculations were performed in AMBER20, and the binding free energy can be computed as follows:

$$\begin{aligned}< \Delta {G_{bind}} > & \\ & = < \Delta {E_{MM}} > + < \Delta {G_{solvation}} > - T < \Delta {S_{MM}}> \\ \end{aligned}$$

(1)

Where < ΔGbind > refers to the calculated average free energy, and < ΔEMM > refers to the average molecular mechanical energy.

$$\begin{aligned}< \Delta {E_{MM}} > & = < \Delta {E_{bond}} > + < \Delta {E_{angle}} > + < \Delta {E_{tors}} > \\ & + < \Delta {E_{vdw}} > + < \Delta {E_{elec}} > \\ \end{aligned} $$

(2)

$$\:<\Delta\:{G}_{solvation}>=<\Delta\:{G}_{GB}>+<\Delta\:{G}_{SA}>$$

(3)

<ΔGsolvation > refers to the desolvation free energy upon ligand binding. The polar contribution of desolvation (< ΔG_GB>) was calculated based on the Generalized Born (GB) model, with the igb parameter set to 2. The dielectric constants for solute and solvent were set to 1 and 80, respectively. The nonpolar contribution of desolvation (< ΔG_SA>) was determined by the solvent accessible surface area (SASA) using the LCPO method [49], with the formula ΔG_SA=0.0072 × ΔSASA. The normal-mode analysis [50] was utilized to estimate the conformational entropy contribution (−T < ΔS>) upon ligand binding, using 500 snapshots extracted from the MM/GBSA calculations.

Construction of pharmacophore modeling

AncPhore [51] is utilized to identify the key pharmacophore features of the interactions between receptors and ligands in the four systems. As a plug-in of PYMOL, AncPhore is implemented in C/C++, and automatically recognizes the essential pharmacophore characteristics for the receptor-ligand complex. The pharmacophore model structure is derived from the average structure of the last 200ns of equilibrium trajectory extracted from the MD simulation. The key pharmacophore features are then established and imported into PYMOL for visualization and mapping.

Results and discussion

Initial poses of inhibitors binding to NLRP3

As shown in Fig. S1, we obtained the bound modes of four inhibitors—MCC950(A), NP3-146(B), (3)(C) and YQ-128(D) —with NLRP3 using molecular docking techniques. MCC950 and NP3-146 share structural similarities, which manifest in their comparable binding interactions with NLRP3. Specifically, both inhibitors establish hydrophobic contacts with a series of the residues including Ala71, Gly73, Arg195, Pro196, Val197, Phe254, ILE255, Thr283, Tyr287, Thr368, Ile418, Arg 422, Phe423, Gln468 and Glu473. Additionally, they form hydrogen bond with Gly73 and Arg422.

However, subtle conformational differences between MCC950 and NP3-146 result in a distinct hydrogen bonding mode; MCC950 forms an extra hydrogen bond with Gln468 that is absent in NP3-146. In contrast, the binding mode of compound (3) with NLRP3 differs significantly from that of MCC950, despite remaining within the confines of the preset binding pocket. As depicted in Fig. S1, compound (3) exclusively forms a hydrogen bond with the residue Arg422, highlighting a shift in the key interaction points compared to MCC950 and NP3-146. The binding mode of YQ-128 with NLRP3 exhibits some similarities to that of compound (3), yet it also displays unique features. YQ-128 establishes hydrogen bonds with residues Ile214, Tyr476, and Asp506, indicating a distinct interaction profile compared to the other inhibitors. These preliminary differences in binding modes were identified through molecular docking, providing a snapshot of the potential interactions between the inhibitors and NLRP3. To delve deeper into these interactions and uncover more nuanced differences, we plan to conduct comprehensive molecular dynamics simulations on these four complex systems. This approach will allow us to gain a more detailed understanding of the dynamic nature of the inhibitor-NLRP3 interactions and potentially reveal new insights into the mechanism of action of these inhibitors.

Stability of MD simulations

The accuracy and reliability of the docking poses of four inhibitors with NLRP3 were assessed through all-atom explicit solvent MD simulations, and a total of 1µs of simulation trajectories was collected for analysis. To assess the stability of the four complexes, we monitored the root-mean-square deviation (RMSD) values computed on the backbone atoms of the protein and the heavy atoms of the ligand throughout the simulation process. As illustrated in Fig. S2, after approximately 500 ns of simulation, both the heavy atoms of the inhibitor and the residues within 5 Å of the binding pocket exhibited relatively small conformational changes, with RMSD values remaining below 1 Å. This indicated that each complex had reached a stable equilibrium state. However, it is worth noting that the residues within 5 Å of the ligand compound (3) and the heavy atoms of compound (3) displayed greater fluctuations compared to the other complexes. This is attributed to the inherent flexibility of compound (3). Despite these fluctuations, the RMSD values remained within a range of ± 1 Å, suggesting that the system had still achieved a dynamic equilibrium. Therefore, we can conclude that the docking poses of all four inhibitors with NLRP3 are reliable and stable under the conditions of our MD simulations.

Binding free energies calculations

The binding free energy of the four complexes was calculated using the MM/GBSA method to investigate the distinct mechanisms by which these inhibitors interact to NLRP3 and to quantify the contributions of various components to the binding free energy(ΔG_bind). The contributions of ΔG_bind are summarized in Table 1.

Table 1 Binding free energy predicted by MM/GBSA method

As evident from Table 1, both van der Waals interactions (ΔE_vdw) and electrostatic interactions (ΔE_ele) play crucial roles in facilitating the binding of the inhibitors to NLRP3. The van der Waals interactions contribute negatively to the binding free energy for all inhibitors, indicating their favorable energetic contribution to binding. Similarly, the electrostatic interactions also contribute negatively for MCC950 and NP3-146, but less so for (3) and YQ128, suggesting a relatively weaker electrostatic attraction in these latter two complexes. Notably, there is no significant variation in the calculated nonpolar component of the binding free energy (ΔG_{nonpolar = ΔEvdw + ΔGSA}) across the complexes. This suggests that the nonpolar interactions, which are primarily driven by van der Waals forces and surface area burial, contribute similarly to the binding affinity of all inhibitors. However, marked differences are observed in the polar component (ΔG_{polar = ΔEele + ΔGGB}), primarily stemming from differences in ΔE_ele and the polar solvation energy (ΔG_GB). The polar component contributes positively to the binding free energy for all inhibitors, but the magnitude of this contribution varies significantly. For MCC950, the polar component is relatively small and negative, indicating a favorable electrostatic interaction that stabilizes the complex. In contrast, for the other inhibitors, the polar component is larger and positive, suggesting a less favorable electrostatic interaction. These disparities result in substantial variations in the final binding free energy values. Specifically, the ΔG_bind for the NLRP3/MCC950 complex is markedly lower than that for the other three complexes. This finding indicates that MCC950 exhibits significantly higher binding affinity for NLRP3 compared to the other inhibitors, which aligns well with the experimental data on inhibitory activity. As the most potent NLRP3 inhibitor among those studied, MCC950’s binding to NLRP3 is characterized by a notable difference in the electrostatic component. The favorable electrostatic interaction in the MCC950 complex is likely due to a combination of favorable charge-charge interactions and a favorable polar solvation energy. Additionally, the van der Waals contribution to the binding free energy also constitutes a significant portion for MCC950, further stabilizing the complex.

On the other hand, the entropy contribution (-TΔS) remains relatively unchanged across the inhibitors. This suggests that the changes in binding affinity are primarily driven by changes in the enthalpic components of the binding free energy, rather than by changes in the entropy of binding. In the subsequent analysis, we will delve into the reasons behind these differences by examining the energy decomposition and binding modes of all inhibitors with NLRP3. This will provide further insights into the distinct mechanisms by which these inhibitors interact with NLRP3 and how their binding affinities are modulated by various energetic contributions.

Binding mechanism of inhibitors in the NLRP3

To conduct a more precise assessment of the distinct binding modes of the four inhibitors with NLRP3, we derived the average structure from the MD equilibrium trajectory and tracked the hydrogen bond occupancy between the receptor and ligand throughout this trajectory. As illustrated in Figs. 2 and 3 and a consistent binding feature of the inhibitors observed for all inhibitors with NLRP3, when compared to the initial structure, is the stable persistence of a hydrogen bond between residue 422 and the inhibitor throughout the entire MD simulation (Table 2), which aligns with the initial binding conformations. Specifically, in the MCC950/NLRP3 complex, hydrogen bonds between residue Ala72 and MCC950 exhibit occupancy rates of approximately 39% and 25%, depending on whether Ala72 acts as the acceptor or donor, respectively. Additionally, several hydrogen bonds with lower occupancy rates were detected between residue Asp506 and MCC950 (~ 34% and 29%). In the NP3-146/NLRP3 complex, stable hydrogen bonds formed by residues Arg422 and Asp506 with NP3-146 were also observed. These findings underscore the significance of hydrogen bonds in enhancing the stability of the interaction between inhibitors and the target protein. In the complex involving compound (3) and NLRP3, hydrogen bonds are formed between residues Arg422 and Val197 with compound (3); however, only the hydrogen bond contribution from residue Arg422 was noted in the YQ128/NLRP3 complex. These observations corroborate the numerical differences in the predicted binding free energies between NLRP3 and the four inhibitors.

Fig. 2

The binding mode of the stable conformations of the four complex systems obtained from the equilibrium trajectories of MD simulations, the green sticks represent key residues in the pocket of NLRP3 and the yellow sticks represent ligands. (A) NLRP3/MCC950 complex; (B) NLRP3/NP3-146 complex; (C) NLRP3/(3) complex; (D) NLRP3/YQ128 complex

Fig. 3

Residues and water that form hydrophobic and hydrogen bond with the inhibitor at the binding site predicted by the LigPlus. (A) NLRP3/MCC950 complex; (B) NLRP3/NP3-146 complex; (C) NLRP3/(3) complex; (D) NLRP3/YQ128 complex

Table 2 Summary of hydrogen bond between the ligand and protein

Per-residue energy contribution analysis of inhibitors binding to NLRP3

We employed the MMPBSA method to calculate the enthalpy change component of the binding free energy between the four inhibitors and NLRP3, and further decomposed this energy into contributions from individual residues. This approach allows us to assess the importance of each residue in mediating the interaction between NLRP3 and the inhibitors based on the energy contribution of a single residue. In Fig. 4, we present the residues whose absolute energy contribution to the binding of the four inhibitors exceeds 0.5 kcal/mol. The energy decomposition profiles for the NLRP3/MCC950 and NLRP3/NP3-146 complexes exhibit similarities, suggesting that MCC950 and NP3-146 bind to NLRP3 in comparable modes. Key residues involved in the interactions of MCC950 and NP3-146 with NLRP3 include Ala71, Ala72, Arg195, Ile255, Val258, Thr283, Thr368, Ile418, Phe419, Arg422, Gln468, Leu472, Glu473, Met505, and Asp506. Notably, the energy contributions of Glu473 and Asp506 are positive, indicating that these residues contribute unfavorably to the binding of MCC950 and NP3-146 to NLRP3. Among these, Arg422 contributes the most significantly to the binding of both inhibitors.In the compound (3)/NLRP3 complex, the residues Arg195 and Asp506 exhibit positive energy contributions, which are unfavorable for the binding of compound (3) to NLRP3. Similarly, in the YQ-128/NLRP3 complex, the residues Lys76, Arg195, Glu473, and Asp506 contribute positively, indicating an unfavorable interaction between YQ-128 and NLRP3. Furthermore, Fig. 2 reveals that the binding mode of YQ-128 with NLRP3 differs markedly from that of other three compounds, which aligns with our initial analysis of the binding modes of the four inhibitors with NLRP3.

Fig. 4

Per-residue interaction decomposition of the binding free energies for (A) NLRP3/MCC950, (B) NLRP3/NP3-146, (C) NLRP3/(3) and (D) NLRP3/YQ128 complexes

Common features shared by NLRP3-inhibitors recognition

The pharmacophore models for the interactions between the four inhibitors and NLRP3 were generated using the Pymol plug-in AncPhore, as illustrated in Fig. 5. For MCC950 (Fig. 5A), the dominant interaction is hydrophobic, involving the residues Met252, Ile255, Val258 and Tyr287. Additionally, the oxygen and hydrogen atoms in the side chain of residue Ala72 form hydrogen bonds with MCC950, acting as donors and acceptors respectively. The nitrogen atom in the side chain of residue Arg422 also forms a hydrogen bond with MCC950 as a donor, while the side chain of residue Asp506 forms a hydrogen bond with MCC950 as an acceptor. Due to structural differences, the pharmacophore model for NP3-146 and NLRP3 deviates significantly from that of MCC950. Specifically, NP3-146 exhibits hydrophobic interactions with residues Met252, Ile255, Thr283, and Thr368. Furthermore, residues Arg422 and Asp506 form hydrogen bonds with NP3-146, acting as hydrogen bond receptors and donors, respectively. The pharmacophore models for compound (3) and YQ-128 with NLRP3 primarily involve hydrophobic interactions. For instance, compound (3) interacts hydrophobically with residues Pro196, Val197, and Tyr476, and the residue Pro196 also forms a hydrogen bond with compound (3) as a receptor. Similarly, YQ-128 interacts hydrophobically with residues Ala71, Ala72, Glu473, and Ile255. The findings obtained from the pharmacophore model analysis are in agreement with the previous binding mode analysis, providing further validation of the interaction modes between the inhibitors and NLRP3.

Fig. 5

The pharmacophore model generated by the protocol AncPhore in Pymol. The purple sticks represent the key residues in the pocket of NLRP3, and the purple sticks represents the studied inhibitors MCC950 (in figure A), NP3-146 (in figure B), (3) (in figure C) and YQ-128 (in figure D)

Conclusion

In this study, we employed a comprehensive computational approach to elucidate the binding modes of NLRP3 with its representative inhibitors. The total binding free energy of MCC950, as predicted by the MM/GBSA method, was found to be lower than that of the other three inhibitors. This finding is in accordance with experimental observations, which indicate that MCC950 exhibits higher binding affinity compared to the other inhibitors. A detailed decomposition of the residue binding free energies revealed that specific residues located within the binding pocket, namely Ala71, Ala72, Arg195, Ile255, Val258, Thr283, Thr368, Ile418, Phe419, Arg422, Gln468, Leu472, Glu473, Met505, and Asp506, play a crucial role in determining the high bioactivity of NLRP3 inhibitors. By examining the conformations extracted from the equilibrium MD trajectories, we were able to identify key pharmacophore features of the inhibitors that interact with NLRP3. Notably, the primary interaction between the inhibitors and NLRP3 is hydrophobic, and the hydrogen bond formed between the residue Arg422 and the inhibitors is vital for enhancing the activity of the ligands. Our simulation results provide valuable insights into the mechanism of action of NLRP3 inhibitors and their interaction with NLRP3. These findings will undoubtedly facilitate the future development of highly selective inhibitors for NLRP3.