Introduction
Imidazoline and α-adrenergic receptor ligands are mainly imidazoline and guanidine derivatives with different pharmacological effects and target tissues [1]. The pharmacological effects of these ligands are generally the results of interaction with three types of imidazoline receptors (I1-IR, I2-IR, and I3-IR) and α2-adrenoreceptors (α2-AR) [2,3]. Clonidine is a well-known first-generation centrally acting antihypertensive drug that lowers blood pressure by stimulating α2-AR and I1-IR. On the other hand, moxonidine and rilmenidine, which belong to the second generation of antihypertensives, act mainly at the I1-IR and produce fewer α2-adrenoceptor-mediated side effects [3-5]. Idazoxan and its analogs have a high affinity for the I2-IR whose involvement in different psychiatric disorders, opiate withdrawal, pain, and Parkinson’s and Alzheimer’s diseases has been investigated [6-8]. Efaroxan is able to induce insulin secretion from pancreatic β-cells by activating I3-imidazoline receptors (I3-IR) [9].
Considering the chemical structures of the above-mentioned ligands, the design of new compounds as drug candidates for the treatment of various cardiovascular or neurological diseases will be directed to compounds that contain one or more basic centers in their chemical scaffold and have the ability to ionize. For ionizable compounds, the extent of ionization together with lipophilicity is of great importance for their pharmacokinetic properties, such as absorption, distribution, metabolism and elimination (ADME), as well as for the drug-receptor interaction [10].
In addition, the extent of ionisation of compounds at different pH values has a strong effect on their behaviour in different separation systems, such as chromatographic and electrophoretic systems. As the degree of ionisation of basic compounds increases, their mobility in capillary electrophoresis increases, while the retention of unionised form of the analyte can be 10–20 times greater than that of the corresponding ionised form for a given mobile phase composition. Currently, high-performance liquid chromatography (HPLC) and capillary electrophoresis (CE) are very useful techniques for the rapid and efficient separation as well as examination of the individual imidazoline and α-adrenergic receptor ligands [11-14]. Therefore, accurate quantitative relationships between chromatographic retention or electrophoretic mobility and pH, which can be established through quantitative structure-retention relationships (QSRR) and quantitative structure-mobility relationships (QSMR) modelling can be very useful for predicting the retention behaviour and electrophoretic mobility of the new compounds, understanding the separation mechanism in a given analytical system and identifying important molecular descriptors [15]. In addition, the processes of drug action are believed to share many similarities with the processes underlying chromatographic separations, as the similar elementary intermolecular interactions are important for the behaviour of chemical compounds in both biological and chromategraphic environments [16]. Therefore, the chromatographic retention constant log Kw has been successfully used to evaluate the lipophilicity of compounds [17] as well as in correlation studies with various physicochemical and biological properties to establish appropriate QSPR (quantitative structure-property relationship) and QSAR (quantitative structure-activity relationship) models [18,19]. Chromatographic and electrophoretic data have also been successfully used for the pharmacological classification of compounds and prediction of their ADME properties [20]. HPLC was used to determine the effective permeability of 40 compounds (18 IRs/a-ARs ligands and 22 CNS drugs) using the parallel artificial membrane permeability assay (PAMPA) to predict brain penetration [21]. The electrophoretic behaviour of a compound is indicative of its acid-base properties. i.e., pKa values and molecular mass, since the charge-to-mass ratio of compounds largely determines their electrophoretic mobility. Thus, all these models, QSPR, QSRR, QSAR and QSMR, can provide us with valuable information about the activity, physicochemical and ADME properties of drug candidates and could, therefore, significantly reduce the number of drug candidates synthesised and tested in the early phase of drug discovery.
The main objective of this work was to investigate and compare the mechanism of chromatographic retention and electrophoretic mobility of 29 compounds with main or side effects to alpha adrenergic/imidazoline receptors. The study was performed at different ionization levels, including physiological pH, and applied to QSRR and QSMR models using multiple linear regression (MLR) and partial least squares regression (PLS). Rigorous validation criteria as proposed by Golbraikh and Tropsha [22,23] and Roy [24], were also applied as an essential part of QSRR/QSMR model development. The developed QSRR/QSMR models can be used for the fast and reliable prediction of the chromatographic and electrophoretic behavior of structurally related imidazoline and alpha adrenergic receptor ligands. In addition, the usefulness of chromatographic and electrophoretic data for predicting permeation across the blood-brain barrier (BBB) was also examined.
Experimental
Reagents and chemicals
Methanol (HPLC grade, 99.8 %) and glacial acetic acid (≥99.7 %) were purchased from J.T. Baker (Deventer, Netherlands). Ammonium acetate, sodium hydroxide, sodium dihydrogen phosphate, and disodium hydrogen phosphate of ACS reagent grade were purchased from Merck (Darmstadt, Germany). Boric acid (ACS reagent, 99.5 to 100.5 %) was purchased from Sigma-Aldrich (St. Louis, MO, USA). Ammonium hydroxide solution 30 % (analytical grade) was purchased from Carlo Erba (Milan, Italy). Aqueous solutions of background electrolytes and buffers for the preparation of mobile phases were prepared with water (HPLC grade).
Clonidine hydrochloride (≥98 %), moxonidine hydrochloride (≥98 %), guanfacine hydrochloride (≥98 %), brimonidine hydrochloride (≥98 %), efaroxan hydrochloride (≥98 %), idazoxan hydrochloride (≥98 %), rilmenidine hemifumarate (≥98 %), harmane (98 %), harmine hydrochloride (≥98 %), tizanidine hydrochloride (≥98 %), triamterene (≥99 %), clopamide (≥98 %), indapamide (≥99 %), naphazoline hydrochloride (≥98 %), xylometazoline hydrochloride (≥99 %), tetrahydrozoline hydrochloride (≥98 %), oxymetazoline hydrochloride (≥99 %), ephedrine hydrochloride (99 %), pseudoephedrine hydrochloride (99 %), maprotiline hydrochloride (>99 %), tamsulosin hydrochloride (≥98 %), mianserin hydrochloride (≥98 %), carvedilol (≥98 %), clozapine (≥98 %) and olanzapine (≥98 %) were purchased from Sigma–Aldrich (St. Louis, MO, USA). Tramazoline hydrochloride (≥98 %) and doxazosin mesilat (≥98 %) were obtained from Zdravlje-Actavis (Leskovac, Serbia). Amiloride hydrochloride (≥98 %) was kindly donated by Galenika (Belgrade, Serbia). Phenilephrine hydrochloride (98 %) was obtained from Ivančić i sinovi (Belgrade, Serbia).
HPLC conditions
Solutions for HPLC analysis
The working solutions were prepared by dissolving the substances in methanol to obtain a concentration of 0.7 mg/mL for rilmenidine and 0.1 mg/mL for the other compounds.
HPLC equipment and working conditions
The HPLC working conditions are described in our previously published work [20]. HPLC analysis was performed using an Agilent Technologies 1200 (Wilmington, DE, USA) liquid chromatography system equipped with an online degasser, a binary pump, a column oven and a photodiode array detector. The data was recorded and analyzed using Agilent's ChemStation software (Chemstation for LC 3D system (Rev. B. 02. 01-SR2 [260]).
Separation was performed on the XTerra® RP18 column, 4.6×100 mm, particle size 3.5 μm (Waters Corporation, Milford, MA, USA) using three different constant ionic strength buffer solutions (I = 25 mmol/L) prepared at pH values 4.4, 7.4 and 9.1, as the aqueous component in the mobile phase while methanol added as an organic modifier. All measurements were performed at 25°C with a constant flow rate of 0.8 mL min-1 and UV detection in the range of 200-280 nm. The injection volume was 20 μL. Retention times were determined in the isocratic elution mode using at least six mobile phases of methanol/buffer (pH 4.4, 7.4 or 9.1) with a methanol concentration of 75-2 % depending on the retention properties of compounds. The retention factor k was calculated for each mobile phase composition according toEquation (1):
where tr is the retention time of the tested compound and t0 is the dead volume measured with KNO3 as a non-retentive marker. The values corresponding to 100 % of the buffered eluent (log kw) were determined by extrapolation according toEquation (2):
where log kw is the intercept, S is the slope of the regression line, and φ is the volume fraction of the organic modifier in the mobile phase. Log kw of the analyte at three different experimental conditions, methanol/buffer pH 4.4, methanol/buffer pH 7.4 and methanol/buffer pH 9.1 (log kw4.4, log kw7.4 and log kw9.1, respectively) were calculated for all investigated compounds and used for the QSRR analysis.
Electrophoretic conditions
Solutions for CE analysis
Stock solutions of harman, harmine, clopamide, indapamide, rilmenidine, mianserin, doxazosin, carvedilol, clozapine, olanzapine and triamterene were prepared in methanol (c = 1 mg/mL) and diluted with water to a final concentration of 5.8 μg/mL for triamterene, 60 μg/mL for rilmenidine and 30 μg/mL for the other substances. Brimonidine was dissolved in 0.1 % formic acid, while the remaining substances were dissolved in water. Acetone 2 vol.% was used as electroosmotic flow (EOF) marker. Depending on UV responsiveness and solubility, samples were prepared at different concentrations, from 5.8 to 60 μg/mL.
CE equipment and working conditions
The CE working conditions are described in our previously published paper [20]. All experiments were performed using the SpectraPhoresis 500-capillary electrophoresis system (Spectra Physics Analytical, USA) with UV detector. Data were recorded and analyzed using ChromQuest software version 4.0 (Thermo Finnigan, USA).
The separations were performed with an uncoated fused capillary (total length 31.5 cm, 50 μm id, effective length 24 cm, Polymicro Technologies, USA) applying a voltage of 11 kV and at 25 °C. The samples were injected hydrodynamically and detection was performed in the 200-280 nm range. The new capillary was rinsed stepwise with 0.1 M NaOH (15 min), water (10 min) and a running buffer (10 min). Finally, a 10 min high voltage was applied through the background electrolyte-filled capillary to equalize the inner surface, stabilize the electroosmotic flow and ensure good reproducibility of injections from run to run. Background solutions with constant ionic strength (I = 25 mmol/L) were prepared at three different pH values by mixing the corresponding amounts of CH3COOH/NaOH, NaH2PO4/NaOH and H3BO3/NaOH for pH 4.4, 7.4 and 9.1, respectively. All compounds were injected in three replicates at each pH (4.4; 7.4 and 9.1). Between runs, the capillary was rinsed with the background electrolyte for 1 minute. The effective electrophoretic mobility, μeff, of the analyte at three different pH values 4.4; 7.4 and 9.1 (μeff4.4, μeff7.4 and μeff9.1, respectively) is expressed in cm2 / V s and was calculated using theEquation (3):
where Lt is the total length of the capillary (cm), Leff is the effective length of the capillary (cm), U is the applied voltage (V), tm is the migration time of the analyte (s) and teof is the migration time of the electroneutral EOF marker acetone (s). The calculated μeff values (μeff4.4, μeff7.4 and μeff9.1) were used as dependent variables in QSMR analysis.
Theory/calculation
Computational methods
Clonidine, moxonidine, tizanidine, brimonidine, rilmenidine, and tramazoline are cyclic guanidines or amidine structures that can exist in two main tautomeric forms (amino and imino). The stability of their amino and imino tautomers was investigated at the B3LYP/6–31G(d, p) level of density functional theory-(DFT) [25,26] using the Gaussian 98 program [27]. Based on the Self Consistent Field Energy calculated for neutral amino and imino tautomers of the compounds, it was shown that clonidine, moxonidine, tizanidine, brimonidine and tramazoline exist as more stable imino tautomers, while the predominant tautomeric form of rilmenidine is the amino form. The calculation of pKa and the selection of the dominant molecule/cation species under experimental conditions (pH 4.4, 7.4 and 9.1) were performed for 29 α-adrenergic and imidazoline receptor ligands using the Marvin 5.5.1.0 ChemAxon program [28].
Structural optimization (selected tautomers and molecules/cations species at three different pH values) was performed using the B3LYP/3–21(d,p) levels of the DFT in the Gaussian 98 program [27]. Molecular descriptors were calculated for all investigated compounds in the programs Dragon [29], Gaussian 98 (B3LYP/3-21G(d,p) basis set) [27], Marvin 5.5.1.0 ChemAxon [28], and Chem3D Ultra 7.0.0 [30]. In addition, descriptors such as hardness (η), global softness, chemical potential (μ), electronegativity (χ) and electrophilicity index (ω), which have been extensively used to interpret various aspects of chemical bonding and reaction mechanism, were included in the set of calculated properties [31]. The partial least squares regression analysis started with 3445 computed descriptors. The descriptors with constant values were excluded from the modelling, as were descriptors with an intercorrelation greater than 0.996. For a pair of correlated descriptors, the one with the better correlation to Y was retained, while the other was excluded from the analysis. In the case of the MLR analysis, the cut-off value for the intercorrelation was set at 0.9 [32].
Quantitative structure-retention relationship and quantitative structure-mobility relationship modelling QSRR and QSMR modelling studies were performed for 29 imidazoline and α-adrenergic receptors ligands using PLS and stepwise MLR statistical methods at pH values of 4.4, 7.4 and 9.1 and on both HPLC and CE systems. An overview of the dataset was performed for each pH value and each chromatographic or electrophoretic system was examined in order to identify and exclude possible outliers before modelling. For this purpose, the scatter plots in the Simca P 12+ program were used [33], in which the pool of calculated molecular descriptors and the response variables are organised as an X-matrix. The scatter plot t1 vs. t2 is a window in X-space, showing how the X observations relate to each other, where t1 and t2 are scores (one vector for each model dimension) representing new variables computed as linear combinations of X. The score t1 (first component) explains the largest variation of the X space, followed by t2 and so on. For a two-dimensional score plot, the tolerance ellipse is presented based on Hotelling's T2. Observations that lie outside the ellipse are outliers. The t-u plot (t1 vs. u1) shows the relationship between X and Y. In addition, the degree of fit (a good fit corresponds to a small scatter around the straight line), indications of curvature and outliers can also be seen. The u score (one vector for each model dimension) is new variable summarizing Y, so as to maximize the correlation with the scores t. Based on the plots (t1 vs. t2 and t1 vs. u1), the remaining compounds were distributed among the training and test sets such that each chemical group (i.e., guanidine, 2-aminoimidazoline, 2-arylmethylimidazoline, phenylethylamine etc.) has a representative in the test set and the log kw4.4, log kw7.4, log kw9.1, μeff4.4, μeff7.4 and μeff9.1 of the compounds in the test set are homogeneously distributed over the entire range of values of log kw4.4, log kw 7.4, log kw9.1, μeff4.4, μeff7.4 and μeff9.1, respectively. The selected training and test sets contained at least 20 compounds used for model building and at least 5 compounds used for external validation. In order to perform a meaningful comparison between the different methods (MLR and PLS) used to develop the QSRR models, the same training and test set was used for each pH value. Therefore, six independent training and test sets were created, one for each dependent variable (log kw4.4, log kw7.4, log kw9.1, μeff4.4, μeff7.4 and μeff9.1). The predictive power and robustness of the models were estimated using the cross-validated parameter R2 (Q2), which was calculated fromEquation (4):
InEquation (4), PRESS is a parameter calculated using the leave-one-out (LOO) approach where each compound from the training set was deleted once and a new model was created with the remaining compounds and used to predict the Y-value of the deleted compound. The procedure was repeated until all compounds had been deleted once. For all created models, the squared sum of the differences between observed and LOO-predicted values (e(i)) (PRESS) was calculated according toEquation (5):
The models with Q2 ≥ 0.5 can be considered to have good predictive capability [34]. The root mean squared error of estimation (RMSEE) was calculated to characterize the predictive ability of the models for the compounds in the training set,Equations (6):
where n is the number of compounds in the training set, while Yobs(training) and Ypred(training) denote the experimental and predicted values for the compounds in the training set.
In addition to internal validation, the predictive power of all created models was estimated by external validation using the following parameters inEquations (7) and(8):
where RMSEP is the root mean square error of prediction, n is the number of compounds in the test set, while Yobs(test) and Ypred(test) are the experimental and predicted values for the compounds in the test set.
where Ypred(test) and Yobs(test) are the predicted and observed values of the dependent variables of the compounds in the test set, respectively, and Y̅training indicates the mean values of the dependent variables of the compounds in the training set. For a predictive QSRR/QSMR model, the R2pred value should be higher than 0.5 [35].
The rigorous external validation parameters proposed by Roy (
) [24] and Golbraikh and Tropsha [22,23] were also applied to all created models and were used to select the most reliable models for the prediction of chromatographic and electrophoretic behaviour at each pH value.
According to Ojha and Roy [35], 
is validation metric calculated according toEquation (9) using observed (y-axis) and predicted (x-axis) values:
where r2 is the determination coefficient for the least squares regression line (with intercept) correlating observed (y-axis) and predicted (x-axis) values and 
is the determination coefficient for the least squares regression line (without intercept) correlating observed (y-axis) and predicted (x-axis) values, 
is validation metric calculated according toEquation (10) using observed (x-axis) and predicted (y-axis) values
where 
is the determination coefficient for the least squares regression line (without intercept) correlating observed (x-axis) and predicted (y-axis) values. 
is average of 
and 
, while 
is the absolute difference between 
and 
.
Models were considered acceptable if they met all of the following conditions:
the values of
and 
should be close to each other and 
should be greater than 0.5;
should be less than 0.2 [24].
According to criteria proposed by Golbraikh and Tropsha, following conditions should be achieved: (i) Q2 > 0.5; (ii) R2 > 0.6, (iii) [(R2 - R02) / R2] < 0.1 or [(R2 - R'02) / R2] < 0.1; (iv) 0.85 ≤ k ≤ 1.15 or 0.85 ≤ k' ≤ 1.15 where Q2 is the cross-validated correlation coefficient calculated for the training set, but all other criteria are calculated for the test set [22,23]. R02 is the coefficient of determination (predicted versus observed values) for regressions through the origin. R'02 is the coefficient of determination (observed versus predicted values) for regressions through the origin. K and K′ are slopes of regression lines through the origin.
Multiple linear regression
Stepwise multiple linear regression was performed with the program STATISTICA [36]. A stepwise multiple linear regression procedure was used to select the molecular descriptors with the greatest influence on the dependent variables. The models were created stepwise according to the following procedure: identification of an initial model, then iterative „stepping“ in which descriptors were added or removed from the model based on the previously established stepping criteria (F to enter and F to remove) according to their statistical significance (F-test ), and finally termination of the search when no more stepping was possible or when a specified maximum number of steps was reached. Stepping criteria were set to a narrow range, and the predictive abilities of the models were tested by leave-one-out cross-validation and external test set prediction methods to avoid obtaining an overfitting model.
Partial least squares regression
The PLS analysis was carried out using the SIMCA P+ 12.0 (soft independent modelling of class analogy) program [33]. The PLS approach is able to analyse data with noisy, collinear and even incomplete independent variables. In PLS modelling, the variable importance in the projection (VIP) parameter was used to assess the significance of the independent variables in the course of PLS model building. The VIP parameter represents a summary of the importance of each variable (X) for both the Y and X matrices. Molecular descriptors with VIP values greater than 1 are most important in explaining the regression model, those with 1.0 > VIP > 0.5 have a moderate influence, while the independent variables with VIP values less than 0.5 are not relevant for the model [33,37]. Therefore, during model generation, the descriptors with the lowest VIP values were successively removed from the model, while statistical parameters such as the squares of the multiple correlation coefficients R2, Q2 (a cross-validated version of R2), the root mean square error of estimation for the training set, and the root mean square error of prediction for the test set were calculated for each new PLS model and compared with the previous one. The procedure was repeated until the best models for the investigated experimental conditions (pH 4.4, 7.4 and 9.1) and analytical systems (HPLC and CE) were formed.
The response permutation test (Y scrambling), as a measure of model overfitting, was used to examine the statistical significance of the R2 and Q2 values and overfitting due to the chance correlation [37]. In this test, the Y-variables were randomly re-ordered 100 times while the X-matrix was remained unchanged. The PLS model was fitted to the permuted data and the new parameters R2, Q2 and VIP were calculated. All model selection steps were repeated with the scrambled Y response data. Regression lines were fitted through the R2 and Q2 values to obtain two separate intercepts. The values of the obtained intercepts for valid QSRR and QSMR models should be below 0.05 for the Q2 intercept and not above 0.4 for the R2 intercept [37].
Results and discussion
The influence of the degree of ionization (under acidic, neutral and basic conditions) of 29 imidazoline and/or α-adrenergic receptor ligands (Figure S1, Supplementary material) on their retention in reversed-phase liquid chromatography and electrophoretic mobility in capillary electrophoresis was investigated. A pH value of 4.4 was chosen for the work in an acidic medium at which the majority of the tested compounds reached the highest degree of ionization, pH = pKa - 2 (based on calculated pKa values in [28] and experimentally determined pKa values found in the literature,Table S1, Supplementary material). For operation in a neutral medium a pH value of 7.4 was used, which simulates physiological conditions, while a pH of 9.1 was chosen as a compromise where 9 compounds move together with the electroosmotic flow marker and the effective mobility is zero such as e.g. brimonidine (experimental pKa determined by potentiometry 7.50), clozapine (experimental pKa determined by HPLC 7.719 ± 0.015), doxazosine (experimental pKa determined by voltammetry 6.89 ± 0.57), guanfacine (calculated basic pKa 8.65), harman (experimental pKa determined by HPLC 7.21), mianserin (experimental pKa determined by potentiometry 7.40), moxonidine (experimental pKa determined by HPLC 7.84), tizanidine (calculated pKa 7.49) and triamterene (experimental pKa determined by pH spectrophotometric method 7.16) (Tables S1 and S4, Supplementary material); about 1/3 of the compounds reach their mean mobility (pH close to the pKa value, amiloride (experimental pKa determined by potentiometry 8.72 ± 0.05), idazoxan (experimental pKa determined by HPLC 9.04), phenylephrine (experimental pKa determined by spectrophotometry 9.17), tamsulosin (calculated pKa 9.28), clopamide (calculated acid pKa 8.85), indapamide (experimental acid pKa determined by potentiometry 8.8 ± 0.2) (Tables S1 and S4, Supplementary material) and the remaining compounds with pKa values around and above 10 are still ionized to a high percentage, such as e.g. efaroxan (experimental pKa determined by HPLC 10.02), ephedrine (experimental pKa determined by spectrophotometry 9.65), maprotiline(experimental pKa determined by potentiometry 10.45 ± 0.02), naphazoline (experimental pKa determined by pH spectrophotometric method 10.81), oxymetazoline (experimental pKa determined by pH spectrophotometric method 10.62), tetrahyrozoline (experimental pKa determined by pH titration 10.51 ± 0.05), tramazoline (experimental pKa determined by pH titration 10.66 ± 0.05)(Table S1 and S4, Supplementary material).
Two different linear statistical methods (MLR and PLS) were applied and compared for each experimental condition and examined systems in order to identify the most reliable models enabling accurate prediction of dependent variables.
Quantitative structure-retention relationship modelling
Determination of log kw values
The retention behaviour in HPLC system was monitored by calculating extrapolated retention factors log kw corresponding to pure buffer as mobile phase. The log kw values are considered to be better lipophilicity indices compared to the isocratic log k, as their values are of the same order of magnitude as the log P / log D octanol-water [38]. The values of log kw (intercepts) and slopes (S) at all investigated pH values: 4.4, 7.4 and 9.1 are shown inTable S2 (Supplementary material). In an acidic environment (pH 4.4), the obtained log kw values are in the range from 0.09 (phenilephrine) to 3.75 (carvedilol). Under the given conditions, only two compounds (clopamide and indapamide) are present in neutral (molecular) form, while the majority of the other compounds are completely ionized (over 99 %) [28]. This is also the reason for the lower log kw values of the investigated compounds at pH 4.4 compared to pH 7.4 and 9.1, where the compounds are present in ionized form to a lesser extent or are predominantly non-ionized (pH 9.1). Therefore, as the pH increases, the lipophilicity of the compounds and, consequently, their log kw values increase (Table S2, Supplementary material).
Quantitative structure-retention relationship models
As a result of the QSRR studies applied to the log kw data obtained at three pH values (log kw4.4, log kw7.4 and log kw9.1), six models were created: MLR-QSRR (log kw4.4), PLS-QSRR (log kw4.4), MLR-QSRR (log kw7.4), PLS-QSRR (log kw7.4), MLR-QSRR (log kw9.1), and PLS-QSRR (log kw9.1). The score plot of the first two components (t1 vs. t2) at each pH shows a uniform data distribution in all four quartiles. Tamsulosin is in group with oxymetazolin, clopamide, indapamide, carvedilol, doxazosin and oxymetazoline. These compounds are characterized by the highest descriptor values with tamsulosin lying outside or on the Hoteling T2 ellipse. Thus, tamsulosin was identified as an outlier in all chromatographic models.
For the generation of QSRR models at pH 4.4, 21 compounds were used in the training set, while the remaining 7 compounds (clopamide, clozapine, ephedrine, harman, mianserin, naphazoline, tizanidine) were used in the test set.
Using the stepwise MLR method, 4 descriptors with the greatest influence on the dependent variable log kw4.4 were selected,Equation (11):
while statistically, the most significant PLS-QSRR model was formed with 6 descriptors,Equation (12):
Better regression and validation parameters (Tables 1 and2) were obtained for the PLS-QSRR(log kw4.4) model (Q2 = 0.865, RMSEP = 0.325, F= 57.524, p = 1.52E-08) compared to the MLR-QSRR(log kw4.4) model (Q2 = 0.764, RMSEP = 0.496, F = 23.502, p = 1.57E-06), so that the PLS-QSRR (log kw4.4) model was selected as optimal for the prediction of retention behaviour at pH 4.4.
The QSRR models at pH 7.4 were created with a training set of 21 compounds, while the remaining 7 compounds (clonidine, ephedrine, maprotiline, olanzapine, tizanidine, triamterene, xylometazoline) were used in the test set. The best MLR-QSRR model was obtained with 3 descriptors,Equation (13):
The PLS analysis resulted in the model with five important descriptors,Equation (14):
Although the MLR-QSRR (log kw7.4) model has statistically more significant values for F, p and R2 than the PLS-QSRR (log kw7.4) model, better validation parameters (Q2 = 0.883; R2test=0.854; RMSEP = 0.378) for the PLS model compared to the MLR model ((Q2 = 0.846; R2test = 0.829; RMSEP = 0.446) (Tables 1 and2) were allocateed this regression model as more reliable for predicting the retention behaviour of related compounds at pH 7.4.
At pH 9.1, the test set consisted of 7 compounds (clozapine, efaroxan, ephedrine, maprotiline, tizanidine, triamterene, xylometazoline) and a training set of 21 compounds.
The most significant MLR model describing the retention behaviour at pH 9.1 consists of 3 molecular descriptors,Equation (15):
The statistically best PLS model was formed with 5 descriptors,Equation (16):
The analysis of the calculated statistical parameters (Tables 1 and2) obtained for the created models at pH 9.1, suggests that the MLR-QSRR (log kw9.1) model can be selected as an optimal and reliable model for predicting the retention behaviour of α-adrenergic and imidazoline receptor ligands in an alkaline environment.
In addition, it should be emphasised that all created MLR-QSRR and PLS-QSRR models meet the external validation criteria proposed by Golbraikh, Tropsha and Roy (Table 2) [22-24].
Plots of observed vs. predicted values for compounds in the training and test sets in selected QSRR models are depicted inFigure 1. The values of molecular descriptors in selected QSRR models are listed inTable S3.
Quantitative structure-mobility relationship modelling
Determination of effective electrophoretic mobility
The effective electrophoretic mobility of the investigated compounds was determined in the presence of acetone as a neutral marker, and the results obtained under the examined conditions (pH 4.4, 7.4 and 9.1) are presented inTable S4 (Supplementary Material). At pH 4.4, 27 compounds that move faster than the neutral marker are in the form of cations and reach their maximum effective mobility, while the effective mobility of clopamide and indapamide, which are in neutral form under the studied conditions, is zero. The range of the obtained μeff4.4 values is between 0 (indapamide, clopamide) and 28.19×10-5 cm2 / V s (harmane). At a physiological pH of 7.4, compounds with pKa values around 7.4 have reached their intermediate mobility, indapamide and clopamide co-migrate with the EOF marker (μeff = 0), and compounds with a pKa value above 10 (Table S1, Supplementary material), a high percentage of which are in the form of cations under the given conditions, have maintained an effective mobility similar to that in an acidic environment. The resulting range of effective mobility is between 0 and 24.18×10-5 cm2 / V s (naphazoline). In an alkaline environment at pH 9.1, negatively charged ions of clopamide and indapamide migrate to the cathode more slowly than EOF markers and exhibit negative values of effective mobility. Uncharged compounds migrate together with the EOF marker and their mobility is zero, while the derivatives of 2-methylimidazoline with pKa values above 10 maintain a high degree of mobility (μeff = 25.05×10-5 cm2 / V s, tetrahydrozoline) (Table S1 andTable S4).
Quantitative structure-mobility relationship models
As a result of the QSMR studies applied to the μeff data obtained at three pH values (μeff4.4, μeff7.4 and μeff9.1), six models were created: MLR-QSMR (μeff4.4), PLS-QSMR (μeff4.4), MLR-QSMR (μeff7.4), PLS-QSMR (μeff7.4), MLR-QSMR (μeff9.1), and PLS-QSMR (μeff9.1).
The QSMR models, which describe the mobility of the investigated compounds at a pH of 4.4, were created with 20 compounds in the training set. Three compounds, one with the highest electrophoretic mobility (olanzapine μeff4.4 = 33.632 ± 0.036 10-5 cm2 / V s) and two with the lowest electrophoretic mobility (clopamide μeff4.4 = 0.000 ± 0.000 and indapamide μeff4.4 = 0.000 ± 0.000) clearly showed a significant deviation from the regression line on the t1 vs. u1 scater plot and were therefore identified as outliers. The remaining 6 compounds (amiloride, carvedilol, ephedrine, mianserin, tizanidine and xylometazoline) were used for external validation.
The most significant MLR-QSMR model at pH 4.4 was formed with 4 descriptors,Equation (17):
The PLS analysis also selected the model with the four descriptors that showed the greatest influence on the dependent variable μeff 4.4,Equation (18)
According to the better statistical parameters of the PLS-QSMR model (Q2 = 0.931, F = 115.508, p = 1.28×10-10, R2test = 0.897 and RMSEP = 0.984) compared to the MLR-QSMR model (Q2 = 0.880, F = 64.338, p = 2.92×10-9, R2test = 0.876 and RMSEP = 1.063) (Tables 3 and4), the PLS-QSMR model can be selected as the optimal and reliable model for predicting the electrophoretic mobility of related compounds at pH 4.4.
QSMR models at pH 7.4 were formed with 20 compounds in the training set, while 5 compounds (ephedrine, olanzapine, oxymetazoline, tamsulosin and tizanidine) were retained in the test set. Considering the fact that the investigated data sets are very different in their acid-base properties, the resulting differences in effective electrophoretic mobility are significant (Table S1, Supplementary material).
Clopamide and indapamide, which carry a weakly acidic sulfonamide group, are still predominantly unionized at pH 7.4 and move together with the electroosmotic flow, resulting in an electrophoretic mobility of zero.
A very low effective electrophoretic mobility is also obtained for triamterene (μeff7.4 = 1.140±0.053××10-5 cm2 / V s) as a weak base and guanfacine (μeff7.4 = 5.083 ± 0.025×10-5 cm2 / V s). The Y-values of these compounds in relation to the calculated X-values showed a significant deviation in the scatter plot t1 vs. u1 and were identified as outliers. The most significant MLR-QSMR model is composed of 4 descriptors,Equation (19):
Using PLS analysis, a model was constructed with 5 molecular descriptors that significantly affect the electrophoretic mobility of the studied compounds at pH 7.4,Equation (19):
For the MLR-QSMR model, lower values for RMSEE and RMSEP, slightly better parameters Q2 and F value as well as significantly higher values for R2trening and R2test were obtained compared to the PLS-QSMR model (Tables 3 and4). In addition, the PLS-QSMR model did not meet external validation criteria established by Golbraikh and Tropsha [22] [(R2 - ) / R2 = 0.178 > 0.1; R'02= -0.818; (R2- R'02 ) / R2 = 2.192 > 0.1] (Table 4); and Roy [23,24] [ =0.446 and = -0.155 are not closed; = 0.145 < 0.5; = 0.601 > 0.2] (Table 4) indicating that the MLR-QSMR regression equation can be reliably applied to predict the electrophoretic mobility of related compounds at pH 7.4.
At a pH of 9.1, nine compounds migrate together with EOF and have an effective mobility of 0. Five compounds with an effective mobility of zero remain in the training set, while four compounds move to the test set taking into account that all chemical clusters have a representative in the test set. Since all created models were unable to predict the mobility of harman and mianserin (μeff = 0) and taking into account that these two compounds had the greatest deviation on the t1 vs. u1 plot, among others from the training set with zero mobility, these two compounds were excluded from further analysis.
QSMR models at pH 9.1 were created with 21 compounds in the training set, while 6 compounds (doxazosin, ephedrine, idazoxan, olanzapine, tizanidine, xylometazoline) were used for external validation.
By MLR analysis, four significant molecular descriptors were selected,Equation (21):
While the statistically best PLS model was formed with 5 important molecular parameters,Equation (22):
Appropriate validation criteria were met for both MLR-QSMR/PLS-QSMR models created (Tables 3 and4), with the exception of the value (R2-R'02)/R2 0.141 > 0.1 of the PLS-QSMR (μeff9.1) model (Table 4). Due to the higher values of Q2, F, R2trening, R2test (Q2 = 0.891, F = 58.736, R2traning = 0.936 and R2test =0.879) for the MLR model than for the PLS model (Q2 = 0.882, F = 29.464, R2trening = 0.916 and R2test = 0.873) and the lower RMSEE and RMSEP (RMSEE = 2.640 and RMSEP = 3.794) for the MLR model than for the PLS model (RMSEE=3.033 and RMSEP = 4.290), the MLR-QSMR model: MLR-QSMR (μeff9.1) = ƒ(SM1_Dz (Z); SM15_EA (dm); HOMO [B3LYP/3--21G (d, p)]; ISH) was selected as the optimal model for predicting the effective mobility of related compounds at a given pH of 9.1.
Plots of observed vs. predicted values for compounds in the training and test sets of the selected QSMR models are depicted inFigure 2.
It is interesting to note that inFigure 2B, which shows the correlation between the calculated and predicted values of the QSMR model at pH 7.4, two clusters of data can be observed. The upper cluster includes compounds whose experimentally determined effective mobility is greater than 19×10-5 cm2 / V s and experimental pKa value is greater than 8, and the lower cluster includes compounds with experimental effective mobility of less than 15×10-5 cm2 / V s and pKa value less than 8 (Table S1 and S4, Supplementary material). The selected MLR-QSMR model could be used for a rough prediction of the pKa values of new imidazoline analogs.
The values of the molecular descriptors in the selected QSMR models are listed inTable S5 (Supplementary Material).
Interpretation and comparison of selected quantitative structure-retention relationship and quantitative structure-mobility relationship models
Based on the presented statistical results (Tables 1-4) obtained after internal and rigorous external validation, optimal QSRR and QSMR models were selected for each investigated pH value. The classes of descriptors in the selected models and their descriptions are presented inTables 5 and6.
The descriptors selected for the optimal PLS-QSRR (log kw4.4) model belong to three different descriptor classes (Table 5), such as 2D matrix-based descriptors (H_Dz(p), VE2_Dz(p)), Burden eigenvalues (SpMax5_Bh(e), SpMax5_Bh(i)) and Edge adjacency indices (SM02_EA(ri), SM04_EA(ri)).
The coefficient plot (Figure 3) shows that all descriptors have a positive influence, with the exception of VE2_Dz(p), which has a negative influence on the log kw4.4 value.
2D matrix-based descriptors are topological indices computed by applying a set of basic algebraic operators to different graph-theoretical matrices representing an H-depleted molecular graph [29]. Depending on which matrix they are derived from (detour matrix Dt, Barysz matrices (Dz(w)), etc.), they can be divided into different sub-blocks. H_Dz(p) (Harary-like index from Barysz matrix weighted by polarizability) and VE2_Dz(p) (average coefficient of the last eigenvector from the Barysz matrix weighted by polarizability) are derived from the Barysz matrices. Barysz matrices (Dz(w)) are weighted distance matrices that take into account the presence of heteroatoms and multiple bonds in the molecule. Polarizability is an important atomic property included in both the VE2_Dz(p) and H_Dz(p) descriptors, which differ in their influence on log kw4.4. Higher values of VE2_Dz(p) are associated with less lipophilic compounds in the studied group that have lower log kw4.4 values (ephedrine, rilmenidine), while lower values are seen for the structures of carvedilol, clozapine and maprotiline with the highest log kw4.4 values. The opposite effect can be observed for H_Dz(p). SpMax5_Bh(e) and SpMax5_Bh(i) are the largest eigenvalue n. 5 of the Burden matrix weighted by Sanderson electronegativity and ionization potential, respectively. Burden eigenvalues descriptors were originally proposed for searching for chemical similarity/diversity in large databases. They are calculated from the Burden matrix as an ordered sequence of the largest positive and smallest negative eigenvalues that reflect relevant aspects of molecular structure and can be useful for similarity searches [29,39].
SM02_EA(ri) and SM04_EA(ri) (spectral moments of order 2 and 4 from edge adjacency mat. weighted by resonance integral) belong to the Edge adjacency indices descriptors. These indices are based on the edge adjacency matrix of a graph. They encode information about the connectivity between graph edges, taking into account bonds between non-hydrogen atom pairs. Dragon calculates the Edge adjacency matrices using the following bond properties: edge degree (ed), dipole moment (dm), conventional bond order (bo) and parameters related to the resonance integral (ri). The spectral moment represents the linear combination of different structural fragments in the graph [29,40,41].
In contrast to chromatographic behaviour, electrophoretic mobility in an acidic medium is primarily determined by 2D autocorrelation descriptors (ATS6s, ATS7s, ATSC6m) (Table 6), which describe how the observed property (in this case, the intrinsic state in ATS6s and ATS7s and the mass in ATSC6m) is distributed along the topological structure.
This class of descriptors represents the homogeneity of the molecular structure and can be weighted by the intrinsic state (s), the Sanderson electronegativity (e), the atomic mass (m) and the van der Waals volume (v) [42]. H0_Dz(e) (Hosoya-like index (log function) from the Barysz matrix, weighted by the Sanderson electronegativity) is 2D-matrix-based descriptor and represents a measure of molecular size and branching [43]. From the coefficient plot (Figure 4), It can be concluded that all selected descriptors have a negative influence on the μeff4.4 values.
At a physiological pH of 7.4, retention behaviour is a function of 2D matrix-based descriptors and molecular properties (Table 5).
Selected 2D matrix-based descriptors are derived from the detour matrix (Ho_Dt, SM3_Dt, VE2_Dt) and the distance/detour matrix (H_D/Dt), in contrast to the 2D matrix-based descriptors at pH 4.4. Ho_Dt is a Hosoya-like index (log function) from the detour matrix and represents a molecular descriptor that depends on the size of the molecule as well as its voluminosity [43-45]. Molecules with a simple structure, such as ephedrine and phenilephrine have the lowest value of Ho_Dt, while complex structures, such as those found in doxazosin and tetracyclic (mianserin) and tricyclic compounds (carvedilol, clozapine) have the highest value of this descriptor. H_D/Dt (Harary-like index from distance/detour matrix) is also a parameter sensitive to steric effects and whose value increases with increasing molecule size and branching [46]. Both descriptors are positively correlated with the log kw7.4 data (Figure 5.), which means that more lipophilic compounds are characterised by a higher value of this descriptor.
The logarithm of distribution coefficient at pH 7.4, log D7.4 [28], is a measure of the lipophilicity of compounds at pH 7.4 and indicates the importance of considering the presence of ionic species when evaluating the lipophilicity of compounds at pH 7.4. The selected model shows that calculations of all proposed descriptors allow a more reliable prediction of compound lipophilicity and that the design of compounds with the desired properties is influenced by the presence of groups responsible for steric effects and the size of the molecule.
In capillary electrophoresis, the electrophoretic mobility of investigated compounds at pH 7.4 can be discussed on the basis of selected 2D autocorrelations (GGI6), GETAWEY and 3D-MoRSE descriptors (Table 6). In contrast to the PLS-QSMR (μeff4.4) model, in which the class of 2D autocorrelation descriptors dominates, the 3D descriptors appear to be more important in the selected MLR-QSMR (μeff7.4) model.
GETAWAY descriptors combine the 3D-molecular geometry provided by the molecular influence matrix and atomic relatedness through molecular topology with chemical information by using different atomic weights (atomic mass, polarizability, van der Waals volume and electronegativity) together with unit weights [47]. GETAWAYs (HATS2i, R7u+) are geometric descriptors that are sensitive to molecular branching and cyclicity, encode information about the effective position of substituents and fragments in molecular space, and may be suitable to describe differences in congeneric molecular series [29].Equation 19 shows that both descriptors (HATS2i, R7u+) have a negative influence on the μeff7.4 values.
3D-MoRSE descriptors (Mor28s) represent the 3D structure of molecules based on electron diffraction [48]. Mor28s is a signal 28 / weighted by I-state and, in contrast to the GETAWAY descriptors, has a positive influence on the μeff7.4 parameter (Equation 19).
All molecular descriptors selected with the MLR-QSRR (log kw9.1) model (Equation (15)) have a positive influence on log kw9.1 and belong to constitutional indices, Burden eigenvalues and molecular properties classes of descriptors (Table 5).
Constitutional descriptors such as the selected nBM (number of multiple bonds), which counts a molecule's double, triple and aromatic bonds, are the simplest and most commonly used descriptors. They represent the chemical composition of a compound without information about the molecular geometry or the connectivity of the atoms [29]. SpMin6_Bh(p) is the smallest eigenvalue n.6 of the Burden matrix weighted by polarizability and belongs to the class of Burden eigenvalue descriptors. The lipophilicity of compounds at pH 9.1 (log D9.1) [28] appears again as an important molecular property considering all forms of the molecule (ionised and unionised), and its high positive correlation with log kw9.1 suggests that more lipophilic compounds have a higher affinity for the apolar stationary phase and greater retention in HPLC.
In capillary electrophoresis at pH 9.1, the descriptors selected using the MLR method (Equation (21)) belong to the 2D matrix-based descriptors, GETAWAY, Edge adjacency indices and quantum chemical descriptors (Table 6).
SM15_EA(dm) is the spectral moment of order 15 from the edge adjacency matrix, weighted by the dipole moment, and is the only descriptor in a model that shows a positive correlation with μeff9.1. In contrast, toSM15_EA(dm), SM1_Dz(Z) (spectral moment of order 1 from the Barysz matrix weighted by atomic number) belongs to the 2D matrix-based descriptors. ISH (standardised information content on leverage equality) is a class of GATAWEY descriptors and mainly encodes information about molecular symmetry. If all atoms have different leverage values, the molecule has no symmetry and ISH=1, otherwise a theoretically perfectly symmetric molecule has ISH=0. In the case of the tested compounds, all ISH values are above 0.82 and show a negative correlation with electrophoretic mobility, suggesting that compounds with a higher degree of symmetry have lower mobility under the given experimental conditions of background electrolyte pH 9.1. This descriptor also provides information about molecular entropy, so it can be useful in modelling physicochemical properties related to entropy and symmetry [47]. The highest occupied molecular orbital (HOMO) energy calculated using the B3LYP/3-21G(d,p) method [27] is negatively correlated with μeff9.1 (Equation (21)), implying that compounds with a lower HOMO energy have greater mobility because they interact less with positively charged buffer cations against the wall of the capillary.
Conclusions
Evaluation of the chromatographic and electrophoretic behaviour of a series of 29 imidazoline and alpha adrenergic receptors ligands at pH 4.4, 7.4, and 9.1 using QSRR and QSMR models allows several conclusions to be drawn. When comparing the statistical PLS and MLR methods used for model building, it is not possible to give preference to only one method, as the optimal models were formed by both. All selected QSRR models: PLS-QSRR (log kw4.4), PLS-QSRR (log kw7.4), MLR-QSRR (log kw9.1) and QSMR models: PLS-QSMR (μeff4.4), MLR-QSMR (μeff7.4), MLR-QSMR (μeff9.1) displayed high accuracy in retention/migration prediction for internal and external test set compounds. Based on the established QSRR and QSMR models, it can be seen that different classes of descriptors appear in the models created at different pH values and for different analytical systems, such as HPLC and CE. Chemical composition of compounds, lipophilicity, and voluminosity have the strongest influence on the retention mechanism, while molecular size, mass, complexity, charge and electronic properties are valuable structural determinants for the electrophoretic mobility of the examined ligands. The created QSRR and QSMR models are a very useful predictive tools for retention and electrophoretic behaviour of new imidazoline/alpha adrenergic receptors ligands at three pH values, which could be used not only during the RP-HPLC/CE method development but also in the initial phase of design of novel imidazoline and alpha adrenergic receptors ligands with optimized physicochemical properties, such as lipophilicity (log kw values) and charge-to-mass characteristics (μeff values).