Background

Sepsis is a life-threatening condition triggered by the systemic response to infection; it disrupts host homeostasis, increasing the risk of mortality (Shi et al. 2021; Singer et al. 2016). Sepsis cases have increased sharply worldwide in recent years, affecting over 30 million patients and causing nearly six million deaths annually (Liu et al. 2024). NLRP3 functions as a multifaceted intracellular sensor that integrates different danger signals by detecting perturbations in cellular homeostasis. When triggered, the mechanism forms an inflammasome structure that processes and releases the key inflammatory cytokines IL-1β and IL-18 (Park et al. 2024). This cascade is mediated by NLRP3 inflammasome assembly and activation, which drives systemic inflammatory response syndrome (SIRS) and subsequently exacerbates tissue injury and multiorgan dysfunction (McCutcheon et al. 2025). The NLRP3 inflammasome serves as a sensor, integrating immune detection, regulatory processes, and inflammatory responses to mount and manage the defense mechanisms in the body. Consequently, pharmaceutical research has focused on targeting the NLRP3 inflammasome pathway as a promising approach to treat sepsis (Ma et al., 2023). The compelling preclinical efficacy of specific, high-affinity NLRP3 inhibitors, most notably MCC950, in various sepsis models provides robust pharmacological validation for this target (Yin et al. 2022). No clinical trials for NLRP3-targeting small-molecule drugs have been initiated, nor have any such molecules gained clinical approval (Dhruv et al. 2020). To systematically investigate the shared pharmacological properties of NLRP3-targeting small molecules and facilitate the development of novel therapeutic agents, in this study, we focused on NLRP3 inhibitors at both the preclinical and clinical stages of development. We collected a dataset of 34 NLRP3 inhibitors based on: (1) Clinical Relevance, focusing on compounds in advanced preclinical/clinical stages (Table 1); and (2) Structural Diversity, including various chemotypes (e.g., parthenolide) and both successful/discontinued candidates to enable a comprehensive mechanistic and pharmacological evaluation. Hence, new NLRP3 inhibitors need to be discovered to accelerate the development of sepsis drugs. In this study, we identified potential therapeutic agents for inhibiting NLRP3 using in silico techniques and cell models.

Table 1 NLRP3 inhibitors in clinical studies
Full size table

Materials and methods

Experimental procedures

This method allows for efficient and focused screening of a wide range of compounds to identify those with robust binding affinity and stability (Fig. 1).

Fig. 1
figure 1

Flowchart of the steps used to identify novel potential NLRP3 inhibitors

Full size image

Dataset preparation

Promising small-molecule inhibitors of NLRP3 were selected as ligands for the protein based on a literature review. Ligands were drawn using ChemDraw 19.0 and converted into *pdb format with Chem3D. The downloaded structure was checked for missing residues. Water molecules were removed, and hydrogen molecules were introduced to the small molecule ligand. All ligand geometries were further optimized using the semi-empirical PM6 method in Gaussian 16 to ensure accurate energy-minimized conformations and validated structural stability. Prior to docking, each ligand’s protonation state was determined at physiological pH 7.4 using the MarvinSketch pKa calculator to confirm the appropriate ionization state for biological relevance. Hydrogen molecules were added to the small molecules, and the ligands were selected.

Protein preparations

The NLRP3 protein structure as a receptor, characterized by its high-resolution architecture and structural integrity, encompasses critical functional domains-NACHT, LRR, and PYD-essential for studying ligand binding and molecular interactions. Thus, Thus, the three-dimensional crystal structures of NLRP3 (6NPY) were obtained from the Protein Data Bank (https://www.rcsb.org/) (accessed 14 April 2024) and formatted as *pdb. Water molecules were removed, and hydrogen molecules were introduced to the NLRP3 protein receptor. The preservation of secondary structure conformations before and after protein manipulation was confirmed using Ramachandran plots (Fig. S1). The receptor structure was optimized using AutoDock 4.2.7 and saved in *pdbqt format (Shahrami et al. 2021). Additionally, the electrostatic potential energy of NLRP3 was visualized using PyMOL.

Grid generation and molecular docking

The NLRP3 receptor and various small molecules available on the market underwent grid generation followed by molecular docking using automated procedures (spacing [angstrom]: 1.000; grid box: number of points in the X dimension: 85; number of points in the Y dimension: 75; number of points in the Z dimension: 60). The genetic algorithm and docking parameters of the search parameters in docking mode were adopted as default parameters. The binding energy was calculated using Lamarckian GA. The binding affinities were calculated using the AutoDock software. Along with hydrogen bonding, other binding forces at play between the ligands and the receptor, such as salt bridges, hydrophobic interactions, metal complexes, and halogen bonds, were considered.

To ensure statistical robustness and sampling diversity, each ligand was subjected to 50 independent docking simulations, generating 20 poses per run. The top 20 conformations were clustered based on RMSD (< 2.0 Å), and representative poses were selected from the lowest-energy clusters for further analysis.

The best docking pose for each ligand was chosen based on multiple criteria, including (i) minimum binding free energy, (ii) optimal number and geometry of hydrogen bonds, (iii) presence of key hydrophobic or electrostatic interactions within the NLRP3 active site, and (iv) consistency of poses within the dominant RMSD cluster.

To validate the docking protocol, redocking was performed using the co-crystallized ligand ADP from the NLRP3 crystal structure (PDB ID: 6NPY). The redocking RMSD between experimental and predicted conformations was 1.62 Å, confirming the reliability and accuracy of the docking parameters and scoring function.

The files in the *pdbqt format were converted into the *pdb format using OpenBabel 2.4.1, whereas the protein–ligand interaction profiler tool was used to elucidate the non-covalent interactions between molecules (https://plip-tool.biotec.tu-dresden.de/plip-web/plip/index) (Adasme et al. 2021). Visual inspection of the docking results involving active compounds and protein targets was performed using PyMOL software (version 2.2, https://PyMOL.org/2/) (Zou et al. 2023).

ADME/T and proprietary drug prediction

The drug-likeness of all compounds was assessed using ADMETlab 2.0 to evaluate their absorption, distribution, metabolism, excretion, and toxicological (ADME/T) properties. Prior to analysis, all ligand structures were standardized and energy-minimized using the PM6 semi-empirical method to ensure geometric consistency with docking outputs. The results of the analysis are detailed in Tables 4 and 5, and were generated using the ADMETlab 2.0 website (https://admetmesh.scbdd.com/) (Sankirtha et al. 2024).

ADME/T modeling was performed using both physicochemical descriptors and predictive models integrated within ADMETlab 2.0, including Lipinski’s Rule of Five, Veber’s rule, and BOILED-Egg models for permeability prediction. For each compound, the following parameters were calculated: Caco-2 permeability, 30% oral bioavailability (F30%), plasma protein binding (PPB), volume of distribution (Vd), cytochrome P450 metabolism (CYP2C9, CYP3A4, and CYP2D6), elimination half-life (t1/2), hepatotoxicity (H-HT), drug-induced liver injury (DILI), Ames mutagenicity, and maximum recommended daily dose (MRDD). Compounds meeting these criteria were considered pharmacologically favorable and selected for cytotoxicity assays and further DFT, MD correlation analysis.

The inclusion criteria for compound selection were:

  • 1. Docking score ≤ –6.0 kcal/mol (strong binding affinity),

  • 2. PPB < 90% (adequate free fraction for target distribution),

  • 3. Moderate or low hepatotoxicity (H-HT < 0.5), and

  • 4. Low predicted DILI probability (< 0.7).

Cytotoxicity

The THP-1 cells (Shanghai Yuchi Biotech Co., China) were cultured in RPMI-1640 medium (Corning, USA) supplemented with 10% fetal bovine serum (FBS) (Corning, USA), 5% amphotericin B-gentamycin solution (Coolaber, China), and β-mercaptoethanol (MCM, USA) in a humidified atmosphere of 5% CO2 at 37 °C. THP-1 cells (6 × 105 cells/mL) were differentiated with 30 ng/mL phorbol-12-myristate-13-acetate (PMA) (MCE, USA) for 48 h and then incubated for 24 h in complete PMA-free medium (Peng et al. 2024). After incubation with MCC950, ZYIL1, selnoflast, parthenolide, and glyburide (MCM, USA) for 24 h, the THP-1 macrophages were incubated with WST-8 (CCK-8) solution (10 μL) for 4 h. The absorbance was read at 450 nm.

$$\text{Cell Viability}\left(\text{\%}\right)=\frac{({\text{OD}}_{\text{treated}}-{\text{OD}}_{\text{blank}})}{({\text{OD}}_{\text{control}}-{\text{OD}}_{\text{blank}})}\times 100\text{\%}$$

Quantitative PCR

Total RNA was isolated from lung tissues with TRIzol reagent (Invitrogen, California, United States) and reverse-transcribed with Hifair® II 1st Strand cDNA Synthesis SuperMix for qPCR (Yeasen, China) following the manufacturer's instructions. Next, qPCR analysis was conducted using SYBR Green Master Mix (Yeasen, China), with GAPDH used as a control gene (Wang et al. 2023). All gene primer sequences were synthesized by Tsingke Biotechnology Co., Ltd. (Hubei, China) as follows:

Density functional theory (DFT)

A DFT study was performed using the Gaussian16 (Revision C.01) (https://gaussian.com/dl/g16c01.enw) program under the B3LYP-D3(BJ)/def2-SVP level-of-theory for optimization and frequency calculations (Smith et al. 2016). Prior to geometry optimization, a full Potential Energy Surface (PES) scan was performed to ensure that the lowest-energy conformer was selected for each ligand, thereby confirming the global minimum energy configuration.

The polarizable continuum model (PCM) with integral equation formalism (IEF-PCM) was adopted to simulate an aqueous dispersion medium corresponding to physiological conditions (dielectric constant ε = 78.5), ensuring that solvent polarization effects were accounted for instead of vacuum-phase approximation (Lu et al., 2024).

The electrostatic potential analysis was conducted by the Multiwfn (development version 3.8) program, with Visual Molecular Dynamics (version 1.9.4) used for better visualization (Humphrey et al. 1996). All optimized structures were verified through frequency analysis to confirm the absence of imaginary frequencies, validating that each optimized geometry represented a true energy minimum.

Molecular dynamics (MD)

MD simulations of the ligand-receptor docked complex were performed using GROMACS (version 2022.5). The protein–ligand complex was centered in a cubic box with a minimum distance of 1.2 nm from the box edge and then solvated with TIP3P water molecules. Na+ and Cl– ions were added for neutralization to achieve a NaCl concentration of 0.15 M. Energy minimization was performed using the steepest descent method with 500 steps and a convergence threshold of 100 kJ/mol/nm. Protein and ligand position restraints were applied during solvent relaxation. Non-bonded interactions were managed with a Verlet cutoff scheme at 0.9 nm for van der Waals and Coulomb cutoff distances. Long-range electrostatics were computed using the PME method. Subsequently, equilibration was conducted under NVT conditions for 125 ps with a 1 fs timestep. The temperature of the system was maintained at 310.15 K using a V-rescale thermostat with a coupling constant of 1.0 ps, which was applied separately to the protein and solvent. The initial velocities were assigned according to the target temperature, and the non-bonded interaction parameters were kept consistent with those in the minimization stage (Vagadia et al. 2016). Position restraints on the protein and ligand were retained during this equilibration phase. Finally, a production MD simulation under NPT conditions was performed for 200 ns with a 2 fs timestep after removing the position restraints. The system temperature was maintained at 310.15 K, while the pressure was controlled isotropically at 1 bar using a Parrinello–Rahman barostat with a coupling constant of 5.0 ps. The other non-bonded interaction parameters (cutoff scheme and distances) remained unchanged. The linear momentum was removed, and an angular momentum correction was applied to the solute to ensure the stability and physical accuracy of the simulation. The free energy of protein–ligand interactions was evaluated via the Poisson-Boltzmann surface area, and the binding energy was calculated using MMPBSA in the gmmpbsa package (Khan et al. 2020; Spoel et al. 2005).

$$RMSD=\sqrt{\frac{1}{N}{\sum }_{i=1}^{N}|{{r}_{final}\left(i\right)-{r}_{initial}(i)|}^{2}}$$
(1)

The RMSD at a given time point can be determined using Eq. (1), where rfinal (i) represents the final position of atom i, rinitial (i) denotes its starting coordinates, and N represents the total atomic count.

$$\text{Rg}=\sqrt{\frac{1}{\text{N}}}\sum_{\text{i}=1}^{\text{N}}{|\text{r}\left(\text{i}\right)-{\text{r}}_{\text{center}}|}^{2}$$
(2)

The Rg can be computed using Eq. (2), where r(i) corresponds to the coordinates of atom i, rcenter is the center of mass, and N represents the number of atoms in the molecule.

$$\text{SASA}=\sum (\frac{\text{R}}{\sqrt{{\text{R}}^{2}-{\text{Z}}_{\text{i}}^{2}}}\times \text{D}\times {\text{L}}_{\text{i}})$$
(3)

The SASA value at any point is determined using Eq. (3). Here, R represents the atom's radius, Li represents the arc length across a specific section i, and Zi represents the perpendicular distance of that section i from the center of the sphere.

$$\begin{aligned} & \Delta \text{ TOTAL}\left({\Delta \text{ G}}_{\text{MMGBSA}}\right)=\Delta \text{ GGAS}+\Delta \text{ GSOLV} \\ &=\Delta \text{ VDWAALS}+{\Delta \text{ E}}_{\text{EL}}+{\Delta \text{ E}}_{\text{GB}}+{\Delta \text{ E}}_{\text{SURF}} \end{aligned}$$
(4)

Results

Structural analysis of NLRP3

After receiving signals from infection or "self-injury", NLRP3 activates inflammatory caspases. NACHT-, leucine-rich-repeat- (LRR), and NLRP3 detect various stimuli, ranging from bacterial toxins and extracellular ATP to particulate substances. The NLRP3 protein is a member of the nucleotide-binding domain (NBD) and LRR-containing protein (NLR) families. The NACHT domain comprises helical domains 1 and 2, winged helix domains (WHD), and helical domain 2 (HD2) (Andreeva et al. 2021; Wang et al. 2021a, b). Analysis of the surface electrostatic potential energy distribution revealed two hydrophobic pockets between the NBD region and specific amino acid residues in the HD1 and LRR regions with high electrophilic potentials (Fig. 2A-C). The reliability of the docking approach was confirmed by initially docking the NLRP3 ligand ADP from the PDB database. The docking results revealed that the ligand ADP forms hydrogen bond interactions with ARG165, THR167, GLY227, LYS230, and ILE232 (Fig. 2D) (Sayed et al. 2023).

Fig. 2
figure 2

Structure of NLRP3 (PDB: 6NPY). A Regional distribution of the NLRP3 structure (blue: PYD and LRR; purple: NBD; green: HD1; orange: WHD; pink: HD2). B Electrostatic potential energy distribution of the NLRP3 protein. C Functional domains of NLRP3. D NLRP3 binds to ADP amino acid residues

Full size image

Molecular docking interactions of NLRP3 small molecules with the NLRP3 protein

To prioritize NLRP3 inhibitors with high clinical translational potential, we performed structure-based virtual screening. This computational strategy used molecular docking (AutoDock 4.2.7) to evaluate the binding affinities of 34 clinically relevant NLRP3 candidates documented in recent inflammasome-focused pharmacological studies. Molecular docking revealed that all compounds bound to the NLRP3 protein, with docking scores ranging from –7.23 to –2.82 (Table 2). Hydrogen bonds and hydrophobic interactions were observed between all compounds and proteins. The docking diagrams of the top 5 compounds are presented in Fig S2. Hydrogen bonds were predominantly concentrated in the NACHT domain, including GLU638 and ARG697. Only a few compounds interact with the LRR region. The hydrophobic interactions were located in the polybasic linker, NACHT domain, and LRR, such as ILE123, LYS124, and LYS163. Additional interaction forces included salt bridges and hydrophobic interactions. Some of the compounds presented metal complexes and halogen bonds. The salt bridge interaction occurred mainly in the NBD and WHD domain regions of NACHT.

Table 2 Sequences of the primers used for RT-qPCR
Full size table

To correlate the docking outcomes with pharmacokinetic performance, the top 10 compounds with the highest binding affinities were analyzed for their ADME/T properties. Compounds with both strong docking energies (≤ –6.0 kcal/mol) and favorable ADME/T profiles—specifically moderate plasma protein binding, low hepatotoxicity, and short half-life—were considered the most pharmacologically viable. Among these, ZYIL1 demonstrated an optimal combination of strong binding affinity (–7.23 kcal/mol), moderate PPB (63.35%), and low DILI probability (0.575), making it the top-ranked candidate. In contrast, MCC950, despite a strong docking score (–6.76 kcal/mol), exhibited high PPB (> 95%) and a shorter distribution profile, suggesting lower tissue permeability.

To establish a mechanistic link between electronic properties and pharmacological behavior, DFT calculations were used to evaluate the molecular polarizability and HOMO–LUMO energy gaps of the leading candidates. ZYIL1 exhibited the smallest energy gap (4.54 eV), indicating high electronic responsiveness and stronger interaction potential within the NLRP3 active site. This quantum-level finding directly supports its superior docking affinity and enhanced ADME/T performance, confirming a structure–activity relationship (SAR) in which improved electron delocalization contributes to more stable ligand–receptor binding.

Overall, these results demonstrate a consistent pattern: compounds with higher docking scores and better ADME/T characteristics exhibited favorable DFT parameters, reflecting complementary molecular stability and reactivity. This integrative correlation confirms that ZYIL1 combines optimal binding efficiency, physicochemical stability, and pharmacokinetic suitability—features that collectively define it as the most promising NLRP3 inhibitor among the screened molecules.

ADME/T analysis

We prioritized the top 10 compounds based on their highest docking scores for subsequent analysis. In the CaCo-2 model, most compounds exhibited poor permeability, except for RRx-001 and parthenolide. The parameter F30% is important for oral absorption, indicating poor oral bioavailability when the value is 0.9–1 (Tang et al. 2023). Intravenous administration is the primary route for drug delivery in sepsis patients. The evaluation of intestinal permeability (Caco-2 model) and oral bioavailability (F30%) is less relevant for these compounds (Waterbeemd et al. 2001). A high PPB indicates a low probability of drug escape from blood vessels, resulting in a limited distribution of the drug. NT-0249, 4b, EmlenoflastIZD-174, and glyburide exhibited high plasma protein binding (PPB > 90%), which may limit their distribution to target tissues because of the reduced free drug concentration. Potential drug Vd values fall within the range of 0.04–20.00 L/kg. All compounds presented Vd values between 0.325 and 2.071 L/kg in our study. In sepsis, a greater volume of distribution (Vd) is associated with improved drug penetration into infection sites, as it reflects greater tissue distribution and access to target tissues (Craig et al., 1998). The compounds were ranked according to their Vd values as follows: selnoflast, GDC-2394, parthenolide, ZYIL1, RRx-001, NT-0249, EmlenoflastIZD-174, MCC950, and glyburide, 4b. The compounds were predicted to be metabolic substrates for CYP2C9. Additionally, the treatment of sepsis requires rapid achievement of therapeutic drug concentrations and sustained maintenance of drug levels within the target therapeutic ranges (Grande et al., 2015). Therefore, compounds with relatively short half-lives (t1/2) are preferred for sepsis treatment. Pharmacokinetic analysis predicted that the t1/2 values of most candidate compounds fall within the range of 0.2 h to 0.8 h (Tables 3 and S2).

Table 3 Molecular docking interactions of NLRP3 small-molecule compounds
Full size table

Since patients with sepsis receive intravenous drug administration, F30% and Caco-2 permeability are not primary considerations (Takano et al. 2025). The top 10 ranked compounds aligned with the practical requirements for sepsis treatment. A PPB value exceeding 90% (MCC950, 4b, EmlenoflastIZD-174, and glyburide) may limit tissue distribution due to a decrease in free drug concentration. In contrast, ZYIL1, with a PPB of 63.35% and a Vd of 0.885 L/kg, not only facilitates distribution to target tissues but also promotes effective drug concentrations at the site of infection. Liver injury is common in sepsis, but developing targeted and efficient therapy remains a critical challenge (Tandon et al. 2021). Therefore, we focused on the effects of drug hepatotoxicity and drug-induced liver injury. The top five compounds with low hepatotoxicity were MCC950, 4b, glyburide, RRx-001, EmlenoflastIZD-174, and ZYIL1. The five compounds with the lowest risk of drug-induced liver injury (DILI) were identified as parthenolide, ZYIL1, selnoflast, MCC950, and GDC-2394. This finding indicates their wide safety margin and suitability for high-dose administration. Additionally, all compounds presented low hERG blocker and AMES toxicity, except RRx-001, which presented high AMES toxicity. Most of the compounds presented MRDDs (Table 3). Additional results are provided in Table S3.

Cell viability

Molecular docking and ADME/T analysis were conducted to predict the results. Considering the pharmacological requirements for sepsis treatment, we selected compounds with moderate t1/2 values, relatively high Vd values, and low risks of DILI and H-HT for cytotoxicity (Evans et al. 2021). ZYIL1, MCC950, selnoflast, parthenolide, and glyburide were selected for assessing their inhibitory potency and calculation of their IC50 values (Table 4). The five compounds exhibited inhibitory activities, with IC50 values of 5.74 ± 0.42 µM, 7.12 ± 0.36 μM, 11.59 ± 0.82 μM, 10.79 ± 0.55 μM, and 14.37 ± 0.76 μM, respectively. Cytotoxicity measurements were utilized to ascertain the effective concentration of the anti-inflammatory agent for subsequent assays. At the same concentration of 1.5 μmol/L, ZYIL1 significantly suppressed inflammatory cytokine expression compared to the other four compounds (Fig. 3). The preliminary result highlights the inhibitory ability of ZYIL1 and MCC950. These findings preliminarily suggested that ZYIL1 and MCC950 have superior inhibitory activity against the activation of the NLRP3 inflammasome.

Table 4 ADME/T values of selected NLRP3 small molecule compounds
Full size table
Fig. 3
figure 3

THP-1 cells were subjected to 6 h of 1 μg/mL LPS treatment, and the IL-1β (A), NLRP3 (B), IL-6 (C), and TNF-α (D) mRNA levels were measured via qPCR (n = 3)

Full size image

DFT analysis

To further investigate the potent inhibitory activity of ZYIL1 and MCC950, DFT calculations and MD simulations were performed to predict their binding mechanisms and interaction stability with the NLRP3 inflammasome. For a constant external electric field strength, the polarizability (α) is directly proportional to the size of the molecule and its volume. The volume of larger molecules results in greater deformation and higher polarizability. A higher polarizability indicates stronger molecular activity (Tandon et al. 2021). The electrostatic potential-colored molecular surface of MCC950 and ZYIL1 is illustrated in Fig. 4A. The oxygen atoms in the sidechain of MCC950 generate a negative region (colored from blue to white) within the hydrophobic area formed by hydrocarbon groups, whereas the sidechain of ZYIL1 does not have this characteristic. Compared to MCC950, the greater polarizability of ZYIL1 enhances its binding affinity to NLRP3. The MMPBSA results also support this result (Table 5). ZYIL1 shows a greater negative electrostatic potential surrounding its sulfonyl group, potentially enhancing electrostatic interactions with positively charged residues such as LYS124 and LYS163 in the NLRP3 binding pocket. This observation aligns with the docking results, revealing several hydrogen bonds and hydrophobic interactions with these residues (Table 2). The DFT results suggested a smaller HOMO–LUMO gap of 4.54 eV for ZYIL1 than for MCC950, which was 4.93 eV (Fig. 4B). Despite the difference in the sidechains, the LUMOs of both molecules are located on the π-conjugated system linked to the sulfone group, changing from –1.27 eV to –1.32 eV. The π-orbital of the benzyl ring acts as the HOMO in MCC950 with an energy of –6.20 eV, whereas the n-orbital (lone pair) of the amine nitrogen atom replaces that with a higher level of –5.86 eV, resulting in a smaller frontier orbital gap in ZYIL1. This result indicates that compared to MCC950, ZYIL1 is more polarizable to bind to the target of NLRP3 and therefore contributes to the greater activity of ZYIL1. To elucidate the structural basis of the binding affinity and stability observed in the docking and MD simulations, we correlated the electrostatic potential (ESP) profiles of ZYIL1 and MCC950 with their interaction patterns with NLRP3.

Fig. 4
figure 4

A Electrostatic potential-colored van der Waals surface of MCC950 and ZYIL1; B Molecular structure and HOMO–LUMO orbital isosurface (isovalue = 0.05 a.u.) of MCC950 and ZYIL1

Full size image
Table 5 The IC50 values of the five NLRP3 small-molecule compounds
Full size table

Molecular dynamics (MD) analysis

We investigated the stability of the interactions between NLRP3 and the small-molecule inhibitors MCC950 and ZYIL1 by calculating the RMSD, Rg, and SASA. Due to the binding of MCC950 and ZYIL1 to NLRP3, the ligand–protein complexes exhibited high stability within the range of 0.7 Å (Fig. 5A). Overall, ZYIL1 demonstrated greater stability than MCC950. Under external force, the complexes formed between MCC950, ZYIL1, and NLRP3 were 3.42 ± 0.02 nm and 3.41 ± 0.02 nm, respectively. A fluctuation was detected in MCC950 at 20 ns, followed by stabilization at 535 nm2. ZYIL1 also gradually decreased at 40 ns, followed by stabilization at 535 nm2 (Fig. 5B–C; Fig. S3). Rg and SASA were not influenced by the two small-molecule inhibitors of NLRP3. Figure 5D shows the change in hydrogen bonding between the protein and water, which exhibited a complementary relationship. To better understand the interaction between ligands and proteins, as well as to study the relative affinity of ligand–protein binding sites and gain insights into the contribution of ligands to specific residues, we used binding free energy calculations and residue contributions. ZYIL1 had greater binding energy (–49.46 kcal/mol) than MCC950 (–48.91 kcal/mol) (Table 6). To simplify the analysis, the residue contributions from MCC950 and ZYIL1 are presented in Fig. 3E and 3F. Since the binding energies of the five positively charged residues GLU174, GLU178, GLU367, GLU627, and ASN654 are all positive, they are not conducive to the stability of the complex.

Fig. 5
figure 5

Molecular dynamics analysis of ligand-receptor complexes. A RMSD reflects the stability of ligand-receptor complexes. C SASA reflects the size of the solvent-accessible surface area of responsive complexes. D Time evolution of the number of hydrogen bonds formed between the protein and ligands. E and F Free energy contribution reflects the specific energy contribution of key amino acid residues to ligand binding (blue: MCC950; red: ZYIL1)

Full size image
Table 6 MMPBSA binding free energies (kcal/mol) for MCC950 and ZYIL1
Full size table

Discussion

The NLRP3 inflammasome plays a key role in the pathogenesis of sepsis by amplifying inflammatory responses through the activation of caspase-1 and the subsequent release of pro-inflammatory cytokines such as IL-1β and IL-18 (Prakash et al. 2023; Vigneron et al. 2023). Although they are involved in septic inflammation, no specific NLRP3 inhibitors have been approved for clinical use. These findings highlight the urgent need to develop novel therapeutic agents targeting this pathway.

To ensure that our docking predictions are reliable, we performed redocking validation using the co-crystallized ligand ADP. After a thorough literature review, 34 compounds were selected for potential docking with the NLRP3 molecule and subsequently ranked based on their binding energy. The top three compounds identified all fell within the sulfonamide structures. In the context of sepsis treatment, swiftly attaining therapeutic efficacy while mitigating the potential for sepsis-associated hepatotoxicity (H-HT) and DILI to the greatest extent possible is imperative. A thorough evaluation encompasses the pharmacokinetic/pharmacodynamic (ADME/T) profiles of medications, their hepatotoxic potential, and the likelihood of DILI, all of which are important considerations for sepsis therapeutic agents. Our analysis focused on the top 10 compounds (ZYIL1, MCC950, selnoflast, RRx-001, 4b, NT-0249, parthenolide, elenoflastIZD-174, GDC-2394, and glyburide), which presented the most elevated molecular docking binding energies.

We demonstrated concordance with established experimental and clinical data. The high PPB (> 90%) predicted for MCC950 aligned with its known pharmacokinetic profile, which reported extensive protein binding and a long in vivo half-life in preclinical models (Rebecca C et al., 2015). Similarly, the high risk of DILI predicted for glyburide (probability: 0.975) was consistent with its documented hepatotoxic potential in pharmacovigilance databases and clinical case reports (Wang X et al., 2025). However, some discrepancies were noted and may reflect the limitations of predictive algorithms or species-specific differences. For example, while our model identified selnoflast as a CYP2C9 substrate, clinical phase I data suggest that its primary metabolism is mediated by CYP3A4 (Barbara et al. 2023). Despite these isolated variances, the overall strong correlation between our predictions and the literature for reference compounds enhances confidence in the projected favorable profiles of our top candidates, particularly ZYIL1, which has a balanced combination of moderate PPB, a suitable volume of distribution, and low predicted DILI risk—a profile meriting further experimental validation. ZYIL1 exhibited greater inhibitory activity in vitro. Functional assessments of the top five compounds revealed that ZYIL1 also demonstrated the most potent suppression of key inflammatory mediators associated with NLRP3 activation. The integration of DFT and MD simulation findings provided additional support for this observation. DFT analysis showed that ZYIL1 possessed a smaller HOMO–LUMO gap (4.54 eV) and higher polarizability than MCC950, indicating greater electronic flexibility and stronger electrostatic complementarity with positively charged residues (such as LYS124 and LYS163) in the NLRP3 binding site.

These quantum–mechanical properties were consistent with the molecular dynamics results, where ZYIL1 maintained lower RMSD fluctuations and a more stable complex configuration throughout 200 ns simulations. MM-PBSA energy calculations confirmed this stability, with ZYIL1 exhibiting a slightly higher binding free energy (–49.46 kcal/mol) compared to MCC950 (–48.91 kcal/mol). Persistent hydrogen-bond interactions between ZYIL1 and residues GLU228B, LYS124B, and LYS163B were observed, remaining intact for most of the simulation time.

This direct correspondence between DFT-derived electrostatic potential maps and MD-observed hydrogen-bond persistence indicates that ZYIL1’s enhanced electronic polarizability translates into superior conformational stability within the NLRP3 active site. Together, these results confirm that the quantum-level reactivity and dynamic binding behavior of ZYIL1 jointly contribute to its high inhibitory potency against NLRP3.

Although ZYIL1 is under clinical investigation for cryopyrin-associated periodic syndrome (CAPS), its potential application in sepsis remains unexplored [15,18]. Our results highlight its promising profile as a suitable candidate for treating sepsis, warranting further mechanistic and in vivo studies to comprehensively elucidate its mode of action and therapeutic efficacy in sepsis models.

Conclusion

The management of sepsis relies on conventional strategies, including antibiotics, fluid resuscitation, and organ support. In this study, we presented a comprehensive and innovative structure-based virtual screening strategy to systematically evaluate the efficacy of NLRP3 small-molecule inhibitors for sepsis treatment. By integrating molecular docking, ADME/T profiling, biochemical activity, and DFT and MD simulations, we thoroughly characterized the binding affinities, pharmacokinetic properties, and mechanistic interactions of 34 clinically relevant NLRP3 inhibitors. Using this multifaceted approach, we not only identified ZYIL1 and MCC950 as promising candidates with high binding stability and favorable drug-like properties but also obtained deeper insights into the molecular basis of NLRP3 inhibition through ESP analysis and binding free energy calculations. The greater performance of ZYIL1 highlights its potential as a novel therapeutic agent. We conducted a preliminary evaluation of NLRP3 inhibitors through an integrative computational and experimental approach. We established a foundational framework for the rational development of NLRP3-targeted therapeutic agents in sepsis, a field that lacks clinically approved inhibitors. To offer a more comprehensive and reliable evaluation of the potential of the molecule, in vitro and in vivo assays need to be conducted to experimentally confirm our computational findings.