QSAR Modeling and Biological Testing of Some 15‐LOX Inhibitors in a series of Homo‐ and Heterocyclic Compounds

1. Introduction

Lipoxygenases (EC 1.13.11.12, LOX) are oxidoreductases with iron or manganese as a cofactor and are the most important enzymes in biological systems. They are found in mammals, plants, fish, mosses, bacteria, yeasts, corals, algae, and fungi [1,2,3,4,5,6,7,8,9,10,11,12,13]. LOXs catalyze the oxidation of free and esterified polyunsaturated fatty acids (PUFAs) containing one or more (1Z, 4Z)-penta-1,4-diene systems into hydroperoxides [1,14], which are then metabolized into various signaling compounds such as leukotrienes and lipoxins in animals [15,16], prostaglandin-like molecules in corals [8], volatile substances of green leaves [17,18,19], jasmonic acids in plants [20,21], and lactones in microorganisms [22]. The LOX superfamily is classified based on their regiospecificity into 5-, 8-, 9-, 10-, 12-, 13-, 15-, fusion-, mini- and Mn-LOX. The specificity of LOXs affects the number of the carbon atom in the PUFA molecule which is attacked by an oxygen molecule followed by the formation of hydroperoxides [2]. Subsequently, the hydroperoxides formed are converted into ketones, aldehydes and alcohols by other enzymes [3]. Of the 10 LOX types mentioned above, 5-, 8-, 9-, 10-, 12-, 13-, 15- LOXs are the so-called classic LOXs. Among these classical LOXs, 9- and 13-LOXs are the most important enzymes in plants while 5-, 12- and 15-LOX are predominantly found in animals [2,23]. The catalytic center of all the six families of classical LOXs (5-, 8-, 9-, 10-, 12-, 13-, and 15-LOX) contains "non-heme" iron Fe(III) in activated state [1]. The mechanism of PUFA oxidation under the action of LOXs has been extensively studied and described in the scientific literature, e.g. in [24,25,26].

The relevance of the search for 15-LOX inhibitors is due to the pathophysiological effect of the products of oxidative metabolism of PUFAs under the action of this enzyme on the organism of animals and humans. Thus, arachidonic acid metabolites formed with participation of 15-LOX in various types of cells and organs of animals and humans are involved in the development of many diseases, including atherosclerosis, hypertension, diabetes, obesity and neurodegenerative disorders [27,28,29]. In addition, 15 LOX inhibitors have been shown to be more efficient than COX inhibitors in suppressing the growth of epithelial cancer cells, regardless of the expression status of each enzyme [28,30].

To date, several types of inhibitors of the catalytic activity of LOX are known: antioxidants and free radical acceptors (a); chelating agents (b); non-competitive redox inhibitors (c); and combinations of chelators and reducing agents (d) [31]. All of these are used in biochemistry and in the food industry to regulate the PUFA metabolism under the action of LOX isoforms. However, most of these bioactive agents are 5-LOX inhibitors. Compared to 5-LOX inhibitors, the 15-LOX inhibitors are covered much less in scientific literature. Therefore, the search for 15-LOX inhibitors is an urgent task for the medical and pharmaceutical chemistry. It should be noted that soybean 15-LOX (15-sLOX) has been used as a model for 5-hLOX due to their similarity in structure and mechanism of action [32,33].

Currently, QSAR/QSPR methods are actively used in the development of lead compounds and drugs based on them. These methods make it possible to solve a wide range of important tasks: 1) unbiased search for potential bioactive compounds in virtual librari es; 2) extension of the scope of application of active components of known drugs; 3) evaluation of the side effects and toxicity of potential drugs and new bioactive compounds; 4) molecular design of new potentially low-toxicity bioactive compounds based on the structures of known drugs and hit compounds.

QSAR/QSPR methods allow for modeling of biological activity and physicochemical properties based on a relatively small number of structures in the training sets, which is their major advantage. The concept of these methods is based on the fundamental postulate of organic chemistry that states that the properties of chemical compounds are determined by their structures. The essence of QSAR/QSPR methods lies in the correct selection of various physicochemical and structural descriptors which allow an unbiased description of each structure of the training sets, as well as mathematical and statistical methods for building regression relationships based on the most informative descriptors. The application of QSAR approaches at the preclinical phase can significantly reduce the time and material costs in the targeted development of new potential drugs, as various researchers repeatedly noted [34,35,36,37,38,39,40,41,42,43,44,45,46]. The use of QSAR/QSPR methods to solve the above problems is widely reported in scientific literature, including research articles, scientific reviews and monographs. Every year, the well-known software packages implementing these methods are improved and new ones are developed [47,48,49,50,51,52,53,54,55,56].

Among the many QSAR/QSPR methods whose classification is based on the choice of certain types of descriptors and machine learning methods for constructing mathematical equations [34,35,36,37,38,39,40,41,42,43,44,45,46,57], 2D- and 3D-QSAR methods are the most popular. The calculation of descriptors in these methods is based on structural formulae and molecular diagrams of chemical compounds (2D-QSAR) and on potentially bioactive conformations of chemical structures (3D-QSAR) [34,36,58,59]. The choice of one of these methods depends mainly on the objectives of a particular study. 2D-QSAR methods, together with 3D to 6D-QSAR methods, can be used to model biological activity where there is a need to find compounds that can either increase or decrease the enzymatic specificity of certain proteins. The lack of reliable crystallographic data on potentially biologically active conformations of organic compounds in the active centres of proteins, such as in the case of 15-LOX inhibitors, increases the demand for 2D-QSAR methods over all other methods for modeling enzyme specificity. It should be understood that the application of 2D-QSAR methods is not limited to modeling biological activity. These methods can be successfully applied to model the physicochemical properties of organic compounds, as noted in our previous work and in the work of others.

GUSAR 2019 (General Unrestricted Structure Activity Relationships) is among the programs that allow the calculation of physicochemical and structural descriptors with subsequent selection of the most significant ones and construction of consensus QSAR/QSPR models based on them. Previous versions of this program are known such as GUSAR 2013 and GUSAR 2011 [47,48,53,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74]. The demo version of the latter can be found at way2drug.com. Regardless of the version, this program has proved its worth in modeling different types of biological activity and a number of physicochemical properties (lipophilicity, biotransformation factor, antioxidant activity) in a number of heterogeneous organic compounds, as reported both by its developers and by us in previous works [60,61,62,63,67,68,69,70,71,72,73,74].

The purpose of the present work was to study the quantitative structure–activity relationship of 15-LOX inhibitors in a series of homo- and heterocyclic compounds with common structural formulae I–XVI (Figure 1) [75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90] using the GUSAR 2019 program and to build statistically significant valid prediction models for the pIC₅₀ parameter designed to search for new potentially efficient 15-LOX inhibitors in virtual libraries and databases.

2. Results and Discussion

2.1. Prediction of the Numerical Values of the IC₅₀ Parameter Using the GUSAR 2019 Program

As a result of QSAR modeling based on the consensus approach implemented in the GUSAR 2019 software, eighteen consensus models M1–M18 were generated. All of these models are designed to predict the numerical values of the pIC₅₀ parameter for LOX inhibitors. The difference between these models lies in the choice of different types of descriptors and the number of partial regression relationships constructed from them. The descriptive power characteristics of these consensus models are shown in Table 1. They were automatically calculated in GUSAR 2019 software based on the comparison of experimental values of the pIC₅₀ parameter for LOX inhibitors and its values predicted by these models. In our previous works [60,61,62,63,67,68,69,70,71,72,73,74], we explained that the coefficients of determination R², Q²_LMO, values of standard deviation SD and F criterion presented in Table 1 are averages calculated taking into account all partial regression models included in each of the Mi consensus models.

The data presented in Table 1 allow us to $\overline{R_{TrSi}^2}$ conclude that all the QSAR consensus models M1–M18 constructed by us feature acceptable stability, since the A parameter for them is smaller than the value of 0.3 allowed in scientific literature [44]. Consensus models M1–M3 and M10–M12 have the highest stability, as they provide small A values (A < 0.1). The SCR method was used to select the descriptors for their construction. The consensus models M7–M9 and M16–M18, which were constructed using the Both method for descriptor selection, show an acceptable level of stability (A ⩽ 0.2). The consensus models M5–M6 and M13–M15 have the lowest stability, with the exception of M4. The selection of descriptors in the construction of these models was performed using the RBF-SCR method. Thus, the data in Table 1 allow us to conclude that, in all cases, the nature of the descriptors on which the consensus models M1–M18 were built did not play a decisive role in their stability. At the same time, the method of descriptor selection was a significant factor that affected the stability of these models.

It should be noted that the statistical characteristics (average R² ($\overline{R^2}$), average Q² ($\overline{Q^2}$), average F ($\overline{F}$), average Q²_LMO($\overline{Q_{LMO}^2}$), average SD ($\overline{SD}$)) whose calculation is provided by GUSAR 2019 do not allow a detailed assessment of the descriptive ability of QSAR models, while assessment of the predictive ability of QSAR models is not performed by this program at all. Therefore, metrics based on different types of R² coefficients of determination (R², R²₀, average R²_m ($\overline{R_m^2}$), Q²_F1, Q²_F2, average R²_mTSi ($\overline{{R_m^2}_{TSi}}$), CCC) were additionally used to objectively assess the descriptive and predictive ability of the M1–M18 QSAR models. In addition, metrics designed to evaluate the prediction errors of pIC₅₀ values (RMSE, MAE, SD) were used to determine the true prediction quality index of the pIC₅₀ parameter based on the M1–M18 models for the test sets TS1–TS2 of compounds [44,45,46,57]. The calculation of all these criteria was performed using the Xternal Validation Plus 1.2 software [91]. The formulae used in this program to calculate all the above criteria are given in Table S1 in the Supplementary Material. The same software was used to check the models for systematic errors.

The process of internal and external validation of the developed QSAR models M1–M18 was based on the structures of training sets TrS1–TrS2 and external test set TS1, respectively. In addition, the predictive ability of the M10–M18 models was evaluated using the structures of the LOX inhibitors contained in the internal test set TS2. In addition, the predictive ability of the M3, M6, M9, M12, M15 and M18 models was evaluated by comparing the experimentally determined (see Section 2.2) and calculated values of the pIC₅₀ parameter for the structures of LOX inhibitors 1–8 contained in the test set TS3.

Table 2 and Table 3 present the threshold (maximum and minimum) values of different types of coefficients of determination and prediction errors of the pIC₅₀ parameter for some of the QSAR models we built. These were calculated using Xternal Validation Plus 1.2 software [91] for 100% and 95% structures of the 15-LOX inhibitors contained in the training sets TrS1–TrS2 and test sets TS1–TS2, respectively. Those QSAR models for which numerical data are not available in Table 2 and Table 3 occupied an intermediate position between the maximum and minimum values for all statistical features listed in this table. The data in Table 2 and Table 3 allow an objective assessment of the range of variability of different statistical criteria characterising the descriptive and predictive ability of the M1–M18 QSAR models. Tables S2–S6 in the Supplementary Material present the full set of criteria calculated using Xternal Validation Plus 1.2 software for the TrS1–TrS2 and TS1–TS2 sets, considering both 100% and 95% of the 15-LOX inhibitor structures contained therein.

The data from Table2 allow us to conclude that all the models M1–M18 showed high predictive power.

Figure 2 how correlations between the experimental and calculated values of the pIC50 parameter for the 15-LOX inhibitors contained in the TrS1 set as an example of the high predictive power of our QSAR models M1–M9.

In this case, according to this Table2, the best descriptive ability was provided by the M4, M13, and M15 QSAR models built using RBF-SCR as the descriptor selection method. The lowest accuracy in reproducing the experimental data contained in the TrS1–TrS2 sets was demonstrated by the M2 and M10 models, in which the descriptors were selected using the SCR method. In particular, the M2 model had the smallest numerical values of different types of coefficients of determination and the largest values of prediction errors of the pIC₅₀ parameter for the structures contained in the TrS1 set. At the same time, the M4 model showed the largest numerical values of different types of coefficients of determination and the smallest error values. This conclusion holds true for both 100% and 95% of the data contained in the TrS1 set.

The conclusions drawn from the analysis of the goodness of fit of the M10–M18 models built on the TrS2 set are less clear. For example, detailed analysis of different types of coefficients of determination and prediction errors shows small differences in the estimation of the descriptive ability of the M10–M18 models depending on the completeness of inclusion of the experimental data (100% and 95% of the data) contained in the TrS2 set. In fact, if 100% of the data in TrS2 are considered, the maximum coefficients of determination R²₀, R²’₀, $\overline{R_m^2}$ , CCC and the minimum ⊗R²_m coefficient are provided by the M13 model. The same model, based on the data in the tables, is characterised by the minimum values of error and standard deviation in the prediction of the pIC₅₀ parameter for the structures of the 15-LOX inhibitors included in the TrS2 set. The maximum value of the R² coefficient is given by the M15 model.

If 95% data of the training set TrS2 is taken into account, the M13 model showed the best descriptive power even on the basis of a smaller set of criteria. In fact, it was characterised by the largest values of coefficients of determination such as R², R²₀, and CCC, along with the minimum value of the ΔR²_m coefficient. The M15 model gave the largest values of the coefficients, R²’₀, mean R²_m in combination with the smallest values of the RMSE, B parameter and standard deviation SD. At the same time, the minimum value of the MAE error is observed for the M14 model based on 95% data of the TrS2 training set.

The data of Table 3 and Tables S4-S6 in the Supplementary Material allow us to conclude that almost all the QSAR models built by us are characterised by moderate predictive ability and absence of systematic errors in predicting the target property for the structures of 15-LOX inhibitors contained in test sets TS1–TS2. The highest numerical values of different types of coefficients of determination in predicting the pIC₅₀ parameter for 100% of the 15-LOX inhibitor structures contained in TS1 were demonstrated by the M14 model. The same model had the smallest value of the RMSE error. The smallest values of MAE error, SD value and B parameter for 100% of the data in the TS1 were provided by the M13 and M17 models, respectively.

The worst results for all the statistical criteria described above were observed for the M9 and M12 models. For example, the M9 model was characterised by the minimum values of the R², R²₀, $\overline{R_m^2}$ and CCC criteria. The M12 model gave the minimum values of the R², R²₀,$\overline{R_m^2}$ , and CCC criteria, while the M12 model provided the minimum values of the R², R²₀, Q²_F1, and Q²_F2 criteria, combined with the maximum value of the RMSEP error. However, the M18 model was characterised by the maximum value of the ⊗R²_m coefficient. The highest value of the MAE error in predicting the pIC₅₀ parameter for 100% of the TS1 structures was demonstrated by the M3 model. The M5 model was gave the maximum values of the SD and B parameters.

The results of evaluating the predictive ability of QSAR models M10–M18 using 100% data of the test set TS2 were more unambiguous. The maximum values of the different types of determination coefficients along with the minimum values of the prediction errors of the pIC₅₀ parameter in this case belonged to the M13 model. The minimum values of the determination coefficients in combination with the maximum values of the prediction errors of the pIC₅₀ parameter for the 15-LOX inhibitors corresponded to the M11 model.

Removal of 5% of the data with the worst prediction results for pIC₅₀ from TS1 and TS2 contributed to a slight increase in the numerical values of different types of coefficients of determination and a decrease in the numerical values of the prediction errors for this parameter. However, the changes in all the statistical criteria were not systematic when the data from both sets were manipulated in this way. In particular, after removing 5% of data from TS1, the highest numerical values of the coefficients of determination R², R²₀, Q²_F1, and Q²_F2 and the $\overline{R_m^2}$ were obtained from the M10 model. The maximum value of the R²'₀ criterion in comparison of the experimental pIC₅₀ values with the predicted ones was observed for the M1 model. The minimum values of RMSE and MAE errors based on 95% of the data of the TS1 set were shown by the M13 and M16 models, respectively. The lowest values of SD and B were shown by the M4 model.

Removing 5% of the data from the TS2 set had almost no effect on the generalizing conclusions. In this case, like in the case of 100% data in the TS2, the M13 model provided the largest numerical values of different types of coefficients of determination and the smallest values of pIC₅₀ prediction errors for the 15-LOX inhibitors, with the exception of SD whose numerical value was the smallest for the M17 model. The minimum numerical values of different types of coefficients of determination combined with the maximum values of the pIC₅₀ prediction errors for both 100% and 95% data in the TS2 were observed for the M11 model.

Based on analysis of the numerical values of different validation criteria presented in Table 1, Table 2 and Table 3 of this section and in Tables S2-S6 in the Supplementary Material, we can conclude that almost all the models showed high descriptive and moderate predictive power, since they met the internal and external validation criteria described in Section 2.3. It should be noted that this condition was met for both 100% and 95% of the data contained in the TrS1-TrS2 and TS1-TS2 sets.

2.2. Experimental Determination of the IC₅₀ Parameter Against 15-LOX for Compounds 1–8

The results of the in vitro analysis of the inhibitory activity of a series of 2H-(benzo)pyran-2-one derivatives 1–8 against 15-LOX are presented in Table 4 and Figure 3.

IC₅₀ values were determined by linear interpolation between the points closest to 50% inhibition (Figure 3).

The inhibitory activity against 15-LOX was experimentally determined for compounds 1–8. The values of semi-efficient inhibitory concentration are in the range of 24–73 μmol/l, which allows us to classify these compounds as moderate inhibitors of the enzyme.

As it can be seen from the plot (Figure 3), a sharp change in inhibitory activity was observed in a narrow concentration range for compounds 7, 8, which can serve as a basis for the assumption that the enzyme is highly sensitive to these inhibitors.

2.3. Evaluation of the predictive ability of the M3, M6, M9, M12, M15, and M18 models based on Compounds 1–8 in the Test Set TS3

Subsequently, the consensus models M3, M6, M9, M12, M15, and M18 were used to predict the numerical values of the pIC₅₀ parameter for compounds 1–8 from the TS3 set. The results of these calculations in comparison with the experimental values of the pIC_{50 exp} parameter and the 2·RMSEP criterion, which corresponded to 100% and 95% of the data from the TS1 and TS2 sets, respectively, for these compounds are shown in Table 5 and Table 6. In selecting the numerical values of the 2·RMSEP criterion corresponding to the M3, M6, M9, M12, M15, and M18 models for subsequent comparison with the ΔpIC₅₀ values, which are in fact equal to the modulus of the difference between the experimental and predicted values of the pIC₅₀ parameter, we focused on the minimum value of this parameter (see Table5). It should be noted that all the tested compounds were within the range of applicability of the M3, M6, M9, M12, M15, and M18 models. These models were selected for the prediction of the pIC₅₀ parameter for the 15-LOX inhibitors 1–8 because each of these consensus models included 320 partial regression relationships, which makes it possible to take into account the structural characteristics of each of the tested compounds in the most objective and complete way. In addition, it was of scientific interest to explore the applicability of these models with satisfactory and worst case predictive performance, judging by the statistical criteria derived from TS1–TS2, to predict the target property for new compounds not included in the modeling.

The data in Table 5 and Table 6 allow us to conclude that almost all numerical values of pIC₅₀ for compounds 1–8 predicted by the M3, M6, M9, M12, M15, and M18 models fall within the 95% confidence interval equal to ± 2·RMSE, i.e. the difference between the predicted and experimentally determined values of the pIC₅₀ parameter for 15-LOX inhibitors 1–8 does not numerically exceed the minimum value of the 2·RMSEP criterion for each of the models (see Table 5 and Table 6). The M3 model showed the highest prediction error by this criterion when the pIC₅₀ parameter was predicted for compound 4, but even in this case the difference between the experimental and theoretically predicted values of the pIC₅₀ parameter fell within the 2·RMSEP range. However, the same model (M3) showed the smallest prediction error of the pIC₅₀ parameter for compound 8. This demonstrates the good predictive ability and correctness of our constructed models, as well as the applicability of the GUSAR2019 program for modeling 15-LOX inhibitors.

Thus, all the QSAR consensus models M1–M18 are characterised by high descriptive and moderate predictive power when comparing experimental and predicted values of pIC₅₀ based on TrS1 and TrS2 training sets structures, external and internal test sets TS1 and TS2, and compounds 1–8. These models can be used for virtual screening of virtual libraries and databases to search for new 15-LOX inhibitors in the series of homo- and heterocyclic compounds with common structural formulae I-XVI.

3. Research Methods

The simulation procedure was performed for the compounds whose formulas are shown in Figure 4.

3.1. The Methodology of the Computational Experiment

QSAR modeling of 15-LOX inhibitors with general structural formulae I–XVI (Figure 4) was performed using the GUSAR 2019 (General Unrestricted Structure Activity Relationships) computer program [53,65,66,92,93,94,95,96,97,98].

The QSAR models were built in several steps that are the basis of the operation of this program and have been described in detail in our previous work [60,61,62,63,67,68,69,70,71,72,73,74]. The stages of QSAR model building are shown schematically in Figure 5.

3.2. Formation of the Training and Test Sets

Training set TrS1 and test set TS1 were generated from the structure set S1. The structure set S1 contained 100 15-LOX inhibitors with their corresponding pIC₅₀ values. Training set TrS2 and test set TS2 were formed based on the structures contained in training set TrS1. Figure 6 shows a scheme that clearly illustrates the strategy of forming the training and test sets TrS1–TrS2, TS1–TS2.

The pIC₅₀ parameter for each compound included in the data set S1 (and the training TrS1–TrS2 and test sets TS1–TS2 generated from it) was calculated as the negative decimal logarithm of its corresponding IC₅₀ value (in mol/l). Numerical IC₅₀ values for the 15-LOX inhibitors being modelled were measured experimentally and are given in [75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90]. The complete list of organic compounds from which the data set S1 was generated with their corresponding pIC₅₀ characteristics is presented in Table S8 in the Supplementary Material.

QSAR models M1–M9 were built using the TrS1 training set, which contained 84 structures of 15-LOX inhibitors with their corresponding values of the pIC₅₀ parameter. To test the predictive ability of the M1-M9 models, a test set TS1 containing 16 15-LOX inhibitors with their corresponding values of the pIC₅₀ parameter was used. Both sets were obtained by partitioning the data set S1, in which all compounds were previously ranked in ascending order of the pIC₅₀ parameter, in 5:1 ratio. The structures were partitioned into the training set TrS1 and the test set TS1 by transferring every sixth compound from the data set S1 to TS1. The remaining 84 structures of the 15-LOX inhibitors were used to form the training set TrS1.

The TrS2 set contained 70 15-LOX inhibitors with their corresponding values of the pIC₅₀ parameter. It was intended for constructing the QSAR models M10–M18. The validity of the QSAR models M10–M18 was tested using the TS2 set. Both sets, TrS2 and TS2, were generated on the basis of TrS1. The same principle was used to generate the training set TrS1 and the test set TS1 from the data set S1. The characteristics of the training sets TrS1, TrS2 and the test sets TS1, TS2 are shown in Table 7 and Table 8 respectively. The data in these tables indicate that the compounds of the training and test sets are fairly evenly distributed over the entire range of pIC₅₀ variation. At the same time, the range of variation of the pIC₅₀ parameter for the 15-LOX inhibitors included in training sets TrS1–TrS2 and test sets TS1–TS2 exceeds the ΔpIC₅₀ value of 3 (i.e. ΔpIC₅₀ > 3), which determines the correctness of the further QSAR modeling process [44]. In addition, as it can be seen from Figure 4, the training sets are characterised by rather a high degree of molecular diversity.

The structures of the compounds in the training and test sets TrS1–TrS2, TS1–TS2 were plotted in Marvin Sketch 23.4 software [99] and then converted to SDF format using Discovery Studio Visualiser software [100].

3.3. Building QSAR Models

The M1-M18 QSAR models were built on the basis of two types of substructural descriptors of atomic neighborhoods: QNA (Quantitative Neighborhoods of Atoms) and MNA (Multilevel Neighborhoods of Atoms) [47,48,53,60,61,62,63,64,65,66,67,68], and three types of whole molecule descriptors (topological length, topological volume, and lipophilicity). These types of descriptors were automatically calculated by the GUSAR 2019 program. At the same time, the QNA and MNA descriptors are unique characteristics of molecules, and their calculation is available in different versions of the GUSAR software (GUSAR 2011, GUSAR 2013, and GUSAR 2019). The ideology of calculating QNA and MNA descriptors was proposed by Professor V.V. Poroikov's research team. It is described in detail in the Appendix and in a number of articles [47,64,65,92,93,94,95,96,97]. The rather complicated mathematical apparatus used for calculating QNA descriptors complicates their physical interpretation. Therefore, they are not explicitly displayed in the section dealing with calculations.

The MNA descriptors are computed using the PASS (Prediction of Activity Spectra for Substances) algorithm [53,66], which predicts approximately 6,400 “biological activities” with an accuracy threshold of an average prediction of at least 95%. These descriptors are generated based on the structural formulae of the chemical compounds without using any pre-compiled list of structural fragments [47,53,66,92]. They are generated as a recursively defined sequence:

• Zero-level MNA descriptor for each atom is the mark A of the atom itself;
• Any next-level MNA descriptor for the atom is the substructure notation A (D1D2 … Di …), where Di is the previous-level MNA descriptor for i–th immediate neighbor of the atom A.

The neighbor descriptors D1 D2 … Di … are arranged in a unique manner. This may be, for example, a lexicographic sequence. The MNA descriptors are generated using an iterative procedure, which results in the formation of structural descriptors that include the first, second, etc. neighborhoods of each atom. The label contains not only information about the type of atom, but also additional information about whether it belongs to a cyclic or acyclic system, etc.

Three methods were used to reduce the descriptor space and select the most significant descriptors:

• Self-consistent regression (SCR) method;
• The method of combining self-consistent regression with radial basis functions (RBF-SCR);
• The Bath method, which combines the simultaneous use of SCR and RBF-SCR methods in a unique way.

All three of these unique methods were also developed by Professor Poroikov's research team and implemented in the GUSAR 2019 program for selecting the most appropriate options. A more detailed description of each method can be found in the Appendix and in the relevant papers [66].

The stability of the models was tested by using a sliding control procedure, with a 20-fold randomized release of 20% of the compounds from the training samples TS1 and TS2. Both of these procedures are automatically implemented in the GUSAR 2019 program [47,48,53,60,61,62,63,64,65,66,67,68].

Each of the eighteen final QSAR M1-M18 models was based on a consensus approach. This approach involves combining several regression equations into one model, which is done automatically based on the similarities between the equations.

Each of the final QSAR models M1-M2 and M4-M5, M7-M8, M10-M11, M13-M14, M16-M17 included 20 partial regression dependences. At the same time, the M1, M4, M7, M10, M13, and M16 models were based on QNA descriptors and three additional descriptors that describe the topological length, topological volume, and lipophilicity of the modeled 15-LOX inhibitors. The M2, M5, M8, M11, M14, and M17 models were built on a similar principle, but they were based on NA descriptors, with automatic addition of the three whole molecule descriptors described above. The M3, M6, M9, M12, M15, and M18 models each included 320 partial regression dependences. At the same time, each of these 320 particular models was built independently of each other based on the three descriptors of the entire molecule described above, with addition of the QNA or MNA descriptors.

Due to the specifics of the calculation process, which is described in detail in the Supplemental Material, QNA and MNA descriptors are not amenable to unambiguous physical interpretation. Therefore, the regression equations based on these descriptors are not explicitly displayed in the GUSAR 2019 program.

3.4. Evaluation of the Descriptive and Predictive Ability of QSAR Models

The descriptive ability of the M1-M18 models was evaluated using several metrics. These included metrics based on the coefficients of determinations R2, R20, average R2m and CCC, as well as metrics estimating errors in predicting pIC₅₀ values (RMS error (RMSE), mean absolute error (MAE), and standard deviation (SD). The parameters of the predictive ability of the M1-M18 models also included metrics based on the coefficients of determination R², R²₀, $\overline{R_m^2}$ , CCC, Q²_F1 and Q²_F2 , as well as metrics estimating errors in predicting pIC₅₀ values (RMS error (RMSE), mean absolute error (MAE), standard deviation (SD)).

These statistical parameters were calculated using the External Validation Plus 1.2 software for 100% and 95% of the data (to account for errors) contained in the training and test samples [100]. The Supplemental Material contains formulas for automatically calculating these criteria in this program. The internal verification of the M1-M6 models was performed using LMO cross-validation (Q2LMO), with 20% of compounds excluded from the training sets.

The threshold values of the validation criteria for the above parameters, based on which the descriptive and predictive ability of the QSAR M1-M18 models was evaluated, are presented in Table 8.

Acceptable values of different types of determination coefficients, based on R2, as well as ranges of variation in MAE and RMSD values and criterion B for assessing the descriptive and predictive ability of QSAR models M1-M18 were calculated taking the recommendations from the leading scientists in the field of QSAR modeling [44,45,46,57] into account.

Thus, the permissible range of variation for the MAE parameter was estimated considering the range of variability in the pIC₅₀ parameter for the compounds in training sets TrS1 and TrS2, using the following formulas:

MAE < 0.1 pIC₅₀ — if the criterion is met, then the models are characterized by a high predictive ability

MAE = [0.1; 0.15} pIC₅₀ —if the criterion is met, then the models are characterized by a moderate predictive ability

MAE > 0.15 pIC₅₀ — if the criterion is met, then the models are characterized by a low predictive ability

The permissible range of variation of criterion B, where B=MAE+3SD, was estimated using the formulas:

B < 0.2 pIC₅₀ — if the criterion is met, then the models are characterized by a high predictive ability

B = [0.2; 0.25] pIC₅₀ — if the criterion is met, then the models are characterized by a moderate predictive ability

B > 0.25 pIC₅₀ — if the criterion is met, then the models are characterized by a low predictive ability

Additionally, the predictive ability of the QSAR M1-M18 models was evaluated by comparing the predicted pIC₅₀ values with the experimental values of the same parameter for the new promising 15-LOX inhibitors 1-8 contained in the TS3 test set (Figure 4). These compounds were missing in the S1 data set and, accordingly, did not participate in the consensus models.

3.5. The Technique of the Biochemical Experiment to Measure Inhibitory Activity

Biological tests of inhibitors 15-LOX 1-8 (Figure 4) were performed at the Federal State Budgetary Educational Institution of Higher Education "Saratov National Research State University named after N.G. Chernyshevsky" under the supervision of Professor O.V. Fedotova [101]. The inhibition of the catalytic activity of 15-lipoxygenase by the substances studied was evaluated according to the procedure reported elsewhere [102,103,104]. Quercetin, a well-known inhibitor of 15-LOX, was employed as the positive control.

For biological tests, Lipoxidase/15-lipoxygenase (Sigma Aldrich, Germany, lyophilizate, 15 million Units, CAS Number: 9029-60-1) was used. The study was carried out in borate buffer (0.2 M, pH 9.0). Linoleic acid (Sigma Aldrich, Germany, purity ≥ 95.0%, CAS Number: 60-33-3) was used as the oxidation substrate. The activity of 15-LOG in the presence of the compounds studied was evaluated by spectrophotometric recording of changes in the concentration of 13-hydroperoxylinoleic acid, the product of oxidative transformation of linoleic acid [105], using a LEKI SS2110UV two-beam scanning spectrophotometer (CJSC LOIP, Russia).

The initial concentration of lipoxygenase in the sample was 167 U/ml. Its final concentration in the sample was 134 µM. In order to study the inhibitory activity of compounds, they were added to the sample as solutions in DMSO.

To conduct the experiment, 1400 µl of a substrate solution was placed in a cell, along with 24 µl of a sample in DMSO (for the prototype) or DMSO (for the control sample). Immediately after adding 76 µl of the enzyme solution, the timer was turned on. After 90 seconds, the optical density (λ = 234 nm) was recorded.

The activity value obtained from the control experiment was used as a baseline, or 100%, of enzyme activity. It was measured with only the enzyme solution and solvent (DMSO) present in the cell, without any active substance.

The percentage of activity inhibition by compounds 1-8 was determined as a relative decrease in the optical density of the solution, using the formula (1):

$$I={\displaystyle \frac{D_c{-D}_t}{D_c}}\times 100\%,$$

(1)

where $D_c$ is the optical density index of the control sample 90 seconds after the start of the reaction; $D_t$ is the optical density index of the sample containing the test compound 90 seconds after the start of the reaction.

The range of working concentrations of compounds 1-8 was 1-100 µM.

The IC₅₀ values of compounds 1-8 were determined using linear interpolation between the points closest to 50% inhibition by means of Microsoft Excel 2016. Enzyme activity experiments were performed in triplicate. The values are expressed as means ±SD. Student’s t-test was employed for determination of statistical significance, using a P value of 0.05 or less as a criterion for significant inhibition.

4. Conclusion

Using the QSAR methodology implemented in the GUSAR 2019 program, a quantitative structure-inhibitory activity relationship has been found for a series of 100 15-LOX inhibitors based on a series of derivatives of phenol, resorcinol, anacardic acid, dimethoxybenzene, alkyl ester of 2-(4-isobutylphenyl)propionic acid, 1,3-diarylprop-2-yn-1-one, proline, pyrrole, alkyl ester of 2-methylfuran-3-carbonyl acid, ketoprofen, naphthalene, 1,4-di-N-oxide of quinoxaline, isoflavone, 4-hydroxy-2-(phenylmethyl)benzofuran and coumarin with general structural formulas I–XIV. These compounds had various degrees of inhibition of the catalytic activity of 15-LOX. The variation range of the pIC₅₀ parameter was 3.873 (pIC₅₀=3.873). Based on a combination of MNA and QNA descriptors with three whole molecule descriptors including topological length, topological volume, and lipophilicity, eighteen statistically significant and valid consensus models (M1-M18) were generated.

All models reproduced the experimental data contained in the training samples with a high degree of accuracy. Cross-validation with 20-fold deletion of 20% of the data from the training samples also showed good results. The reliability of the prediction of the pIC₅₀ parameter based on the evaluation of this parameter for compounds of two test samples and ten subsequently experimentally studied compounds showed a moderate predictive ability of the QSAR M1-M18 models.

The satisfactory match of the theoretically calculated pIC₅₀^pred values with the experimental pIC₅₀^exp values for compounds of test sets TS1-TS2 and compounds 1-8 opens up prospects for application of M1-M18 models in virtual screening of virtual libraries and databases in the search for new potentially efficient 15-LOG inhibitors in these sources.