Structure-odor relationship in pyrazines and derivatives: A physicochemical study using 3D-QSPR, HQSPR, Monte Carlo, molecular docking, ADME-Tox and molecular dynamics

2.1 Selection of the data set

In this new study, a series of 78 selected pyrazine derivatives with odor threshold (t) values were extracted from reference, these molecules were studied to build the CoMSIA, HQSPR, Topomer CoMFA, and monte Carlo models, 59 molecules were selected to propose the quantitative model (training set), and 19 molecules were used to test the performance and quality of the proposed model (test set). Table 1 shows the chemical structures of the studied molecules, and the property values of the experimental odor thresholds in log(1/t) (Belhassan et al., 2019).

Table 1 Chemical structures of pyrazine compounds and their experimental odor threshold properties.

N°	R1	R2	R3	R4	Log(1/t)	N°	R1	R2	R3	R4	Log(1/t)
M1	H	H	H	H	3.523	M40	C5H11	SC2H5	H	H	9.000
M2	CH3	H	H	H	4.523	M41	C8H17	SC2H5	H	H	8.699
M3	C2H5	H	H	H	5.398	M42	C10H21	SC2H5	H	H	6.921
M4	C3H7	H	H	H	6.523	M43	H	SC6H5	H	H	6.398
M5	C4H9	H	H	H	6.398	M44	CH3	SC6H5	H	H	6.523
M6	C5H11	H	H	H	8.301	M45	C3H7	SC6H5	H	H	7.046
M7*	C6H13	H	H	H	6.699	M46	C5H11	SC6H5	H	H	8.000
M8	C7H15	H	H	H	7.000	M47	C8H17	SC6H5	H	H	7.097
M9*	C8H17	H	H	H	6.398	M48*	C10H21	SC6H5	H	H	6.523
M10*	C10H21	H	H	H	5.955	M49	H	OCH3	H	H	6.398
M11	CH3	H	H	CH3	6.398	M50*	CH3	OCH3	H	H	8.155
M12*	CH3	H	H	OCH3	7.770	M51	C2H5	OCH3	H	H	8.000
M13*	CH3	H	H	OC2H5	8.301	M52*	C3H7	OCH3	H	H	9.921
M14	CH3	H	H	SCH3	7.699	M53	C4H9	OCH3	H	H	10.301
M15	CH3	H	H	COCH3	6.523	M54*	C5H11	OCH3	H	H	10.699
M16	C2H5	H	H	CH3	7.398	M55*	C6H13	OCH3	H	H	10.155
M17	CH3	H	CH3	H	7.097	M56*	C7H15	OCH3	H	H	10.585
M18	CH3	H	OCH3	H	7.699	M57	C8H17	OCH3	H	H	10.222
M19*	CH3	H	OC2H5	H	7.921	M58	C10H21	OCH3	H	H	7.398
M20	CH3	H	SCH3	H	7.222	M59	(CH2)2CH(CH3)2	OCH3	H	H	11.201
M21*	CH3	H	COCH3	H	6.398	M60	(CH2)3CH = CH2	OCH3	H	H	10.523
M22	C2H5	H	CH3	H	7.796	M61	CH(CH3)C2H5	OCH3	H	H	10.398
M23*	CH3	CH3	H	H	6.097	M62	CH2CH(CH3)C3H7	OCH3	H	H	11.097
M24	C2H5	CH3	H	H	6.301	M63	CH2CH(CH3)2	OCH3	H	H	10.347
M25	C3H7	CH3	H	H	7.222	M64	CH2CH(CH3)C2H5	OCH3	H	H	10.921
M26*	CH(CH3)2	CH3	H	H	7.796	M65	(CH2)3CH(CH3)2	OCH3	H	H	11.222
M27	CH3	COCH3	H	H	7.699	M66	CH(CH3)2	OCH3	H	H	10.620
M28	H	SCH3	H	H	6.699	M67	(CH2)3CH = CHCH3 (E)	OCH3	H	H	9.886
M29*	CH3	SCH3	H	H	8.398	M68	(CH2)3CH = CHCH3 (Z)	OCH3	H	H	9.301
M30*	C2H5	SCH3	H	H	7.398	M69	H	OC2H5	H	H	7.097
M31	C3H7	SCH3	H	H	9.000	M70*	CH3	OC2H5	H	H	9.097
M32	C5H11	SCH3	H	H	9.921	M71	C2H5	OC2H5	H	H	7.699
M33	C8H17	SCH3	H	H	9.155	M72	C5H11	OC2H5	H	H	10.097
M34	C10H21	SCH3	H	H	7.699	M73	C8H17	OC2H5	H	H	8.699
M35	CH(CH3)2	SCH3	H	H	10.328	M74*	C10H21	OC2H5	H	H	7.222
M36	H	SC2H5	H	H	6.046	M75	H	OC6H5	H	H	7.523
M37	CH3	SC2H5	H	H	7.155	M76	CH3	OC6H5	H	H	6.699
M38	C2H5	SC2H5	H	H	7.222	M77	C5H11	OC6H5	H	H	7.301
M39	C4H9	SC2H5	H	H	8.398	M78	C10H21	OC6H5	H	H	7.155

Odor threshold (t); Test set (*).

2.2 Modelling and molecular alignment

We used the Sybyl software to draw molecules from the database, and these registered molecules were reduced using molecular mechanics. We picked the Alignment database technique, which has a force field parameter (Tripos), an energy gradient convergence threshold of 0.01KCal/mol, and a charge type of Gasteiger-Hückel with a maximum of 1000 iterations. The structural alignment phase is included in the 3D-QSPR study because alignment is a necessary first step in statistical analyses of flexible structures in 3D-QSPR models (Vyas et al., 2017).

The alignment databases were chosen for their ability to produce satisfactory results, and the M65 compound was chosen as the model molecule against which the other compounds were matched (Fig. 1).

Close

Fig. 1

Database alignment.

Export to PPT

2.3 Monte Carlo regression

The method for optimal hybrid descriptors is based on SMILES structures and molecular graphical representations, both of which are based on correlation weight optimization. The Monte Carlo regression model was created using the CORAL application (https://www.insilico.eu/Coral). CORAL works out as follows:

The molecular structure can be represented using SMILES and/or a molecular graph. The two molecular structure representations stated above have an effect on optimal hybrid descriptors. Maximizing the correlation weight of the SMILES attributes with the correlation weight of the chart invariants yields optimal hybrid descriptors. A parameter prediction model is fitted with the best DCW hybrid descriptor (T, N): $$\mathit{PropertyBiological}=C0+C1\times DCW(T,N)$$

The optimum descriptor of the molecular feature function taken from SMILES is DCW (T, N). T and N are the best threshold and number of epochs determined during construction, respectively (Liman et al., 2022).

T is the threshold for classifying molecular qualities (structural features) into two groups: rare and frequent. If a structural feature is widespread in energizing substances, it can be linked to property; nevertheless, unusual features have less value and are viewed as noise. The weights of correlation (CW) of active structural attributes (ASs) are used to construct the ideal descriptor: for rare attributes, the CW is set to zero, signaling that they are not included in the model.

The N represents the number of times a DCW of a molecular characteristic is changed. The DCW values are obtained during construction using Monte Carlo optimization. The goal of the Monte Carlo approach is to maximize the correlation coefficient between the ideal descriptor and the parameter. Optimization is done with the formation and invisible formation ensembles to avoid over-formation (Toropov et al., 2020).

In this work, optimization was performed by balancing the correlations of the TF1 and TF2 target functions in the CORAL 2017 software to build the QSPR models, which are used to obtain the best results. In order to avoid over-training, the correlation balance method, which requires a high correlation between training and validation packages, was proposed.

2.4 Hologram-QSPR (HQSPR)

The Tripos HQSPR method is a computational approach for predicting molecular properties based on 2D molecular fragments. The method takes advantage of the fact that molecular properties can be correlated with the molecular structure, and that the molecular structure can be represented as a series of 2D molecular fragments.

In the first step of the HQSPR method, 2D molecular structures are divided into all possible linear and branched fragments, and each unique fragment is assigned a specific integer using a cyclic redundancy check (CRC) algorithm. These integers are then used to create a molecular hologram that serves as a structural descriptor, containing topological and compositional information about the molecular structure (Vyas et al., 2017).

In the second step, the molecular holograms are correlated with the corresponding biological properties using partial least squares (PLS) analysis. A drop-out cross-validation (LOOCV) is performed to identify the optimal number of explanatory variables or components that result in an optimal model. The final step involves deriving a mathematical regression equation that links the values or components of the molecular holograms with the corresponding biological properties.

The Tripos HQSPR method is a powerful tool for predicting molecular properties and has potential applications in fields such as drug discovery, materials science, and chemical engineering. By allowing for the automated analysis of large datasets, the Tripos HQSPR method can help speed up the process of discovering new drugs and materials, as well as improving our understanding of molecular structure–property relationships (Tabti et al., 2022c). $${\mathit{BP}}_i=\mathrm{C}\mathrm{s}\mathrm{t}+\sum _{J=1}^L{X}_{\mathit{ij}}{C}_j$$ where ${\mathit{BP}}_i$ is the biological property of the i^th compound, x_ij is the occupancy value of the molecular hologram of the i^th compound at position or bin j, C_j is the coefficient for bin j derived from the PLS analysis, and L is the hologram length. One of the drawbacks of HQSPR is the fragment collision problem, which occurs during the fragment hashing process (Fig. 2).

Close

Fig. 2

Development of the predictive equation obtained by partial least squares analysis of the HQSPR model.

Export to PPT

Although hashing reduces the length of the hologram, it results in different fragments in the same bin. The hologram length, a user-definable parameter, controls the number of bins in the hologram, and changing the hologram length can result in a change in the bin occupancy pattern. The program provides 12 default lengths but after several trials, we chose the first four which proved to give good predictive models on different data sets.

2.5 CoMFA and CoMSIA analysis

CoMSIA was created to address some of the shortcomings of CoMFA. In CoMSIA, molecular similarity indices derived from the modified SEAL similarity criteria are used as descriptors to account for steric, electrostatic, hydrophobic, and hydrogen bonding properties all at the same time. Individual indices are calculated by comparing the similarity of each molecule in the dataset to a common probe atom (with radius 1A, character + 1, and hydrophobicity + 1) located at the intersections of a surrounding grid/network. The mutual distance between the probe atom and the atoms in the aligned data set molecules is also considered when calculating similarity at all grid points (Wang et al., 2021).

Gaussian-like functions are used to characterize the distance dependence and calculate the molecular properties. Because the underlying Gaussian-like functional forms are “smooth” and free of singularities, their slopes are not as steep as those of the Coulomb and Lennard-Jones functions in CoMFA, thus eliminating the need for arbitrary cutoff limits. CoMSIA and CoMFA are included in the Sybyl software from Tripos Inc.

2.6 Topomer CoMFA

The CoMFA topomer method was developed to overcome the alignment issue in CoMFA; this technique takes into account not only the steric and electrostatic descriptors, but also the molecular topomer descriptors (Tabti et al., 2022a). A topomer is a molecular fragment with a single internal geometry or “Pose.” To produce better models, the two main fragmentation strategies (R1, R2, R3 and R4 – Common Core) were used (Fig 3). The calculation of steric and electrostatic fields is similar to the CoMFA approach.

Close

Fig. 3

Cutting mode in Topomer CoMFA.

Export to PPT

2.7 External validation plus

Validation is the most important stage in constructing quantitative structure–activity relationship (QSPR) models. It ensures the dependability of the final QSPR model as well as the acceptance of each phase in the model generation process. External validation (using an independent test set) is commonly employed to assess the predictive accuracy of a QSPR model (Ouabane et al., 2022).

The external predictivity of QSAR models is represented by computing and analyzing various validation metrics, which can be broadly classified into two classes or categories, namely, $R^2$ -based metrics such as $R_{\mathit{test}}^2$ , $Q_{\mathit{ext}}^2$ (F1), $Q_{\mathit{ext}}^2$ (F2), average $R_m^2$ , Δ $R_m^2$ , and purely error-based metrics such as root mean square error (RMSE), mean absolute error (MAE), and so on.

Furthermore, before doing external validation on a model, it is necessary to check for the presence of a systematic mistake that violates the basic assumptions of the least square’s regression model. If the model contains significant systematic error (bias), it should be ignored, and any external validation test performed on such a biased model is worthless.

External Validation Plus is a tool that checks for systematic errors in the model (only if the user chooses; optional) and then computes all of the required external validation parameters while evaluating the performance of a QSAR model's actual prediction quality using our recently proposed MAE-based criteria (Yu et al., 2023).

Systematic error is present when any one or more of the following conditions is true:

1. $\mathit{NPE}/NNE>5orNNE/NPE>5$

NPE: Number of positive errors, and NNE: Number of negative errors

1. $\mathit{ABS}(MPE/MNE)>2orABS(MNE/MPE)>2$

MPE: Mean of positive error, MNE: Mean of negative errors, and ABS: Absolute value

1. $\mathit{AAE}-ABS(AE)<0.5\times AAE$

AAE: Average of absolute prediction errors, and AE: Average of predictions errors

1. $R^2$  greater than 0.5 for residuals sorted on experimental response values

In the external Validation tool, the above threshold values are set as default and user can change the threshold values as required.

If any of the conditions listed above are met, the program will alert you to the presence of systematic error in the output file and recommend that you discard the model (only if the user chooses to check for systematic error). Furthermore, examining the external validity parameter is pointless if a systematic inaccuracy exists.

2.8 In silico pharmacokinetic ADME-Tox study

The ADME/Tox examination is a crucial step in the drug discovery process as it determines the pharmacokinetic properties of a drug candidate, including its absorption, distribution, metabolism, excretion, and toxicity. However, this process can be lengthy and costly, and not all drug candidates meet the required criteria (Ugale et al., 2017).

To overcome these limitations, there are now various online and offline tools available that can help predict ADME/Tox properties in silico. The online pkCSM tool and the SwissADMET online server are two examples of these tools. These tools use computer models and algorithms to estimate the pharmacokinetic properties of new compounds without the need for extensive experimental testing (Hajji et al., 2021a).

By using these online tools, researchers can quickly and cost-effectively predict the ADME/Tox properties of drug candidates, allowing them to make informed decisions about which candidates to pursue in the drug discovery process. This can contribute to speeding up the drug discovery process and increasing the success rate of developing new drugs.

2.9 Molecular docking study

The mechanism of action between olfactory receptors and odorant molecules is a crucial step in the development of new odorant molecules for use in various pharmaceutical, cosmetic and agro-food fields. In this study, molecular docking simulation is a useful tool to determine this binding between pyrazine derivatives and the active site of the protein responsible for the transmission of the odorant molecules and we have chosen a receptor (PDB ID: 1DZK) to make comparisons with the reference ligand 2-isobutyl-3-metoxypyrazine, this odorant molecule noted M63 in the data base (Fig. 4).

Close

Fig. 4

Representation of the crystalline structure 1DZK.pdb, chain A in green (ligand: 2-methoxy-3-(2-methylprop-2-en-1-yl) pyrazine) and chain B in blue (ligand: 2-methoxy-3-(2-methylpropyl) pyrazine).

Export to PPT

In this study, we selected the most and least active compounds, which the reference ligand methoxypyrazine in the database on the one hand, and on the other hand the new compounds proposed from the molecule M65 with the highest odor threshold. All steps were performed by various available software such as CB-Dock (Ahfid et al., 2023), AM-Dock, Autodock-Vina (Mahfuz et al., 2022), Pyrex (En-nahli et al., 2022) and Sybyl. The molecular docking results were visualized by Discovery Studio (Hajji et al., 2021b).

The 3D structure of the receptor was chosen from the crystal structure (PDB ID: 1DZK) with a resolution of 1.48 Å, obtained from the RCSB protein database (https://www.rcsb.org). The structure of porcine odorant binding protein (pOBP) is a monomer of 157 amino acid residues, purified in abundance from porcine nasal mucosa. These residues were minimized by Sybyl software, then we added non-polar hydrogen atoms to the protein chain, Gasteiger-Huckel type charges were chosen, the convergence energy gradient was set to 0.01 kcal/mol maintaining 1000 times of optimization as the maximum limit cycle number.

The grid parameters are: positions P1: (x = 12.03, y = 4.71, z = 24.10), P2: (x = 14.65, y = 4.28, z = 25.12), dimensions (x = 30, y = 30, z = 30), energy range of order 10 and completeness equal to 8. We also used the re-Docking method to validate the Docking procedure, this technique based on the superposition of the M63 molecule after and the reference ligand were extracted before the molecular Docking, the RSMD value was displayed by the Discovery studio visualization software. Docking is reliable if the root means square deviation (RMSD) value is low (typically (1.5–2) Å depending on ligand size) (Shi et al., 2022).

2.10 Molecular dynamics simulations

To study the stability of complex (ligand–protein) interactions, molecular dynamic simulation (MD) calculations were performed on the best molecular docking poses. Using the CHARMM force field parameters, generate protein input files using the c solution generator and ligand topology using the CHARMM general force field Param-Chem server. The building blocks of the CHARMM-GUI solution consist of five steps. In the first step, the coordinates of the protein–ligand complex are read by the tool. The second step is to solve the protein–ligand complex and determine the shape and size of the system. ${\mathit{Na}}^{+}$ and ${\mathit{Cl}}^{-}$ ions are added at this stage to neutralize the system. Periodic boundary conditions (PBC) are defined in the third step, which allows a large system to be approximated using a single cell which is then replicated in all directions. The simulation is only for atoms inside the PBC box. Bad contacts are removed at this stage by performing a short minimization. The fourth and fifth steps are to balance the system and production. Balancing is carried out in two steps: NVT and NPT, ensuring the system reaches the desired temperature and pressure (Koubi et al., 2023).

The input files for balancing and production are then downloaded and the desired changes are made, including the number of steps in the MD simulation, the frequency of trajectory recording and energy calculation, etc. GROMACS 2020.2 was used for both equilibration and production of all MD simulations. All complexes were first solvated in a cubic box of TIP3P water, then ${\mathit{Na}}^{+}$ and ${\mathit{Cl}}^{-}$ ions were added to neutralize the net atomic charge of the whole system by randomly replacing the water molecules. Periodic boundary conditions (PBC) were imposed taking into account the shape and size of the system. The unrelated interactions were treated with a cut-off distance of 12 Å and the neighbors search list was buffered with the Verlet cut-off scheme, and the long-range electrostatic interactions were treated with the Ewald particle mesh method (SME) (Khaldan et al., 2022).

The CHARMM36 force field was applied to the protein–ligand complex. Prior to the production simulation, system energy minimization was performed using the gradient descent algorithm (5000 steps) (Balupuri et al., 2020). The complex was then balanced to stabilize its temperature and pressure by subjecting it to an NVT and NPT assembly and simulating it for 125 ps at a temperature of 300,15 K using 400 Kj. ${\mathit{mol}}^{-1}$ . ${\mathit{nm}}^{-2}$ and 40 Kj. ${\mathit{mol}}^{-1}$ . ${\mathit{nm}}^{-2}$ on the main axis and side chains respectively. Finally, the complex is subjected to a production simulation of 100 ns in NPT assembly at 300.15 K and 1 bar. To maintain the temperature, a Nose-Hoover thermostat was used and to maintain the pressure, a Parrinello-Rahman barostat was used. The LINCS algorithm was used to constrain hydrogen bonds using inputs provided by CHARMM-GUI. The 300 K V-rescale thermostat with a 1 ps coupling constant was used. The trajectories were recorded every 2 ps. Simulations of 100 ns in NPT assembly were performed for the production stage (Mahfuz et al., 2022).

3.1 Results obtained by Monte-Carlo approach

QSPR model based on SMILES is a type of quantitative structure–property relationship model that uses Simplified Molecular Input Line Entry System (SMILES) notation to represent molecular structures. SMILES strings are used as input features to predict various physical and chemical properties of the molecules. The aim of such models is to establish a correlation between molecular structure and property, allowing for the prediction of properties for new molecules based on their SMILES representation. A total of 78 molecules were developed from two target function types, TF1 (without IIC) and TF2 (with IIC), using equations (TF1) and (TF2) respectively, to select consistent statistical performance. $$Pred\mathit{p}\mathit{t}=-5.7746958(\pm 0.1095872)+0.1383070(\pm 0.0011208)\ast DCW(1,10)(TF1)$$ $$Pred\mathit{p}\mathit{t}=-0.0333622(\pm 0.1840351)+0.0751853(\pm 0.0018147)\ast DCW(8,10)(TF2)$$

The results of the statistical parameters obtained for the QSPR models generated from the odor threshold of odorant molecules based on SMILES descriptors for four sets randomly (+): Training set, (-): Invisible training, (#): Calibration and (*): test set are recorded in Table 2 and 3.

The results indicate that the (TF2) model gives better results than the (TF1) model, which means that the addition of IIC (TF2) to the model weight improves its prediction. The comparison of experimental log(1/t) values to the calculated values for both models, presented in Table 4, clearly shows the statistical reliability of the TF2 models and that they meet the criteria established by external validation.

3.1.1

3.1.1 Y-randomization test

Y-randomization is a tool used in the validation of QSPR/QSAR models that compares the performance of the original model in the data description ( $R^2$ ) to models produced for permuted (randomly mixed) response, depending on the original model's descriptor directory and creation procedure. The validation method's results are presented in Table 5.

Table 5 Y-randomization test results for both TF1 and TF2 models.

	Training		Invisible training		Calibration
	TF2	TF1	TF2	TF1	TF2	TF1
N	20	20	20	20	19	19
Origin	0.9569	0.7840	0.9721	0.7616	0.9085	0.8173
1	0.0086	0.0357	0.0016	0.3507	0.0681	0.0001
2	0.1002	0.0310	0.0070	0.0022	0.1032	0.1039
3	0.0279	0.0165	0.1187	0.0321	0.1032	0.1378
4	0.1570	0.1145	0.0128	0.0912	0.0022	0.0268
5	0.0805	0.0456	0.2302	0.0036	0.0528	0.0088
6	0.0184	0.0282	0.0022	0.0003	0.0356	0.1381
7	0.0035	0.0035	0.0766	0.1430	0.0147	0.0125
8	0.0180	0.0293	0.0000	0.0006	0.0926	0.0005
9	0.0355	0.0222	0.0203	0.0011	0.3408	0.0622
10	0.1139	0.0027	0.0024	0.1460	0.0251	0.0145
$R_{\mathit{Average}}^2$	0.0563	0.0329	0.0472	0.0771	0.0940	0.0505
CRp2	0.9283	0.7674	0.9482	0.7220	0.8602	0.7916

The T threshold and N epochs were selected to produce the best statistical indicators for the calibration set. TF1 eliminates the effect of CII on the log (1/t) odor threshold and its predicted values differ from the experimental values, whereas in the case of TF2, which takes into account the influence of IIC on the odor threshold, the experimental odor threshold values are closer to the calculated values (Fig. 5). The (T, Nepoch) values for both TF1 and TF2 models demonstrate that the threshold and number of epochs are not identical.

Close

Fig. 5

Models TF1 and TF2 were used to calculate the olfactory correlation thresholds of experimental odorous compounds.

Export to PPT

3.2 Hologram-QSPR (HQSPR)

HQSPR is a new QSPR technique that avoids many problems associated with conventional 3D-QSPR approaches. Only structures and properties are required as input, no complex process of descriptor selection or 3D molecular alignment is required. HQSPR converts molecules in a data set into numbers of their constituent fragments. The fragment counting patterns from the molecules in the dataset is then linked to the targeted biological activity data using the PLS partial least squares analysis. Both steps, fragment counting, and PLS analysis are very fast. Nevertheless, the method is robust and highly predictive for many data sets. holograms 53, 59, 61, and 71 were chosen with a fragment number varying between 4 and 7 to construct all possible combinations of the Atoms, Bonds, Connections, Hydrogen Atoms, Chirality, and number of donor/acceptor hydrogens that are represented in Table 7.

A total of 63 models obtained using HQSPR shows that the Bonds and Connection fragments are explanatory fragments (comparison between the following models: 12, 32, 33 and 63) which generated the olfactory threshold prediction on the other hand we find that no effect of Hydrogen fragments Atoms and Chirality as a result of many of the donor and acceptor hydrogens having a negative effect on the statistical values of the prediction of the targeted biological property at cause of absence of donor groups, asymmetric carbons, enantiomers, and diastereomers except that two conformations E (M67) and Z M(68) among the 78 odorant molecules. Model 12 (Bonds/Connections) was chosen to develop the number of atoms connecting for each fragment, the set of results obtained present in Table 8.

The change in the number of fragment connections demonstrates what we previously stated; on the other hand, our choice of starting point, minimum and maximum number of atoms in a fragment yields good results, as shown in the picture below (Fig. 6).

Close

Fig. 6

The influence of the minimum and maximum number of atoms has an impact on the model performance (B/C).

Export to PPT

3.2.1

3.2.1 PLS HQSPR-Hologram

In this work, a PLS technique was used to compare the prediction values of each Hologram (53, 59, 61, 71, and AVG) and to establish the best fragment of the HQSPR predictive models represented in Table 9.

Table 9 Experimental and calculated values by different Hologram models.

	N	pt	59	71	61	53	AVG		N	pt	59	71	61	53	AVG
Training set	M1	3.52	3.76	4.81	4.60	4.39	4.39		M57	10.22	9.06	9.14	9.12	9.14	9.11
	M2	4.52	5.17	5.39	5.12	5.49	5.29		M58	7.40	8.22	8.24	8.36	8.43	8.31
	M3	5.40	5.59	5.84	5.83	5.57	5.71		M59	11.20	11.20	10.90	10.94	10.89	10.98
	M4	6.52	6.06	6.26	6.49	6.08	6.22		M60	10.52	9.69	9.65	9.97	9.66	9.74
	M5	6.40	7.16	7.05	6.99	6.59	6.95		M61	10.40	10.70	11.48	10.89	10.98	11.01
	M6	8.30	7.46	7.39	7.31	7.00	7.29		M62	11.10	11.43	11.24	11.57	11.52	11.44
	M8	7.00	6.77	6.68	6.73	6.22	6.60		M63	10.35	10.27	10.07	10.02	10.24	10.15
	M11	6.40	7.00	5.68	6.03	6.19	6.23		M64	10.92	11.46	11.11	11.10	10.80	11.12
	M14	7.70	8.24	7.54	8.14	7.87	7.95		M65	11.22	10.28	10.55	10.81	11.16	10.70
	M15	6.52	6.41	7.44	7.18	6.98	7.00		M66	10.62	10.09	10.30	10.00	10.13	10.13
	M16	7.40	6.62	7.04	7.15	7.51	7.08		M67	9.89	9.67	9.53	10.17	9.29	9.66
	M17	7.10	7.28	6.67	6.31	6.85	6.78		M68	9.30	9.67	9.53	10.17	9.29	9.66
	M18	7.70	7.55	7.60	7.59	7.67	7.60		M69	7.10	6.11	6.08	5.88	5.83	5.98
	M20	7.22	7.55	7.60	7.59	7.67	7.60		M71	7.70	7.67	7.57	7.71	7.76	7.68
	M22	7.80	7.66	7.85	7.57	7.90	7.74		M72	10.10	9.51	9.36	9.22	9.44	9.38
	M24	6.30	6.45	6.40	6.22	6.99	6.51		M73	8.70	8.41	8.20	8.26	8.31	8.30
	M25	7.22	6.79	6.79	6.84	7.26	6.92		M75	7.52	6.91	6.46	6.55	6.85	6.69
	M27	7.70	8.01	7.86	7.97	7.50	7.83		M76	6.70	6.67	6.31	6.75	6.46	6.55
	M28	6.70	6.54	6.76	6.76	6.57	6.66		M77	7.30	8.19	8.47	8.13	8.16	8.24
	M31	9.00	8.93	9.19	9.09	9.42	9.16		M78	7.16	6.23	6.42	6.41	6.33	6.35
	M32	9.92	10.17	10.30	10.07	10.26	10.20	Test set	M7	6.70	7.19	7.12	7.11	6.58	7.00
	M33	9.16	9.06	9.14	9.12	9.14	9.11		M9	6.40	6.35	6.23	6.35	5.87	6.20
	M34	7.70	8.22	8.24	8.36	8.43	8.31		M10	5.96	5.50	5.34	5.59	5.17	5.40
	M35	10.33	10.09	10.30	10.00	10.13	10.13		M12	7.77	8.24	7.54	8.14	7.87	7.95
	M36	6.05	6.11	6.08	5.88	5.83	5.98		M13	8.30	7.80	6.96	7.35	7.13	7.31
	M37	7.16	7.15	7.04	6.78	6.98	6.99		M19	7.92	7.12	6.92	6.71	6.94	6.92
	M38	7.22	7.67	7.57	7.71	7.76	7.68		M21	6.40	7.03	6.76	7.20	7.23	7.06
	M39	8.40	9.22	9.02	8.90	9.04	9.04		M23	6.10	6.50	5.89	5.30	6.65	6.08
	M40	9.00	9.51	9.36	9.22	9.44	9.38		M26	7.80	7.19	7.57	7.80	8.20	7.69
	M41	8.70	8.41	8.20	8.26	8.31	8.30		M29	8.40	7.88	7.85	7.83	7.70	7.81
	M42	6.92	7.56	7.31	7.50	7.61	7.49		M30	7.40	8.33	8.51	8.57	8.58	8.50
	M43	6.40	6.91	6.46	6.55	6.85	6.69		M48	6.52	6.23	6.42	6.41	6.33	6.35
	M44	6.52	6.67	6.31	6.75	6.46	6.55		M50	8.16	7.88	7.85	7.83	7.70	7.81
	M45	7.05	6.94	7.36	7.15	7.31	7.19		M52	9.92	8.93	9.19	9.09	9.42	9.16
	M46	8.00	8.19	8.47	8.13	8.16	8.24		M54	10.70	10.17	10.30	10.07	10.26	10.20
	M47	7.09	7.08	7.31	7.17	7.03	7.15		M55	10.16	9.91	10.03	9.88	9.84	9.91
	M49	6.40	6.54	6.76	6.76	6.57	6.66		M56	10.59	9.49	9.58	9.50	9.49	9.51
	M51	8.00	8.33	8.51	8.57	8.58	8.50		M70	9.10	7.15	7.04	6.78	6.98	6.99
	M53	10.30	9.88	9.95	9.75	9.86	9.86		M74	7.22	7.56	7.31	7.50	7.61	7.49

Following comparing the predicted values of each Hologram, we can see that Hologram 59 produces excellent results; we will apply this model to continue research investigation, and we will show the correlation between experimental and calculated values in the image below (Fig. 7).

Close

Fig. 7

Residual graphs between experimental and predicted for HQSPR models.

Export to PPT

3.2.2

3.2.2 HQSPR contribution map

The QSPR analysis focused on the contributions of biologically owned fragments of odorous molecules for each individual compound. The graphical results of the atomic contribution maps of the M65 compounds with the highest olfactory threshold value are represented in the color-coded structure diagram in Fig. 8.

Close

Fig. 8

The HQSPR model contribution map on the M65 compound, and the coefficient of influence of each contribution.

Export to PPT

The color of each atom implies its contribution to the olfactory threshold of the compound. The green color (green–blue, and yellow) represents a positive effect while the red color (red, red–orange, and orange) reflects the negative contribution to the olfactory threshold, and the intermediate contributions are white.

3.3 CoMSIA statistical results

Indicates Comparative Molecular Similarity Indices Analysis (CoMSIA), a statistical method used in molecular biology and computational chemistry to study the relationship between chemical structures and biological activity. The results summarized in Table 10 possibly reflect the results of different combinations of molecular descriptors used as inputs to the CoMSIA model, and their impact on the accuracy and robustness of the model. The goal of these tests is probably to determine the optimal combination of descriptors that provide the best prediction of biological activity for a given set of molecules.

The combination of different 3D-CoMSIA fields with their results is a technique used to improve the performance and reliability of QSPR/QSAR models. This approach combines multiple molecular descriptors and multiple modeling algorithms to produce a more accurate and robust model. The search for this combination depends on the optimal number of components (N) choice, most of the later works use the standard Sybyl software parameters, but in our work, we studied the optimal number of component variation with respect to other statistical parameters to determine the relationship between them. The results obtained are represented in the Table 11.

In this part, we determine the optimal number of components N of the 3D-QSPR model using the Comparative molecular similarity indices analysis (CoMSIA/SEH) method studied, we choose the parameters $Q^2$ , $R_{\mathit{cv}}^2$ , $R_{\mathit{ncv}}^2$ , $R_{\mathit{bs}}^2$ , SEE and F to find the optimal number of components, It can be seen that the f-test value is proportional to the number of optimal components in a linear way, but the other parameters $Q^2$ , $R_{\mathit{cv}}^2$ , $R_{\mathit{bs}}^2$ and SEE are proportional to the number of optimal components in a logarithmic way, N = 6 maximum threshold of $Q^2$ , N = 10 maximum threshold of ( $R_{\mathit{ncv}}^2$ , and $R_{\mathit{bs}}^2$ ) and N = 13 maximum threshold of SEE, according to the maximum value of $Q^2$ , we choose the optimal number of components N = 6 (Fig. 9).

Close

Fig. 9

Graphical representation of each parameter studied as a function of the number of optimal components.

Export to PPT

Upon determining the best number of components to be N = 6, it was found that the CoMSIA/SEH model can be trusted to give an explanation and prediction of the olfactory threshold of pyrazine derivative odor molecules. The results of the prediction and the comparison between the calculated and observed olfactory threshold are presented in Table 12 and Fig 10.

Close

Fig. 10

Graph of observed and calculated properties for training set and test set process based on the CoMSIA model.

Export to PPT

3.3.1

3.3.1 CoMSIA contour map analysis

The molecule M65 was chosen to display the contribution of steric, electrostatic and hydrophobic fields represented in Fig 11:

Close

Fig. 11

Contribution map of CoMSIA/SEH model on the compound M65.

Export to PPT

- Steric: the steric congestion between the R1 ((CH₂)₃CH (CH₃)₂) and R2 (OCH₃) substitutions of the M65 molecule promotes a steric field of 80% in green and 20% unfavorable in yellow shown in Fig. 11. The volume occupied by the methoxy group in space produces steric genes with branched substitutions more than linear substitutions, this explains that the olfactory threshold of compound M65 (log (1/t) = 11.222) smaller than M58 (log (1/t) = 7.398, R1 (C₁₀H₂₁), R2 (OCH₃)), M54 (log (1/t) = 10.699, R1 (C₅H₁₁) and R2 (OCH₃)). On the other hand, the favorable steric field clutter of hydrogen atoms in positions R3 and R4 is too small. This steric resistance is often used to propose new odor molecules.

- Electrostatic: a favorable electropositive field in blue color around the atoms of nitrogen of the common structure of pyrazines, the non-binding doublet of nitrogen acts as a nucleophile that will be attacked by an electrophilic group, for example reduction of pyrazine. however, we worked on pyrazine derivatives, so proposing new molecules cannot be based on this field. As well as the electron field near the oxygen atom of the methoxy group in position R2, the presence of electron substitution reduces the olfactory threshold of odor molecules. In general, the results demonstrated that the olfactory threshold of molecules containing the oxygen atom in position R2 is lower than the olfactory threshold of molecules containing the sulfur atom because of electronegativity (electron field), as the electronegative field and electronegativity are inversely proportional.

- Hydrophobic: 80% favorable, in yellow, corresponding to a hydrophobic grouping around positions R2 and R4. The + M mesomeric effect of methoxy substitution in R2 makes it possible to delocalize the electrons in the conjugate system P-σ-π, it directs the position para (R4). Which means that the addition of a hydrophobic group in this position lowers the olfactory threshold. On the other hand, 20% unfavorable, in white color, corresponds to a hydrophilic grouping around R1 and R3. Usually, alkyls are + I donor groups, they also orient the para (R3) position. So, the addition of a hydrophilic group in position R3 will increase the olfactory threshold.

3.4 Topomer CoMFA results

The chosen way of cutting the fragments has a considerable influence on the model's quality. In this work, we selected fragments R1, R2, R3, and R4 to investigate the effect of each substitution on the common structure of pyrazine, and we used the compound M65 as a template to display fragment positions. After severing the fragments, we discovered that molecule M1 was removed due to the lack of at least one heteroatom attached to the common pyrazine structure; thus, the spreadsheet must also include the relevant input structures and activity data. Furthermore, the computation of activity prediction on these molecules utilizing the current Topomer CoMFA model. The only format available for saving Topomer CoMFA templates and results is an SYBYL (.tbl) table file, as shown in Table 13.

Analysis of the results obtained by the Topmer CoMFA model shows that the mean coefficient of the steric and electrostatic fields of the hydrogen atom depends on fragments such as the mean coefficients of orders of 0.13, 0.09, 0.44 and 0.41 respectively corresponding to the fragments follow R1, R2, R3, and R4. In this case, the steric clutter between the two adjacent fragments R3 and R4. This also causes to restrict the torsion of the common pyrazine structure. In addition, the comparison of the different substituents in the same fragment R2 revealed that the olfactory threshold of molecules that contains the oxygen atom is smaller than those that contain the sulfur. Their mean steric and electrostatic coefficients were in the order of 2.91, 2.15, 1.14 for substituents -SCH₃, -SC₂H₅, and -SC₆H₅ and 3.19, 2.30, 1.30 for substitutions –OCH₃, -OC₂H₅, -OC₆H₅ respectively. These results are in agreement with those obtained by the CoMSIA/SEH model whereas the oxygen atom more electronegative than sulfur. Although, the mean coefficients of the oxygen atom larger than the sulfur atom, for this the olfactory threshold of molecules that contains the oxygen atom smaller than the substitution that contains the sulfur atom.

By comparing the predicted values and the experimental, we can see that this model produces excellent results; we will apply this model to further research, and we will show the correlation between the experimental and calculated values in the Fig 12.

Close

Fig. 12

correlation of experimental and predicted values of the Topomer ComFA model and their residues.

Export to PPT

3.4.1

3.4.1 Topomer CoMFA contour map analysis

The molecule M65 was chosen to display the contribution of the steric and electrostatic fields represented in Fig. 13.

Close

Fig. 13

Contour Map Analysis of the Topomer CoMFA Model (A, B, C, and D correspond to the electrostatic fields R1, R2, R3, and R4, respectively) and (E, F, G, and H correspond to the Electrostatic field R1, R2, R3, and R4 respectively).

Export to PPT

The generation of contour maps obtained by the Topomer CoMFA model is based on the R1, R2, R3, and R4 fragments, the Topomer CoMFA model shows that it is a significant detection model, Fig 13 represents the contour maps generated around each fragment.

- Electrostatic field: the contribution of the electrostatic field of 80% favorable in blue color around fragment R4, Fig. 13 (D) and 20% unfavorable in red color around fragments R2 and R3, Fig. 13 (B and C) indicate that the addition of an electropositive group in fragment R4 and of an electropositive group in fragments R2 and R3 can be decreases the olfactory threshold of the odorous molecules.

- Steric field: The 80% favorable green steric contribution around the R1, R2, R3, and R4 fragments, represented in Fig. 13 (E, F, G, and H), indicates that adding a bulky group to these fragments can reduce the odor threshold of odorant molecules.

3.5 Validation of prediction models

External Validation Plus 1.2 is a tool that computes all of the essential external validation parameters https://teqip.jdvu.ac.in/QSAR_Tools/, as well as judging the performance of prediction quality of a QSPR model using MAE-based criteria. It can also check for the presence of systematic mistakes in the model.

All calculated external validation parameters (shown in Table 14), along with some basic information about the data structure, are necessary for assessing the quality of QSAR model predictions. We have selected an option regarding the output file will also indicate whether systematic error is present or not, designated as “present” or “absent” based on the 5 conditions mentioned before. It summarizes the information again, all validation parameters are necessary to evaluate the performance of the quality predictive of each QSPR model. Based on these results, it can be observed that the HQSPR and Monte-Carlo (TF2) models are good predictive models.