Permeability Data of Organosulfur Garlic Compounds Estimated by Immobilized Artificial Membrane Chromatography: Correlation Across Several Biological Barriers

Among healthy vegetables, those of the genus Allium stand out. Antioxidant and anti-inflammatory properties have been associated with these vegetables, attributed mainly to organosulfur compounds (OSCs). In turn, they are linked to a protective effect counteracting cardiovascular disease development. Now, to really ensure the bioactive efficacy of the said compounds once consumed, it is necessary to previously evaluate the ADME (absorption, distribution, metabolism, and excretion) profile. Alternatively, in vitro and in silico methods attempt to avoid or reduce experimental animals’ use and provide preliminary information on drugs’ ability to overcome the various biological barriers inherent in the ADME process. In this sense, in silico methods serve to provide primary information on drugs’ bioavailability mechanisms. High-performance liquid chromatography (HPLC) using a stationary phase composed of phospholipids, the so-called immobilized artificial membrane (IAM), has been widely recognized as a valuable alternative method to extract and quantify information about the structure and physicochemical properties of organic compounds which are extensively used in studies of quantitative structure–activity relationships (QSARs). In the present study, the chromatographic capacity factors (log k’ (IAM)) for 28 OSCs were determined by IAM-HPLC. In order to evaluate the ability of the IAM phase in assessing lipophilicity of the compounds under study, several quantitative structure–retention relationships (QSRRs) were derived from exploring fundamental intermolecular interactions that govern the retention of compounds under study on IAM phases. As expected, the hydrophobic factors are of prime importance for the IAM retention of these compounds. However, the molecular flexibility and specific polar interactions expressed by several electronic descriptors (relative negative charge, RNCG, and Mulliken electronegativity) are also involved. We also evaluated the IAM phase ability to assess several ADME parameters for the OSCs under study obtained using the SwissADME web tool integrated into the SwissDrugDesign workspace and the PreADMET web tool. The human gastrointestinal absorption (HIA), blood–brain barrier (BBB) permeation, and skin permeability were investigated through QSAR modeling, using several chemometric approaches. The ADME properties under study are strongly dependent on hydrophobic factors as expressed by log k’(IAM), which provide evidence for the great potential of the IAM phases in the development of QSAR models.

Among healthy vegetables, those of the genus Allium stand out. Antioxidant and antiinflammatory properties have been associated with these vegetables, attributed mainly to organosulfur compounds (OSCs). In turn, they are linked to a protective effect counteracting cardiovascular disease development. Now, to really ensure the bioactive efficacy of the said compounds once consumed, it is necessary to previously evaluate the ADME (absorption, distribution, metabolism, and excretion) profile. Alternatively, in vitro and in silico methods attempt to avoid or reduce experimental animals' use and provide preliminary information on drugs' ability to overcome the various biological barriers inherent in the ADME process. In this sense, in silico methods serve to provide primary information on drugs' bioavailability mechanisms. High-performance liquid chromatography (HPLC) using a stationary phase composed of phospholipids, the so-called immobilized artificial membrane (IAM), has been widely recognized as a valuable alternative method to extract and quantify information about the structure and physicochemical properties of organic compounds which are extensively used in studies of quantitative structure-activity relationships (QSARs). In the present study, the chromatographic capacity factors (log k' (IAM)) for 28 OSCs were determined by IAM-HPLC. In order to evaluate the ability of the IAM phase in assessing lipophilicity of the compounds under study, several quantitative structure-retention relationships (QSRRs) were derived from exploring fundamental intermolecular interactions that govern the retention of compounds under study on IAM phases. As expected, the hydrophobic factors are of prime importance for the IAM retention of these compounds. However, the molecular flexibility and specific polar interactions expressed by several electronic descriptors (relative negative charge, RNCG, and Mulliken electronegativity) are also involved. We also evaluated the IAM phase ability to assess several ADME parameters for the OSCs under study obtained using the SwissADME web tool integrated into the SwissDrugDesign workspace and the PreADMET web tool. The human gastrointestinal absorption (HIA), blood-brain barrier

INTRODUCTION
The current guidelines recommend diets rich in fruits and vegetables as an important factor in reducing the risk of developing chronic diseases (Yu et al., 2018). Among healthy vegetables, those of the genus Allium stand out. Particularly, abundant literature supports garlic and its by-products' healthpromoting effects associated with bioactive compounds present in their matrix (Corzo-Martinez et al., 2007;Santhosha et al., 2013;Ramirez et al., 2021). Amid these bioactivities, cardiovascular protection is one of the most well-known beneficial effects related to garlic or garlic-based products' intake (Banerjee and Maulik, 2002;Santhosha et al., 2013;Schwingshackl et al., 2016). In this sense, evidence has shown that organosulfur compounds (OSCs) present in different garlic matrices interfere against inflammation, oxidative stress, obesogenic effects, and mitochondrial dysfunction (Quesada et al., 2020) by multiple actions, such as modifying signaling pathways that trigger chronic diseases. Now, if we focus on the bioactive efficacy, we must first ensure that the phytochemical compounds manage to reach the target sites within the body where they will exert their biological effect, so it is necessary to previously evaluate in vivo the ADME (absorption, distribution, metabolism, and excretion) profile of the said compounds (Rein et al., 2013). Considering the previously mentioned biological steps, cellular permeability of bioactive compounds is one of the key physicochemical parameters to be considered when selecting potential therapeutic candidates (Gozalbes et al., 2011). To estimate cellular permeability there have been several absorption models that are proposed, including cell-based models, using the human colorectal adenocarcinoma (Caco-2) cell line and the Madin-Darby Canine Kidney (MDCK) cell line and tissuebased models, such as in situ rat intestinal perfusion or the Ussing chamber system employing rat or human intestinal tissues (Berben et al., 2018). Although these models are useful in terms of permeability prediction, they have several drawbacks, including high-cost and complex experimental procedures as well as ethical issues associated with animal-based assays (Grumetto et al., 2015a).
A recent work from our research group has determined that intestinal absorption proved to be a limiting step for OSC availability (Torres-Palazzolo et al., 2018). However, it is necessary to delve deeper into the structural basis of the permeability stage mechanisms based on a molecular level (Butina et al., 2002). Alternatively, computational approaches (in silico) serve to provide primary information on drugs' bioavailability mechanisms (Yen et al., 2005). In this context, alternative methods, such as chromatographic ones, employing immobilized artificial membrane (IAM) chromatographic columns, have been widely recognized as a valuable tool to model processes in the biophase since the components of both chromatographic and biological systems are comparable. This is supported by the fact that some biological processes, such as passive cellular diffusion, have much in common with the processes that take place in some chromatographic separations (Escuder-Gilabert et al., 2003). To delve deeper into the physicochemical processes that are involved in cellular permeability, the information resulting from chromatographic data using membrane-like systems, such as stationary phases containing phosphatidylcholine or similar lipid-like ligands, is used to develop statistical models called quantitative structure-activity relationship (QSAR) models, using several chemometric tools.
Based on the above-mentioned models, the chromatographic capacity factors (log k' (IAM)) for 28 OSCs were determined by IAM-HPLC and various quantitative structure-retention relationships (QSRRs) were derived to explore the fundamental intermolecular interactions that govern the retention of compounds under study on the IAM phases. Thus, the partial least squares (PLS) approach was used to model the IAM chromatographic data determined for these compounds. The molecular characterization was carried out by calculating numerous nonempirical descriptors, which were subsequently used to construct QSRR models.
On the other hand, we also attempted to evaluate the IAM phase ability to assess several ADME parameters for the OSCs under study using the SwissADME web tool integrated into the SwissDrugDesign workspace and the PreADMET web tool. The human gastrointestinal absorption (HIA), blood-brain barrier (BBB) permeation, and skin permeability were investigated through QSAR modeling, using several chemometric approaches.
E-Z Ajoene and vinyldithiins compounds were synthesized following a previously reported procedure by our research group (Ramirez et al., 2017). Then, the synthesized OSCs were concentrated under reduced pressure and characterized by UV-spectroscopy and GC-MS analysis.
Potassium phosphate monobasic and dibasic (anhydrous) were of PA grade and were purchased from Biopack. Acetonitrile (ACN) was of chromatographic grade and was purchased from Merck (United States). Ultrapure water (18 MΩcm) was obtained from a Milli-Q water purification system (Millipore, France).

Chromatographic Conditions and Instrumentation
Chromatographic analysis was performed using a Shimadzu LC 20A chromatograph coupled to a diode array detector, operated at a range of 210-260 nm, depending on the analyzed compound. A biomimetic commercially available column was used in all experiments: IAM. PC.MG12 µm (150 × 4.6 mm) was purchased from Regis Chemical Company (Morton Grove, IL).
OSC separation conditions under IAM-HPLC were adapted from Luco et al. (2003) and Barbato et al. (1998), with modifications. The chromatographic mode consisted in an isocratic elution used as the mobile phase ACN/0.035 M phosphate buffer −pH 6.8− (30:70 (v/v)) at 0.5 ml min −1 . The column dead time (T 0 ) was estimated by ascorbic acid retention time which was measured at 210 nm. The obtained retention data were used to calculate the chromatographic parameter corresponding to the capacity factor k' as follows: k' (T r -T 0 )/T 0 . All capacity factors given represent the mean of 2-4 determinations of each sample solution.
For this study, we considered 28 phytochemical compounds containing sulfur-based groups. Supplementary Table S1 shows their name and chemical structure. Most of these compounds are garlic-derived substances with bioactive functions. However, we also included seven phytochemicals-sulforaphane, sulforaphene, I3C, allyl ITC, phenethyl ITC, phenyl ITC, and benzyl ITC-that correspond to glucosinolates degradation products, such as isothiocyanates and indole-3-carbinol. We chose them because of their bioactive functionality and their similar metabolism within the plants and human body. Each analytical standard stock solution was prepared in the mobile phase at concentrations ranging from 100 to 800 µg ml −1 , depending on the solubility of the compound and their DAD detector response. Peak identification in samples was carried out by comparing retention times with reference standards.

Molecular Descriptors and ADME Properties
To characterize each OSC, a set of molecular descriptors were calculated using different software packages, including Dragon vs. 3.0 and the Online Chemical Modeling Environment (OCHEM) web platform. As for indicators of the molecular size and structural descriptors, the following were considered: molar volume (Vm), molecular weight (MW), rotatable bonds, H-donors, H-acceptors, and TPSA (topological polar surface area). To explain lipophilicity effects, several calculated partition coefficients were obtained: log D, iLOGP, XLOGP3, WLOGP, MLOGP, Silicos-IT Log P, and Consensus Log P. These parameters were obtained from different online web tools, such as SwissADME (developed and maintained by the Molecular Modeling Group of the Swiss Institute of Bioinformatics) and Chemicalize © 1998-2021.
Another group of structural descriptors included quantum chemical indexes obtained by the HyperChem package (release 7.5 for Windows). Three-dimensional molecular structures were built using two different software: the MM + molecular mechanics potential energy function employing the HyperChem package and on the other side the CORINA Classic online service. In a follow-up procedure, complete optimization of the geometrical parameters was carried out by using the PM3 method, implemented in the standard version of MOPAC 6.0. The following indices obtained from molecular orbital calculations were considered: total energy (E total ), heat of formation (ΔHf), energy of highest occupied molecular orbital (HOMO), energy of lowest unoccupied molecular orbital (LUMO), dipole moment (µ), absolute total charge (Qtotal), the most positive and the most negative absolute charges (q p,max , q n,max ), and the positive and negative relative charge (RNCG and RPCG). Furthermore, other descriptors obtained from Dragon software were considered such as Mulliken electronegativity, LDip, and several WHIM descriptors. Multivariate analyses including principal components analysis (PCA) and a stepwise multiple regression procedure, based on the algorithms forward selection and backward elimination, were used for the inclusion or rejection of descriptors in the screened models.
Finally, to assess permeability parameters for the OSCs under study, several calculated properties were obtained from the SwissADME web tool and the PreADMET web-based application ( © 2005-2017 BMDRC. | Designed by Y.-M. Kang). The human gastrointestinal absorption (HIA), blood-brain barrier (BBB) permeation, and skin permeability (Log K P ) were investigated.
To ensure the reliability of the in silico data about HIA, BBB passage, and Kp skin calculated from the SwissADME and PreADMET web platforms, the predictor space of each platform was analyzed. By doing so, the applicability domain (AD) of each QSAR model which serves to calculate the prediction values offered in the web tools can be tested. In silico prediction systems feature the capability to evaluate if a query molecule is a part of the developed model's AD or not and, therefore, gave an improved sense of reliability at the individual prediction level (Kar et al., 2018). For this, a series of physicochemical descriptors were calculated, both for each of the compounds that make up the predictor space associated with each permeation response (241 compounds from the data set to build the HIA permeation model (Zhao et al., 2001), 260 compounds from the data set were used to calculate BBB passage (Daina and Zoete, 2016), and 96 compounds from the data set were used to calculate Kp skin , as well as for our 28 OSCs. From these data, we then performed a PCA, and it was evaluated whether our OSCs belong to the predictive space generated by the compounds used in such platforms.

Statistical Methods
A multivariate analysis was used primarily. Multiple regression analysis (MLR) and partial least squares projections (PLS) were the selected methods to build different statistical models between IAM chromatographic parameters and the structural descriptors.
The determination of the significant number of PLS components was made by cross-validation. MLR analysis was performed using MINITAB ® Release 17, whereas PLS analysis was carried out using SIMCA-P 7.01 software obtained from Umetri AB, Umea, Sweden. To avoid overestimations or difficulties in the interpretation of the resulting models, pairs of variables with r ≥ 0.75 were classified as intercorrelating ones, and only one of these was included in the screened model. The predictive ability of the models were evaluated by the squared correlation coefficient (R 2 ) and the Fisher F statistic (F). In PLS models, we also evaluated the cross-validation coefficient (Q 2 ), which is based on the prediction error sum of squares (PRESS). Considering the QSAR models to assess permeability parameters, we used TableCurve 2D ® to build a nonlinear curve-fitting model for HIA vs. log k'IAM, MINITAB ® Release 17 to carry out a discriminant analysis to assess BBB permeation, while skin permeation was evaluated using a PLS analysis (SIMCA-P 7.01 software). Principal component analysis to check OSCs belonging to the predictor space of each permeation response, according to the different web platforms, was carried out using MINITAB ® v.17.

OSC Lipophilicity Determination Using IAM and Physicochemical Parameters
First, we evaluated the IAM phase ability in assessing OSC lipophilicity. For this, we studied to what extent the log k'IAM values were related to the calculated values of the partition coefficient log P, obtained from the SwissADME web tool. Table 1 shows the chromatographic data and the corresponding physicochemical parameters of the compounds under study. In this case, we built a MLR model based on iLOGP (implicit Log P). This parameter is a calculated value corresponding to the n-octanol/water partition coefficient descriptor based on the GB/SA approach. It has been proven to be an efficient physicochemical parameter showing a strong linear correlation between the computed solvation free energy in implicit solvents and the experimental log P o/w on more than 17,500 molecules (Daina et al., 2014). The regression equation using the whole set of compounds was as follows: In Eq. 1 and the following equations, n is the number of compounds, s is the standard deviation, R 2 is the squared correlation coefficient, R 2 (cv) is the squared cross-validation coefficient, and F is the Fisher F statistic. Values in parentheses correspond to the standard deviations and p-values of coefficients, and the term W represents the standardized regression coefficient.
The model statistical quality is suitable as seen in the quality of fit and predictive ability, expressed by R 2 and R2 cv coefficients, respectively. However, some aspects should be considered to further improve the statistics of Equation 1. The relationship between the log k'(IAM) and iLOGP values is depicted in Figure 1. There, it is possible to distinguish different subgroups of compounds, marked by different colors, corresponding to nonionizable OSCs (mostly sulfide-based compounds), ionizable OSCs (ACSOs and cysteine-containing compounds), sulfoxide-based compounds (allicin, sulforaphane/ ene, and E-Z-ajoene), and two special cases (I3C and AMS). This suggests a divergence between the different compounds that may be due to the sulfoxide groups' different interactions, both with phosphatidyl choline in the IAM phase as well as with octanol, when we consider the octanol-water partition. These interactions will be explained in more detail in the next section, and they will be discussed accordingly.

QSRR Models to Assess OSC Interactions With the IAM Phase
Several QSRR equations were built to evaluate the fundamental intermolecular interactions that govern the retention of the different OSCs on the IAM phases. As shown in the previous section, lipophilicity is an essential factor to be considered in the quantitative structure-retention relationship (QSRR) models. Nevertheless, since lipophilicity can be described by the contribution of two components, namely molecular size and molecular polarity, models that incorporated both molecular properties were examined. The best regression equation that included the iLOGP parameter was as follows: which is fundamental to reach a correct physicochemical interpretation. An interesting point to highlight in this equation is an improvement in the statistical parameters when compared to Eq. 1, which suggests not only a better predictive ability but also, as hypothesized, that specific interactions regarding OSC behavior in the IAM phase precisely emerge from the other terms that are added to the model, besides iLOGP.
As an alternative way to evaluate the reliability and robustness of the QSRR model, we carried out a different statistical approach using the same descriptors included in the MLR model. Accordingly, we performed the analysis of molecular descriptors using partial least squares regression. The statistical significance of the QSRR-PLS derivative models was evaluated by means of the variance of the matrices X and Y (R2X and R2Y), the square correlation coefficient (R 2 ), standard deviation (RMSS), and the statistical F. The predictive ability was evaluated by the cross-validation coefficient (Q) which is based on the prediction error sum of squares (PRESS). The PRESS statistic is computed as the squared differences between observed and predicted values when the observations are kept out of the derived model (Luco et al., 2003).
The model was built using the total set of OSCs, resulting in two statistically significant components (CPLS), whose parameters are detailed in Table 2. As shown in Figure 2, the agreement between the measured and calculated data from the derived PLS model is very satisfactory. There was also a linear tendency of the residuals when observing the normal probability plot of the said values (Data not shown), which indicates that the error terms were normally distributed.
The only compound that showed an outlier behavior was the AMS, which, while observing the standardized residuals, showed more than two standard deviations with respect to the centralized tendency of the rest of the compounds. Thus, considering this outlier, a second PLS model was built without AMS (compound 14), and as expected the PLS model's-with two components-statistical quality improved (R 2 0.953; Q 0.763; and F 240.76).

IAM Phase Ability to Assess OSC Human Intestinal Absorption
To evaluate the effectiveness of the IAM phase in the prediction of human intestinal absorption percentage (HIA%), we built several QSAR models based on MLR analysis. Since there is a lack of information regarding organosulfur compounds' HIA in-vivo values, we used calculated values obtained from the PreADMET web application. First, we assessed the reliability of PreADMET data in relation to our OSCs. As shown by Kar et al. (2018), to analyze how much a certain model's prediction for a query compound can be trusted, in our case for the target OSCs, it is necessary to evaluate the applicability domain (AD) of that particular QSAR model. In this sense, it is important to recognize that QSARs are associated with restrictions in terms of the grouping of chemical structures, physicochemical properties, and mechanisms of action for which the models can produce FIGURE 1 | Relationship between log k' (IAM) and iLOGP for the 28 compounds shown in Table 1. Key: Nonionizable organosulfur compounds Ionizable organosulfur compounds Compounds with sulfoxides Unusual observations. We also calculated these parameters corresponding to our OSCs set. We then carried out a PCA, since principal components analysis has showed to be a valid approach to assess the AD space of a QSAR model (Kar et al., 2018). Results are shown in the Supplementary Material section. Shortly, it was possible to confirm OSCs belonging to the original data predictor space that supports PreADMET web server values. The Mahalanobis distance plot (Supplementary Figure S2) serves to probe the previous approach, where it can be seen that OSCs showed acceptable Mahalanobis distance. The Mahalanobis distance  Frontiers in Chemistry | www.frontiersin.org September 2021 | Volume 9 | Article 690707 7 measures the distance from each point in multivariate space to the overall mean or centroid, utilizing the covariance structure of the data. All these results serve to assure that HIA predicted values of the OSCs under study are actually reliable to build the following QSAR models (Minitab v.17; Minitab Inc.; State College, PA 2010).
Focusing on our QSAR models, considering a linear estimation, a poor statistical correlation was found for the target OSCs (S 6.04354; R 2 (adj) 46.8%; PRESS 371.285; F 6.27). Hence, based on the previous results and OSC disposition when plotting log k' (IAM) vs. HIA (Figure 3), a nonlinear estimation was performed to model the data. In this case, employing Even though a reasonable estimation was achieved to model OSC behavior regarding HIA, notably, a deviation in the fit was observed due to the behavior of ionizable OSCs (ACSOs, cysteine containing compounds, and I3C). In this sense, previous reports have treated ionizable compounds separately based on their lipophilicity values, and the occurrence of polar and/or electrostatic intermolecular interaction forces (Grumetto et al., 2015b;Russo et al., 2017). So, for better adjusting estimations, these compounds were separately modeled. In this case, a simple regression analysis and a fitted line plot were carried out for the seven ionizable compounds. As stated by Russo et al. (2017), the polar surface area is a critical factor to model HIA when dealing with ionizable compounds, since electrostatic and polar interaction forces affect the analyte-phospholipid partition. Therefore, the total polar surface area (named TPSA_1 to differentiate from the acronym TPSA-topological polar surface area-was considered in the model as follows. HIA ionizable OSCs (%) 100(4.182) − 0.143(0.023)TPSA_1 (0.000) (0.000) R 2 0.884; R 2 (cv) 0.837; s 3.087; n 7; F 38.20; W TPSA_1 −0.940 (4) Figure 4 shows the relationship between HIA % and total polar surface area (TPSA_1) for the following compounds: alliin, isoalliin, methiin, S-allyl-L-cysteine, S-methyl-L-cysteine, S-propyl-L-cysteine, and indole-3-carbinol.

Blood-Brain Barrier Permeation QSAR Modeling
Given that IAM chromatography has been successfully applied to drug permeability predictions across the blood-brain barrier (Luco and Marchevsky, 2006), a QSAR model based on IAM phase parameters was built to assess the BB barrier permeation of the target OSCs. Acknowledging that there is no BBB-empirical data reported for OSCs, the BBB-QSAR model was built using estimated data from the SwissADME web tool. The reliability of SwissADME BBB permeation data with regard to our OSCs was evaluated following the steps previously discussed in IAM Phase Ability to Assess OSCs Human Intestinal Absorption. Results from a PCA including several physico-chemical parameters   Figure 6. The relationship between actual and calculated log KP skin permeability values using the PLS approach.
Frontiers in Chemistry | www.frontiersin.org September 2021 | Volume 9 | Article 690707 that our OSC data set actually belongs to the BBB-permeation predictor space generated from the 260 molecules used to developed BBB-permeation QSPR model by SwissADME (Daina and Zoete, 2016); therefore, OSC BBB-permeation data from SwissADME could be used to build the corresponding QSAR model. BBB permeation data from the SwissADME web tool shows which compound could permeate the BBB (marked as Yes value, in Table 1) and which one could not (marked as No value, in Supplementary Table S1). Since BBB permeation was a categorical variable, a discriminant analysis using crossvalidation was carried out. The grouping variable was BBB permeation, while the predictors were log k' (IAM) and TPSA. The discriminant ability was assessed by the correct classifications percentage, and the discriminant function quality was evaluated using the Wilks parameter, λ, which was obtained by a multivariate analysis of variance that tests the equality of group means for the variable in the discriminant model. The standardized discriminating function was The DF equation is highly significant as p < 0.00001, and amongst the 28 observations used to fit the model, 28 (or 100.0%) were correctly classified. The correct classification rate was evaluated using cross-validation, and 100% was obtained (28 of 28). Figure 5 plots the relationship between log k' (IAM) vs. TPSA, considering OSC distribution using BBB permeation as the grouping classification variable. Therefore, it is possible to observe two different groups: BBB permeating organosulfur compounds (group 1) and BBB notpermeating organosulfur compounds (group 0).

IAM Phase Ability to Assess OSC Skin Permeability
To study the influence of OSC chemical structures on skin permeability a PLS analysis was carried out. For this purpose, several calculated values obtained from the SwissADME web tool were used, including log K P , a skin permeability coefficient resulting from a quantitative structure-property relationship model built by Potts and Guy (1992). As in the previous sections, we checked whether our OSC data set belonged to the data predictor space used to build K P predictions values in the SwissADME web tool. These results are shown in the Supplementary Material section.
As an indicator of polarity and size, and given that these properties showed influence in other skin permeability in silico models (Guth et al., 2014), the TPSA descriptor was also FIGURE 5 | Scatterplot showing the relationship between log k' (IAM) vs. TPSA, considering OSC distribution using BBB permeation as the grouping classification variable. BBB permeating organosulfur compounds (group 1), and BBB not-permeating organosulfur compounds (group 0). considered in this model as a predictor. The other descriptors used here were Tm_1 (WHIM descriptor) and TE1 (Charge descriptor), which resulted from a stepwise regression analysis (SRA) and VIP (variable importance in the projection) parameter. Statistical significance of the QSRR-PLS derivatives models was evaluated by means of variance of the matrices X and Y (R2X and R2Y), the square correlation coefficient (R 2 ), standard deviation (RMSS), and the statistical F. The predictive ability was evaluated by the cross-validation coefficient (Q). The complete set of OSCs was used, resulting in two statistically significant components (CPLS), whose parameters are detailed in Table 3. The obtained regression coefficients are shown in Table 4. Graphically, a high agreement between observed log K P skin permeability values and the corresponding calculated ones derived from the PLS-1 model can be seen in Figure 6. There was also an approximately linear tendency of the residuals, when observing the normal probability plot of the said values (data not shown), which indicates that the error terms were normally distributed, except for isoalliin, which showed more than 2 standard deviations from the central tendency of the remaining OSCs.

DISCUSSIONS
The main purposes of the present study were, first, to build different QSRR models to obtain some insight into the mechanism of IAM retention for the 28 phytochemicals with sulfur-containing groups under study. Secondly, the present study aimed to establish whether, and to what extent, the IAM-chromatographic measures for the target OSCs, serve to assess ADME parameters. Given the chemical nature of the phospholipids that make up the IAM columns' stationary phase, it was only logical to, first, study lipophilicity parameters associated with OSCs. It is well-known that passive absorption of orally consumed substances is related to drug lipophilicity, commonly expressed as the logarithm of the n-octanol/water partition coefficient, log P. And this partition phenomenon is a driving force for the membrane passage (Grumetto et al., 2015b). So, in order to evaluate IAM phase ability in assessing OSC lipophilicity, we studied the relation between log k'IAM values and calculated partition coefficient log P values, obtained from the SwissADME web tool (iLOGP). The regression model's obtained results showed that even though there was a suitable fitting between the iLOGP and log k'IAM variables (both physicochemical parameters carry similar information regarding lipophilicity), they do not encode the same information regarding specific interactions, as shown below. This fact agrees with the other reports (Grumetto et al., 2015a;Russo et al., 2017) which proposed that when a compound is ionizable and/or presents a polar surface area ≠ 0, the interactions with the two partition phases, i.e., n-octanol and phospholipids, are differently affected by electrostatic and polar interaction forces. Thus, it may be concluded that the balance between polar and hydrophobic effects is not the same toward the n-octanol and IAM phases. Figure 1 supports this conclusion, where these "extra" intermolecular interactions become indirectly visible because of the observed divergence between the different classes of compounds. In view of this, it seems reasonable to consider a common interaction mechanism between the compounds containing sulfoxide groups, in which some of them, in turn, present charged groups and different pKa values (see Table 2) and the phospholipid constituents of the IAM stationary phase. This fact explains the modest fit observed between the log k' (IAM) and iLOGP parameters, indicating the octanol partition system's inadequacy to mimic the complex FIGURE 6 | Relationship between actual and calculated log K P skin permeability values using the PLS approach.
Frontiers in Chemistry | www.frontiersin.org September 2021 | Volume 9 | Article 690707 interactions between the IAM phase and the ionized forms of the OSCs under study. At this point, it is important to highlight that we did not address log D evaluation because, at the experimental pH used during OSC chromatographic analysis, the log D and log P values were the same. All these results gave us the guidelines to further explore other molecular descriptors to explain OSC-retention behavior in the IAM phase adequately. From a stepwise regression analysis and using the VIP (variable importance in the projection) parameter, different descriptors were selected, and an MLR, as well as a PLS model, were built based on the best subset of structural parameters. From the analysis of the results obtained, and as previously discussed, the specific polar interactions were better modeled by two electronic descriptors, the Mulliken electronegativity and the relative negative charge (RNCG_1). These results are in agreement with the findings by other authors, who have shown that the polar/electrostatic molecular properties affect the phospholipid-analyte partition (Grumetto et al., 2015b;Russo et al., 2017). The presence of rotatable bonds descriptor in the model requires a separate discussion. In the most general sense, it may be considered as conveying information about molecular flexibility, and interestingly, this property appears to be also involved in the IAM retention process. The molecular flexibility plays an important role when analyzing orally administered drugs, since increased rotatable bonds count (meaning a greater molecular flexibility) has been reported to negatively affect the permeation rate (Li et al., 2003). In this way, our results agree with the latter, since the negative sign of the rotatable bonds variable in the developed models suggests that the higher the molecular flexibility, the weaker the interaction with the IAM phase, and consequently, a shorter retention time could be observed. However, it must be considered that molecular rigidity is a much more complex issue than the simple counting of rotatable bonds. Thus, this molecular property's appearance in the models suggests the importance of conformational constraints in the OSC-IAM partition.
Once the molecular factors that govern OSCs' retention in the IAM column were elucidated, we assessed the predictive performance of IAM retention data on some permeability parameters obtained from the SwissADME web tool, and the PreADMET web-based application. Accordingly, the human gastrointestinal absorption percentage (HIA%), blood-brain barrier (BBB) permeation, and skin permeability (Log K P ) were investigated. It should be highlighted that there are no experimental permeability data for these OSCs in the literature. To assure the reliability of these predicted values, we first checked our OSC data set belonging to each predictor space that supports the corresponding HIA, BBB passage, and K P skin values obtained from the web servers.
Regarding HIA (%), it has been previously shown by other authors that the nonlinear approaches outperform linear models in explaining HIA% data variance, since linear models tend to overestimate the HIA % for low permeable drugs, but this bias is not so marked in nonlinear models (Talevi et al., 2011). Given that an MLR model based on HIA (%) vs. log k' (IAM) showed a poor statistical correlation between the target compounds, and based on OSC disposition when plotting log k' (IAM) vs. HIA (%) (Figure 3), a nonlinear estimation was performed to model the data. As proposed by Talevi et al., the nonlinear estimation explained better the OSC-IAM partition. Moreover, Yen et al., 2005 have discussed that when molecules contain one or more ionizable groups, their lipophilicity changes with respect to pH resulting in a nonlinear relationship with gastrointestinal absorption. However, considering the results from the QSRR models and based on the OSC scatterplot of HIA (%) vs log k' IAM, it is possible to see that log k' (IAM) describes better the neutral compounds' behavior (Ermondi et al., 2018). But it is not possible to achieve a successful prediction of HIA permeability using only log k' IAM for all 28 OSCs.
A separate modeling of ionizable compounds shows that the total polar surface area (TPSA_1) explains better the behavior of sulfoxides-and cysteine-containing compounds, whose HIA permeation is reduced compared to the other OSCs. This descriptor has been already proven as a fundamental molecular parameter to describe the electrostatic/polar interaction in the absorption process (Grumetto et al., 2015b;Russo et al., 2017). Singh and Singh (2010) have found similar results with regard to HIA estimation of ACSOs, which showed moderate permeation capabilities, and they also discussed the importance of polar surface interactions in this matter. In short, membrane partitioning, during passive passage, has been exposed as a twostep process. In the first step, lipophilicity impulses the compounds into the phospholipid bilayer. On the other side, in the second step, the drug is transferred into the phospholipid bilayer's interior, and the main interactions between the bilayer and polar parts of a solute are related to hydrogen bonding and polarity. Therefore, the permeation rate depends mainly on simple molecular descriptors such as hydrogen bond capacity, lipophilicity, and size and charge of the molecule. In this sense, the polar surface area (PSA), which is generally assumed to be related to hydrogen bonding capacity, becomes an excellent predictor of passive membrane permeability, especially for ionizable and/or polar compounds (Artursson and Bergstrom, 2003).
On another note, organosulfur blood-brain barrier permeation (OSC-BBB) was studied using a QSAR model based on the obtained IAM retention parameters. In this case, a linear discriminant analysis (LDA) was used to fit the classification of OSCs as BBB(+) or BBB(-) penetrators of the BBB. The use of the log k' (IAM) parameter was necessary to obtain a reasonable classificatory model, but it is not enough to fully explain the OSC-BBB permeation. Thus, the combination of log k' (IAM) and TPSA yielded the best OSC-BBB model classificatory. TPSA has also been related to BBB permeation models in several investigations (Luco, 1999;Luco and Marchevsky, 2006;Janicka et al., 2020) showing that the polar surface area is one of the most important factors determining the passage of molecules through the BBB.
Finally, the IAM phase ability to assess OSC skin permeability, expressed as the permeability coefficient Log K P skin permeability , was evaluated by a PLS model. As found in the previous models presented here, lipophilicity and polar interactions would appear to affect the skin permeation ability of the target OSCs, also in agreement with other authors (Barbato, 2006). But it is noteworthy that the other two molecular features seem to influence the skin biopartitioning phenomenon. Tm_1 (T total size index/weighted by mass) is a WHIM-weighted holistic invariant molecular-descriptor, that accounts for 3D molecular size and atom distribution (Todeschini and Gramatica, 1997), and TE1 is a topographic electronic descriptor. These descriptors, therefore, show the influence of molecular shape, size, and charge in OSC permeation through the skin.
In relation to ADME data about these compounds, it has been argued that allicin, one of the major bioactive garlic-related compounds, can easily permeate cell membranes of phospholipids bilayers, where it carries out its activity intracellularly and interacts with the SH group. Even though this has high permeability, allicin is rapidly metabolized and transformed into other active compounds. Other OSCs such as SAC, diallyl disulfide, and diallyl trisulfide have also shown high cellular permeability, although their permeability passage varies from each compound and tissue (Rahman, 2007). On another note, vinyldithiins' pharmacokinetic properties were investigated. 1,3-Vinyldithiin seems to be less lipophilic and was rapidly eliminated from serum, the kidney, and fat tissue, whereas 1,2-vinyldithiin is more lipophilic and showed a tendency to accumulate in fat tissues (Egen-Schwind et al., 1992). Torres-Palazzolo et al. (2018) showed that allicin, DAS, and 2VD would be able to permeate less than 5% through the intestine. Something that has been extensively discussed is the rapid reaction of OSCs with blood components, resulting in different metabolites with possible bioactive effects in the organism (Lawson and Gardner, 2005). To summarize, something to highlight from the models built here in relation to the different permeation capacities of the OSCs under study is the possibility of delving into the molecular mechanisms that govern these properties and ends up affecting the variations in OSC permeation rates in different tissues. We could observe the influence of the hydrophobic factors that govern OSC membrane passage, as well as the high effect of simple molecular descriptors such as hydrogen bond capacity, size, and charge of the molecules, particularly considering the polar surface area-related parameters, which greatly affect the ionizable compounds' passage. It is necessary though, to delve deeper in organosulfur compounds' ADME properties, for example, by comparing our results with cellbased and tissue-based models to validate the QSAR models proposed here for future uses.

CONCLUSION
To assure OSC bioefficacy, a first evaluation must be conducted in terms of the pharmacokinetic potential of the said phytochemical agents. So far, our results proved that an in silico-based QSAR modeling was a convenient approach to assess different OSC permeability parameters. All this information serves to give an insight into the mechanisms implied during the OSC-membrane interaction. In brief, the evidence presented here suggests that in general lines, organosulfur compounds with sulfide groups showed a higher chromatographic retention using the IAM column, and accordingly, a higher permeation capability, which was mainly associated with lipophilic/polar interactions. On the contrary, ionizable OSCs showed a lower permeation capability as a result of different membrane-analyte interactions. Therefore, it may be concluded that the ADME properties under study are strongly dependent on hydrophobic factors as expressed by log k' (IAM), which provides evidence for the great potential of the IAM phases in the development of QSAR models.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding authors.