Virtual screening of flavonoids as potential RIPK1 inhibitors for neurodegeneration therapy

Ligand library

We downloaded a flavonoid ligand library of 4,858 compounds previously prepared by Jiménez-Avalos et al. (2021) (File S1). Using the Open Babel version 3.1.1, we created sdf files from SMILES (O’Boyle et al., 2011). Next, pdbqt files with 3D coordinates were generated in a multithreaded parallelized ligand preparation pipeline (Samdani & Vetrivel, 2018) leveraging Open Babel (O’Boyle et al., 2011) and GNU parallel (Tange, 2018). Here, the MMFF94 forcefield was used and energy minimization was done with the steepest descent algorithm.

Protein 3D structure

We downloaded a 3D structure of the human RIPK1 from the RCSB Protein Data Bank (PDB ID: 6NYH). In a previous study, we used this pdb file for molecular docking (Bepari et al., 2022). The x-ray crystal structure has a resolution of 2.10 Å. However, there are many unmodelled residues and several partially modeled residues. Therefore, we generated a template-based homology model of the RIPK1 kinase domain through the SWISS-MODEL webserver (https://swissmodel.expasy.org/) (Waterhouse et al., 2018). Here, we used the target sequence retrieved from UniProt (https://www.uniprot.org/uniprotkb/Q13546/entry#sequences), and 6NYH was fed as the template. Next, we refined the pdb file for docking in UCSF Chimera version 1.17.1 by removing the ligand and adding hydrogen atoms and charges. Finally, a pdbqt file of the macromolecule was generated using the AutoDockTools (Morris et al., 2009).

Molecular docking

We adopted the docking protocols from a previous study (Bepari et al., 2022). We utilized a 26*26*26 Å gridbox, centered at the XYZ coordinates −17.094, 13.895, and −37.761, corresponding to the center of the reference ligand L8D in 6NYH. We employed AutoDock Vina (Trott & Olson, 2010) for molecular docking, which was parallelized using POAP (Samdani & Vetrivel, 2018). The number of individual runs (exhaustiveness) for each docking simulation was eight. Following the virtual screening, POAP provided us with a ranked list of ligands based on their binding energy (kcal/mol). The pdb files of the top complexes were used for further analysis and MD simulations.

Pharmacokinetic screening

We predicted the pharmacokinetic parameters of the top 20 compounds using the SwissADME website (Daina, Michielin & Zoete, 2017). A list of canonical SMILES was provided as the input, and the prediction results were downloaded as a CSV file.

Molecular dynamics simulations

For MD simulations, we used the RIPK1 complexes with the top hits from docking simulations to generate input files through CHARMM-GUI (Jo et al., 2008). We initiated the solution builder tool (Jo et al., 2008) to prepare the protein solution system from the pdb file of the complex. For the ligand, the hetero chain residues were read from a mol2 file previously prepared using Open Babel (O’Boyle et al., 2011), and CHARMM processes were invoked. We adopted the default parameters with some minor modifications. Briefly, an octahedral waterbox with an edge distance of 10 Å was specified to fit the protein. Forty-seven potassium (K⁺) and 40 chloride (Cl⁻) ions were placed using the Monte-Carlo method to neutralize the system. The number of solvent molecules ranged from 14,203 to 14,238.

Energy minimization was set to perform a maximum of 50,000 steps with the steepest descent algorithm. Before proceeding, the system’s maximum force (Fmax) needed to be below 1,000 kJ mol⁻¹ nm⁻². The system was then equilibrated for 100 ps at 300 K using the V-rescale thermostat (NVT: constant number of particles, volume, and temperature). Subsequent simulation steps were run with the leap-frog integrator. Applying the lincs constraint algorithm, the force constants for backbone atoms (BB) and side-chain atoms (SC) were set to 400.0 and 40.0 kJ/mol/nm^2, respectively. Placing the same constraints, 100-ps NPT equilibration was performed at 300 K and one bar using the V-rescale thermostat and the C-rescale pressure coupling (NPT: constant number of particles, pressure, and temperature). The equilibrated system was then subjected to the 500-ns production MD runs (two replicates for each system).

We mostly used the gmx analysis tools of Gromacs to analyze the trajectories. For root-mean-square deviation (RMSD) and root-mean-square fluctuation (RMSF) calculations, we used reference structures drawn from average structures from the trajectories spanning 200 to 500 ns and fitted to the protein backbone. To determine clusters of conformation, we used the gmx cluster tool with an RMSD cutoff of 0.15 nm. Clustering was done for the last 400 ns, and every 500th frame was analyzed. We computed the solvent accessible surface area (SASA) using the gmx sasa tool. For plotting RMSD, the radius of gyration (Rg), and SASA over time, we used the simple moving average of 100 data points.

We have uploaded the topology files of the best ligand for both replicates as raw data (File S2). The initial and final pdb structures for the complexes of the best ligands are available as raw data (File S3). The trajectories for both replicates are available in Figshare (DOI: 10.6084/m9.figshare.24455332).

Binding free energy

To calculate binding free energies for ligands, we employed the Molecular mechanics/Poisson–Boltzmann surface area (MMPBSA) method (Miller et al., 2012) using gmx_MMPBSA v1.5.0 (Valdés-Tresanco et al., 2021). We extracted 1,000 frames from the last ten ns of the trajectory for the MMPBSA calculations. The binding free energy of the ligand (ΔTOTAL) was calculated by subtracting the free energies of the receptor and ligand from the complex. For each system, the total free energy was the sum of the gas phase energy (GGAS) and the solvation energy (GSOLV). GGAS was contributed mainly from van der Waals interactions (VDWAALS) and electrostatic interactions (EEL), and GSOLV was the sum of polar (EPB) and nonpolar (ENPOLAR) interactions (Miller et al., 2012).

Data visualization

Data visualization and statistical analysis were performed using open-source Python packages such as pandas, matplotlib, seaborn, and pingouin in a Jupyter Notebook (Kluyver et al., 2016).

Homology modeling of RIPK1

The human RIPK1 protein has 671 amino acids and a mass of 75,931 Da, where the kinase domain resides at the N-terminal (amino acids position 17–289). We selected the x-ray crystallographic structure 6NYH (Fig. 1A) as the template for homology modeling. The 3D structure of the kinase domain we obtained from SWISS-MODEL spanned the amino acid range 8–294, and the missing amino acid residues of the template were successfully modeled (Fig. 1B). The co-crystallized ligand L8D, a potent RIPK1 inhibitor, served as a reference for our docking and molecular dynamics simulations. The inhibitor L8D was buried in a hydrophobic pocket of RIPK1 (Fig. 1C) (Xie et al., 2013).

Virtual screening of flavonoids by molecular docking

We first redocked the reference ligand L8D to the RIPK1 homology model. The docking pose revealed that our AutoDock Vina docking protocol can reliably imitate the binding conformation of the inhibitor observed in the crystal structure 6NYH (Fig. 1D).

Next, we moved to the initial screening of the flavonoid ligand library with docking simulations. The curated library consisted of 4,858 compounds: 119 from the CHEMBL (https://www.ebi.ac.uk/chembl/), four from the DrugBank (https://go.drugbank.com), and the remaining 4,735 from the Metabolomics (http://metabolomics.jp/wiki/Main_Page) database. For Vina docking, we limited the grid box to the L8D binding pocket (Fig. 1C), which is also the binding pocket for necrostatins, a group of compounds that prevent necroptosis mainly through RIPK1 inhibition (Xie et al., 2013).

We short-listed the top 20 compounds based on the Vina docking scores (kcal/mol), with more negative values indicating stronger binding (Table 1). Scores ranged from −11.7 to −10.6 kcal/mol, and molecular weights ranged from 360 to 724. All the top hits are available in the Metabolomics database: five anthochlors (FL1), four flavanones (FL2), five flavones (FL3), two dihydroflavonols (FL4), two flavonols (FL5), and two isoflavonoids (FLI).

Pharmacokinetic evaluation of top flavonoids

We predicted physicochemical descriptors and pharmacokinetic properties of the top 20 hits from molecular docking (Table S1). Oral bioavailability and blood-brain-barrier permeation are illustrated in Fig. 2A using the BOILED-Egg model (Daina & Zoete, 2016). Twelve molecules showed optimal parameters for gastrointestinal absorption. In contrast, only three molecules, Nitiducarpin (ID: FLID1CNP0005, PubChem CID: 630656), Pinocembrin 7-O-benzoate (ID: FL2FA9NS0008, PubChem CID: 42607888), and Paratocarpin J (ID: FL2FAANP0013, PubChem CID: 15200532), were candidates for blood-brain-barrier penetration. P-glycoprotein (Pgp) transporters translocate substrate out of the cell and can compromise drug actions. SwissADME (Daina, Michielin & Zoete, 2017) revealed that Paratocarpin J is a Pgp substrate while Nitiducarpin and Pinocembrin 7-O-benzoate are not. In the SwissADME output, the physicochemical space for oral bioavailability is further illustrated in radar charts (Fig. 2B), which consolidated six molecular descriptors: lipophilicity (LIPO) as a function of Log P_o/w (XLOGP3), molecular weight (SIZE), polarity (POLAR) as a function of topological polar surface area (TPSA), insolubility (INSOLU) as a function of Log S (ESOL), insaturation (INSATU) as a function of the fraction of sp3 hybridized carbon atoms (Fraction Csp3), and flexibility (FLEX) as a function of the number of rotatable bonds. The radar charts revealed that the three best hits satisfied nearly all the criteria for the above parameters except that Pinocembrin 7-O-benzoate has a higher insaturation than allowed (Fig. 2, Table S1). At this stage, we considered these compounds, Nitiducarpin, Pinocembrin 7-O-benzoate, and Paratocarpin J, the best ligands targeting the human RIPK1 protein.

Molecular dynamics simulations for top ligands

To validate RIPK1 interactions with the top three ligands, we performed replicated 500-ns MD simulations using Gromacs. As described previously, the co-crystallized ligand L8D in 6NYH was a positive control for the inhibition of RIPK1 (Bepari et al., 2022). RIPK1 complexes with the ligands were placed in an octahedral box, which we solvated with water and then neutralized the system with KCl. Using the steepest descent algorithm, we energy-minimized the protein-ligand complexes (Fig. S1). For the RIPK1-L8D complex, the steepest descents converged to Fmax < 1,000 within 264 steps. The complex with Nitiducarpin also converged very quickly, within 316 steps. The Pinocembrin 7-O-benzoate and Paratocarpin J complexes took 1,004 and 299 steps, respectively, for energy minimization (Fig. S1). Once energy minimization succeeded, the protein complexes in solutions underwent sequential equilibration simulations using NVT and NPT ensembles (Fig. S1). With negligible variations, the MD systems equilibrated at the reference temperature of 300 K and the reference pressure of 1 bar within 100 ps (Fig. S1).

Next, we performed 500-ns production MD runs at NPT. Afterward, the simulation trajectories were analyzed for the stability of the RIPK1 complexes and protein-ligand interactions. The RMSD of the protein backbone atoms (Fig. 2A) is a powerful predictor of protein stability. A low and stable RMSD indicates stable protein conformations within the complex, while a high RMSD points to conformational instability. The reference ligand L8D (also known as GNE684) is a potent pharmacological inhibitor of human RIPK1 in vitro (Patel et al., 2020). A low RMSD of the protein backbone in the RIPK1-L8D complex (Fig. 3A) thus strongly aligned with the experimental data (Patel et al., 2020) and validated our MD simulation protocol. The mean protein RMSD was 2.048 Å for L8D, and the RMSD converged very quickly (Fig. 3A). Interestingly, the RMSD for the Nitiducarpin complex converged within the first 100 ns and had the lowest value (1.713 ± 0.388 Å) among the studied complexes (Fig. 3A).

The protein backbone in the Pinocembrin 7-O-benzoate complex also showed considerable stability with a mean RMSD of 2.080 Å. On the other hand, the backbone RMSD was the highest in the Paratocarpin J complex, and we observed convergence only after 400 ns (Fig. 2A).

We next calculated the ligand RMSD (Fig. 3B), a robust metric for validating ligand binding stability following molecular docking simulations (Guterres & Im, 2020). Again, the reference ligand L8D exhibited minimal fluctuations from the protein backbone, with a mean RMSD of 0.972 Å (Fig. 3B), indicating the stability of the ligand within the binding pocket. Nitiducarpin briefly exited the binding cavity at around 210 ns in one of the replicates (Replicate-2) (Fig. 3B). Nonetheless, Nitiducarpin had the lowest ligand RMSD (1.191 ± 0.498 Å) and achieved the fastest convergence among the three hits in our study (Fig. 2B).

Although Pinocembrin 7-O-benzoate exhibited low fluctuations in one simulation, the RMSD varied considerably in the replicated simulation. Paratocarpin J displayed higher structural fluctuations early in the MD simulations. The RMSD converged below the mean at 50 ns in Replicate-1 and at 280 ns in Replicate-2 (Fig. 3A). Taken together, the protein backbone and the ligand RMSD values convincingly demonstrate the stable binding of Nitiducarpin in the inhibitor-binding cavity of RIPK1.

We then calculated the clustering of RIPK1 complexes by evaluating the RMSD values (Fig. 4A). Structures were grouped if their RMSD values did not exceed 0.15 nm. Additionally, we derived the average 3D structures for the clusters (Fig. 4B). Notably, the RIPK1-Nitiducarpin exhibited a high degree of stability, displaying predominantly a single conformation throughout the entire simulation period (Fig. 4). For the reference ligand L8D, we identified ten and three clusters in Replicate-1 and Replicate-2, respectively. In one replication, the RIPK1 complex with the ligand Pinocembrin 7-O-benzoate underwent a few initial conformational changes and then stabilized at the 4^th cluster for the remainder of the simulation. The same complex transitioned through 16 clusters in another replication, with rapid conformational changes occurring predominantly within the last 140 ns. Conversely, the Paratocarpin J complex did not appear to converge. Instead, it formed a rapidly evolving series of protein-ligand conformations throughout the simulation period, as illustrated in Fig. 4. Remarkably, conformational changes were evident primarily at the activation loop of RIPK1 (residues 156–196) (Fig. 4B, orange).

For the top complexes, we estimated local fluctuations of the RIPK1 C-alpha atoms by evaluating the per residue RMSF values (Fig. 5A). Generally, loose ends and loop regions of proteins undergo higher conformational changes. The RIPK1 P-loop corresponding to the residues 24–31 displayed RMSF values in the range of 4 to 6 Å for all complexes. However, the maximum conformational rearrangements were evident for the activation loop residues (amino acids 156–196), the Paratocarpin J complex generating the highest peaks in both replicates (Fig. 5A). Residues (amino acids 198–208) adjacent to the activation loop also showed higher fluctuations, although there were no remarkable differences among the top complexes.

The compactness of a protein, often estimated by the Rg from MD simulation trajectories, infers protein stability in a biological system. A lower Rg value suggests more compactness and increased protein stability. For the top ligands, the mean Rg values ranged from the minimum of 2.034 nm for Nitiducarpin to the maximum of 2.053 nm for Paratocarpin J (Fig. 5B). Notably, the standard deviation was also the lowest (±0.015 nm) for RIPK1 complexed with Nitiducarpin, denoting a more stable binding of this ligand.

The SASA of proteins is a robust estimator of protein folding and stability and has been employed in predicting protein-ligand binding and pathogenicity of protein mutations (Núñez, Venhorst & Kruse, 2010; Wei, Brooks & Frank, 2017; Savojardo et al., 2021). In MD simulations of protein-ligand complexes, a significant increase in SASA can be translated into a departure of the ligand from the binding pocket. We calculated protein SASA values from the trajectories using the gmx sasa tool (Fig. 6). Interestingly, RIPK1 complexed with Nitiducarpin exhibited the most stable SASA values (Mean: 159.475 nm², SD: ±2.395 nm²) while the variation was maximum for Paratocarpin J (Fig. 6A). Nevertheless, per residue SASA values did not display any remarkable changes among the complexes (Fig. 6B).

We estimated the binding free energies of the top ligands using the MMPBSA method (Fig. 7). The gas phase energy (GGAS) was contributed mainly by van der Waals interactions (VDWAALS) and electrostatic interactions (EEL). The reference ligand L8D displayed the strongest VDWAALS and EEL. However, the differences among the top three ligands were not substantial (Fig. 7). All ligands showed small and indistinguishable ENPOLAR values. Variations in GSOLV values were influenced by EPB values, with Nitiducarpin having the lowest GSOLV and L8D exhibiting the highest (Fig. 7). Total enthalpies were −28.857, −24.624, −29.699, and −34.741 kcal/mol for Nitiducarpin, Pinocembrin 7-O-benzoate, Paratocarpin J, and L8D, respectively (Fig. 7). To decipher the roles of binding residues, we also analyzed per residue total energy decomposition using the MMPBSA method (Table 2). Residues within 4 Å of the ligand were selected for analysis. Val31, Ile43, Met92, Ala155, Leu157, and Leu159 were found to play roles in ligand binding. Notably, Met92 and Leu157 made the most significant contributions (Table 2).

Combining the above RMSD, RMSF, and SASA analyses for the MD simulation trajectories of the top protein-ligand complexes (Figs. 3–6), Nitiducarpin demonstrated the most favorable binding to RIPK1. To further reveal the stability of the RIPK1-Nitiducarpin complex, we generated 2D protein-ligand interaction diagrams from snapshots retrieved from the simulation trajectory at 100 ns intervals (Fig. 8). We observed a clear predominance of nonpolar interactions of Nitiducarpin (Fig. 8), as expected from the target hydrophobic binding cavity in our molecular docking step. A stable ligand binding was governed primarily through VDWAALS and Alkyl/Pi-Alkyl interactions (Fig. 8).

Similar to canonical kinases, the kinase activity of RIPK1 largely depends on the catalytic triad residues Lys45, Glu63, and Asp156, P-loop residues 24–31, and the catalytic loop residues 136–143 (Xie et al., 2013). We observed that Nitiducarpin formed stable Alkyl/Pi-Alkyl interaction with Lys45 and VDWAALS with Asp156, the two residues of the kinase catalytic triads (Fig. 8). An interaction with the P-loop residue Val31 was also evident at all time points. Other significant interacting residues include Met67, Val76, Leu78, Leu90, Met92, Leu157, and Phe162 (Fig. 8). Overall, Nitiducarpin is predicted to interact with many key residues governing a stable binding of the reference ligand L8D to RIPK1 (Fig. 8) (Xie et al., 2013).

To confirm that Nitiducarpin could indeed form stable interactions with RIPK1, we calculated the distances of the two key residues, Met92 and Leu157, from the ligand (Fig. 9). The analysis revealed that Nitiducarpin stayed within a mean distance of 4.060 Å and 4.062 Å from Met92, and Leu157, respectively. We observed brief departures from Met92 on several occasions. Nevertheless, stable interactions were maintained with Leu157 throughout the simulation period (Fig. 9).

In summary, we screened 4,858 flavonoids using molecular docking simulations. We shortlisted 20 ligands based on binding affinities, which were then subjected to ADME screening, returning three orally bioavailable and blood-brain barrier penetrant compounds: Nitiducarpin, Pinocembrin 7-O-benzoate, and Paratocarpin J. We validated binding interactions using replicated molecular dynamics simulations, compared the results with the strong RIPK1 inhibitor L8D, and identified Nitiducarpin as the most potent ligand for RIPK1. To be noted, Nitiducarpin (PubChem CID: 630656) is an isoflavonoid of the pterocarpan subclass with the chemical formula C₂₆H₂₆O₅ and a molecular weight of 418.48 g/mol.