Identification of Potential COX-2 Inhibitors for the Treatment of Inflammatory Diseases Using Molecular Modeling Approaches

Non-steroidal anti-inflammatory drugs are inhibitors of cyclooxygenase-2 (COX-2) that were developed in order to avoid the side effects of non-selective inhibitors of COX-1. Thus, the present study aims to identify new selective chemical entities for the COX-2 enzyme via molecular modeling approaches. The best pharmacophore model was used to identify compounds within the ZINC database. The molecular properties were determined and selected with Pearson’s correlation for the construction of quantitative structure–activity relationship (QSAR) models to predict the biological activities of the compounds obtained with virtual screening. The pharmacokinetic/toxicological profiles of the compounds were determined, as well as the binding modes through molecular docking compared to commercial compounds (rofecoxib and celecoxib). The QSAR analysis showed a fit with R = 0.9617, R2 = 0.9250, standard error of estimate (SEE) = 0.2238, and F = 46.2739, with the tetra-parametric regression model. After the analysis, only three promising inhibitors were selected, Z-964, Z-627, and Z-814, with their predicted pIC50 (−log IC50) values, Z-814 = 7.9484, Z-627 = 9.3458, and Z-964 = 9.5272. All candidates inhibitors complied with Lipinski’s rule of five, which predicts a good oral availability and can be used in in vitro and in vivo tests in the zebrafish model in order to confirm the obtained in silico data.


Introduction
Inflammatory processes stand out for their pathophysiological aspect, as they are caused by pathogenic microorganisms, such as viruses, bacteria, fungi, and parasites that invade host cells to

Introduction
Inflammatory processes stand out for their pathophysiological aspect, as they are caused by pathogenic microorganisms, such as viruses, bacteria, fungi, and parasites that invade host cells to reproduce, resulting in a complex and heterogeneous group of diseases that cause morbidity and mortality to those affected by them. At the same time, the need for specific attention to protect the integrity of human organisms from harmful or exogenous agents is emphasized [1].
Inflammatory mediation is regulated by the action of neutrophils, mast cells, eosinophils, macrophages, dendritic cells, and epithelial cells. It is a complex process involving vasodilation, chemotaxis, and increased permeability. Sensors called Toll-like receptors (TLR) recognize the products of pathogens, such as endotoxins and bacterial DNA that are located in the plasma membrane and in the endosomes, and they are capable of detecting intra and extracellular microorganisms [2].
When an injury occurs, platelets release complement proteins, where mast cells degranulate releasing histamine (vasodilation) and serotonin (cell permeability and diapedesis), where neutrophils are activated and migrate to the site of action induced by chemokines. Neutrophils phagocytose pathogenic organisms releasing mediators and attracting macrophages, increasing the release of pro-inflammatory mediators (prostaglandins and leukotrienes) and cytokines (interleukin 1 (IL-1), interleukin 6 (IL-6) and tumor necrosis factor (TNFα)). When cells are activated, the arachidonic acid (AA) of the cell membrane is converted by enzymes for the synthesis of prostaglandins and leukotrienes.
In one of the stages of the inflammatory process, prostaglandins, also called eicosanoids, are synthesized by triggering various stimuli that activate cell membrane receptors, when coupled with a regulatory protein that results in the activation of phospholipase A2 or with an increase in the concentration of Ca +2 . This type of enzyme hydrolyzes membrane phospholipids, consequently releasing the cascade of arachidonic acid, which is a substrate for the synthesis of physiopatological and inductive prostaglandins; see Figure 1 below [3].  Prostaglandins-endoperoxide synthase (PTGS) are known as cyclooxygenases responsible for the synthesis of prostaglandins that are percussive biologically active molecules. Figure 1 shows that the conversion of AA into signaling molecules takes place in 2 moments. (i) The enzyme cyclooxygenase-1 planning of drugs, with already countless cases of success involving or using computer simulations citing as an example the main factors: losartan, atorvastatin, and celecoxib.
Mathematical analyses accompany in silico studies in order to enable the reduction of costs and time to obtain positive results, observing the molecular structures and their possible affinities with the therapeutic target, using the quantitative structure-activity relationship (QSAR). This method aims to build parametric models for predicting inhibitory activity (IC 50 ) correlated with dependent variables such as physical-chemical, biological, and toxicological properties [10].
In these analyses, reduction filters are applied to the models used to predict inhibitory activity, correlating the molecule structures with their activity and toxicological potential, and comparing them with a positive control. In 2016, Brick et al. [10] applied the QSAR analysis to identify new antimalarial inhibitors from 1H-Imidazol-2-IL-Pyrimidine-4,6-Diamines, with reducing filters to eliminate descriptors that did not show correlation or information relevant to the process of statistical and toxicological analysis, beginning the screening with 107 compounds from the ZINC database and ending with four more promising compounds.
In this context, one of the main types of studies in progress by the scientific community is the in silico and in vivo studies of inhibitors for the inflammatory processes to search for new selective molecules. In parallel, the QSAR study (quantitative structure-activity relationships) uses multiparametric models that interrelate biological activity with the physical-chemical properties of selected molecules in order to predict their inhibitory capacity against the inflammatory mechanism [10]. Therefore, the objective of this work is the virtual screening of analogs of rofecoxib (Vioox ® ), based on pharmacophore and QSAR analysis, understanding approaches and molecular modeling techniques through free software that is easily accessible by the scientific community, in parallel with the prediction of pharmacokinetic properties and toxicity that show the possible effectiveness of the selected structures, according to the methodological scheme presented in Figure 2 (see more details in the Materials and Methods section).
Molecules 2020, 25, x FOR PEER REVIEW 4 of 32 The process from a biological target project to a new drug discovery can take an average of 10 years or more, and the computational chemistry comes with an excellent direction in the rational planning of drugs, with already countless cases of success involving or using computer simulations citing as an example the main factors: losartan, atorvastatin, and celecoxib.
Mathematical analyses accompany in silico studies in order to enable the reduction of costs and time to obtain positive results, observing the molecular structures and their possible affinities with the therapeutic target, using the quantitative structure-activity relationship (QSAR). This method aims to build parametric models for predicting inhibitory activity (IC50) correlated with dependent variables such as physical-chemical, biological, and toxicological properties [10].
In these analyses, reduction filters are applied to the models used to predict inhibitory activity, correlating the molecule structures with their activity and toxicological potential, and comparing them with a positive control. In 2016, Brick et al. [10] applied the QSAR analysis to identify new antimalarial inhibitors from 1H-Imidazol-2-IL-Pyrimidine-4,6-Diamines, with reducing filters to eliminate descriptors that did not show correlation or information relevant to the process of statistical and toxicological analysis, beginning the screening with 107 compounds from the ZINC database and ending with four more promising compounds.
In this context, one of the main types of studies in progress by the scientific community is the in silico and in vivo studies of inhibitors for the inflammatory processes to search for new selective molecules. In parallel, the QSAR study (quantitative structure-activity relationships) uses multiparametric models that interrelate biological activity with the physical-chemical properties of selected molecules in order to predict their inhibitory capacity against the inflammatory mechanism [10]. Therefore, the objective of this work is the virtual screening of analogs of rofecoxib (Vioox ® ), based on pharmacophore and QSAR analysis, understanding approaches and molecular modeling techniques through free software that is easily accessible by the scientific community, in parallel with the prediction of pharmacokinetic properties and toxicity that show the possible effectiveness of the selected structures, according to the methodological scheme presented in Figure 2 (see more details in the Materials and Methods section).

Molecular Optimization and QSAR Analysis
The structure of rofecoxib was selected as a pivot given its potential to mitigate the gastrointestinal effects compared to other selective drugs. Although this drug has the unwanted effect of acute myocardial infarction, the objective is to detect essential pharmacophore characteristics through virtual screening so that the selected promising structures have the same effectiveness. On the other hand, this is one of the only structures that presents the complexed structure in the Protein Data Bank (PDB) (5KIR, https://www.rcsb.org/) for the Homo sapiens organism, bringing the results of an ideality in front of the human organism.
The 32 molecules (Rofecoxib as a pivot) for the analysis were selected in the BindingDB database (https://www.bindingdb.org/bind/index.jsp) obeying an increasing order of IC 50 , with specific activity related to COX-2 and the Homo sapiens organism, in addition to not repeating inhibitory activity values, which could impact false-positive results through a straight line adjustment facilitated by of statistical analysis.
The molecular optimization values are shown in Table S1. The overlapping process was carried out by selecting molecules with the lowest energy value (PM3), since the optimization of molecular structures aims to bring the real structure closer to the energy minimum conformation, and with the observed experimental data, the optimization time quantification aims to elucidate the computational cost, as it is an expensive and time-consuming process [11].
Later, they were submitted to the PharmaGist software (https://bioinfo3d.cs.tau.ac.il/PharmaGist) for the extraction of physicochemical properties and construction of structure-activity relationship models (QSAR). The characteristics were analyzed with the aid of the Statistica ® software, where the most relevant ones were used to predict the inhibitory activity as a function of the pIC 50 value to decrease statistical inconsistencies. This software is capable of predicting the relationship between the inhibitor structure and its inhibitory activity, with a Pearson correlation cut-off of 0.4, obtaining a training set with n = 20 structures (methodology adopted by Santos, Cruz and Santos) [12][13][14]. Table 1 shows the selected descriptors. The atoms (ATM) characteristic presented the best correlation among all descriptors, with a value of 0.7651, allowing inferring that the number of atoms significantly interferes in the pIC 50 responses of the selected molecules. However, it must be noted that the selected regression model is tetra-parametric, so the prediction analysis must take into account the contribution of each descriptor in the process of prediction of the inhibitory activity value, as is the case of aromatic characteristics (ARO) with a p value of 0.7358 and acceptors (ACC) with a p value of 0.6399, which also contributes to the prediction of the inhibitory activity.
This result can also be accompanied by the analysis of hierarchical cluster analysis (HCA) ( Figure S1) performed with the aid of the Minitab ® Trial software, allowing observation of the similarity between the physical-chemical characteristics and the inhibitory activity of the respective molecules, corroborating with the data obtained by Pearson's correlation. The characteristics of ATM and ARO show greater proximity to the predicted pIC 50 . The descriptor ACC is inversely proportional to the predicted pIC 50 value, which indicates that the presence of hydrophilic groups capable of establishing hydrogen-bonding interactions can increase the inhibitory potential of the selected structures.
The ATM characteristic may not be essential when analyzed individually; however, we observed that the greater the number of atoms present in a structure, the greater its volume and topological polar surface area (TPSA), which are both characteristics that are essential for good oral absorption of the medication in the body, consequently obeying the Lipinski rule. On the other hand, it is not the only relevant characteristic for the prediction of the values of inhibitory activity and must be corroborated with the other characteristics provided by the statistical analysis [15]. Figure S2 represents the PCA analysis for the selected molecules. It correlates its characteristics with the inhibitory activity; the compounds with the lowest activity are in red, and those with the best activity are in blue. Molecule 11 is displaced from the others because it presents values of hydrogen donors equal to 5 In addition, the values of the number of atoms provide a better forecast: molecule 11 (37 atoms) shows one of the lowest values for ATM, as well as 15, 18, 19, and 20. It is observed that for the most active molecules, the ATM characteristics are relevant to the value of inhibitory activity, shifting them to the most active side. All the most active structures have four aromatic groups, except for molecule 4 which has only 2; however, its inhibitory activity is accentuated by the number of atoms in its structure (50).
In parallel, it must be understood that the predicted activity depends on the correlation between the four selected characteristics and their relative weight, and that the objective of the preliminary QSAR analysis is to investigate the most relevant characteristics among the data presented in the sampling of the structures already reported with inhibitory activity that are selective for COX-2. Figure S3 shows the analysis of HCA for the selected inhibitors. The HCA analysis gathers in hierarchical groups by similarity; the most active are represented in blue, and the least active ones are in red. It is observed that the data from the cluster followed the data obtained previously via QSAR and principal component analysis (PCA) analyses.
The tree-like dendrogram is seeking the structural similarities and the response to the inhibitory activity. It is noticed that 2 to 9 have four aromatic groups, with the exception of 4, which has only two, but its activity is enhanced by the number of atoms of the structure. Figure 3 shows the structure of the eight most active molecules. The common group observed in these inhibitors is pyrrole, in addition to the methylsulfone group, with the exception of 4. The presence of the methylsulfone group resembles the structure of rofecoxib and differs from the structures of the other coxibs (celecoxib and etoricoxibe), which have a sulfonamide functional group that is responsible for their toxicological characteristics. Pyrrole gives the appropriate lipophilic character to the molecule, which can help the molecule enter the COX-2 active channel [16].  Table 2 shows the regression data of the descriptors used to verify the best model for predicting the inhibitory activity. A combination has been used to evaluate the statistical parameters and select the parametric prediction equation according to the best fit. It is observed that the physical-chemical parameters ATM, ACC, and ARO are significant for the final calculated pIC50 values. The best statistical parameters were obtained for the parametric tri and tetra models with R 2 values of 0.9599 and 0.9617, and variance ratios of 62.6373 and 46.2719, respectively. It is emphasized that the greater the number of equated scores, the greater the quality of the predictor model, although the values are statistically close [10]. Note that the parametric tetra prediction values were better adjusted with a correlation R = 0.9617 and standard error of estimate (SEE) = 0.2238, with a notable predictive capacity. Such alignment can be compared with the residual values found during the validation step of the equation, with Δ4 differences close to 0.2, which demonstrates the ability to predict the values of inhibitory  Table 2 shows the regression data of the descriptors used to verify the best model for predicting the inhibitory activity. A combination has been used to evaluate the statistical parameters and select the parametric prediction equation according to the best fit. It is observed that the physical-chemical parameters ATM, ACC, and ARO are significant for the final calculated pIC 50 values. The best statistical parameters were obtained for the parametric tri and tetra models with R 2 values of 0.9599 and 0.9617, and variance ratios of 62.6373 and 46.2719, respectively. It is emphasized that the greater the number of equated scores, the greater the quality of the predictor model, although the values are statistically close [10]. Note that the parametric tetra prediction values were better adjusted with a correlation R = 0.9617 and standard error of estimate (SEE) = 0.2238, with a notable predictive capacity. Such alignment can be compared with the residual values found during the validation step of the equation, with ∆4 differences close to 0.2, which demonstrates the ability to predict the values of inhibitory activity. Table 3 analyzes the predicted models for the molecules, allowing inference of the difference (∆ = residual) between the pIC 50 values found in an experimental way and the statistical prediction values for a defined parametric model and the equation was determined considering the highest statistical correlation.

Prediction Equations
Mono pIC 50  The best results in question were for 10 (∆4 = 0.0075) and 12 (∆4 = 0.0020) inhibitors, although 2 (∆4 = 0.4764), 16 (∆4 = −0.2779), and 18 (∆4 = −0.2625) present residues greater than 0.2 in relation to the experimental data obtained. The margin of error (SEE = 0.2238) allows us to infer that the two may be within the desired perspective of residue, mainly, less than 1. The internal validation set demonstrated the detection of anomalous samples, which were excluded from the test set because they reduced the statistical correlation of the applied parametric method, with residues greater than 0.4, increasing the estimated error and deviating the correlation of the predicted values with the experimental data, which is justified by its exclusion from the test set initially. However, the reinclusion of these aims to determine the predictive capacity of the model, and results with residuals less than 1 are significant.  Table 4 shows the predicted values for the external validation set, applying the equations according to Table 4. For the validation test, a set between 20% and 30% of the total of the original set were used in the predictive model in order to prove its robustness [10]. It is observed that the values have good predictive quality for the molecules selected as an external set, with greater proximity for 26, 27, and 30, with residues close to 0.1. This shows that the model has a significant correlation between the descriptors used.

Virtual Screening and Analysis of Pharmacokinetic and Toxicological Properties
After selecting the best inhibitors, these were used on the Protox II and Molinspiration servers to select the reduction filters; the compounds were directed to the ZincPharmer database through the Pharmit web server (http://pharmit.csb.pitt.edu/search.html) shown in Table S2, with the maximum and minimum values being selected in order to limit the promising structures to those within the preapplied characteristics through the QSAR analysis, with a maximum limit of 2000 structures to be selected. The filters are applied on the online server, as well as the pharmacophore coordinates (below) elucidated for possible data reproduction and comparison with statistical analysis.
The pharmacophore structure obtained through the PharmaGist server is demonstrated, aligning the similarities of the twenty selected molecules (Table 5). It is observed that the characteristics of the pharmacophore follow the data obtained through statistical analysis, presenting two aromatic groups, two hydrophilic groups, and hydrogen receptors. Such pharmacophore characteristics are essential when compared to the central process molecule, which has two ARO groups and four ACC groups, allowing the tracking of molecules with physical and chemical  Table 4 shows the predicted values for the external validation set, applying the equations according to Table 4. For the validation test, a set between 20% and 30% of the total of the original set were used in the predictive model in order to prove its robustness [10]. It is observed that the values have good predictive quality for the molecules selected as an external set, with greater proximity for 26, 27, and 30, with residues close to 0.1. This shows that the model has a significant correlation between the descriptors used.

Virtual Screening and Analysis of Pharmacokinetic and Toxicological Properties
After selecting the best inhibitors, these were used on the Protox II and Molinspiration servers to select the reduction filters; the compounds were directed to the ZincPharmer database through the Pharmit web server (http://pharmit.csb.pitt.edu/search.html) shown in Table S2, with the maximum and minimum values being selected in order to limit the promising structures to those within the pre-applied characteristics through the QSAR analysis, with a maximum limit of 2000 structures to be selected. The filters are applied on the online server, as well as the pharmacophore coordinates (below) elucidated for possible data reproduction and comparison with statistical analysis.
The pharmacophore structure obtained through the PharmaGist server is demonstrated, aligning the similarities of the twenty selected molecules (Table 5). It is observed that the characteristics of the pharmacophore follow the data obtained through statistical analysis, presenting two aromatic groups, two hydrophilic groups, and hydrogen receptors. Such pharmacophore characteristics are essential when compared to the central process molecule, which has two ARO groups and four ACC groups, allowing the tracking of molecules with physical and chemical characteristics closer to rofecoxib.  On the other hand, in studies by Chakraborty, Sengupta, and Roy (2004) [17], linear multiple regression (RML) analyses were used to deduce statistically acceptable equations. The variation ratios were 0.675 for COX-1 and 0.842 for COX-2, observing three important pharmacophores groups: methyl sulfonyl portion, central phenyl ring, and terminal phenyl ring. These are relevant when compared to their affinity with the lipophilic channel present in the active sites of the enzymes, corroborating with the data obtained in the present study.
After the application of the reduction and selection filters of the 2000 compounds, they were submitted to similarity to Tanimoto to find out which ones are closer to the characteristics of the pivot molecule used in the process (Rofecoxib). The fifty eight (58) molecules were obtained with a Tanimoto index greater than 0.35 (see Table 6), which is a value that is considered reasonable for the application of toxicological and pharmacological prediction studies in silico, with three promising molecules being selected during these tests as reported below, subsequently applying the molecular coupling and molecular dynamics tests [18]. Table 6. Similarity studies of molecules using the Tanimoto Index.

Zinc
Tanimoto Index (Ti) > 0. 25 Table 7 shows the best results of the toxicological tests applied to the three inhibitors selected through the Tanimoto index and tests performed through the online server PreADMET (https://preadmet.bmdrc.kr/adme/) in order to screen those who present better absorption, distribution, and metabolism values besides limiting the possibilities of mutagenicity through toxicological tests. It is observed that molecules showed high LD50 values, with the exception of Z-627, but it presents good values for the absorption and distribution tests, contributing to its selection in the molecular docking tests. Carcinogenicity tests for rats and mice demonstrate a possibility of mutation for all the selected inhibitors; however, when compared to the control compound, it also observed this important side effect, and accordingly, it did not prevent the selection of these molecules. At the same time, the Ames test was used as a cut-off parameter between the most promising and those that would be excluded from the subsequent steps, where those that showed a positive result were eliminated from the process. This test assesses the possibility of mutagenicity of chemical compounds in media with a low histidine concentration, which allows the strains of Salmonella typhimurium to change and return to a prototypical state, which directly influences the carcinogenic response. On the other hand, in studies by Chakraborty, Sengupta, and Roy (2004) [17], linear multiple regression (RML) analyses were used to deduce statistically acceptable equations. The variation ratios were 0.675 for COX-1 and 0.842 for COX-2, observing three important pharmacophores groups: methyl sulfonyl portion, central phenyl ring, and terminal phenyl ring. These are relevant when compared to their affinity with the lipophilic channel present in the active sites of the enzymes, corroborating with the data obtained in the present study.

Pharmacophore
After the application of the reduction and selection filters of the 2000 compounds, they were submitted to similarity to Tanimoto to find out which ones are closer to the characteristics of the pivot molecule used in the process (Rofecoxib). The fifty eight (58) molecules were obtained with a Tanimoto index greater than 0.35 (see Table 6), which is a value that is considered reasonable for the application of toxicological and pharmacological prediction studies in silico, with three promising molecules being selected during these tests as reported below, subsequently applying the molecular coupling and molecular dynamics tests [18].  Table 7 shows the best results of the toxicological tests applied to the three inhibitors selected through the Tanimoto index and tests performed through the online server PreADMET (https: //preadmet.bmdrc.kr/adme/) in order to screen those who present better absorption, distribution, and metabolism values besides limiting the possibilities of mutagenicity through toxicological tests. It is observed that molecules showed high LD 50 values, with the exception of Z-627, but it presents good values for the absorption and distribution tests, contributing to its selection in the molecular docking tests. Carcinogenicity tests for rats and mice demonstrate a possibility of mutation for all the selected inhibitors; however, when compared to the control compound, it also observed this important side effect, and accordingly, it did not prevent the selection of these molecules. At the same time, the Ames test was used as a cut-off parameter between the most promising and those that would be excluded from the subsequent steps, where those that showed a positive result were eliminated from the process. This test assesses the possibility of mutagenicity of chemical compounds in media with a low histidine concentration, which allows the strains of Salmonella typhimurium to change and return to a prototypical state, which directly influences the carcinogenic response. The pharmacokinetic data for distribution are shown in Table 8. The plasma protein binding values (PPB) refer to the degree of binding of the inhibitors with the proteins present in the blood and C brain /C blood represents the permeability of the blood-brain barrier. Compounds with C brain /C blood values less than 1 do not have activity on the central nervous system (CNS). It is observed that Z-964 shows 100% of binding with the plasma proteins, inferring the possibility of its bioaccumulation and a consequent increase in its half-life within the organism, since the unbound portion is metabolized, consequently is excreted, and the bound part is slowly released in order to maintain the balance of the medium. In parallel, Z-627 showed an 85% binding, indicating that 15% of the fraction will not be bound, which increases the efficiency of diffusion and penetration into the cell membranes [19]. All selected compounds have no activity on the central nervous system, as they show values below one. In silico values for absorption are shown in Table 9. The selected drug candidates showed high values for intestinal absorption (HIA> 94%), being one of the most important absorption, distribution, metabolism and excretion (ADME) properties [20]. The drug molecules are transported from the gastro-enteric tract to the blood circle and permeate the gastro-enteric membrane by various mechanisms, and among them, the activity of the P-glycoprotein must be taken into account. This P-glycoprotein is a common transporter in the intestinal penetration of drugs, inferring in the hypothesis that the inhibitors Z-964 and Z-814, because they present an in vitro inhibition of P-gp, decrease the efflux process through the passive permeability of the inhibitors, which is mediated by this protein. However, they have considerable absorption values when compared with those of the other molecules screened in this study.
The P MDKC permeability value is significant for the Z-814 inhibitor (28.3061 nm/sec), being higher than for the control compound. Values greater than two indicate a significant medication efflux. The compound Z-964 showed a low permeability MDCK (0.0517 nm/sec) and Z-627 approached the ideal (1.4352 nm/sec) [21,22]. In parallel, the Caco-2 permeability assay measures the flow rate of a compound through Caco-2 cell monolayers to predict the in vitro drug absorption, where values greater than two present drug efflux. Inhibitors Z-814 (PCaco-2 = 12.4185 nm/sec) and Z-964 (PCaco-2 = 42.9100 nm/sec) have higher values than rofecoxib (PCaco-2 = 2.7291 nm/sec), with the Z-627 inhibitor (PCaco-2 = 0.6460 nm/sec) having a lower value; however, it is not a P-gp inhibitor, which can significantly interfere with intestinal absorption [22,23]. Table 10 demonstrates the predicted data for the biological activity of the selected inhibitors and compares the results against the selected controls rofecoxib and celecoxib, which were obtained from the PASS server (http://www.pharmaexpert.ru/passonline/). Celecoxib is used in everyday clinical practice, being part of the set of external validation and molecular docking of this research. The three selected inhibitors showed Pa > Pi values, indicating the possibility of activity in relation to the reported biological activities, mainly in terms of anti-inflammatory responses. Z-627 has the best values for anti-arthritic (Pa = 0.985) and anti-inflammatory (Pa = 0.852) activities higher than controls (Z-627 = 0.852, rofecoxib = 0.828, celecoxib = 0.663). All the candidate inhibitors have the possibility of activity against COX-2, although it is below the value of the reference compounds. Nonetheless, it serves as a reference for possible activities that they may present during the in vivo tests to be performed. In addition, a prediction of adverse effects that they may have on the organisms was performed, verifying that all of them present a propensity of activity similar to the other selected control compounds (celecoxib and rofecoxib) in the case of extrapyramidal effects. However, they have a lower propensity for the emergence of gastrointestinal problems, such as ulcers, which is the main focus in the development of new selective anti-inflammatory drugs, and the Z-964 structure did not present the possibility of presenting such an adverse effect. Table 11 shows the physical-chemical data of the selected inhibitors. Knowing the possibility that the structures present adverse cardiotoxic risks, the results were compared with the molecule still marketed in celecoxib, now considering its performance. However, it is observed that it presents risks of myocardial infarction and heart failure, which are analyzed through the Metatox (http://way2drug.com/mg2/gen_meta_all.php) and the hERG study (Human ether-a-go-go), through the PreADMET server, which refers to the blocking of the potassium channel, and that may cause cardiac collateral damage; see Table 11 below.
In view of the foregoing, this fact still did not allow its withdrawal from the market, as follow-up and adequate dosage reduce side effects and toxicological risks, which is a valid narrative for every drug currently sold; therefore, its cost-benefit must be evaluated. For now, it is seen that the molecules present a risk similar to or below the molecule withdrawn from the market (rofecoxb) and that which is still on the market (celecoxib), which does not preclude the possibility of being evaluated as candidates for specific COX-2 inhibitors. It is observed that the structures Z-814 and Z-627 present low and medium hERG risk, respectively, being better or equal to the molecules already commercialized, which makes their application as candidates for inhibitors of cyclooxygenase-2 possible. The Z-964 structure remains in the study due to the good results of bioavailability, in order to evaluate its experimental response in another study, as well as the others.
All candidate inhibitors present physical and chemical data within the acceptable range, showing no violation (Nv) or violating Lipinski's rule of five. This rule says that drugs with good oral bioavailability must obey four physicochemical parameters: molecular weight (MW) ≤ 500 g/mol, octanol/water partition coefficient (log P) ≤ 5, the number of hydrogen-bond donor groups (nHD) ≤ 5, and the number of hydrogen-bond acceptor groups (nHA) ≤ 10, see Table 12.  Table 12 shows that they have good absorption or permeability [24]. The pIC 50 values (nM) were predicted for the selected molecules, see Table 13, according to the equations of Table 4, demonstrating acceptable values. The three selected inhibitors with Tanimoto index are shown in Figure 5.  Table 12 shows that they have good absorption or permeability [24]. The pIC50 values (nM) were predicted for the selected molecules, see Table 13, according to the equations of Table 4, demonstrating acceptable values. The three selected inhibitors with Tanimoto index are shown in Figure 5.   Figure 6 shows the poses calculated in relation to the deposited PDB complexes, with the deviation of the mean square root (RMSD) calculated at 0.91 Å for rofecoxib (RCX; 5KIR PDB code), 0.63 Å for celecoxib (CEL; 3LN1 PDB code) and 0.71 Å for indomethacin (IMS; 2YOE PDB code). Such a methodology provides alignment values for a maximum of 2 Å for the study of molecular docking, and accordingly, it validates the protocols used [25,26].   Figure 6 shows the poses calculated in relation to the deposited PDB complexes, with the deviation of the mean square root (RMSD) calculated at 0.91 Å for rofecoxib (RCX; 5KIR PDB code), 0.63 Å for celecoxib (CEL; 3LN1 PDB code) and 0.71 Å for indomethacin (IMS; 2YOE PDB code). Such a methodology provides alignment values for a maximum of 2 Å for the study of molecular docking, and accordingly, it validates the protocols used [25,26].

Molecular Docking
Molecules 2020, 25, x FOR PEER REVIEW 14 of 32 Table 12 shows that they have good absorption or permeability [24]. The pIC50 values (nM) were predicted for the selected molecules, see Table 13, according to the equations of Table 4, demonstrating acceptable values. The three selected inhibitors with Tanimoto index are shown in Figure 5.   Figure 6 shows the poses calculated in relation to the deposited PDB complexes, with the deviation of the mean square root (RMSD) calculated at 0.91 Å for rofecoxib (RCX; 5KIR PDB code), 0.63 Å for celecoxib (CEL; 3LN1 PDB code) and 0.71 Å for indomethacin (IMS; 2YOE PDB code). Such a methodology provides alignment values for a maximum of 2 Å for the study of molecular docking, and accordingly, it validates the protocols used [25,26].    It is known that COX-1 and COX-2 have practically identical tertiary structures; however, the main difference between both is the replacement of Ile 434 , His 513 , and Val 434 residues in COX-1 by Val 434 , Arg 513 , and Val 523 in COX-2, respectively. This allows an increase of approximately 25% of the active site that consists of a more accessible pocket with Arg 513 as a fundamental bonding site [16,27]. Figure 7a shows the main interactions of rofecoxib, within the pocket that provides a selective inhibition.

Molecular Docking
interactions with the amino acids of the hydrophobic region of the β leaf, Ser 530 , and Val 523 for the former, and Ser 530 , Phe 518 , Val 523 , and Leu 358 for the latter (Figure 7 and Table S3). The lipophilic channel of the enzyme is also constrained by the presence of Tyr 355 and Arg 513 on the enzyme surface, with the additional hydrogen-bond interaction between the Ala 527 and Val 523 phenolic group and the oxygen sulfone atoms of the structures.
On the other hand, in COX-2, there are some interactions that allow greater accessibility in the lipophilic channel in this isoform than in COX-1, which can be observed, indicating a greater ease of interaction with the active site of COX-2 via Phe 518 . This effect can also be translated by the negative Gibbs free energy required for the interaction to occur ( Figure S6). The inhibitors selected may have an equivalent affinity in relation to the selected control compounds. The interaction with Ser 353 in Z-814 demonstrates the possibility of a binding activity associated with low IC50 values [28][29][30]. In Figure 8 and Table S4, it is possible to verify the interactions of the selected inhibitors with the reference drug (CEL) against Mus musculus. Hydrophobic interactions are observed with Val 509 , It is observed that the selected molecules Z-814, Z-964, and Z-627 present a similarity of interactions with the amino acids of the hydrophobic region of the β leaf, Ser 530 , and Val 523 for the former, and Ser 530 , Phe 518 , Val 523 , and Leu 358 for the latter (Figure 7 and Table S3). The lipophilic channel of the enzyme is also constrained by the presence of Tyr 355 and Arg 513 on the enzyme surface, with the additional hydrogen-bond interaction between the Ala 527 and Val 523 phenolic group and the oxygen sulfone atoms of the structures.
On the other hand, in COX-2, there are some interactions that allow greater accessibility in the lipophilic channel in this isoform than in COX-1, which can be observed, indicating a greater ease of interaction with the active site of COX-2 via Phe 518 . This effect can also be translated by the negative Gibbs free energy required for the interaction to occur ( Figure S6). The inhibitors selected may have an equivalent affinity in relation to the selected control compounds. The interaction with Ser 353 in Z-814 demonstrates the possibility of a binding activity associated with low IC 50 values [28][29][30].
In Figure 8 and Table S4, it is possible to verify the interactions of the selected inhibitors with the reference drug (CEL) against Mus musculus. Hydrophobic interactions are observed with Val 509 , Phe 504 , Gly 512 , Ser 339 , and Leu 338 , and hydrogen-bond interactions are observed with Gln 178 , Phe 504 , and Ser 339 for CEL. In parallel, we can observe the interactions for the selected inhibitors, where Z-627 shows interactions with the hydrophobic residues Ser 339 and Val 509 as well as the control, and in addition, it presents a hydrogen-bond with Ala 513 , showing selectivity [31,32].
On the other hand, the molecule Z-964 shows greater interactions with Ala 513 , Ser 339 , and Val 509 in the lipophilic region present in the β-leaf of the enzyme and the hydrophilic residue Leu 338 , which links to the sulfonic group of both inhibitors (sulfonamide for CEL and methylsulfone for Z-964). The Z-814 molecule showed a lower affinity than the others, but it showed relative selectivity when it comes to the amino acid residues that are part of the interactions (Ser 339 , Val 509 , and Phe 504 (fluorine)) in the hydrophilic region of the molecule, which allows for interactivity in parallel with the CEL molecule. Furthermore, the data corroborate the QSAR analyses carried out when dealing with the connections with hydrogen acceptors, which are mainly influenced by the electronegativity of the selected structures. These interactions have already been observed in other studies, corroborating with the affinity data shown in Figure S5, in which Z-967 shows an energy of 10.00 kcal/mol and Z-964 shows an energy of 9.50 kcal/mol. These data are considered the most important ones [30,31]. Docking studies corroborate the preliminary QSAR results, as they consider that the presence of aromatic groups can influence the inhibitory activity of such molecules; nevertheless, chemical On the other hand, the molecule Z-964 shows greater interactions with Ala 513 , Ser 339 , and Val 509 in the lipophilic region present in the β-leaf of the enzyme and the hydrophilic residue Leu 338 , which links to the sulfonic group of both inhibitors (sulfonamide for CEL and methylsulfone for Z-964). The Z-814 molecule showed a lower affinity than the others, but it showed relative selectivity when it comes to the amino acid residues that are part of the interactions (Ser 339 , Val 509 , and Phe 504 (fluorine)) in the hydrophilic region of the molecule, which allows for interactivity in parallel with the CEL molecule. Furthermore, the data corroborate the QSAR analyses carried out when dealing with the connections with hydrogen acceptors, which are mainly influenced by the electronegativity of the selected structures. These interactions have already been observed in other studies, corroborating with the affinity data shown in Figure S5, in which Z-967 shows an energy of 10.00 kcal/mol and Z-964 shows an energy of 9.50 kcal/mol. These data are considered the most important ones [30,31].
Docking studies corroborate the preliminary QSAR results, as they consider that the presence of aromatic groups can influence the inhibitory activity of such molecules; nevertheless, chemical changes are necessary in order to decrease in the cytotoxic effect of the inhibitors when compared with the reference drugs, such as the replacement of the sulfonamide by a methylsulfone group (rofecoxib analogs) [28]. The QSAR analysis demonstrates a structure-activity relationship, as is the case with the characteristics ARO, ACC, and DONN, being closely linked with the possibilities of their interaction with the active enzyme site. Lipophilicity deals with an intrinsic relationship of the possibility of permeation and good oral availability, which was previously reported with obedience to the rule of five by Lipinski (logP ≤ 5) interacting with the side pocket of the enzyme [30,31].
The three structures were subjected to molecular coupling tests (Figure 9) to assess the possibility of being selective for COX-1 as well, which would determine the possibility of the appearance of undesirable adverse effects, such as gastrointestinal problems. They demonstrate a lower affinity possibility to the COX-2 enzyme as previously reported, with low bond energies changes are necessary in order to decrease in the cytotoxic effect of the inhibitors when compared with the reference drugs, such as the replacement of the sulfonamide by a methylsulfone group (rofecoxib analogs) [28]. The QSAR analysis demonstrates a structure-activity relationship, as is the case with the characteristics ARO, ACC, and DONN, being closely linked with the possibilities of their interaction with the active enzyme site. Lipophilicity deals with an intrinsic relationship of the possibility of permeation and good oral availability, which was previously reported with obedience to the rule of five by Lipinski (logP ≤ 5) interacting with the side pocket of the enzyme [30,31]. The three structures were subjected to molecular coupling tests (Figure 9) to assess the possibility of being selective for COX-1 as well, which would determine the possibility of the appearance of undesirable adverse effects, such as gastrointestinal problems. They demonstrate a lower affinity possibility to the COX-2 enzyme as previously reported, with low bond energies (Z-627 = −8.    It is observed that the activity ratio on COX-1 is lower when compared to COX-2, suggesting that the structures will not be highly selective for isoform 1, emphasizing that all the NSAIDs already prescribed have a small selectivity for this, which decreases the possibility of side effects [7]. On the other hand, the perspective of the advent of adverse effects can be compared in parallel with the prediction of biological PASS activity (Table 10), indicating a probability of few gastrointestinal effects. Furthermore, it is noted that selectivity in relation to COX-2 is given by the substitution of valine for an isoleucine in COX-1 in position 523, which in this case interacts with the phenolic ring of the selected structures. In addition, most inhibitors selective for COX-2 suggest not having free carboxylate groups, which contributes to this low affinity to isoform 1, and the high affinity is expressed by the interaction as the amino acid residue Arg 120 [33].
In this case, the Z-627 structure presents this interaction relationship with the Arg 120 residue in the pyrroline portion of the structure, showing a possible structural rigidity and suggesting a possibility of structural modification in order to further limit the relationship estimate. However, the addition of the methylene group in the residue from Ile 523 indicates that interactions are restricted in access to the COX-1 side pocket, directly impacting the time-dependent competitive inhibition process in relation to COX-2 [34,35].

Structure-Activity Relationship of the Most Promising Molecules
The selected compounds ( Figure 5) show a similar structure to that of the pivotal compound, having in their structures the methylsulfone group, showing no cytotoxic effect in relation to the sulfonamide group (Z-627 and Z-964). Small substituents are the best, because they influence the volume of the molecule and possible van der Waals interactions with COX-2, which is a fact observed in docking studies. The introduction of fluorinated groups may show more significant activity.
According to Hayashi et al. 2012 [36], substituted analogs by acceptable hydrogen-bond groups potentiate the inhibitor activity. On the other hand, substitutions in the isoindoline nucleus can contribute to the inhibitor-enzyme stabilization, further demonstrating the fundamental role that the electrostatic and dipole-dipole interactions can play [37,38].
At the same time, the endocyclic nitrogen atoms included in five-or six-membered cycles such as pyrrole, pyridine, and pyrimidine, among others, may produce an increase in selectivity. The fiveamino group in the isoindoline ring may favor the inhibitory activity of Z-627 and, moreover, possible hydrogen-bond interactions through the methylsulfone group. The inhibition mechanism depends on the prostaglandin biosynthesis by means of arachidonic acid (AA), estimating that AA fits into the channel cavity surrounded by amino acid residues with aromatic, aliphatic, and phenolic groups that establish several interactions. Therefore, competitive or selective inhibitors bind to Val 523 in COX-2, interfering with the arachidonic acid cascade and preventing the peroxidase action, as well as the formation of It is observed that the activity ratio on COX-1 is lower when compared to COX-2, suggesting that the structures will not be highly selective for isoform 1, emphasizing that all the NSAIDs already prescribed have a small selectivity for this, which decreases the possibility of side effects [7]. On the other hand, the perspective of the advent of adverse effects can be compared in parallel with the prediction of biological PASS activity (Table 10), indicating a probability of few gastrointestinal effects. Furthermore, it is noted that selectivity in relation to COX-2 is given by the substitution of valine for an isoleucine in COX-1 in position 523, which in this case interacts with the phenolic ring of the selected structures. In addition, most inhibitors selective for COX-2 suggest not having free carboxylate groups, which contributes to this low affinity to isoform 1, and the high affinity is expressed by the interaction as the amino acid residue Arg 120 [33].
In this case, the Z-627 structure presents this interaction relationship with the Arg 120 residue in the pyrroline portion of the structure, showing a possible structural rigidity and suggesting a possibility of structural modification in order to further limit the relationship estimate. However, the addition of the methylene group in the residue from Ile 523 indicates that interactions are restricted in access to the COX-1 side pocket, directly impacting the time-dependent competitive inhibition process in relation to COX-2 [34,35].

Structure-Activity Relationship of the Most Promising Molecules
The selected compounds ( Figure 5) show a similar structure to that of the pivotal compound, having in their structures the methylsulfone group, showing no cytotoxic effect in relation to the sulfonamide group (Z-627 and Z-964). Small substituents are the best, because they influence the volume of the molecule and possible van der Waals interactions with COX-2, which is a fact observed in docking studies. The introduction of fluorinated groups may show more significant activity.
According to Hayashi et al. 2012 [36], substituted analogs by acceptable hydrogen-bond groups potentiate the inhibitor activity. On the other hand, substitutions in the isoindoline nucleus can contribute to the inhibitor-enzyme stabilization, further demonstrating the fundamental role that the electrostatic and dipole-dipole interactions can play [37,38].
At the same time, the endocyclic nitrogen atoms included in five-or six-membered cycles such as pyrrole, pyridine, and pyrimidine, among others, may produce an increase in selectivity. The five-amino group in the isoindoline ring may favor the inhibitory activity of Z-627 and, moreover, possible hydrogen-bond interactions through the methylsulfone group. The inhibition mechanism depends on the prostaglandin biosynthesis by means of arachidonic acid (AA), estimating that AA fits into the channel cavity surrounded by amino acid residues with aromatic, aliphatic, and phenolic groups that establish several interactions. Therefore, competitive or selective inhibitors bind to Val 523 in COX-2, interfering with the arachidonic acid cascade and preventing the peroxidase action, as well as the formation of prostaglandins or thromboxanes (pro-inflammatory eicosanoids). In parallel with the studies carried out by Hayashi et al. 2012 [36], the best inhibitors Z-627 and Z-814 have two hydrogen-bond donors, as well as low values of TPSA and MW. For the three selected inhibitors we proposed theoretical synthetic routes-Supplementary Materials Figures S9-S11.

Molecular Dynamics Results and Affinity Energy
The studies of molecular dynamics simulations were carried out to understand more deeply the modes of interaction of the selected compounds with the target proteins. The results obtained through molecular dynamics simulations have served as support for the detailed evaluation of conformations over time observed in drug-receptor complexes [39][40][41]. Thus, understanding that the dynamics and changes in the movement of a protein are closely related to its biological function allows us to understand that the observation of these phenomena is extremely important. In this way, we carried out the investigation of the protein structure during the 100 ns of molecular dynamics simulations using the methods of root mean square deviations (RMSD) and root mean square fluctuations (RMSF). To plot the RMSD of the ligands, all the heavy atoms of the molecules were used, while to plot the RMSD and RMSF of the protein backbone, the Cα carbon atoms were used. In Figure 11, the graphs of the compounds that were bound to COX-2 of Homo sapiens are plotted, while in Figure 12, the RMSD plot of the complexes established with COX-2 of Mus musculus is displayed.
Molecules 2020, 25, x FOR PEER REVIEW 19 of 32 prostaglandins or thromboxanes (pro-inflammatory eicosanoids). In parallel with the studies carried out by Hayashi et al. 2012 [36], the best inhibitors Z-627 and Z-814 have two hydrogen-bond donors, as well as low values of TPSA and MW. For the three selected inhibitors we proposed theoretical synthetic routes -supplementary material from Figures S9, S10, and S11.

Molecular Dynamics Results and Affinity Energy
The studies of molecular dynamics simulations were carried out to understand more deeply the modes of interaction of the selected compounds with the target proteins. The results obtained through molecular dynamics simulations have served as support for the detailed evaluation of conformations over time observed in drug-receptor complexes [39][40][41]. Thus, understanding that the dynamics and changes in the movement of a protein are closely related to its biological function allows us to understand that the observation of these phenomena is extremely important. In this way, we carried out the investigation of the protein structure during the 100 ns of molecular dynamics simulations using the methods of root mean square deviations (RMSD) and root mean square fluctuations (RMSF). To plot the RMSD of the ligands, all the heavy atoms of the molecules were used, while to plot the RMSD and RMSF of the protein backbone, the Cα carbon atoms were used. In Figure 11, the graphs of the compounds that were bound to COX-2 of Homo sapiens are plotted, while in Figure 12, the RMSD plot of the complexes established with COX-2 of Mus musculus is displayed. Along the trajectories of MD simulations, COX-2 showed differences in the RMSDs of the complexes. The maximum Plator rising by the backbone RMSD was 3 Å, which was visualized in the The evaluations of the regions of the protein that obtained the greatest fluctuations along the trajectories of molecular dynamics were performed using the RMSF plot (see Figure 13). In general, the RMSF graphs showed a similar profile, even in the regions that suffered the greatest fluctuations. The greatest fluctuations were observed in the N-terminal portion of the protein. This region is exposed to the solvent, being formed by alpha helices and beta leaves that are connected by long loop regions. Structurally, loops are the most flexible regions of the protein, so a region that exhibits many loops has a tendency to be flexible, as was observed in the RMSF plots displayed. Although this region is close to the active site of the ligands, its flexibility did not compromise the binding of the compounds, since all the ligands showed energy of favorable affinity with the protein, according to the molecular mechanics/generalized born surface area (MM-GBSA) results obtained. In addition, the fluctuation of this region did not affect the conformational stability of the ligand along the trajectory of molecular dynamics, as the RMSD graphics of the ligands Along the trajectories of MD simulations, COX-2 showed differences in the RMSDs of the complexes. The maximum Plator rising by the backbone RMSD was 3 Å, which was visualized in the COX-2-Z814 system, and the smallest fluctuation was observed in the COX-2-rofecoxib system. Despite the fluctuations displayed, this did not impair the interaction with the complexed ligands. It is important to note that the RMSD of the ligands showed low fluctuations and had a low RMSD value; this means that the ligands did not undergo drastic conformational changes after settling at the protein binding site.
Similar phenomena were observed for the complexes established between Mus musculus COX-2 and ligands. The backbone RMSD Plator was approximately 3 Å, and the ligands also remained in equilibrium throughout the 100 ns simulations, as observed in the RMSD graphs with small fluctuations.
The evaluations of the regions of the protein that obtained the greatest fluctuations along the trajectories of molecular dynamics were performed using the RMSF plot (see Figure 13). In general, the RMSF graphs showed a similar profile, even in the regions that suffered the greatest fluctuations. The greatest fluctuations were observed in the N-terminal portion of the protein.
This region is exposed to the solvent, being formed by alpha helices and beta leaves that are connected by long loop regions. Structurally, loops are the most flexible regions of the protein, so a region that exhibits many loops has a tendency to be flexible, as was observed in the RMSF plots displayed. Although this region is close to the active site of the ligands, its flexibility did not compromise the binding of the compounds, since all the ligands showed energy of favorable affinity with the protein, according to the molecular mechanics/generalized born surface area (MM-GBSA) results obtained. In addition, the fluctuation of this region did not affect the conformational stability of the ligand along the trajectory of molecular dynamics, as the RMSD graphics of the ligands demonstrated that they remained stable along the trajectories without showing drastic changes in the RMSD plot.  In addition to structural analysis of the protein and ligands using RMSD and RMSF, we also evaluated whether the compounds are capable of interacting favorably with molecular targets. For this, we use the MM-GBSA method. The results obtained are summarized in Table 14. All ligands have been shown to be able to interact favorably with COX-2. The selected compounds showed great affinity with COX-2 when we compared their values of affinity energy with the value obtained for the positive control of protein in the human body and Mus musculus. In addition to structural analysis of the protein and ligands using RMSD and RMSF, we also evaluated whether the compounds are capable of interacting favorably with molecular targets. For this, we use the MM-GBSA method. The results obtained are summarized in Table 14. All ligands have been shown to be able to interact favorably with COX-2. The selected compounds showed great affinity with COX-2 when we compared their values of affinity energy with the value obtained for the positive control of protein in the human body and Mus musculus.

Optimization of Molecular Structures and Determination of Pharmacophore Characteristics
The selected inhibitors were pre-optimized by means of Molecular Mechanics (MM+), followed by calculations of Austin Model 1 (AM1) and PM3 in the Hyper Chem 7 ® program (Table 2), with the lowest energy value used as a parameter of choosing the best model to carry out the construction of the pharmacophore hypothesis. Subsequently, the input was made to the PharmaGist Web Server 15 to determine the following characteristics: atoms (ATM), spatial characteristics (SF), characteristics (F), aromatic (ARO), hydrophobic (HYD), acceptor (ACC), and donor of hydrogen (DONN). The initial set presented 25 molecules, which were aligned according to the similarity with the selected pivot molecule, allowing the generation of pharmacophore models with the aid of the Discovery Studio ® v. 4.0 program, following the methodology developed by us [10,12,14,[57][58][59].

QSAR and PCA/HCA Studies
The inhibitory activity values were transformed into pIC 50 (−log (IC 50 )) in order to reduce the inconsistencies of the data obtained in an experimental way and homogenize the dataset, following the adopted methodological proposal [10,57,59]. In parallel, the importance of each pharmacophore descriptor was attributed-atoms (ATM), spatial characteristics (SF), characteristics (F), aromatic (ARO), hydrophobic (HYD), acceptor (ACC) and hydrogen donor (DONN); these were used for prediction in order to assess notoriety regarding the response to the pIC 50 value through the Pearson correlation (p), using the software Statistica 7.0 ® and Minitab 19 ® , adapting the methodology adopted by Santos et al. 2015 and Ferreira et al. 2019 [12,59]. Pearson's coefficient (Equation (1)) measures the degree of linearity between two variables, assuming a value between +1 and −1. If one variable tends to increase while the others decrease, the value is negative. On the other hand, if both increase, the coefficient is positive. Moreover, x is the sample mean for the first variable; s x is the standard deviation for the first sample; y is the sample mean for the second variable; s y is the standard deviation for the second sample; and n is the column length.
The best pharmacophore descriptors were obtained considering the statistical quality relations of multiple linear regression (MLR), such as correlation coefficient (R), correlation coefficient squared (R 2 ), explained variance (adjusted R 2 ), standard error of estimate (SEE), and variance ratio (F), and they were transformed into parametric models for predicting the inhibitory activity at pIC 50 values. The combinations were obtained using four parameters indicated by Pearson's correlation without repetition [12,59], according to Equation (2), where C = number of combinations, p = model type (p 0 and p = 4), and n = number of variables (n = 4).
For the prediction of the best model, in the internal validation stage, the random correlations between the descriptors and the inhibitory activity were measured to normalize the data obtained, applying the technique of detecting anomalous samples (outliers), in order to obtain a homogeneous set. This subset is considered as internal validation, for analysis of the prediction capacity of the selected model, comparing the data obtained during the two validations (internal and external; Figure S8). Principal component analysis (PCA) together with hierarchical cluster analysis (HCA) were applied in order to verify whether the model obtained corresponds to the degree of similarity, using Pearson's squared distance as a measurement parameter in the latter [60,61]. For the respective analyzes, Minitab v. 19 ® trial version was used.

Virtual Screening and Selection of Inhibitor Compounds
After selecting the best model via QSAR analysis, the selected molecules were superimposed to form a pharmacophore model. After inputing the pharmacophore, the search was performed within the ZINC database, selecting the 2000 most similar molecules, using the partition coefficient (log P), surface area (TPSA), number of atoms (Natoms), Molar Mass (MW), hydrogen acceptors (nHA), hydrogen donors (nHD), number of violations (Nv), number of revolutions (Nrotb), and volume, which were values as filters determined via Protox II servers (http://tox.charite.de/protox_II/) and molinspiration (https://www.molinspiration.com/cgi-bin/properties). The RMSD (Equation (3)) value was used as a reference parameter, which is the measure of the average distance between the atoms of the overlapping inhibitors, given in Angstroms, representing the quantitative similarity relationship between them. The lower the RMSD value, the better the model is compared to the target structure. δ 2 i is the distance between atom i of any reference structure or the average position of N equivalent atoms.
Then, the Tanimoto test was performed via the BindingDB server. The similarity was determined according to the chemotype of the compounds screened with the pivotal molecule of the selection process to reduce and optimize the selection of compounds for determining pharmacokinetic and toxicological characteristics, using a cut-off index of 0.3, applying Equation (4) below [22].
where M 11 -total number of attributes where A and B have a value of 1; M 01 -total attributes where A is 0 and B is 1; M 10 -total attributes where A is 1 and B is 0; M 11 -total attributes where A and B have a value of 0.

Prediction of Toxicological and Pharmacokinetic Properties
Pharmacokinetic and toxicological studies were applied to inhibitors extracted from Pharmit via the ZincPharmer server. PreADMET v. 2.0 (https://preadmet.bmdrc.kr/) was used, which is an application based on a database on the web that is used for the prediction of ADME data (Absorption, Distribution, Metabolism, Excretion) with the following being selected: blood-brain barrier (BBB) penetration, in vitro permeability in Caco2 cells, human intestinal absorption (HIA), in vitro permeability of MDKC cells, in vitro P-glycoprotein inhibition, plasma protein binding (PPB), and Toxicological for Ames_Test, Carcinogenicity for Rats and Mice. The LD 50 values were determined via Protox II servers (http://tox. charite.de/protox_II/) as well as toxicity class.

Prediction of Biological Activity of Selected Inhibitors
Activity predictions were made using the online PASS server (http://www.pharmaexpert.ru/ passonline), which allows you to predict the biological effects of compounds based on their formula using multilevel atom neighbors (VMA) descriptors, suggesting that the inhibitor's activity is expressed in terms of its chemical structure. Molecules with activities reported for anti-inflammatory and cyclooxygenase inhibitor effects were selected [25,61].

Molecular Docking
For this step, only the molecules with the best pharmacokinetic, toxicological, and biological parameters were selected for the study of molecular docking, in order to evaluate the interactions with selected inhibitors and the respective targets through the measurement of free energy interaction with amino acid residues and binding affinity. The crystallographic poses were extracted from the The selected inhibitors and proteins were prepared with the aid of Discovery Studio ® 4.0 software, and the evaluation of the complexes with the ligand was evaluated using the AutoDock 4.2/Vina 1.1.2 software and the PyRx graphical interface version 0.8.30 (https://pyrx.sourceforge.io), with the standard exhaustiveness parameter of the software being the best conformation obtained through the analysis of the RMSD value. The validation protocol was based on the determination of the x, y, and z coordinates according to the average region of the active site; these values are observed in Table 16. The energy function score was used to evaluate the free binding energy (∆G) of the interaction of the receptors with the ligands [25]. The calculation of binding affinity (∆G) was also performed in order to compare the actual data obtained and the values predicted in silico, which was the same methodology adopted by Santos et al., 2020 [14], according to Equation (5). ∆G = −RT ln K i (5) where R (gas constant) is 1.987.10 −3 kcal·mol −1 ·K −1 , the temperature is 310 K for rofecoxib/celecoxib, and K i is 310.10 −9 M for rofecoxib and 340.10 −9 M for celecoxib [28,32,52].

Molecular Dynamics Protocol
The initial structure for the system was obtained from molecular docking methods. The restrained electrostatic potential (RESP) protocol with the HF/6-31G* basis sets was applied to obtain the partial atomic charges of the atoms of each ligand [62][63][64][65]. The parameters of the ligand were constructed with the Antechamber module [66] using General Amber Force Field (GAFF) [67].
The amino acid protonation state was characterized using the PDB2PQR server [68]. The systems were built with the tLEaP module of the Amber 16 package [69][70][71]. The force field used to describe the protein in all simulations was ff14SB [72]. The protein-ligand system was solvated in an octahedron periodic box containing water molecules in the TIP3P model [73]. The partial charges were neutralized by adding counter-ions.
Energy minimization occurred in four stages. First, the water molecules and ions were optimized using 2000 cycles of the steepest descent and 3000 cycles of conjugate gradient. Then, the position of receptor-ligand hydrogen atoms was optimized using 4000 steps of the steepest descent algorithm and 3000 steps of the conjugate gradient. At the third stage, hydrogen atoms, water molecules, and ions were further optimized using 2500 steps of the steepest descent algorithm and 3500 steps of the conjugate gradient. Finally, all atoms were minimized using 3000 steps of the steepest descent algorithm and three steps of the conjugate gradient.
Molecular dynamics simulations were performed at a constant volume by heating the systems up to 298 K. This heating was performed in five steps for a duration of 1 ns. After 100 ns, production runs were performed for each system.
The Particle Mesh Ewald method [74] was used for the calculation of the electrostatic interactions, and the bonds involving hydrogen atoms were restricted with the SHAKE algorithm-Restriction algorithm used to ensure that the distance between points of mass is maintained [75]. The temperature control was performed with the Langevin thermostat [76] within a collision frequency of 2 ps −1 .
The free energy of bonding (∆G bind ) is the summation of the interaction energy of the gas phase among the protein-ligand (∆E MM ), desolvation free energy (∆G solv ), and system entropy (−T∆S). ∆E MM is the result of the sum of internal energy (∆E internal , sum of the energies of connection, angles and dihedral) electrostatic contributions (∆E ele ), and the van der Waals term (∆E vdW ). ∆G solv is the sum of the polar (∆G GB ) and non-polar (∆G NP ) contributions. ∆G SA was determined from the solvent accessible surface area (SASA) estimated by the linear combination of pairwise overlaps (LCPO) algorithm.

Conclusions
After the pharmacophore-based virtual screening, the QSAR analysis demonstrated a good line fit with R 2 = 0.96 and an equation with four main prediction parameters for pIC 50 , ATM, ARO, ACC, and DON, where the ARO, ACC, and DON report the relationship with the three new and promising compounds selected and the pivot structure (rofecoxib). The development of the predetermined multiple linear regression model predetermined the pIC 50 values for the selected compounds Z-814 = 7.9484, Z-627 = 9.3458, and Z-964 = 9.5272. In database searches to evaluate possible applications that may have already been carried out, these substances are not used in specific biological activities (https://scifinder.cas.org/ and https://zinc.docking.org/).
The analyzes of toxicological prediction and bioavailability confirm the possibility of significant activity of the structures with a reduction of possible undesirable effects, of which Z-627 was considered the most promising in view of all the tests applied via ADME analysis, without consequences for the CNS; this was corroborated with the main compounds selected. All selected compounds have the methyl sulfone group, unlike coxibs, which have the sulfonamide group. These three molecules do not present toxicological risks; they comply with the Lipinski rule of five, which provides for good oral availability, and PASS provides for a specific activity with a high probability of showing promising anti-inflammatory activity, in addition to dim side effects in relation to the compound's selected controls. Molecular coupling tests demonstrate strong energy affinity with isoform 2 and low activity with isoform 1 through relationship analysis, which induces a possibility of minor side effects. Finally, zebrafish larvae should be analyzed to assess anti-inflammatory activity in the treatment of inflammatory disorders to confirm in silico results.
Supplementary Materials: The following are available online, Table S1: Energy values of the optimized molecules, Table S2: Filters applied according to the properties of the selected molecules, Table S3: Distance of Interactions for the structures for the PDB 5KIR, Table S4: Distance of Interactions for the structures for the PDB 3LN1, Table S5: Distance of Interactions for the structures for the PDB 2OYE, Figure S1: Dendrogram representing clustering of pharmacophores, Figure S2: Analysis of the main components for the sorted molecules. Scores (a) and Loading Graph (b), Figure S3: Dendrogram of selected molecules. More active (blue) and less active ones (red), Figure S4: Binding affinity results of compounds, including Vioxx bound (COX-2 -Homo sapiens), Figure S5: Binding affinity results of compounds, including celecoxib (COX-2 Mus musculus), Figure S6: Binding affinity results of compounds, including Indomethacin (COX-1 Ovis aries), Figure S7: Structures used in the molecular modeling, Figure S8: Structures used in the external validation set, Figure S9: Theoretical synthetic route for the preparation of compound A (Z-814), Figure S10: Theoretical synthetic route for the preparation of compound B (Z-964), Figure S11: Theoretical synthetic route for the preparation of compound C (Z-627).