Computational Identification of Potential Anti-Inflammatory Natural Compounds Targeting the p38 Mitogen-Activated Protein Kinase (MAPK): Implications for COVID-19-Induced Cytokine Storm

Severely ill coronavirus disease 2019 (COVID-19) patients show elevated concentrations of pro-inflammatory cytokines, a situation commonly known as a cytokine storm. The p38 MAPK receptor is considered a plausible therapeutic target because of its involvement in the platelet activation processes leading to inflammation. This study aimed to identify potential natural product-derived inhibitory molecules against the p38α MAPK receptor to mitigate the eliciting of pro-inflammatory cytokines using computational techniques. The 3D X-ray structure of the receptor with PDB ID 3ZS5 was energy minimized using GROMACS and used for molecular docking via AutoDock Vina. The molecular docking was validated with an acceptable area under the curve (AUC) of 0.704, which was computed from the receiver operating characteristic (ROC) curve. A compendium of 38,271 natural products originating from Africa and China together with eleven known p38 MAPK inhibitors were screened against the receptor. Four potential lead compounds ZINC1691180, ZINC5519433, ZINC4520996 and ZINC5733756 were identified. The compounds formed strong intermolecular bonds with critical residues Val38, Ala51, Lys53, Thr106, Leu108, Met109 and Phe169. Additionally, they exhibited appreciably low binding energies which were corroborated via molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) calculations. The compounds were also predicted to have plausible pharmacological profiles with insignificant toxicity. The molecules were also predicted to be anti-inflammatory, kinase inhibitors, antiviral, platelet aggregation inhibitors, and immunosuppressive, with probable activity (Pa) greater than probable inactivity (Pi). ZINC5733756 is structurally similar to estradiol with a Tanimoto coefficient value of 0.73, which exhibits anti-inflammatory activity by targeting the activation of Nrf2. Similarly, ZINC1691180 has been reported to elicit anti-inflammatory activity in vitro. The compounds may serve as scaffolds for the design of potential biotherapeutic molecules against the cytokine storm associated with COVID-19.


Introduction
Coronavirus disease 2019  is caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) [1,2]. It is a positive-sense single-stranded RNA virus [3] Atom (AA) force field was used to generate the protein topology and position restrain files [50,51]. Periodic Boundary Conditions (PBC) were applied to the structure with the protein centered 1 nm from the edge of a cubic box to monitor the movement of all particles and avoid edge effects on the surface atoms [52]. The system was solvated with SPC water [53,54] and neutralized with six Na atoms. The structure was energy minimized using the steepest descent algorithm at 50,000 steps.

Validation of Docking Protocol
The easyROC [58] was used to generate the receiver operating characteristic (ROC) curve. The ROC curve was used to validate the ability of Autodock Vina to discriminate between active and decoy compounds. The SMILES of the eleven inhibitors served as inputs for the generation of 50 decoys each via Directory of Useful Decoys (DUD-E) [59]. Each decoy has similar physicochemical properties to a known inhibitor but is chemically distinct [59]. A total of 561 compounds comprising 11 actives and 550 decoys were screened against the p38 MAPK structure. The 11 actives comprised the known inhibitors of p38 MAPK.

Virtual Screening of Ligands
AutoDock Vina [60] was used to screen the library against the p38α MAPK protein structure. The pre-filtered library was imported into OpenBabel [61] and minimized using the Universal Force Field (Uff) for 200 steps and optimized using the conjugate gradient. A grid box of dimensions (78.35, 47.84, 49.86) Å and center (9.97, 30.65, 20.29) Å was used for docking with a default exhaustiveness of eight.

Pharmacological Profiling
The pharmacokinetic and physicochemical profiles of the compounds were predicted via SwissADME [62] with the Simplified Molecular Input Line-Entry System (SMILES) of the compounds as inputs.

Elucidation of the Protein-Ligand Interactions
LigPlot+ [63] was used to generate the 2D protein-ligand interactions with default settings. The best poses of the hits were saved in ".pdb" file formats and then visualized using PyMOL.

Prediction of Biological Activities of Hit Compounds
The biological activities of the hits were predicted using the Bayesian-based Prediction of Activity Spectra for Substances (PASS) [64]. The SMILES files of the compounds were used to perform structural similarity searches for antiviral and anti-inflammatory compounds with a DrugBank [65] similarity threshold of 0.7.

Molecular Dynamics Simulation of Protein-Ligand Complexes
MD simulations of the protein-ligand complexes were performed using GROMACS 2020 version [49]. The 3ZS5 protein topologies were generated using the GROMOS96 43a1 force field [66] in GROMACS, whereas the ligand files were generated via the PRODRG [67]. Each complex was solvated with water molecules in a cubic box of size 5.0 nm and neutralized with six Na ions. Energy minimization of each complex was conducted for 50,000 steps using the steepest descent algorithm. The ligands were restrained before the canonical ensemble, and then the isothermal-isobaric ensemble. Equilibration of each complex was performed for 100 ps apiece and final MD simulation was conducted for 100 ns with time steps of 2 fs under PME. The free binding energies were computed using G_mmpbsa [68]. The binding free energy contribution per residue was calculated using MM-PBSA and the output plots generated with R.

Protein Structure Retrieval and Analysis
The PDB database contains several solved p38α MAPK structures. The p38α is the most expressed isoform of the p38 MAPK [69]. The selected protein structure with PDB ID: 3ZS5 was solved using x-ray diffraction at a high resolution of 1.6 Å. This ensured the detailed representation of the water-mediated hydrogen-bonding network [45]. 3ZS5 has also being used in a previous computational study [70] and it is bound to the inhibitor SB 203580 [45]. The structure is a monomer consisting of 362 residues and is made up of two domains comprising N-and C-terminal domains. The N-terminal domain is composed of largely β-sheets, while the C-terminal domain is essentially helical [71]. The catalytic ATP binding site is found between the two domains ( Figure 1). This site is characterized by a peptide flip, resulting in the presentation of two hydrogen-bond donors in the ATP pocket, a clear distinction from the acceptor-donor-acceptor pattern normally present in the ATP pocket of other protein kinases [45]. One large binding cavity with a surface area of 1040.309 Å 2 and volume 1365.552 Å 3 was predicted via CASTp 3.0 [72].

Prediction of Biological Activities of Hit Compounds
The biological activities of the hits were predicted using the Bayesian-based Prediction of Activity Spectra for Substances (PASS) [64]. The SMILES files of the compounds were used to perform structural similarity searches for antiviral and anti-inflammatory compounds with a DrugBank [65] similarity threshold of 0.7.

Molecular Dynamics Simulation of Protein-Ligand Complexes
MD simulations of the protein-ligand complexes were performed using GROMACS 2020 version [49]. The 3ZS5 protein topologies were generated using the GROMOS96 43a1 force field [66] in GROMACS, whereas the ligand files were generated via the PRODRG [67]. Each complex was solvated with water molecules in a cubic box of size 5.0 nm and neutralized with six Na ions. Energy minimization of each complex was conducted for 50,000 steps using the steepest descent algorithm. The ligands were restrained before the canonical ensemble, and then the isothermal-isobaric ensemble. Equilibration of each complex was performed for 100 ps apiece and final MD simulation was conducted for 100 ns with time steps of 2 fs under PME. The free binding energies were computed using G_mmpbsa [68]. The binding free energy contribution per residue was calculated using MM-PBSA and the output plots generated with R.

Protein Structure Retrieval and Analysis
The PDB database contains several solved p38α MAPK structures. The p38α is the most expressed isoform of the p38 MAPK [69]. The selected protein structure with PDB ID: 3ZS5 was solved using x-ray diffraction at a high resolution of 1.6 Å. This ensured the detailed representation of the water-mediated hydrogen-bonding network [45]. 3ZS5 has also being used in a previous computational study [70] and it is bound to the inhibitor SB 203580 [45]. The structure is a monomer consisting of 362 residues and is made up of two domains comprising N-and C-terminal domains. The N-terminal domain is composed of largely β-sheets, while the C-terminal domain is essentially helical [71]. The catalytic ATP binding site is found between the two domains ( Figure 1). This site is characterized by a peptide flip, resulting in the presentation of two hydrogen-bond donors in the ATP pocket, a clear distinction from the acceptor-donor-acceptor pattern normally present in the ATP pocket of other protein kinases [45]. One large binding cavity with a surface area of 1040.309 Å 2 and volume 1365.552 Å 3 was predicted via CASTp 3.0 [72].

Docking Protocol Validation
The docking protocol of AutoDock Vina [60] was validated using the ROC curve. Eleven inhibitors with their respective decoys were screened against the p38 MAPK structure to generate a ROC curve. The ROC curve evaluates the performance of AutoDock Vina to classify inhibitors from decoys [73]. An area under the curve (AUC) value of 1 is considered perfect, while an AUC below 0.5 signifies poor classification [74]. Additionally, an AUC value of 0.7 to 0.8 is considered acceptable, 0.8 to 0.9 is very good, and above 0.9 is excellent [75,76]. The obtained AUC computed from the ROC curve was 0.704 ( Figure 2), which is considered acceptable [75]. This indicates that AutoDock Vina reasonably distinguishes between active and inactive compounds of p38 MAPK. AutoDock Vina is more accurate than AutoDock 4 in selecting anti-SARS-CoV-2 M pro compounds [77]. Moreover, AutoDock Vina was employed successfully to screen compounds against p38 MAPK, which were found to inhibit p38 MAPK in vitro [70].
binding pocket is colored green. The structure is a monomer with a binding cavity of 1365.552 Å 3 .

Docking Protocol Validation
The docking protocol of AutoDock Vina [60] was validated using the ROC curve. Eleven inhibitors with their respective decoys were screened against the p38 MAPK structure to generate a ROC curve. The ROC curve evaluates the performance of AutoDock Vina to classify inhibitors from decoys [73]. An area under the curve (AUC) value of 1 is considered perfect, while an AUC below 0.5 signifies poor classification [74]. Additionally, an AUC value of 0.7 to 0.8 is considered acceptable, 0.8 to 0.9 is very good, and above 0.9 is excellent [75,76]. The obtained AUC computed from the ROC curve was 0.704 (Figure 2), which is considered acceptable [75]. This indicates that AutoDock Vina reasonably distinguishes between active and inactive compounds of p38 MAPK. AutoDock Vina is more accurate than AutoDock 4 in selecting anti-SARS-CoV-2 M pro compounds [77]. Moreover, AutoDock Vina was employed successfully to screen compounds against p38 MAPK, which were found to inhibit p38 MAPK in vitro [70].

Pre-Filtering of Library and Molecular Docking Studies
A total of 7282 out of the 38,271 natural products were obtained after pre-filtering. These compounds had molecular weights between 250 g/mol and 350 g/mol, and were predicted to be non-toxic. The compounds were screened against the ATP-binding pocket of p38 MAPK and were ranked based on their binding energies. AutoDock Vina combines both empirical scoring functions and knowledge-based potentials to compute the binding energies of compounds [60]. ZINC38321631, a compound annotated in the TCM database, Figure 2. Evaluating the performance of the virtual screening via ROC curve. Binding energies from the screening of inhibitors and decoys against the p38 MAPK receptor were used to generate the curve. An AUC of 0.704 was obtained from the ROC curve.

Pre-Filtering of Library and Molecular Docking Studies
A total of 7282 out of the 38,271 natural products were obtained after pre-filtering. These compounds had molecular weights between 250 g/mol and 350 g/mol, and were predicted to be non-toxic. The compounds were screened against the ATP-binding pocket of p38 MAPK and were ranked based on their binding energies. AutoDock Vina combines both empirical scoring functions and knowledge-based potentials to compute the binding energies of compounds [60]. ZINC38321631, a compound annotated in the TCM database, had the lowest binding energy of −12.4 kcal/mol, whereas ZINC5734567 had the highest binding energy of −2.5 kcal/mol. Among the known inhibitors, SB 202190 had the least binding energy of −11.0 kcal/mol, followed by SB 203580 with −10.9 kcal/mol. Losmapimod that is currently undergoing clinical trial had binding energy of −9.3 kcal/mol, whilst that of Dilmapimod was −7.9 kcal/mol. In this study, a threshold of −11.0 kcal/mol was used for selecting hits. Selected compounds had lower binding energies relative to the known inhibitors. A previous study used a threshold of −10.0 kcal/mol for p38 MAPK [78]. Another docking study on the potential inhibitors of p38 MAPK identified a promising compound with a binding energy of −11.1 kcal/mol [79]. A total of 42 compounds and two inhibitors (SB 202190 and SB 203580) were selected for downstream analysis.

Pharmacological Profiling of Hit Compounds
Pharmacological profiling studies are important aspects of the drug development pipeline. The human intestinal absorption, permeability glycoprotein (P-gp) binding, and cytochrome P450 3A4 (CYP450) inhibition were considered [80]. The human intestinal absorption measures the likelihood of absorption for orally administered drugs into the bloodstream. Two of the hits, namely ZINC85550217 and ZINC85543655 were predicted to have low GI absorption. Additionally, 16 out of the 42 compounds including ZINC38321631 were predicted as P-gp substrates with likely decrease in drug bioavailability [81,82]. CYP3A4 is responsible for the metabolism of about 50% of all drugs [83]. The two standard inhibitors (SB 202190 and SB 203580) were predicted as CYP3A4 inhibitors. A total of 17 hits were predicted to be CYP3A4 inhibitors suggesting they may interfere with the metabolism of other drugs. After eliminating compounds predicted to have poor pharmacological profiles, 18 compounds were selected as promising hits (Supplementary Table S1) with none of them violating Lipinski's rule of five for evaluating drug-likeness.

Biological Activity Prediction
The prediction of activity spectra for substances (PASS) [64] was employed to predict relevant biological activities of the hits. PASS utilizes structural-activity relationship to predict the potential biological activities and mechanisms. When probable activity (Pa) is greater than the probable inactivity (Pi) for a particular compound activity and Pa > 0.3, it is essential to test in vitro the predicted activity [57,88]. The biological activities that were considered included anti-inflammation, kinase inhibition, antiviral, platelet aggregation inhibition and immunosuppression. PASS predictions are usually corroborated experimentally [89,90]. ZINC5519433 was predicted as a platelet aggregation inhibitor (Pa = 0.509 and Pi = 0.008), and antiviral (rhinovirus) with Pa = 0.470 and Pi = 0.037. Platelet activation accelerates inflammation and can lead to disseminated intravascular coagulation, a condition observed in critically ill COVID-19 patients [91,92]. ZINC5733756 was predicted as antiviral (influenza)

The Rationale for the Selection of Compounds
A lead compound is a molecule that is likely to be of therapeutic relevance and may be modified to improve potency, selectivity, and pharmacological profiles [94]. After screening a library of 7282 pre-filtered compounds against the p38 MAPK structure, 42 compounds with binding energies ≤ −11.0 kcal/mol were selected for pharmacological profiling. A total of 18 compounds were selected from this lot as promising hits due to predicted good pharmacological profiles. The protein-ligand interactions of the promising hits were evaluated to identify residues critical for binding. Further biological activity predictions were performed to identify compounds with relevant biological activities and mechanisms of action. The multistage techniques informed the selection of five compounds, namely ZINC5519433, ZINC5733756, ZINC95486106, ZINC1691180, and ZINC4520996 ( Table 2). The selected compounds were predicted as having plausible pharmacological profiles including high gastrointestinal absorption, non-Pgp substrates, and non-CYP3A4 inhibitors. Assessment of their predicted biological activities showed all five to have anti-inflammatory and antiviral propensities, with Pa > Pi and Pa > 0.3.

Molecular Dynamics Simulation of Selected Compounds
A 100 ns MD simulation was performed on eight structures consisting of the unbound protein, the complexes of the five selected compounds, and the two known p38 MAPK inhibitors SB 203580 and SB 202190. This was aimed at probing the structural stability and conformational changes in physiological conditions [108]. The parameters evaluated were the root mean square deviation (RMSD), the radius of gyration (Rg), and the conformational changes with time.
The RMSD is a metric used to assess the stability of a protein structure [109]. An RMSD plot of the p38 MAPK-ligand complexes that tend to converge is indicative of a stable and well-equilibrated system [110]. Considering the p38 MAPK-ZINC5519433 RMSD plot, a steady rise to 0.38 nm was observed during the 1 ns of simulation. It then descended and maintained an average RMSD of 0.3 nm (Figure 4a). Moreover, the RMSD of the p38 MAPK-ZINC5733756 complex gradually rose to 0.42 nm around 15 ns and quickly dropped to an average RMSD of 0.3 nm around 17 ns (Figure 4a). It maintained the average RMSD of 0.3 nm till about 62 ns, where a sharp rise to an average of 0.45nm was observed till the end of the 100 ns simulation period (Figure 4a). The p38 MAPK-ZINC4520996 complex experienced a rise in RMSD to 0.45 nm around 13 ns and maintained stability till the end of the simulation time (Figure 4a). An initial rise to 0.40 nm during the first 6 ns was observed for the p38 MAPK-ZINC1691180 complex (Figure 4a).
ZINC5519433 had binding energy of −11.6 kcal/mol (Table 1). It formed hydrophobic interactions with residues including Lys53, Thr106 and Phe169, which have been reported in other studies to be critical in ligand binding (Table 1). It is 71% structurally similar to macelignan via DrugBank [95]; a natural compound that attenuates the activation of p38 MAPK and has exhibited anti-inflammatory activity [96]. Previous experimental studies showed ZINC5519433 as a kinase inhibitor and selective to c-Jun N-terminal kinases (JNK) [97]. This corroborates our studies on the potential lead-likeness of ZINC5519433. A further lead optimization may likely increase its selectivity towards p38 MAPK [98].
ZINC1691180 and ZINC4520996 were both predicted to have a binding energy of −11.6 kcal/mol (Table 1). Both formed hydrogen bonding with Tyr35, with varying lengths and hydrophobic contacts with some critical residues (Table 1). ZINC1691180 has been reported to evoke anti-inflammatory properties via regulating the production of inflammatory factor TNFα in mice [93]. TNFα mediates the activation of p38 MAPK [99] and obstruction of this pathway with ZINC1691180 could be useful in severe COVID-19 cases. ZINC5733756 had binding energy of −11.1 kcal/mol due to strong intermolecular interaction with residues including Thr106 and Phe169 (Table 1). Moreover, ZINC5733756 is 73% structurally similar to estradiol, a hormone that exerts anti-inflammatory activities by targeting the activation of Nrf2 [100].
Furthermore, ZINC95486106 formed interactions with critical residues namely Lys53, Leu108, and Phe169 contributing to its binding energy of −12.1 kcal/mol. It was also predicted to be GPCR and Beta-adrenergic receptor kinase inhibitors via PASS with Pa > 0.8.
The inhibitory constants (Ki) of the selected compounds and the known inhibitors were calculated (Supplementary Table S2) using a previously described method [57,101]. All the selected compounds had lower predicted Ki values than all the known inhibitors used herein (Supplementary Table S2), indicating their inhibitory potential [102]. For the selected compounds, ZINC95486106 was predicted to have the lowest  Table S2) [103][104][105][106][107]. ARRY-797 with the lowest IC 50 value of 4.5 nM had binding energy of −9.5 kcal/mol, which is higher than that of SB 20358 (−10.9 kcal/mol) with the highest IC 50 value 41.2 nM.

Molecular Dynamics Simulation of Selected Compounds
A 100 ns MD simulation was performed on eight structures consisting of the unbound protein, the complexes of the five selected compounds, and the two known p38 MAPK inhibitors SB 203580 and SB 202190. This was aimed at probing the structural stability and conformational changes in physiological conditions [108]. The parameters evaluated were the root mean square deviation (RMSD), the radius of gyration (Rg), and the conformational changes with time.
The RMSD is a metric used to assess the stability of a protein structure [109]. An RMSD plot of the p38 MAPK-ligand complexes that tend to converge is indicative of a stable and well-equilibrated system [110]. Considering the p38 MAPK-ZINC5519433 RMSD plot, a steady rise to 0.38 nm was observed during the 1 ns of simulation. It then descended and maintained an average RMSD of 0.3 nm (Figure 4a). Moreover, the RMSD of the p38 MAPK-ZINC5733756 complex gradually rose to 0.42 nm around 15 ns and quickly dropped to an average RMSD of 0.3 nm around 17 ns (Figure 4a). It maintained the average RMSD of 0.3 nm till about 62 ns, where a sharp rise to an average of 0.45nm was observed till the end of the 100 ns simulation period (Figure 4a). The p38 MAPK-ZINC4520996 complex experienced a rise in RMSD to 0.45 nm around 13 ns and maintained stability till the end of the simulation time (Figure 4a). An initial rise to 0.40 nm during the first 6 ns was observed for the p38 MAPK-ZINC1691180 complex (Figure 4a). The RMSD dropped to about 0.32 nm around 20 ns and was maintained till the 100 ns period (Figure 4a). Additionally, the RMSD plot for p38 MAPK-ZINC95486106 revealed stability after 15 ns, with an average RMSD of 0.33 nm (Figure 4a). Considering the known inhibitors, average RMSDs of 0.35 nm and 0.38 nm were observed for the SB 203580 and SB 202190-p38 MAPK complexes, respectively. The p38 MAPK complexes of ZINC5519433, ZINC95486106, ZINC1691180, ZINC4520996, SB 203580 and SB 202190 had lower average RMSDs relative to the unbound protein (Figure 4a). This may suggest stability induced by ligand binding.
In a pharmacophore-guided study to identify p38 MAPK inhibitors, a 10 ns MD simulation was performed on the most potent molecule (compound 48) and the protein complex [111]. It was observed that the RMSD ranged between 0.5 to 2.5 Å (0.05 to 0.25 nm) [111]. Other MD simulation studies on p38 MAPK-ligand complexes reported RMSD values up to about 0.6 nm [112,113] consistent with the results presented herein. The Rg was measured to evaluate the compactness of the protein systems during the MD simulation [114]. The Rg of p38-MAPK-ZINC5519433 fell from 2.25 nm to 2.20 nm during the first 5 ns and it further rose and peaked at 2.26 nm around 9 ns. It then fluc-tuated over the remaining time with an average Rg of 2.20 nm (Figure 4b). For the p38 MAPK-ZINC5733756 complex, an initial Rg of 2.2 nm was obtained and it then fluctuated around 2.15 nm until 60 ns, where it dropped to 2.02 nm. The p38 MAPK-ZINC4520996, p38 MAPK-ZINC1691180, p38 MAPK-ZINC95486106, p38 MAPK-SB 203580, and p38 MAPK-SB 202190 complexes were observed to be stable over the simulation time with average Rg of 2.12 nm, 2.13 nm, 2.13 nm, 2.14 nm, and 2.15 nm, respectively. These were close to the average Rg of 2.12 nm exhibited by the unbound protein. A stably folded protein maintains a reasonable steady radius of gyration throughout the simulation [114].   The distance between the centre of mass of a protein and an inhibitor influences the molecular interactions, which in turn play a crucial role in the structural stability of the protein-ligand complexes [115]. Snapshots at 25 ns intervals (time step = 0, 25, 50, 75, and 100 ns) were generated to elucidate the position and binding modes of the ligands during the 100 ns MD simulation period ( Figure 5 and Supplementary Figure S1). In the snapshot analysis, all the ligands were observed to stably bind in the ATP binding pocket of the p38 MAPK protein throughout the 100 ns simulation time. Moreover, the number of hydrogen bonds over the simulation time with a distance cut-off of 0.35 nm was also generated ( Figure 6). The 2D protein-ligand interaction profiles of the snapshot frames of the protein-ligand complexes at time intervals of 25 ns were also generated to study the interactions during the 100 ns simulation period ( Supplementary Figures S2 and S3). For the p38 MAPK-ZINC1691180 complex, the hydrogen bond of length 3.17 Å with Tyr35 before MD (Supplementary Figure S2a) was lost at 25 ns and was not formed again throughout the simulation time. However, at 75 ns, ZINC1691180 formed a new hydrogen bond with Ser32 of bond length 2.96 Å, which was lost at 100 ns ( Figure S3a). For the p38 MAPK-ZINC4520996 complex, the hydrogen bonding with Tyr35 pre-MD was not present at 25 ns ( Supplementary Figures S2b and S3b). However, at 50 ns, ZINC4520996 formed a new hydrogen bond with Asp112 of bond length 3.01 Å (Supplementary Figure S3b).
The hydrogen bond with Asp112 was lost at 75 ns and re-formed at 100 ns with a bond length of 3.2 Å ( Figure S3b). For the p38 MAPK-ZINC5519433 complex, the ligand was observed to form a hydrogen bond at 25 ns with Gly170 of bond length of 3.15 Å (Supplementary Figure S3d). However, at 50, 75, and 100 ns, there were no hydrogen bonds observed. For the protein-ZINC95486106 complex, the hydrogen bond with Met109 (bond length of 3.03 Å) was lost throughout the simulation time ( Supplementary Figures S2c and S3c). However, at 100 ns, two hydrogen bonds were observed with Ser32 of a bond length of 2.73 Å and Lys15 with a bond length of 3.06 Å (Supplementary Figure S3c). For the p38 MAPK-ZINC5733756 complex, ZINC5733756 did not form any hydrogen bond with the protein before MD (Supplementary Figure S2e). At 25 ns, ZINC5733756 was observed to form one hydrogen bond with Asp112 of length 2.86 Å. At 50 ns, it formed two hydrogen bonds with the p38 MAPK protein composed of one with Asp112 (bond length of 2.71 Å) and the other with Asn115 (bond length of 2.70 Å). Hydrogen bonds with Asp112 and Asn115 of lengths 2.65 and 3.14 Å were observed at 75 ns, respectively (Supplementary Figure S3e). At 100 ns, ZINC5733756 maintained only one hydrogen bond of length 2.70 Å with Asp112 (Supplementary Figure S3e). The multiple hydrogen bonds formed indicate the strong interactions between ZINC5733756 and the p38 MAPK protein throughout the simulation period and could influence the activity of the ligand [116].

Evaluation of Selected Compounds via MM-PBSA Calculations
MM-PBSA technique was used to estimate the free binding energies of the seven complexes. This technique addresses limitations with current scoring functions [117,118]. The binding free energy is a more reliable method of evaluating docking studies [77,117,119]. The other parameters computed included the van der Waal, electrostatic, polar, and non-polar solvation energies (Table 3) [120,121]. The selected molecules comprising ZINC5519433, ZINC95486106, ZINC5733756, ZINC1691180 and ZINC4520996 had predicted average free binding energies of −185.122, 30.620, −61.726, −146.008, and −151.561 kJ/mol, respectively. SB 202190 and SB 203580 were predicted to have plausible free binding energies of −166.369 kJ/mol and −236.175 kJ/mol, respectively (Table 3). Since ZINC95486106 was the only selected compound with positive binding free energy, it was eliminated as a potential lead. Optimization of the ligand structure in future studies can improve its affinity to the p38 MAPK protein. The poor binding free energy could largely be due to the high electrostatic energy relative to the other compounds (Table 4) and hence may limit its lead likeness [57]. In all seven compounds assessed, the van der Waal forces contributed favorably to the free binding energies [29,120,121]. However, the polar solvation energy contributed large positive energies to binding in all the complexes [122].   and (e) ZINC5733756-p38 MAPK complexes. The ligands are represented as spheres and the protein as cartoons. All the ligands were observed to bind stably in the ATP binding pocket. Figure 6. The total number of hydrogen bonds formed between the selected compounds and the protein structure. ZINC1691180, ZINC4520996, ZINC95486106, ZINC5519433, ZINC5733756, SB 202190, and SB 203580-p38 MAPK complexes are represented as black, red, green, blue, yellow, brown and ash, respectively.

Evaluation of Selected Compounds via MM-PBSA Calculations
MM-PBSA technique was used to estimate the free binding energies of the seven complexes. This technique addresses limitations with current scoring functions [117,118]. The binding free energy is a more reliable method of evaluating docking studies [77,117,119]. The other parameters computed included the van der Waal, electrostatic, polar, and non-polar solvation energies (Table 3) [120,121]. The selected molecules comprising ZINC5519433, ZINC95486106, ZINC5733756, ZINC1691180 and ZINC4520996 had predicted average free binding energies of −185.122, 30.620, −61.726, −146.008, and −151.561 kJ/mol, respectively. SB 202190 and SB 203580 were predicted to have plausible free binding energies of −166.369 kJ/mol and −236.175 kJ/mol, respectively (Table 3). Since ZINC95486106 was the only selected compound with positive binding free energy, it was eliminated as a potential lead. Optimization of the ligand structure in future studies can improve its affinity to the p38 MAPK protein. The poor binding free energy could largely be due to the high electrostatic energy relative to the other compounds (Table 4) and hence may limit its lead likeness [57]. In all seven compounds assessed, the van der Waal forces contributed favorably to the free binding energies [29,120,121]. However, the polar solvation energy contributed large positive energies to binding in all the complexes [122].
Additionally, a per-residue decomposition of the binding energy was performed to gain useful insight into important interactions of key residues (Figure 7 and Supplementary Figure S4). Residues contributing binding free energies greater than 5 kJ/mol or less than −5 kJ/mol are worthy of consideration as critical for binding [57]. For the p38 MAPK-ZINC5519433 complex, critical residues Phe169 and Thr106 contributed favorable energies of −9.2470 kJ/mol and −6.1224 kJ/mol, to ligand binding, respectively (Table 4 and Figure 7). Asp168 was found to contribute unfavorable energy of 7.1894 kJ/mol (Figure 7).  Additionally, a per-residue decomposition of the binding energy was performed to gain useful insight into important interactions of key residues (Figure 7 and Supplementary Figure S4). Residues contributing binding free energies greater than 5 kJ/mol or less than −5 kJ/mol are worthy of consideration as critical for binding [57]. For the p38 MAPK-ZINC5519433 complex, critical residues Phe169 and Thr106 contributed favorable energies of −9.2470 kJ/mol and −6.1224 kJ/mol, to ligand binding, respectively (Table 4 and Figure 7). Asp168 was found to contribute unfavorable energy of 7.1894 kJ/mol (Figure 7). Moreover, only Phe169 contributed binding free energies of ±5.0 kJ/mol to ZINC5733756 binding (Table 4 and Supplementary Figure S4d). From the energy decomposition plot, Tyr30, Tyr38, Leu167, and Phe169 were predicted as critical residues that contributed favorably to ZINC1691180 binding (Table 4 and Supplementary Figure S4a). However, Lys53 and Asp112 contributed adversely towards ZINC1691180 binding with energies of 8.4199 and 9.7176 kJ/mol, respectively. For ZINC4520996 complex, Lys53, Leu167, Asp168, and Phe169 contributed energies of ±5.0 kJ/mol to binding (Table 4 and Supplementary Figure S4b). Altogether, Lys53 was found to contribute high repulsive positive energies in the p38 MAPK-ZINC5519433, ZINC1691180, ZINC4520996, SB 203580 and SB 202190 complexes (Table 4), unfavorable for ligand binding in the ATP binding pocket. Phe169 contributed favorable energies <−5.0 kJ/mol in all complexes examined except ZINC95486106, corroborating its critical role in ligand binding as suggested in previous studies [87,110].
Therefore, this study sought to identify potential p38 MAPK inhibitors which could be explored as useful in ameliorating the cytokine storms in severely ill COVID-19 patients. Pharmacoinformatics-based approaches including molecular docking and dynamics studies predicted four potential lead compounds with good binding affinity to the p38 MAPK protein and negligible toxicity. The predicted biological activities of the selected compounds showed that they possessed anti-inflammatory, kinase inhibition, antiviral, platelet aggregation inhibition, and immunosuppression activities. Dexamethasone, which is an anti-inflammatory compound is used as a treatment option in COVID-19 patients [26]. Moreover, Losmapimod, which is a p38 MAPK inhibitor, is currently undergoing clinical trials for use as a treatment option in severely ill COVID-19 patients (https://clinicaltrials.gov/ct2/show/NCT04511819 (accessed on 8 April 2021)). Herein, ZINC5519433 was predicted as a platelet aggregation inhibitor and anti-rhinovirus. ZINC5733756 was also predicted as an anti-inflammatory, anti-influenza, and anti-rhinovirus. Moreover, ZINC1691180 and ZINC4520996 were predicted to be anti-inflammatory and antiviral. Interestingly, ZINC1691180 has previously been shown to elicit an anti-inflammatory response in mice [93]. The effects of the predicted compounds as potential inhibitors of p38 MAPK to possibly reduce the cytokine storm can be investigated experimentally. The plethora of predicted activities and mechanisms serve as clues that warrant further evaluation of these compounds. The molecules can be optimized as well as used for fragment-based de novo drug design of potential novel biotherapeutics. These predicted compounds can help fuel the pace of searching for effective antivirals with anti-inflammatory activity.

Conclusions
The study used multistage computational techniques to predict four natural productderived molecules comprising ZINC5733756, ZINC1691180, ZINC4520996 and ZINC5519433, which have the potential to inhibit the p38 MAPK. The compounds were selected based on high binding affinity and plausible binding mechanisms with the p38 MAPK protein structure obtained through molecular dynamics simulations including MM-PBSA. Key residues Tyr38, Ala51 and Phe169 were corroborated as critical for binding, which could guide the design and selection of future p38 MAPK inhibitors. Furthermore, these drug-like compounds were predicted as anti-inflammatory, kinase inhibitors, antiviral, platelet aggregation inhibitors, and immunosuppressive. Since the activation of the p38 MAPK leads to hyperinflammation in severely ill COVID-19 patients, the potential of these molecules to attenuate the cytokine storm can be explored. Further in vitro and in vivo evaluations of the suggested molecules could be undertaken to corroborate the predicted inhibitory activity. The study can serve as a clue for the design of novel p8α MAPK selective inhibitors which may have therapeutic implications in COVID-19-induced cytokine storm.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/biom11050653/s1, Table S1: Predicted pharmacokinetics of the 18 selected hits and two inhibitors. The profiles used were gastrointestinal (GI) absorption, cytochrome P450 (CYP) 3A4 and permeability glycoprotein (P-gp) substrates. Also included is the number of violations of Lipinski's rule of five, Table S2: Predicted inhibitory constants (Ki) and binding energy of the selected compounds and those of the known inhibitors. Also, included are the experimentally determined inhibition constants of known inhibitors., Figure S1: