Pharmacophore as Fundamental Base of the NSAIDs Design: Molecular Docking and SAR/QSAR Analysis of Reported Experimental Data of Cyclooxygenase 1 and 2

1. Introduction

The inflammatory process goes through distinct phases, and the activation of phospholipases is essential for generating mediators for the free radical pathway, such as prostaglandins, lipoxins, thromboxanes, and leukotrienes, including platelet-activating factor (PAF). [1,2]

In contemporary medical practice, Nonsteroidal Anti-Inflammatory Drugs (NSAIDs) represent a fundamental pillar in the treatment of various painful and inflammatory conditions. Since their discovery and subsequent commercialization in the mid-20th century, these drugs have revolutionized the management of chronic and acute pain, as well as inflammatory diseases such as rheumatoid arthritis and osteoarthritis. However, its widespread use is not exempt from critical scrutiny due to its potential adverse effects, particularly at the gastrointestinal and cardiovascular levels, which have motivated a constant search for safer, more effective formulations. [3,4]

The chemical structure of NSAIDs plays a crucial role in their pharmacological activity and selectivity. These compounds usually consist of a central nucleus, frequently an aromatic ring, to which a carboxylic acid or an enol is attached. This basic structure facilitates binding to specific active sites on COX, allowing interference with enzymatic activity, and provides it with physicochemical and molecular properties shared by COX. Over the decades, researchers have explored various modifications to this structure to improve therapeutic efficacy and reduce side effects. [5,6]

Structure-Activity Relationship (SAR) studies have been instrumental in the present contribution, which describes the development and optimization of NSAIDs. These studies aim to identify how subtle modifications in molecular structure can influence the potency, actions, and pharmacokinetics of these drugs. [7,8,9,10,11] For example, it has been shown that the addition of substituent groups at specific positions on the aromatic ring can improve the slight tendency toward COX-2 without significantly affecting COX-1. Therefore, in this work, was identified by molecular docking and SAR/QSAR analyses the necessary structural moieties for the recognition and inhibitory activity on both COXs, were for the recognition being the aryl-carboxylate with captodative properties (primarily shown in salycilates and fenamates) the guide and accessory moieties, which led to establishment the pharmacophore structure for the COX non-selective NSAIDs.

Therefore, it is important to highlight that the present study identified through SAR/QSAR analysis and molecular docking the structural portions necessary for recognition and inhibitory activity on both COX, highlighting a carboxylate on an aryl group with captodative properties (mainly evident in salicylates and fenamates), which led to the establishment of the pharmacophore for non-selective COX NSAIDs.

2. Results and Discussion

2.1. Fenamates and Salicylates

SAR/QSAR analysis is fundamentally valuable for elucidating in greater detail the mechanisms of action of both novel and known molecules. Furthermore, it serves as a critical tool for establishing pharmacophoric groups, which function as foundational elements in drug design, by identifying the physicochemical and molecular rationales supporting such proposals. In this study, we specifically focus on cyclooxygenase inhibitors.

One approach to demonstrate the involvement of specific physicochemical and/or molecular properties in biological activity (COX inhibitory activity) is to construct simple correlations between biological activity and individual properties. These correlations for each COX inhibitor be helpful in identifying the structural requirements essential for biological activity. Consequently, we generated diagrams comparing physicochemical properties across distinct groups of NSAIDs with known experimental activity - specifically examining the relationship between molar refractivity (MR) and the partition coefficient of structurally related salicylate and fenamate derivatives (Figure 1A).

For the salicylate family, structural requirements for COX-1 and COX-2 recognition indicate optimal ranges of logP between 1.0 and 3.5 and molecular size between 33.0 and 102.0. cm³/mol. It is important to emphasize that these drugs are considered non-selective COX inhibitors [12,13], as confirmed in this study, where the average pIC₅₀ values for both enzymes showed no significant difference in vitro and in silico studies (Figure S1). This demonstrates that these physicochemical properties are essential for protein recognition.

For the fenamate family, no significant differences were observed compared to salicylates (Figure S1). The mean inhibitory activities against COX-1 and COX-2 were remarkably similar and did no differ significantly. However, a slight preference for COX-1 affinity wasnoted, suggesting that the molecular requirements for logP and MR properties are more favourable for COX-1 than COX-2. The optimal physicochemical ranges were logP 3.3-4.8 (indicating greater lipophilicity than salicylates) and MR 66.0-75.0 cm³/mol (a narrow mid-range), implying that molecular size has less influence on selectivity than lipophilicity.

The core structures of salicylates and fenamates differ primarily by the isosteric replacement of -OH with -NH₂ [salicylic acid (logP = 1.20, MR = 33.9 cm³/mol) vs. anthranilic acid (logP = 0.79, MR = 36.9 cm³/mol), respectively], conferring similar physicochemical properties. The observed differences in activity and selectivity stem from their respective accessory moieties.

Molecular docking revealed that fenamates form ion-ion interactions between Arg120 and the carboxylate group, except for flufenamic acid (FFA), which interacts with Ala527. Salicylates (sulfasalazine/SSZ and diflunisal/DF) showed ionic interactions with Arg120, while acetylsalicylic acid (ASA) and salicylic acid (SA) formed hydrogen bonds with serine. [14,15,16]

Like the previous analysis, dipole moment (μ) and lipophilicity were evaluated using logP and logμ as key physicochemical descriptors (Figure 1B). The results indicate that active salicylates typically exhibit logP values between 1.0–3.5 and dipole moments (μ) ranging from 0.80–1.4 Debye. In contrast, fenamates display a distinct profile characterized by higher lipophilicity (logP 3.3–4.8) and a narrow dipole moment range (1.0–1.2 Debye).

Notably, the COX-1 inhibitory activity of both drug classes follows a positive parabolic correlation with dipole moment, with optimal compounds (primarily fenamates) clustering in the 1.05–1.20 Debye range – consistent with the logP-logμ diagram. Furthermore, the dipole moment vector is oriented from the accessory ring moiety toward the nitrogen atom of the anthranilic acid group (Figure S2 and S3).

2.2. Figures, Tables, and Schemes

All figures and tables should be cited in the main text as Figure 1, Table 1, etc.

These two drug classes exhibit a negative parabolic correlation with LUMO energy, where specific fenamates (FFA and MCF) demonstrate the highest activity (Figure S4). Conversely, their correlation with TPSA follows a negative linear relationship, with fenamates (MFA, FFA, MCF, NFA) again showing superior activity (Figure S5). Based on these descriptors, we propose the following model (Figure 2):

$${pIC}_{50}=-6.832{log\mu }^2+14.065log\mu +492.539E_{LUMO}^2-116.515E_{LUMO}-10.152logTPSA+20.599$$

n = 8, F = 2.88, Q² = 0.8861

The model indicates that fenamates are more potent inhibitors than salicylates. Fenamates exhibit intermediate dipole moments, low polar surface distribution, and electrophilic character in the carboxylic acid-containing ring. The electron-withdrawing carboxyl group delocalizes a positive charge across this aromatic system.

Most analyzed drugs (DF, MCF, MFA, and NFA – the most active compounds) form ionic interactions between the carboxylate group and Arg120. In contrast, the less active compounds (SSZ, ASA, SA, and FFA) lack this Arg120 interaction. For both fenamates and salicylates, the carboxylate-containing ring interacts with nonpolar amino acids (Leu352, Val523, Gly526, Ala527, etc.), facilitating negative charge transfer to the electrophilic ring, a feature linked to the LUMO orbital.

Therefore, it is interesting that for both enzymes, the substitutions at position 2 with respect to the carboxylate group are essential for the recognition and stabilization of free radicals due to the captodative properties of salicylates and fenamates; the latter being more active since the accessory moieties in this position better stabilize the charge by possessing, and additional aromatic system.

2.2. Acetates

For the acetate-class drugs, correlations were identified between physicochemical properties and binding to both COX-1 and COX-2, consistent with their non-selective inhibition profile observed experimentally and in silico (Figure S6). Regarding COX-1, significant correlations emerged with: partition coefficient (logP), showing a positive parabolic relationship (Figure S7); LUMO energy (E_LUMO), showing a negative parabolic relationship (Figure S8); and molecular volume (V), showing a negative linear relationship (Figure S9). These correlations enabled the development of a QSAR model incorporating all three physicochemical properties influencing COX-1 inhibitory activity (Figure 3A):

$$\mathrm{A}.{pIC}_{50}=0.271{logP}^2-1.6557logP-1678.2937E_{LUMO}^2+192.6651E_{LUMO}-9.2588logV+25.5944$$

n = 7, F = 1.774, Q² = 0.9009

$$\mathrm{B}.{pIC}_{50}=0.54913{logP}^2-3.00186logP-19.03265logMR+14.08181logOv+24.82122$$

n = 7, F = 224.5, Q² = 0.9889

According to the model, the partition coefficient significantly influences inhibitory activity. The exponential component indicates an upward-opening parabola with broad curvature, suggesting low system sensitivity. This is consistent with the negative linear coefficient, which positions the parabola in the second quadrant (logP range: 1.6–4.2) – a range compliant with Lipinski's Rule of Five. [17,18] This positive correlation with biological activity reflects interactions with nonpolar amino acids (Leu, Val, and Ala). Notably, diclofenac (DIC), one of the most active compounds in this group, contains a dichlorinated phenyl ring as its accessory moiety, which enhances lipophilicity.

The LUMO orbital energy also shows a parabolic relationship, but with a downward-opening curve. The most active compounds (DIC, KT, and TOL) feature electron-withdrawing groups on their phenyl rings that delocalize partial positive charge, positioning the LUMO orbital at these sites. This facilitates molecular recognition with electron-donating residues (Leu352 and Trp387) (Figure 4A).

The size descriptor analysis reveals that smaller molecular volumes increase drug-receptor recognition. This trend is observed for DIC, KT, and TOL, which are significantly more compact than indoleacetic acid derivatives.

With the acetates, a positive parabolic relationship was observed between the partition coefficient (logP) and COX-1 inhibitory activity, indicating that lipophilicity is a critical determinant of both target recognition and inhibitory potency. The analysis demonstrates two optimal lipophilicity ranges (logP 1.5-2.0 and 3.5-4.0) that maximize inhibitory activity for key acetate-class drugs (ketorolac, tolmetin, diclofenac, and indomethacin), with intermediate values showing reduced efficacy.

In the same vein, correlations were also found with other physicochemical properties for the acetate family against COX-2. In this case, they were with the partition coefficient (logP), the molar refractivity (MR), and the ovality (Ov). A positive parabolic relationship was observed with the partition coefficient (Figure S10), while negative linear correlations were found with molar refractivity and ovality (Figures S11 and S12). From these correlations, a QSAR model was obtained that involves three descriptors that influence the enzyme (COX-2) inhibitory activity (Figure 3B).

The QSAR model reveals that the partition coefficient (logP) significantly influences inhibitory activity, with a tighter parabolic relationship that for COX-1. This suggests enhanced enzyme recognition, where optimal lipophilicity positively correlates with biological activity through interactions with nonpolar residues (Leu, Val, and Ala). Ketorolac (KT), the most potent compound in this series, demonstrates this principle through its phenyl accessory moiety, which forms closer interactions with Leu352 (5.05 Å) and Trp387 (5.10 Å) in COX-2 versus COX-1 (5.43 Å and 5.21 Å, respectively).

Using molar refractivity as a size criterion, it indicates that less bulky drugs are favoured for recognition by the enzyme; once again, KT, TOL, and DIC are the smallest compared with the indoleacetic acid derivatives.

Finally, the molecular shape also influences biological activity, as the compounds should ideally present an ovoid shape, as in KT, TOL, and DIC. However, when ovality is significantly higher, activity decreases due to reduced affinity, as observed for ETD, SUL, OX, and IND (Figure S13).

2.3. Propionates

For the propionate group, correlations were found with specific physicochemical properties, both for drugs binding to COX-1 and to COX-2, which are not selective for either isoform, as shown both experimentally and in molecular docking studies (Figure S14). For COX-1, the identified correlations were with the acidity constant (pK_a) and molecular volume (V); the pK_a correlation was negative and parabolic (Figure S15), whereas the molecular volume correlation positive and linear (Figure S16). From these correlations, a QSAR model was constructed using two physicochemical properties that influence the enzyme (COX-1) inhibitory activity (Figure 5A).

$$\mathrm{A}.{pIC}_{50}=-44.912{pK}_a^2+426.223{pK}_a-7.201logV-987.923$$

n = 5, F = 1.23, Q² = 0.7802

$$\mathrm{B}.{pIC}_{50}=19.97{pK}_a^2-185.28{pK}_a-300.15E_{HOMO}+384.48$$

n = 5, F = 1.028, Q² = 0.7323

According to the QSAR model, pK_a has a positive influence on biological activity; in general terms, more acidic pK_a values are favoured. Nevertheless, the drugs with the highest inhibitory activity against COX-1 correspond to KP, FEN, and FP (Figure S15). Docking analysis shows that the carboxylate group of KP, FEN, and FP establishes a unique ion–ion interaction with Arg120; whereas in the cases of IBU and NAP, which are the least active, the negative charge of the carboxylate is distributed between Arg120 through an ion–ion interaction and Tyr355 through an ion–dipole interaction (Figure 4B).

The QSAR analysis demonstrates that molecular volume significantly impacts the inhibitory activity of propionate-class drugs against COX-1, with larger molecules exhibiting enhanced potency (Figure S16). This observation is supported by molecular docking studies, which reveal that high-activity compounds (KP, FEN, and FP) possess two aromatic rings that establish extensive hydrophobic interactions with nonpolar residues (Val349, Leu352, Tyr385, Ile523, and Ala527). In contrast, lower-activity compounds such as IBU and NAP exhibit reduced binding efficiency due to structural limitations: IBU contains only a single aromatic ring, while NAP features a fused ring system with restricted conformational flexibility, both resulting in fewer optimal contacts with key hydrophobic residues (Figure 4B).

For COX-2 inhibition, distinct structure-activity relationships emerged. A positive parabolic correlation was observed between biological activity and acidity constant (pK_a), indicating an optimal pK_a range for maximal inhibition (Figure S17). Additionally, a linear correlation with HOMO energy (E_HOMO) was identified, suggesting that electronic properties further modulate inhibitory potency (Figure S18). These relationships were integrated into a robust QSAR model that accounts for the combined influence of these descriptors on COX-2 inhibition (Figure 5B).

This COX-2 QSAR model reveals a compelling relationship between inhibitory activity and two interconnected physicochemical properties: pK_a and HOMO energy (E_HOMO). The observed trends are mechanistically significant: (1) lower pK_a values reflect greater carboxylate ionization at physiological pH, enhancing ionic interactions, and (2) more negative E_HOMO values (Figure S18) indicate stronger nucleophilic character, as the HOMO is localized on the carboxylate group. Consequently, compounds with both low pK_a and negative E_HOMO (FEN, KP, FP) exhibit superior activity compared to less active analogues (IBU, NAP).

Structurally, all drugs interact with Arg120 via ion-ion interactions and with Tyr355 via ion-dipole contacts with the carboxylate. The ionic interaction is governed by carboxylate acidity, while the concurrent nucleophile (HOMO)-electrophile (LUMO of protonated Arg120) interaction further stabilizes binding (Figure 4B). Activity differences within the series arise partly from electronegative substitutions at the 3- and 4-positions of the propionate-bearing ring (Table 1), which fine-tune electronic distribution and thus modulate carboxylate reactivity.

2.4. Pharmacophore Structure

The integrated QSAR and molecular docking studies highlight the critical role of the ionic interaction between Arg120 and the carboxylate group of NSAIDs for initial guidance and recognition by both COX-1 and COX-2, particularly in high-activity compounds. Furthermore, the aryl group containing the carboxylic acid interacts with non-polar amino acids (Table 2), forming a complementary hydrophobic pocket that stabilizes the drug-enzyme complex. These dual ionic interactions at Arg120 and hydrophobic interactions with nonpolar residues are relevant for binding affinity and inhibitory potency. Similarly, the arachidonic acid is oriented from the carboxylic acid to interact with Arg120 and Tyr355 aminoacids. In contrast, the Val349, Ala527, Ser353, Ile523 (in case COX-1) and Val523 (in case COX-2), Ile523 (in case COX-1) interact between 1 and 5 position of arachidonic acid (Figures S19 and S20), which are the same aminoacids that interact with accessory moieties of the NSAIDs, where the most critical aminoacid is Tyr385 that initiates the oxidation via free radicals (Table 2). [2]

The moiety responsible for the anchoring and orientation of NSAIDs on both COXs consists of an aryl group (aromatic generally) which is binding to a carboxylic acid with 0, 1, or 2 carbon atoms of distance of separation (Figure 6). The accessory moieties are characteristic for each drug group, the substitutions in position 2 correspond to the salicylates and fenamates that show relevant captodative properties; while the substitutions in the position 3 is characteristic of the acetates and some propionates; also the substitutions in the position 4 are similar between the acetates, propionates and salicylates (sulfasalazine) drugs (Figure 7), such drugs show decrease in the biological activity on both enzymes.

The core pharmacophore of NSAIDs consists of a carboxylic acid group, connected either directly or via a single carbon bridge, to a 5- or 6-membered aryl ring (Figure 7). This structural arrangement is essential for biological activity, where the carboxylate forms a key ionic interaction with Arg120. At the same time the aromatic ring provides structural rigidity and electron density required for molecular interactions. The accessory moieties of these molecules interact with various amino acids in the COX active site. In position 2, the accessory moieties show plane geometry, while the X-groups are bioisosters between them; such interactions are more specific, particularly aromatic-aromatic contacts with Phe381, Phe518, Tyr385, and Trp387 residues. In positions 3 and 4, the Y and Z substituents establish hydrophobic interactions with nonpolar residues such as Val349, Leu352, and Ala527, which contribute to stabilizing the enzyme-drug complex. This pattern is especially relevant for indoleacetic acid derivatives, where the indole system engages in π-stacking interactions with residues like Tyr385 (Table 2). Therefore, the description in this manuscript show the molecular moieties in the anchoring and orientation with COX-1 and COX-2 enzymes, thus the proposed pharmacophore is applicative for explain the ligand-receptor recognition and the non-selective NSAIDs rational design; also is have established the molecular bases for the search the moieties and stereochemistry necessaries for confer the selectivity on these enzymes, particularly on the COX-2, whose key moiety should be related to the aminoacids that configure the guide molecular portion.

3. Materials and Methods

3.1. Bibliographic Search

By means of an exhaustive bibliographic search of articles published from 1990 to 2023, the median inhibitory concentration (IC₅₀) of COX-1 and COX-2 enzymes was found for each of the Non-steroidal anti-inflammatory drugs (NSAIDs) [21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118] most widely prescribed for the pain, inflammatory and antipyretic diseases, which are the first cause of consulting at the first level of care in the world. The databases consulted were IUPHAR Database, Taylor & Francis, Scopus, PubMed, Google Scholar, SciELO, NCBI, Elsevier, UpToDate, Springer Nature, and Wolters Kluwer. The NSAIDs groups most frequently administered for pain and anti-inflammatory diseases are Salicylates (salicylic acid, acetylsalicylic acid, diflunisal, sulfasalazine), Fenamates (mefenamic acid, flufenamic acid, meclofenamic acid, niflumic acid), Acetates (diclofenac, ketorolac, indomethacin, oxamethacin, sulindac, etodolac, tolmentin), Propionates (ibuprofen, ketoprofen, feniprofen, flurbiprofen, naproxene). [119]

3.2. QSAR Analysis

All the IC₅₀ values obtained from the literature were averaged to express the activity of NSAIDs over COX-1 and COX-2 enzymes. These IC₅₀ values were converted to molar concentration (M) to evaluate the structure-activity relationship of each group. The pIC₅₀ values of each drug were correlated with its respective partition coefficient (logP, representing liposolubility) and molar refractivity (MR, indicating polarity). These two physicochemical properties were selected with the ACD/ChemSketch and CS ChemDraw Pro v.6 software. [120] The experimental pK_a values (an electronic criterion) were taken from the PubChem platform. [121]

Briefly, each NSAID was built with the Gaussian 16 and GaussView 6 program [122,123], and the conformational analysis search was carried out with molecular mechanics (with the Merck molecular force field, or MMFF) after establishing the most relevant torsions for the molecule. Once the conformer set was defined, the most stable conformer was optimized at the DFT level with the B3LYP 6-31G(d,p) basis set. The molecular volumes, HOMO and LUMO energies, and dipole moment were taken from the output file, and then the dipole moment vector and the HOMO and LUMO orbitals were visualized. [124]

The partition coefficient is an important physicochemical property that corresponds to the ratio of the concentrations of a substance in two immiscible phases in equilibrium in a mixture. Usually, the solvents used are water and n-octanol. This property often serves as a criterion to predict the ease of simple diffusion of the analyte across membranes, particularly the blood-brain barrier [125]. The partition coefficient is estimated according to the following equation:

$$P={\displaystyle \frac{{\left[A\right]}_{oct}}{{\left[A\right]}_w}}$$

where P is the partition coefficient, [A]_oct is the analyte concentration in the n-octanol phase, and [A]_wis the analyte concentration in the aqueous phase. The biological activity is relevant when the logP values are greater than 1.0, indicating that the analytes can probably crossmembranes and interact with their molecular targets.

MR, a measure of the total polarization of one mole of a substance, is related to molecular volume because small molecules are more likely to cross membranes and produce biological activity. This parameter is calculated with the Lorentz-Lorenz equation [126]:

$$MR=\left({\displaystyle \frac{MW}{\rho }}\right)\left({\displaystyle \frac{n^2-1}{n^2+2}}\right)$$

where MR is the molar refractivity, MW is the molecular weight, $\rho$ is the density, and n is the refractive index of the molecule in question.

The values of the Topological Polar Superficial Area (TPSA) and molecular ovality were determined using the Molinspiration Platform [127] and Spartan ’14 software respectively, which are topological descriptors related like a measure of the polarity-area and the deviation from sphericity of a molecule. [128]

The apparent partition coefficient is called the distribution coefficient (log D). It is related to pK_a values and is determined by the pH of the system. In the present study, a physiological pH (7.4) was assumed, and log D was calculated from the following equation [129]:

$$\log D=\log P_n-\log (1+{10}^{\mathrm{p}H-p{K}_a})$$

where D is the distribution coefficient, P_n is the partition coefficient for the ionizable form of the analytes in the organic phase, pH is the logarithmic form of the inverse of the H⁺ concentration, and pK_a is the logarithmic form of the inverse of the acidity constant. The log D is the most accurate method for evaluating drug diffusion across membranes, accounting for the medium pH and the ionizable and non-ionizable forms.

The polarity of a molecule helps explain several physicochemical properties, such as chromatographic retention on a polar stationary phase [130]. Among the many descriptors proposed to quantify the effects of polarity, the most obvious and most used parameter is the dipole moment of the molecule. This parameter relates to the number of donors and acceptors in the hydrogen bond. Biological activity is higher for molecules with fewer donors and acceptors.

HOMO (the highest occupied molecular orbital) and LUMO (the lowest unoccupied molecular orbital) are quantum chemical descriptors that play a key role in the control of chemical reactions and the determination of electronic band gaps in solids. In the frontier molecular orbital (FMO) theory of chemical reactivity, the formation of a transition state is due to an interaction between the border orbitals (HOMO and LUMO) of the reactive species. [131,132]

While the HOMO energy is directly linked to the ionization potential and characterizes the susceptibility of the molecule to an attack by electrophiles, the LUMO energy is directly connected to electronic affinity. It characterizes the susceptibility of the molecule to an attack by nucleophiles. According to FMO, the electronic density of the frontier orbitals provides essential information for a detailed characterization of donor-acceptor interactions [123], based on the idea that most chemical reactions take place at the position and in the orientation in which the overlap of HOMO and LUMO of the respective reagents reaches a maximum [131]. Whereas the atomic HOMO density is critical for charge transfer in a donor molecule, the atomic coefficients of LUMO are essential for an acceptor molecule [132]. HOMO-LUMO interactions are relevant for understanding ligand-receptor recognition in biological systems.

Correlations were established between the biological activity of the drugs on COX-1 and COX-2 and their physicochemical and molecular properties mentioned above. Correlations with r values > 0.75 were considered significant.

3.3. Molecular Docking

Preparation of the protein (receptor): Searching and obtaining the three-dimensional structure of the protein in the Protein Data Bank platform, the following crystals were selected: for COX-1 (PDB 3N8Z) and COX-2 (PDB 3PGH). In Autodock Tools 1.5.6, water molecules and other substances were removed from the protein structure with which it was crystallized, leaving only the binding site of the endogenous ligand.

Ligand preparation: The structures of the compounds were drawn in their ionized forms, and the ligand was sent to GaussView 6.0. 16 and Gaussian 16, where the optimization was carried out at a DFT level.

Validation: The RMSD (Root Mean Square Distance) of both docked ligands is within the reliable range of 2 Å, with a value for both of 0.00 Å, verifying that flurbiprofen can interact with the 3N8Z and 3PGH crystal structures in a similarly manner to the pre-existing co-crystallized flurbiprofen (Figures S21 and S22).

Grid generation and docking algorithm: It was decided that targeted molecular docking would be carried out, so a three-dimensional grid was defined with the following dimensions: 66 x 66 x 66, 0.375 Å spacing, center Grid Box: X center: 8.833, y center: -5.250, Z center: 24. 833, for COX-1, while for COX-2: They were used with the following dimensions: 66 x 66 x 66, 0.375 Å spacing, center Grid Box: X center: -19.500, y center: -9.333, Z center: -21.528, which formed the grid around the receptor binding site, where ligand complementarity was assessed. These grids guided the search for the ligand within the search space. The Lamarckian algorithm was applied to explore the ligand’s possible conformations at the binding site.

3.5. Statistical Analysis

Linear regression analysis was performed using the least squares method and one-way ANOVA. The values of the ordinate, slope, and linear correlation coefficient were evaluated using Student's t-test. [133]

Second-order polynomial regression was performed using the least-squares method and one-way ANOVA. The values of the ordinate, slope, and correlation coefficient were evaluated using Student's t-test. Calculations were performed from the equations derived from the second-order polynomial regression model for both positive and negative parabolic curves. [134,135,136]

Data were analyzed using least squares regression and one-way ANOVA. The slopes were adjusted to the origin, and the correlation coefficients were evaluated using Student's t-test, with p < 0.05 considered significant in all cases. The QSAR models were examined by multiple regression analysis using the determinant method and one-way ANOVA. Statistical analyses were performed using Sigma Plot 15.0 [137] and the R framework. The external predictive ability of the QSAR models was assessed by the squared predictive correlation coefficient (Q²). [138]

4. Conclusions

Through multiple QSAR and molecular docking analyses of COX-1 and COX-2, the pharmacophoric features that explain the recognition of non-selective COX inhibitors (NSAIDs) were identified. The pharmacophoric structure show two regions: the guide moiety (carboxylate) that directs the binding mode of the NSAIDs and arachidonic acid to the Arg120 through electrostatic interaction; and the accessory moieties essential to blocking substrate oxidation via free radicals for both enzymes that catalyze the same reaction, the salicylate and fenamate families show a relevant captodative effect by electronic interaction between the carboxylate and an electron-donating group in position 2 on the aromatic ring. The general physicochemical and molecular properties that determine such ligand–enzyme recognition correspond to topological descriptors (dipole moment and topological polar surface area), geometric descriptors (volume, molar refractivity, and ovality), thermodynamic descriptors (partition coefficient), and electronic descriptors (acidity constant, HOMO and LUMO energies). Therefore, the elucidation and identification of the pharmacophore of non-selective COX inhibitors opens the door to a new approach for designing drugs with inhibitory activity on these enzymes and lays the foundation for a more detailed understanding of the physicochemical and molecular properties, as well as the intermolecular interactions, that may be involved in isoform selectivity.