Investigation of structure–activity relationship: In silico studies of [1, 2, 4]triazolo[4, 3-a]pyridine ureas as P38 kinase inhibitors

P38 kinases are the members of serine/threonine kinases family and play a vital role in the progression of inflammation. In the past two decades, numerous p38 kinase inhibitors have been reported, and few of them have failed in clinical trials. Recently, some of the p38 kinase inhibitors have entered in clinical trials for the treatment of Alzheimer’s disease. A potential opportunity exists for medicinal chemistry for the discovery of potent and safe p38 kinase inhibitors. In view of this challenging opportunity, the present manuscript is aimed towards development of a 3D quantitative structure–activity relationship (QSAR) model and the docking and dynamic simulation studies. A statistically robust 3D QSAR model was developed by employing 21 training set molecules which is attributed with appreciating cross-validation coefficient (q2) of 0.0.6269 and conventional correlation coefficient (r2) of 0.8783 respectively. The predicted correlation coefficient (r2pred) was found to be 0.8644 and standard error of 0.3331. The molecular docking analysis of all the p38 kinase inhibitors revealed that the analogs were well docked into the DFG (Asp-Phe-Gly motif) out pocket of p38 kinase and exhibited hydrogen bond interactions with Asp186 and Lys71. Extension of docking studies to the molecular dynamics simulation study informed that the ligand displayed the strong conformational stability within the active site of p38 kinase forming maximum two hydrogen bonds until 100 ns respectively.


Introduction
The p38 MAP kinase pathway is the third major mitogenactivated protein kinases signaling cascade and is considered as a key for inflammatory diseases because of its activation by pro inflammatory cytokines, stress, temperature and high osmotic stress [1], majorly, in inflammation mediated by tumor necrosis factor-α (TNFα) and interleukin-1β (IL-1β) [2]. There are four variants of p38 protein kinase, namely p38α, which is also called p38 kinase, a 38 kDa protein, p38β, p38γ and p38δ. The existence and distribution of these types of p38 kinases differ in the body, p38α, p38β are majorly present in the cells and tissues, p38γ exists only in skeletal muscle cell whereas p38δ is profound in glandular tissue. [3]. p38 MAP kinase pathway gets activated due to a number of extracellular stimulus or intracellular events and activates different cellular functions depending on the cell types. The activation of p38 isoforms can be specifically controlled through different regulators and co activated by various combinations of upstream and downstream regulators [4].
The p38α and β isoforms are expressed in inflammatory cells and are key modulators of TNF-α production, and inhibition of them in cells results in decreased levels of TNF-α expression. Also, administering inhibitors of p38α in animal models of inflammatory disease has proven that such inhibitors are effective in treating inflammatory diseases. Accordingly, the p38 kinase serve an important role in inflammatory processes mediated by IL-1 and TNF-α [5].
Inhibition of P38 kinase has played as an important strategy for anti-inflammatory drug discovery. Numerous p38 kinase inhibitors are under development like Pamapimod, Losmapimod, Doramapimod and VX-745. Doramapimod (BIRB-796) is one of the clinical candidates that is a pyrazolurea based compound which is relatively potent against P38α isoform. Doramapimod had entered clinical trials for the treatment of rheumatoid arthritis; however, it is withdrawn from further development due to CNS-related adverse effects in clinical trials [6,7]. Pamapimod (R-1503, Ro4402257) is a potent and selective inhibitor of p38 MAP kinase and was developed for use in the treatment of rheumatoid arthritis; however, it was discontinued later on because of lack of efficacy and refusal of continued treatment (which included noncompliance and withdrawal of consent) [8]. Losmapimod, a biphenyl amide, is considered to be a novel, oral p38 MAPK inhibitor and acts on both p38 MAPK α and β isoforms (Fig. 1). It is based on biphenyl amide or bicyclic nicotineamide scaffold that is currently undergoing clinical trials for rheumatoid arthritis (RA), asthma, cardiovascular diseases and neuropathic pain. [9]. VX-745 was under development for the potential treatment of rheumatoid arthritis (RA) and is a lead MAP kinase inhibitor and anti-inflammatory candidate, but this was also discontinued after phase II clinical trials because of adverse effects on the nervous system [10].
In spite of more than 20 candidates which are in clinical trials have been reported for p38 MAP kinase inhibitors [11], there is still not a single approved p38 MAP kinase inhibitor as drug in the market; hence, the floor of discovery and research for exploring the new chemical entities in this area looks fascinating and broad.

Data preparation
All the molecular modeling studies have been carried out on Windows 10 operating computer system and Ubuntu 18.04 standalone system. 29 N-pyrazolyl-N′-triazolo[4,3-a]pyridyl thiobenzyl ureas bearing p38 kinase inhibitory activity were collected from literature [20]. The SMILE notations of the compounds were converted into 3D structures on Open Babel.

Computational and alignment of database
The calculation of independent descriptors was done by using desktop-based molecular design tool Vlife MDS 4.6. Merck molecular force field (MMFF) model of force field and batch process minimization was applied by adding partial atomic charges to molecules, and the energy of all the analogs was minimized [21]. The crucial step of the quantitative structure-activity relationship (QSAR) study is the alignment of the analogs; hence, templatebased alignment method was utilized. The database alignment execution was done on the highly potent [1,2,3]triazolo [4,3-a]pyridine analogs of the data set [22]. The by default 3D alignment of the data set in X, Y and Z-axis with three-dimensional cubic lattice of grid space 2.0 Å was created based on the calculation of the 3D descriptors like steric, electrostatic and hydrophobic fields using QSAR module. The addition of Gasteiger-Marsili charges in the data set was achieved with a probe of carbon atom with SP 3 and + 1 charge. The cutoff values selected for electrostatic and steric field were 10 kcal/mol and 30 kcal/mol respectively, and the dielectric constant was 1.0 [23][24][25].

Partial least square analysis
The correlation between multiple independent variable and one dependent variable was established by the statistical method of partial least squares regression (PLS). PLS allows identification of specific structural trends of the compounds which can be related to the observed properties. High values of cross-validated coefficient (q 2 ) and low standard error value (SE) with sufficient number of compounds shall help us to get optimum measures for the development of QSAR model fit with a consistency. The cross-validation coefficient method was performed by removing one molecule from the set, and the model was developed with the remaining molecules. The model is considered to be statistically acceptable if the q 2 value > 0.5 and r 2 value > 0.8 and can be processed for further evaluation [26]. The calculation of q 2 and r 2 value is done with following formulae: RSS is calculated from the data on test set, and PRESS is calculated from data of calculated set.

Predicted correlation coefficient
The external validation was performed on the nine compounds that were grouped in the test set. The validation of the developed QSAR model is done by the calculation of values for the standard deviation (SD) and predictive residual sum of square (PRESS) [27]. SD is the measure of deviation of test set compounds biological activity and the mean of training set compounds biological activity. The press value is the sum of squared deviation of biological activity of the test set and calculated set.

Molecular docking
The interactions of the ligand with the active site amino acids and the required stabilization of the conformations of complex were identified by molecular docking simulation. All the desired molecular docking studies were performed with the help of Glide (Schrodinger Suit) software tool [28]. X-ray crystal structure of p38 kinase protein was collected from the RCSB PDB ID 2YIS. The Maestro software and Ligprep were used for drawing the structure of all the [1,2,4]triazolo[4, 3-a] pyridine compounds and ligand energy optimization respectively. The protein preparation was refined with protein preparation wizard, and ligand docking was performed with Glide (gridbased ligand docking with energetics) [29,30]. Kcal/mol is considered as unit of measurement for the efficiency of interaction of compounds with the macro-molecules and the visualization of respective log file thus generated during the simulation studies was done in a cross platform molecular

Homology modeling
PBD 2YIS when visualized and evaluated has revealed about thirty-three missing amino acids, and the missing amino acid residues have to be added prior to MD simulation.
The Modeller 9.22 command line tool was used to 'fill in' missing amino acid residues by treating the original structure (without the missing residues) as a template and building a comparative model using the full sequence [33]. The most comprehensive catalog of protein sequence and functional annotation, UniProt, helped us to get the amino acid sequence of p38 protein in FASTA format. The secondary structure of the protein was further evaluated with Ramachandran graph [34].

Molecular dynamics simulation
The stability of the inhibitor complex in the active site of the protein under physiological condition was analyzed by molecular dynamics simulation. Groningen Machine for Chemical Simulation (GROMACS 2018.1) software codes were utilized for exploring molecular dynamic simulations for 100 ns [35][36][37]. The snapshot of ligand 21 was considered for the study for which the temperature and pressure were set to 300 K and 1 bar respectively [38]. The energy minimization protocol was applied, and the respective charges were neutralized by addition of ions.
In the equilibration step, each energy minimized structure is simulated for 100 ns with all protein backbone atoms restrained to their initial positions as defined by the crystallographic structure coordinates. The binary files of simulation were analyzed to extract the values of root mean square deviation (RMSD), radius of gyration (Rg), root mean square fluctuation (RMSFs) and hydrogen bond. The three-dimensional visualization of the simulation was aligned and visualized using PyMOL [39].

3D QSAR analysis
Vlife MDS tool was employed for establishing quantitative structure-activity relationship studies. Alignment is an important part of 3D QSAR, and template-based alignment was adopted based on the common scaffold. Calculation of different descriptors is required for current QSAR studies which include electrostatic, steric and hydrophobic fields which were considered to design the QSAR model. The statistical data correlation ship afforded three 3D QSAR models and was as listed in Table 1 As briefed above, the statistical data of the model 3 is best as compared to the rest of two models on the basis of the correlation coefficient and the predictive r 2 of the test set. The value of 86.44% predictive ability hints promising as compared to the rest of two QSAR models. The steric and electrostatic field descriptors helped to draw the 3D contour maps of model 3 which in turn contributed to the 3D QSAR model with its coefficient 792.4230 (S_109), −0.1322 (S_576), 278.2350 (S_1496) for the steric fields and − 0.0475 (E_653) for the electrostatic fields as shown in Fig. 2 and Eq. 1.
Steric fields contribute positively for the P38 kinase inhibitory activity from the QSAR equation generated.
The above equation and the contribution plot afforded a perfect understanding that the steric parameter S_109 and S_1496 contribute positively, while S_576 is a negative contributor of enzyme inhibitory activity. The electrostatic parameter E_653 is also found as a negative contributor for P38 kinase inhibitory activity. The 3D QSAR model 3 has

Homology modeling
Script-based Modeller tool of homology modeling was used to add the missing amino acid in 2YIS, which was deposited in the repository with thirty-three amino acid residues found missing. Structural evaluation was performed on best model based on the dope score among twenty models. The Ramachandran plot of the best homology modeled protein structure depicted that 92.7% of amino acids are in the favored region and 7% were found in additional allowed and generously allowed regions as shown in Fig. 3. The alignment of the protein crystal structure and homology modeled protein structure was done using PyMOL which was suggested at a favorable RMSD of 2.0 Å as shown in Fig. 4.

Molecular docking
Schrodinger's Glide XP docking protocol was employed in order to dwell deep into binding interactions of study p38 kinase inhibitors. The cocrystal ligand was redocked, and best conformation of the docked molecule was aligned over the cocrystal ligand using PYMOL which was in acceptable range as shown in Fig. 5. In the study, p38 kinase inhibitors were found to be type 2 kinase inhibitors and interacted with amino acid residues in DFG out mode.

Molecular dynamic simulation studies
Molecular dynamics simulation (MD) was employed in order to corroborate of molecular docking results. Molecule 21 (most active compound) was selected for MD analysis at 100 ns by using explicit solvent model that helped in detail analysis of structural and dynamic performance of p38 kinase protein and protein-ligand complex respectively. The higher the RMSD value, the lesser the stability with large conformational changes, and the lower the RMSD value, the higher the conformational stability of the macromolecular system. The p38 kinase backbone RMSD revealed 0.15 Å at the primary stage, and a rise in spike was noticed at 20 ns to 0.35 Å till it equilibrated to 0.25 Å at 30 ns. Rise in the RMSD was spotted at 45 ns for 0.43 Å, and a drastic decrease in RMSD was observed till 50 ns. In further course of MD trajectory, a little decline in the RMSD was stated at 63 ns of 0.38 Å, and such cycles of increase and decrease in RMSD were observed from 66 to 80 ns in the range of 0.36-0.4 Å. From 81 to 100 ns protein backbone equilibrated to its stable conformation at 0.35 Å till the end of the MD duration discretely which is clearly highlighted in Fig. 8. The ligand conformational fluctuations of trajectories unearthed a handy similarity with the RMSD of the protein backbone. Conformational alterations in the initial spikes were observed with RMSD of 0.15 Å, which upraised to 0.35 Å at 15 ns and subsequently declined to 0.2 Å in 15-20 ns as rendered in Fig. 8. A petite upsurge in the spikes was detected in 20-30 ns MD simulation time scale, and marginal higher fluctuations were perceived from 31 to 40 ns amid 0.2-0.35 Å. Ligand conformational stability was presented from 42 to 70 ns, and additional spike after 70 ns was comprehended with corresponding RMSD from 0.25 to 0.35 ns, and subsequently, ligand conformation was equilibrated to 0.3 Å ns at the end of 100 ns MD simulation, which hinted at strong conformational stability within the active site of the protein.

Root mean square fluctuation (RMSF)
Another parameter in MD simulation study considered for evaluation is root mean square fluctuation, which tips off the excellence of molecular dynamic simulation by signifying the amino acid residue's average vacillations during the time scale from the initial coordinates. RMSF is the average point or fluctuation of amino acids from the initial coordinates. Figure 9 illustrated the RMSF of p38 kinase amino acids in which α-helical part is displayed in magenta while green parts are β-sheets, whereas the blue color is implied for amino acids in the active site which may have a probable chance of interaction with the inhibitor. Higher fluxes have been perceived for the terminal amino acids and were not in close propinquity with the binding site. However, minor conformational changes which are insignificant could be seen in the active sites and chains. Overall, the amino acid fluctuations which have juxtaposition with inhibitors were within the acceptable range which is less than 2.0 Å.

Radius of gyration (Rg)
The compactness of protein structure, especially tertiary structure, in the MD study can be explained on the basis  of radius of gyration. Radius of gyration explains about the mass distribution within the molecular structure from reference point and is expressed as root mean square. MD trajectory file used in the calculation of Rg of p38 kinase inhibitor 21 complex has been depicted in Fig. 10. Rg at the early part of MD simulation was found to be 2.25 nm and set a direction in downward to be able to attain to 2.2 nm at 10 ns which is followed by rise and declined RG, and at the end it retained a constant RMS of 2.25 nm.

Hydrogen bond analysis
The MD simulation of protein ligand complex also helps to understand the contribution of hydrogen bonding interactions between the protein and ligand. The trajectories revealed that the macromolecule and inhibitor 21 consistently sustained at least six hydrogen bonding interactions throughout the stipulated MD simulation thus suggesting a strong and substantial contribution of hydrogen bonding interactions within the complex as depicted in Fig. 11. The superimposed structure of the ligand throughout the molecular dynamics study which included the initial and post dynamic simulation study has been shown in Fig. 12.

Conclusions
Inhibitors of P38 kinase signaling pathway are considered as important in controlling inflammatory diseases including Alzheimer's disease. Till today, no p38 kinase inhibitors have crossed clinical trial hurdles and are in clinical utility, and hence, an indispensable need of new p38 kinase inhibitor discovery persists among the medicinal chemist. A statistically significant 3D QSAR model was developed using twenty-one training set molecules which resulted in q 2 value of 0.6269, correlation coefficient r 2 of 0.8783 and predicted r 2 of 0.8644. Both the training and test set molecules exhibited interactions in the binding site which were typical type II p38 kinase inhibitors and maximum molecules were hydrogen bonded with Asp186 and Lys67. Also, compound 21 was subjected for 100 ns molecular dynamics simulation studies which indicated stable position of the inhibitor.
Author contribution CA: project execution and writing; KM: compilation of data; BC: execution of molecular modeling studies especially MD; AM: data compilation and manuscript correction; RK: projection conception, planning, data collection and manuscript construction.
Availability of data and materials All the data out from software have been included in the manuscript.

Code availability
We have used maximum GUI-based software except gromacs which is used for molecular dynamic simulation on UBUNTU. Commands have been executed which are available free in public domain as referred. We have not created any new codes during the entire study.

Conflict of interest
The authors declare no competing interests.