**Biomimetic Artificial Membrane Permeability Assay over Franz Cell Apparatus Using BCS Model Drugs**

**Leonardo de Souza Teixeira 1,**† **, Tatiana Vila Chagas 2,**† **, Antonio Alonso 3 , Isabel Gonzalez-Alvarez 4 , Marival Bermejo 4 , James Polli <sup>5</sup> and Kênnia Rocha Rezende 2, \***


Received: 8 September 2020; Accepted: 10 October 2020; Published: 19 October 2020

**Abstract:** A major parameter controlling the extent and rate of oral drug absorption is permeability through the lipid bilayer of intestinal epithelial cells. Here, a biomimetic artificial membrane permeability assay (Franz–PAMPA Pampa) was validated using a Franz cells apparatus. Both high and low permeability drugs (metoprolol and mannitol, respectively) were used as external standards. Biomimetic properties of Franz–PAMPA were also characterized by electron paramagnetic resonance spectroscopy (EPR). Moreover, the permeation profile for eight Biopharmaceutic Classification System (BCS) model drugs cited in the FDA guidance and another six drugs (acyclovir, cimetidine, diclofenac, ibuprofen, piroxicam, and trimethoprim) were measured across Franz–PAMPA. Apparent permeability (Papp) Franz–PAMPA values were correlated with fraction of dose absorbed in humans (Fa%) from the literature. Papp in Caco-2 cells and Corti artificial membrane were likewise compared to Fa% to assess Franz–PAMPA performance. Mannitol and metoprolol Papp values across Franz–PAMPA were lower (3.20 × 10 <sup>−</sup><sup>7</sup> and 1.61 <sup>×</sup> <sup>10</sup> −5 cm/s, respectively) than those obtained across non-impregnated membrane (2.27 × 10 <sup>−</sup><sup>5</sup> and 2.55 <sup>×</sup> <sup>10</sup> −5 cm/s, respectively), confirming lipidic barrier resistivity. Performance of the Franz cell permeation apparatus using an artificial membrane showed acceptable log-linear correlation (R <sup>2</sup> = 0.664) with Fa%, as seen for Papp in Caco-2 cells (R <sup>2</sup> = 0.805). Data support the validation of the Franz–PAMPA method for use during the drug discovery process.

**Keywords:** Franz–PAMPA; BCS drugs; biomimetic membrane; Franz cell; passive drug transport

#### **1. Introduction**

Favorable absorption, distribution, metabolism, and excretion (ADME) of orally administrated drugs are essential for therapeutic activity in vivo. Poor oral bioavailability contributes to a very high failure rate during pre-clinical drug development [1,2]. In this regard, the Biopharmaceutic Classification System (BCS) proposed by Amidon and co-workers [3] have been widely used as an important tool to support early drug development [4–6]. For orally administered drugs, gastrointestinal physiology is a key factor impacting on the rate and extent of drug absorption [7]. Transcellular passive diffusion across membranes is the major route and is governed by several molecular properties such as partition and distribution coefficient, as well as molecular weight [8,9].

Currently, important tools based on physicochemical properties and in vitro assays are used to predict in vivo gastrointestinal absorption [10]. In vitro methodologies include animal [11,12] or human tissues [13], cultured cells [14,15] and artificial membranes [16–18]. The Caco-2 cell monolayers in vitro model is thoroughly studied and generally mimics major transport pathways in the gastrointestinal tract [19]. However, this method is limited by long cell growth and differentiation cycles, risks of microbial contamination, and high implementation costs [19–21].

Cell-free permeation systems using artificial membranes are gaining progressively more interest as an alternative model to cell-based systems that can be simpler, less time consuming, and cost-effective [22,23]. Depending on the composition of the barrier, it can be classified as biomimetic barrier which is constructed from (phospho)lipids or, alternatively, from non-biomimetic barrier containing dialysis membrane [24].

Particularly, there is a growing interest in PAMPA studies with direct comparisons to Caco-2 cells using a consistent number of drugs displaying equally well prediction of in vivo data between them [25]. In this regard, major differences of key components amid cell-free membranes currently used in permeability systems was highlighted by Berben et al. (2018) [23].

Here, a previously validated biomimetic artificial permeability membrane comprising of a microfilter impregnated by a phospholipid solution [5] was mounted on horizontal Franz-cells diffusion chambers (Microette™, Hanson Research) [20]. This new setup approach, herein called Franz–PAMPA (Figure 1), was challenged to assess permeability of BCS model drugs simulating gastrointestinal permeation. Therefore, the aim of this study was to validate this Franz–PAMPA system by evaluating the correlation power between apparent permeability (Papp) for BCS model drugs to their fraction of drug absorbed (Fa%) in humans for rapid and reliable information about passively transported drugs [25,26].

**Figure 1.** Franz cells apparatus (Microette™, Hanson Research) mounted with a previously validated PAMPA membrane from Corti et al. [5,16] for simulating gastrointestinal permeation. (**A**) System control for injection pistons of upper chambers; (**B**) upper and lower diffusion chambers with temperature and stirring control; (**C**) automated module for sampling and collection.

#### **2. Materials and Methods**

#### *2.1. Materials*

Membrane supports were purchased from Millipore® (Mixed Cellulose Esters VCWP 047000; 0.1 µm × 47 mm, white plain, New York, NY, USA). All 19 compounds for permeation studies (acyclovir, amoxicillin, atenolol, caffeine, cimetidine, diclofenac, furosemide, hydrochlorothiazide, ibuprofen, mannitol, metoprolol, naproxen, piroxicam, propranolol, ranitidine hydrochloride, trimethoprim, and verapamil hydrochloride) were of analytical grade and kindly supplied by ICF (Pharmaceutical Sciences Institute, Goiânia, Brazil). All organic solvents were of HPLC grade and solid reagents were of analytical grade.

The spin labels 5-doxyl stearic acid (5-DSA) and 16-doxyl stearic acid (16-DSA) used for electron paramagnetic resonance (EPR) spectroscopy were purchased from Sigma-Aldrich Chem Co. (St. Louis, MO, USA). The spin labels 1-palmitoyl-2-stearoyl-(5-doxyl)-sn-glycero-3-phosphocholine

(5-PC) and 1,2-dipalmitoleoyl-sn-glycero-3-phosphocholin (16-PC) were purchased from Avanti (Avanti Polar Lipids, Inc., Alabaster, AL, USA).

#### *2.2. Methods*

#### 2.2.1. Impregnation of Membrane Support

Membranes were impregnated by immersion for 60 min (22 ± 1 ◦C) with a lipid solution (mixture of phospholipids), as previously reported [5]. Briefly, the lipid phase solution for impregnation was a mixture of 1.7% phospholipids (Lipoid® E 80, Ludwigshafen, Germany), 2.1% cholesterol (Sigma–Aldrich Chemical Co., Milan, Italy), and 96.2% n-octanol (Synth, Diadema, Brazil). Excess lipid was absorbed with cellulose filter paper over 30 min. Next, all impregnated membranes (N = 20) were weighed on a microanalytical scale (Mettler Toledo, mod. XPE56DR, Columbus, OH, USA) and evaluated to check for its accuracy (211.2 mg ± 6.0%). Prior to use, impregnated membranes were protected from moisture atmosphere and refrigerated (−8 ◦C, 24 h). It is worth mentioning that all membranes were stabilized prior to use. Stabilization was confirmed by EPR spectra which did not show any signals of physicochemical degradation: none of membranes showed any difference on <sup>14</sup>N-hyperfine coupling constant value (14.8 G) demonstrating its stability [25]. EPR signals were compared just after 24 h of refrigeration and post-run permeability studies as well as after a month of refrigerated storage time (*data not shown*).

#### 2.2.2. Electronic Paramagnetic Resonance (EPR)

The biomimetic membranes were impregnated, as described above. Spin labeling technique was employed to examine the conformational structure of the membrane using 5-DSA or 16-DSA. EPR was performed using a Bruker ESP 300 spectrometer (Bruker, Rheinstetten, Germany) equipped with an ER 4102 ST resonator. The instrument settings were microwave power of 2 mW; modulation frequency of 100 KHz; modulation amplitude of 1.0 G; magnetic field scan of 100 G; sweep time of 168 s; and a detector time constant of 41 ms. EPR spectral simulations were performed using the nonlinear-least-squares (NLLS) program for an isotropic model. The biomimetic membrane was introduced into flat, quartz EPR cell to perform the EPR measurements at room temperature (~25◦C).

#### 2.2.3. Permeation Studies

Permeation studies were performed using a Franz vertical diffusion cell (MicroettePlus, Hanson Research, CA, USA). Impregnated artificial membranes (Franz–PAMPA) were positioned between upper and lower part of diffusion cells and, the donor (1 mL) and receptor (7 mL) compartments holding phosphate-buffered solution (PBS) pH 7.4 (USP 32). In order to minimize the unstirred water layer (UWL), receptor compartment media was stirred (500 rpm). The temperature was kept constant (37.0 ± 0.5 ◦C). Each drug (n = 3) was added in the donor compartment at a fixed concentration (=10 mg/mL). One milliliter of saturated drug solutions was transferred to the donor compartments and capped to prevent evaporation. The experiments were performed under 'infinite dose' conditions [26,27], except for caffein, metoprolol, propranolol, naproxen, ranitidine, and atenolol (D<sup>0</sup> ≤ 0.01). Individual drug solubility is further shown in results section. Metoprolol was used as a low/high BCS permeability class boundary reference drug for the Franz–PAMPA assay [28].

Samples from permeation studies were collected during 12 h (0.25; 0.5; 1.0; 2.0; 3.0; 4.0; 5.0; 6.0; 10.0, and 12.0 h) and analyzed by HPLC (Shimadzu Class VP; Kyoto, Japan or Agilent 1220, Santa Clara, CA, USA) according to official compendiums (USP 32 or Brazilian Pharmacopeia 4th edition). The sampling volume was immediately replaced with the same volume of fresh PBS prewarmed solution at 37◦ ± 0.5 ◦C. Calibration curves were performed at least at three concentration levels for each drug tested, in a GLP-accredited laboratory (Institute of Pharmaceutical Sciences, Goiânia, Goiás, Brazil). The validated chromatographic conditions used for the drug permeability assay are given in Table 1.


#### **Table 1.**Pharmacopeial methods applied on drug analysis and their respective limit of quantification (LOQs).

#### 2.2.4. Permeability Calculations

The diffusion area (*A*) was calculated from the radius of the Franz cell and was 1.77 cm<sup>2</sup> . Flux through membrane to receptor compartment (*J*; µg/cm<sup>2</sup> /s) was calculated by dividing the amount of drug accumulated in the receptor compartment by *A*. The Fick's first law was derived to calculate flux (*J*) at steady state (Equation (1)):

$$J = dQ \sharp dt^\* A \tag{1}$$

where *dQ* is the amount of drug across the membrane (in moles), *dt* the permeation time (in seconds), and *A* the diffusion area (in cm<sup>2</sup> ). Note that *J* was obtained from the slope of the curve at steady state from typical mean cumulative concentration-time plots (minimum of triplicates), as further shown in results section (Figure 2). Coefficient of variation (CV) of flux for each drug was also measured.

τ **Figure 2.** Experimental (solid line) and best-fit (empty circles) EPR spectra for several spin labels in BAMPA. The isotropic 14N-hyperfine coupling constant, a<sup>0</sup> , showed equal spectra values of 14.8 G, consistent with a spin label in phospholipidic bilayer of eukaryotic cells. The rotational correlation time value, τC, is also showed. 5-DSA or 16-DSA: 5- or 16-doxylstearic acid); 5-PC or 16-PC:5-0r 16-phosphatidylcholine).

The apparent permeability (*Papp*) was calculated normalizing the flux (*J*) over the drug concentration in the donor compartment *C*0, as described by the following Equation (2):

$$Papp = f \pounds\_0 \tag{2}$$

<sup>−</sup> τ

This approximation was used in all cases, even when sink conditions do not hold and donor concentrations change with time, as already described for some experiments [29]. In addition, the following equation was used to account for the fact that in most cases sink conditions were not maintained [30].

$$C\_{\text{nucier},t} = \frac{Q\_{\text{total}}}{V\_{\text{nucier}} + V\_{\text{dona}}} + \left( (\mathbb{C}\_{\text{nucier},t-1} \cdot f) - \frac{Q\_{\text{total}}}{V\_{\text{nucier}} + V\_{\text{dona}}} \right) e^{-P\_{eff} \cdot S \cdot \left( \frac{1}{V\_{\text{nucier}}} + \frac{1}{V\_{\text{dona}}} \right) \cdot \Delta t} \tag{3}$$

τ

−

where *Creceiver*,*<sup>t</sup>* is the drug concentration in the receiver chamber at time *t*, *Qtotal* is the total amount of drug in both chambers, *Vreceiver* and *Vdonor* are the volumes of each chamber, *Creceiver,t*−*<sup>1</sup>* is the drug concentration in receiver chamber at previous time, *f* is the sample replacement dilution factor, *S* is the surface area of the membrane, ∆*t* is the time interval and *Pe*ff is the permeability coefficient. This equation considers a continuous change of the donor and receiver concentrations, and it is valid in either sink or non-sink conditions. The curve-fitting is performed by non-linear regression, by minimization of the sum of squared residuals (*SSR*), where:

$$\text{SSR} = \sum \left[ \mathbb{C}\_{r,i,obs} - \mathbb{C}\_{r,i}(t\_{end,i}) \right]^2 \tag{4}$$

*Cr,i,obs* is the observed receiver concentration at the end of interval *i*, and *Cr,I(tend,i)* is the corresponding concentration at the same time calculated according to Equation (3) [29].

Classification as high permeability was established if the calculated permeability (under sink or non-sink conditions was higher than 0.8\* metoprolol permeability [31].

The in vitro permeability (Papp) of each drug studied was compared to in vivo absorption in humans (Fa%), Papp in Corti artificial membrane [16], and Papp in Caco-2 cells.

#### **3. Results**

#### *3.1. EPR Analysis and Membrane Stability*

The Franz–PAMPA was characterized by EPR spectroscopy of lipid spin labels of doxyl class. The spectra showed a movement consistent with lipid bilayer (Figure 1). Two analogs of stearic acid, 5-DSA and 16-DSA, and two analogs of phosphatidylcholine, 5-PC and 16-PC, having the nitroxide radical positioned at the 5th and 16th carbon atom of the acyl chain, respectively, were used to examine the molecular dynamic at two regions into the bilayer. The EPR spectra of these four spin labels are shown in Figure 2.

The EPR parameter—isotropic <sup>14</sup>N-hyperfine coupling constant, a0—increased with increasing dielectric constant (i.e., solvent polarity) in which the nitroxide radical is dissolved. The measured value of 14.8 G is consistent with a spin label in a membrane [32]. The spin labels 5-DSA and 5-PC with the nitroxide moiety in the region near the polar head group of the bilayer showed more restricted rotational motion relative to their positional isomers 16-DSA and 16-PC, in which the nitroxide radical is more deeply inserted in the hydrophobic core. These results indicate the existence of a gradient of flexibility along the acyl chain, with more restricted motion in the polar region. This pattern is consistent with the properties of lipid bilayers from eukaryotic cells. The rotational motion at the polar interface of the membrane was more restricted for the spin label analog of phosphatidylcholine (5-PC) with <sup>τ</sup><sup>C</sup> of 14.2 <sup>×</sup> <sup>10</sup>−<sup>10</sup> s than for the stearic acid one (5-DSA) whose <sup>τ</sup><sup>C</sup> was of 8.4 <sup>×</sup> <sup>10</sup>−<sup>10</sup> s (Figure 2).

Membrane barriers from similar models such as PAMPA and PVPA have been proven to be stable in a pH range from 2 to 8 [33]. Here, EPR spectra were also recorded before and after permeation studies to check for the integrity of biomimetic membranes. No leaching of barrier-constituents such as phosphatidylcholine and lipids into the donor compartment could be evidenced as none of membranes showed any difference on <sup>14</sup>N-hyperfine coupling constant value (14.8 G) demonstrating its stability [34]. Likewise, using the same chemical composition as Corti (2006) [5], acidic and basic drugs also showed pH-dependent permeability according to the pH partition theory [25,35]. Accordingly, close Person's correlation coefficient was seen (r = 0.7355) to our data from Franz–PAMPA *versus* PAMPA pH 7.4 data from literature.

In this regard, pHs of drug solutions were all measured to assure buffer capacity and drug stability. Some authors correlated membrane flux with the fraction absorbed in human, showing that the flux through the egg lecithin/dodecane membrane correlated better than octanol/water logD values with the fraction absorbed in humans [17]. Later, an in-depth investigation of pH impact on drug Franz–PAMPA

permeability will be necessary to increase the biomimetic and absorption predictive power of this method, although the study of this factor was beyond the scope of this work.

#### *3.2. Membrane Validation and Performance*

Studies here deals with a modified PAMPA method over Franz cell apparatus. The biomimetic membrane (Franz–PAMPA) has been previously described by Corti and coworkers [5] as a modified version from Kansy et al. (1998) [36]. Mannitol and metoprolol were used as a marker for the cutoff point between low and high permeability drugs.

The apparent permeability coefficient (Papp) values found for mannitol and metoprolol, over the lipid impregnated membrane (2.27 <sup>×</sup> <sup>10</sup>−<sup>5</sup> and 2.55 <sup>×</sup> <sup>10</sup>−<sup>5</sup> cm/s, respectively) were higher when compared to the non-impregnated one across Franz–PAMPA (3.20 <sup>×</sup> <sup>10</sup>−<sup>7</sup> and 1.61 <sup>×</sup> <sup>10</sup>−<sup>5</sup> cm/s, respectively), indicating the resistivity of the lipid membrane itself.

Membrane performance was assessed using 14 representative model drugs (Table 2) cited in the FDA BCS guidance [37]. Class I model compounds were caffeine, metoprolol, propranolol and verapamil. Class II model compounds were diclofenac, ibuprofen, naproxen and, piroxicam. Class III model compounds were atenolol, cimetidine, ranitidine, and trimethoprim. Class IV model compounds were acyclovir, furosemide and hydrochlorothiazide. This classification was based in the permeability class indicated in the FDA guidance [37] and, on solubility from literature [28].

Cumulative drug transport through Franz–PAMPA was plotted over 12 h and the Papp was calculated from the slopes obtained from linear regressions (Figure 3, Table 2). Of the 21 compounds studied by Corti and coworkers and of the 14 compounds studied here, there were 11 common compounds tested in both studies: acyclovir, atenolol, caffeine, cimetidine, furosemide, hydrochlorothiazide, metoprolol tartrate, naproxen, propranolol, ranitidine, and trimethoprim. For these drugs Caco-2 Papp values were also surveyed from literature and compared here (Table 2). − − − −

For high permeability drugs (BCS I and II, Table 2), the Papp values showed to be in the range of 4.6–75.2 <sup>×</sup> <sup>10</sup>−<sup>6</sup> cm/s in Franz–PAMPA. For Caco-2 assay, values were narrower (15.8–52.5 <sup>×</sup> <sup>10</sup>−<sup>6</sup> cm/s) and, for Corti membranes they were most narrow (39.7–48.8 <sup>×</sup> <sup>10</sup>−<sup>6</sup> cm/s).

For low permeability drugs (BCS III and IV, Table 2), the Papp coefficient found were consistently much lower than high permeability drugs. Franz–PAMPA, Caco-2, and Corti membrane provided value ranges of 0.2–24.6 <sup>×</sup> <sup>10</sup>−<sup>6</sup> cm/s, 0.1–83.0 <sup>×</sup> <sup>10</sup>−<sup>6</sup> cm/s, and 3.2–45.5 <sup>×</sup> <sup>10</sup>−<sup>6</sup> cm/s, respectively. Permeability of most drugs tested here showed Papp <sup>&</sup>gt;1.0 <sup>×</sup> <sup>10</sup>−<sup>5</sup> cm/s (Figure 4).

**Figure 3.** *Cont.*

**Figure 3.** Cumulative transport (µg/mL, h) of drugs across Franz-PAMPA: (**a**) BCS class I; (**b**) **Figure 3.** Cumulative transport (µg/mL, h) of drugs across Franz-PAMPA: (**a**) BCS class I; (**b**) BCS class II; (**c**) BCS class III; and (**d**) BCS class IV. Permeability was calculated from the linear portion (R<sup>2</sup> ). Data are presented as mean ± SD, n = 3 − − −

−

− −

− −

● **Figure 4.** Permeability values for Franz-PAMPA *versus* D<sup>0</sup> drugs using metoprolol as reference drug. Most drugs (•; 14 out of 19) were in accordance with previous Biopharmaceutic Classification System (BCS) classification. Some of them (; 5 out of 19) disagreed.

■

−

−


**Table 2.**Papp calculated values in Franz-PAMPA and using Non-Sink Artursson equation. BCS classification of studied drugs and literature data to all other parameters.

nC = non classified 1 Yamashita et al., 2000 (17) and Zhu et al., 2002 [38]; 2 Corti and co-workers, 2006 [5]; 3 Di Cagno et al. [22]; 4 Lindenberg et al. 2004 [39], 5 Kasim et al., 2004 [28] (AT) actively transported drugs. (HP) high permeability drug [38]—Data not available for non-sink calculations.

Typically, PAMPA methods are affected by high variability, and therefore, data can be somehow noisy for poorly permeable drugs. Variability is also an issue that impacts permeability for Caco-2 [5] and other in situ [19] and in vivo [24] models. For low permeability drugs (Fa < 80%), Avdeef and coworkers (2003) [6] measured variability for more than 200 different drugs accounting for more than 600 measurements. Papp values close to 10 <sup>×</sup> <sup>10</sup>−<sup>6</sup> cm/s showed variability of around 10%. Such error can increase slightly for higher Papp values but is larger for Papp <sup>&</sup>lt; 0.1 <sup>×</sup> <sup>10</sup>−<sup>6</sup> (60%), with 0.01 <sup>×</sup> <sup>10</sup>−<sup>6</sup> values exhibiting variability of 100% or more. Although, currently BBB (blood brain barrier) [40] or Skin-PAMPA [41] methods can achieve higher precision and reproducibility with some other controlled protocols. Specific adjustments include setting incubation time as low as possible, increasing sensitivity of analytical methods, controlling membrane homogeneity either on the filter or among filters, besides the rationale for compounds dataset amongst others.

Likewise, permeability of small hydrophilic compounds is frequently underestimated in PAMPA since the membrane has hydrophobic nature besides being a cell-free system [42]. For the FDA-listed drugs, PAMPA Papp displayed values ranged from 0.00 to 2.35 <sup>×</sup> <sup>10</sup>−<sup>5</sup> cms−<sup>1</sup> , indicating it was not sensitive enough to discriminate and rank poorly permeable compounds. In contrast, Franz–PAMPA showed values in a wider Papp range of 0.4–68.1 <sup>×</sup> <sup>10</sup>−<sup>6</sup> cm/s. This could be tentatively explained due to the hydrophilic nature of membrane support and pH-dependent characteristics of the drugs [22,24,31]. Moreover, Franz cell stirring clearly reduces the unstirred water layer resistance in the system.

Additionally, variability of Papp values was also addressed by the calculation methods. A more sophisticated analysis is done using Artursson's equation [15] for sink and non-sink conditions as well as checking the impact of extracting a permeability coefficient from data that are not at true steady state and, thus, possibly impacted by dose depletion. Note that for both the sink and the non-sink equation, Papp values showed a particularly good correlation between them (0.8984). Similarly, Papp values obtained by us showed to be very alike to values calculated according to Artursson's non-sink equation (Table 2, Figure 5). The reason is that we used the same systematic procedure, i.e., the best fit method through the linear portion, to calculate all the slopes characterizing an accurate permeability flow, so that the impact from dose depletion is considered not above average. As a result, all drugs got the same BCS classification in both methods.

**Figure 5.** Papp calculations using equations by (**a**) per Artursson *non-sink versus sink* conditions and, (**b**) per Artursson *non-sink* compared to *sink* conditions best fit method of the linear portion.

In this context, Franz–PAMPA profile is mimicking biological permeation in a graphical pattern related to permeation through Caco-2 cells (R<sup>2</sup> = 0.826). Obtained Papp values *versus* fraction of dose absorbed in humans (Fa%) showed log linear correlation (Figure 6), as also described by Zhu et al. [38] when analyzing permeability performance of 93 commercial drugs as for artificial membranes. As expected, Franz–PAMPA also showed a significantly improved log linear correlation (R<sup>2</sup> = 0.6982) when actively transported compounds ranitidine, trimethoprim, and verapamil were not incorporated in the regression analysis. In contrast, the Fa% *versus*. Corti membrane correlation was linear (R<sup>2</sup> = 0.904). Such discrepancy from Franz–PAMPA and Caco-2 reveals that passive permeability of tested drugs

through Corti membrane was greater and better suitable, especially for low and moderate permeability drugs, as discussed elsewhere [39]. PAMPA and Caco-2 technique would be best suited for compounds with medium and high permeabilities. For low permeability compounds, small differences in measured Papp are expected to yield large differences in Fa% values resulting in imprecise measurements.

■ **Figure 6.** Demonstration of method suitability from Franz–PAMPA assay permeability and fraction of dose absorbed in humans (Fa%) compared to Caco-2 cells (o) and Corti membrane (). Actively transported drugs were removed for R<sup>2</sup> calculation. Corti membrane Papp ( ■ rane Papp () ) correlation to %Fa (R<sup>2</sup> = 0.890) was essentially unchanged.

Currently, a promising biomimetic barrier also adapted to Franz diffusion cells Permeapad™— [22] was reported for six drugs concurrent to our model (acyclovir, atenolol, caffeine, ibuprofen, and metoprolol). Even if a satisfactorily comparative analysis was not straightforward, BCS classification of most drugs (4 out of 5) showed to be identical with similar Papp rank order (Table 2).

#### **4. Conclusions**

− − The Franz–PAMPA method provided a permeability pattern similar to those from Caco-2. Methodologically, the advantages of Franz–PAMPA over Caco-2 are the lower costs and simplicity of membrane preparation (e.g., reagents and artificial membrane are commercially available). Furthermore, the method is very versatile and could be transformed in a high-throughput in vitro method to detect and classify compounds absorbed by passive diffusion.

Using metoprolol as a high permeability marker (Papp <sup>=</sup> 1.61 <sup>×</sup> <sup>10</sup>−<sup>5</sup> cm/s; Figure 1), seven drugs were classified as highly permeable (best fit method): atenolol, caffeine, cimetidine, diclofenac, ibuprofen, naproxen, and propranolol (Table 2). Only atenolol and cimetidine were misclassified as highly permeable drugs, relative to their prior literature classification as BCS 3 drugs.

Additionally, 10 out of 17 drugs were classified as low permeability drugs in Franz–PAMPA. Nevertheless, only naproxen, piroxicam, and verapamil (3 out of 10) had their permeability underestimated according to BCS, as they performed as low permeability drugs instead of BCS2 drugs.

Summing up, a potential limitation of our study is that the Papp values were calculated with an equation in which the underlying assumptions are constant donor concentration and sink conditions. In order to account for that, we also did the calculations to estimate permeability values under non-sink conditions. The obtained values are about the same compared with the true values (i.e., assuming donor concentration change and non-sink conditions). Although the relative estimation error does

change across high *versus* low permeability compounds [29], the practical implications for predicting oral fraction absorbed would only be a "shift" to the left on the abscissa. In the case of a direct correlation with Caco-2 values, it would be reflected in a different slope, but it would not change the significance of the regression line. In the case of the use of the apparent permeabilities for classification of compounds, the reference value of metoprolol is also underestimated, so the classification outcome would not be changed [29].

As a final comment, the ability of Franz–PAMPA to classify drugs was good and can be potentially challenged at different pH conditions to predict intestinal permeability of drugs showing passive transport. Eventually, the Franz–PAMPA cell diffusion can be modulated in lipid composition and may be a suitable alternative for studying other biological barriers such as blood–brain barrier, skin, and mucosal barriers as buccal or nasal. The current dataset adds valuable information for future analysis of drug-molecular interactions at the lipid layer and in silico model development. Additionally, all apparatus and supplies experimentally used on Franz–PAMPA are commercially available and affordable to facilitate drug discovery method application.

**Author Contributions:** Conceptualization, funding acquisition, project administration, K.R.R.; T.V.C. and L.d.S.T.; methodology, formal analysis, investigation, data curation, A.A. methodology, investigation, data curation, discussion; I.G.-A., M.B.; J.P. discussion, supervision. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by CNPq, FINEP and FAPEG.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Development of a Hierarchical Support Vector Regression-Based In Silico Model for Caco-2 Permeability**

**Giang Huong Ta 1 , Cin-Syong Jhang 1 , Ching-Feng Weng <sup>2</sup> and Max K. Leong 1, \***


**Abstract:** Drug absorption is one of the critical factors that should be taken into account in the process of drug discovery and development. The human colon carcinoma cell layer (Caco-2) model has been frequently used as a surrogate to preliminarily investigate the intestinal absorption. In this study, a quantitative structure–activity relationship (QSAR) model was generated using the innovative machine learning-based hierarchical support vector regression (HSVR) scheme to depict the exceedingly confounding passive diffusion and transporter-mediated active transport. The HSVR model displayed good agreement with the experimental values of the training samples, test samples, and outlier samples. The predictivity of HSVR was further validated by a mock test and verified by various stringent statistical criteria. Consequently, this HSVR model can be employed to forecast the Caco-2 permeability to assist drug discovery and development.

**Keywords:** intestinal absorption; intestinal permeability; human colon carcinoma cell layer (Caco-2); hierarchical support vector regression (HSVR)

#### **1. Introduction**

Clinically, the majority of drugs are orally administered [1]. Prior to reaching the blood circulation system, the administered pharmaceutical agents have to pass through the intestinal barrier via passive diffusion, active uptake, and/or efflux transport processes [2–4], as illustrated by Figure 10.2 of Proctor et al. [2]. In passive diffusion, drug molecules can permeate the epithelial cell layers through the transcellular pathway, in which they penetrate through the cell membrane, or the paracellular pathway, in which they can cross the epithelial cell layer through the tight junction between cells [5]. The significance of active transporters on intestinal absorption has been detailed elsewhere [6]. Principally, active transport can be modulated by the efflux transporters of the ATP-binding cassette (ABC) family as well the influx transporters of the solute carrier (SLC) family [6], of which the efflux transporters can pump the administrated drugs out of enterocytes, leading to the reduction of the accumulated concentration, whereas the influx can enhance the intestinal uptake, resulting in the increased drug accumulation [7]. Of various active influx and efflux transporters, P-glycoprotein (P-gp), also termed multidrug resistance 1 protein (MDR1/encoded by *ABCB1* gene), breast cancer resistance protein (BCRP/*ABCG2*), organic anion transporting polypeptide 2B1 (OATP2B1/*SLCO2B1*), and peptide transporter 1 (PEPT1/*SLC15A1*) play predominant roles in intestinal absorption [8].

Passive diffusion depends on a number of physicochemical properties, whereas active transport relies on the characteristics of specific binding sites on the transport proteins [9]. The uncharged and modest hydrophobic drugs such as testosterone [10] can permeate through the membrane. Conversely, it is very difficult for highly hydrophobic molecules to get across cells, since they can be adhered to the membrane [5]. On the other hand,

**Citation:** Ta, G.H.; Jhang, C.-S.; Weng, C.-F.; Leong, M.K Development of a Hierarchical Support Vector Regression-Based In Silico Model for Caco-2 Permeability. *Pharmaceutics* **2021**, *13*, 174. https:// doi.org/10.3390/pharmaceutics 13020174

Academic Editor: Maria Isabel Gonzalez-Alvarez Received: 27 December 2020 Accepted: 21 January 2021 Published: 28 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

hydrophilic drugs such as mannitol predominantly pass through the paracellular pathway [10].

Of various drug absorption, distribution, metabolism, elimination, and toxicity (ADME/Tox) properties, drug absorption plays a pivotal role in drug discovery, since they substantially contribute to the earlier preclinical go/no-go decisions for the drug candidates [10,11] to achieve the "fail fast, fail early" paradigm [12]. As such, numerous in vivo and in situ assays have been developed to evaluate the intestinal absorption [13,14]. For instance, the in situ single-pass intestinal perfusion (SPIP) model measures the appearance of the drug in plasma after intravenous and intraintestinal drug administration [13,15]. The drug is orally administrated or directly given into the intestine or stomach in some animal species in in vivo assay [13,14,16].

In addition to in vivo and in situ assays, various in vitro assays have been devised, since they have more advantages such as low cost and time efficiency as compared with their in situ and in vivo counterparts [15]. Of various in vitro assays to evaluate intestinal absorption, human colon carcinoma monolayer cells (Caco-2) [3], parallel artificial membrane permeability (PAMPA) [17,18], and Madin–Darby canine kidney cells (MDCK) [19] are most frequently used. In fact, a comprehensive drug absorption profile should include the Caco-2, MDCK, and PAMPA permeability data to explore drug solubility and bioavailability [20]. Moreover, Caco-2, which can be adopted to evaluate the drug permeability through the cytoplasm (transcellular uptake) or between cells (paracellular uptake) and active transport [6], has become the golden standard for predicting intestinal drug permeability and absorption because of its similarity in morphology and function with human enterocytes [21–23]. The Caco-2 protocol has been clearly described in detail by Hubatsch et al. As compared with the biological membrane, the Caco-2 system still suffers from a range of disadvantages such as high technical complexity, the limitations related to the differences between cell monolayers and intestinal membrane structurally and functionally [24], in addition to its long culture periods (21-24 days) with the significantly extensive costs, contributing to the major concerns in practical applications [21,25].

The Caco-2 permeability is normally expressed by the apparent permeability coefficient (*P*app), in which the drug solution is added to the apical side, viz. the donor compartment, and the *P*app value in the basolateral side, viz. the receiver compartment, is measured [23]

$$P\_{\rm app} = \frac{dQ}{dt} \times \frac{1}{\left(A \times \mathcal{C}\_0\right)}\tag{1}$$

where *dQ*/*dt* is the linear appearance rate of mass in the receiver solution transported during sink conditions, *A* is the membrane surface area, and *C*<sup>0</sup> is the initial concentration at the donor compartment [26]. However, it is not uncommon to observe in vitro permeability variations among different from research groups, because the cultured cells can vary based on culture conditions, passage number, monolayer age, seeding density, and stage of differentiation [27,28], as exemplified by those compounds listed in Table 3 of Lee et al. [29]. Furthermore, Yamashita et al. have found that the different pH values of apical medium and the different solvents can produce different drug absorption values [30]. For instance, the *P*app values of alprenolol are (6.06 ± 0.18) × 10-6 cm/s and (30.0 ± 1.8) × 10-6 cm/s at pH 6.0 and pH 7.4, respectively. More examples of *P*app variations at different pH values can be found in Table 1 of Yamashita et al. [30].

In silico technologies have become an essential component in drug discovery and development according to the fact that they can provide guidance in the early stages in the drug discovery process such as the activity classification (high/moderate/poor) or quantitative predictions [31,32]. As such, a great number of in silico models have been established to predict the ADME/Tox properties [33]. The relationship between biological activity and chemical characteristics can be established by quantitative structure–activity or structure–property relationships (QSAR and QSPR) [34]. Numerous QSAR models have been generated to predict Caco-2 permeability based on a variety of physicochemical and physiological descriptors [35–51]. Nevertheless, the difficulties in developing sound in silico models to predict the intestinal permeability still remain unanswered mainly due to the fact that Caco-2 permeability is a dramatically perplexing process that can take place through numerous non-linear routes (vide supra).

More specifically, the ABC transporters, which are efflux transporters, can reduce the drug absorption, whereas the SLC transporters, which are influx transporters, can enhance the drug uptake, leading to the decrease and/or increase of drug absorption, respectively. In fact, such controversy can establish a paramount barrier in model development. For instance, the number of aromatic rings (*n*Ar) can enhance the compound hydrophobicity [52] and facilitate the passive diffusion consequently. Conversely, *n*Ar is also an important feature for P-gp substrate recognition and modulates the compound efflux correspondingly [53]. Thus, *n*Ar can simultaneously affect the active efflux and passive diffusion.

It is exceedingly difficult, if not nearly impossible, to derive a robust in silico model, which can properly render the complex relationships between the selected descriptors and Caco-2 permeability. However, the hierarchical support vector regression (HSVR) scheme, which is an innovative machine learning-based scheme initially developed by Leong et al. [54], can properly address the complicated and varied dependencies of descriptors that, in turn, can be greatly contributed to its advantageous features of both a local model and a global model, namely wider coverage of applicability domain (AD) and a higher capability of prediction, respectively. When comparing with most theoretical models, which are vulnerable to the outliers that represent mathematic extrapolations, HSVR can still show consistent performance, as demonstrated elsewhere [1,54–57]. Herein, the objective of this study was to develop an in silico model based on the HSVR scheme to predict Caco-2 permeability in conjunction with previously published PAMPA permeability, intestinal absorption, and MDCK efflux in silico models [1,55,57] to facilitate drug discovery and development, since medicinal chemists can employ these models to predict the drug absorption of (virtual) hit compounds as well as drug metabolism and pharmacokinetics (DM/PK) scientists can adopt these models to prioritize the lead compounds.

#### **2. Materials and Methods**

#### *2.1. Data Collection*

The *P*app values were collected from the various sources after a comprehensive literature search [22,23,58–66]. Assay systems were carefully scrutinized to ensure data consistency, since various assay conditions such as pH value and solvent system, for example, can affect the Caco-2 permeability [30]. Only *P*app values, which were measured in the Hank's balanced salt solution (HBSS) buffer and 4-(2-hydroxyethyl)-1-piperazineethanesulfonic acid (HEPES) including ca. 1% dimethylsulfoxide (DMSO) at pH 7.4 were chosen in this study. The average *P*app value was selected to warrant better consistency in case there was more than one *P*app value for a given compound within a near range. Finally, 144 compounds were chosen in this study and their corresponding logarithm *P*app values, simplified molecular input line entry system (SMILES) strings, Chemical Abstracts Service (CAS) registry numbers, and references to the literature are listed in Table S1.

#### *2.2. Molecular Descriptors*

The density functional theory (DFT), Becke 3-parameter Lee–Yang–Parr (B3LYP) method was employed to do full geometry optimization by the Gaussian package (Gaussian, Wallingford, CT, USA) for all recruited samples with the selection of basis set 6-31G (*d,p*). The solvent system was taken into consideration by the polarizable continuum model (PCM) [67,68]. The atomic charges, upon which the dipole moments depend, were calculated by the molecular electrostatic potential (MEP) [69]. The frontier orbitals energies, namely the highest occupied molecular orbital energy (*E*HOMO) and the lowest unoccupied molecular orbital energy (*E*LUMO), molecular dipole (*µ*), as well as the maximum absolute component of *µ* (|*µ*|max) were also recovered from the optimization calculations.

In total, more than 100 descriptors, which feature one-, two-, and three-dimensional ones and can be categorized into a variety of classes consisting of topological descriptors, electronic descriptors, thermodynamic descriptors, structure descriptors, spatial descriptors, and *E*-state indices, were enumerated by Discovery Studio (BIOVIA, San Diego, CA, USA) and E-Dragon (available at the website http://www.vcclab.org/lab/edragon/). The logarithm of the *n*-octanol–water partition coefficient at pH 7.4 (log *P*) was calculated by *XLOGP3* of SwissADME (available at the website http://www.swissadme.ch/index.php). Furthermore, the cross-sectional area (CSA), which has been implicated in membrane permeability [70,71], was calculated using the method modified by Muehlbacher et al. [72]. The collected compounds were divided into 4 ion classes [73], namely zwitterion, base, acid, and neutral ions according to their p*K*<sup>a</sup> values. The neutral ions only have one p*K*<sup>a</sup> value, the zwitterion ions are those whose strongest acidic p*K*<sup>a</sup> values are larger than 7 and the strongest basic ones are smaller than 7, the acidic ions have all their p*K*<sup>a</sup> values smaller than 7, whereas the basic ions have all their p*K*<sup>a</sup> values larger than 7.

#### *2.3. Descriptor Selection*

Descriptor selection was initially executed by removing those descriptors missing more than one molecule or displaying little or no distinction among all molecules. Furthermore, the Spearman's matrix between calculated descriptors was constructed to minimize the chance of spurious correlations, and those descriptors with intercorrelation values of *r* <sup>2</sup> > 0.80 were discarded, since the threshold was proposed by Topliss and Edwards [74]. In this study, a more conservative value of *r* <sup>2</sup> ≥ 0.64 was taken to further ensure the quality of derived models.

Descriptor values can span a wide range due to their diverse nature (vide supra). It is of necessity to transfer descriptors into a more consistent range to decrease the chance of descriptors with broader ranges overriding those with narrower ranges [75]. Accordingly, descriptors were subjected to normalization by centering and scaling

$$\mathfrak{X}\_{ij} = \frac{\mathfrak{x}\_{ij} - \langle \mathfrak{x}\_{j} \rangle}{\sqrt{\sum\_{i=1}^{n} \left( \mathfrak{x}\_{ij} - \langle \mathfrak{x}\_{j} \rangle \right)^{2} / (n - 1)}} \tag{2}$$

where *xij* and *x*ˆ*ij* symbolize the *j*th original and normalized descriptors of the *i*th molecule, respectively; *xj* is the average value of the original *j*th descriptor; and *n* is the number of molecules.

The descriptor selection is of pivotal importance in the performance of QSAR models [76]. Thus, genetic function approximation (GFA) bundled in the QSAR module of Discovery Studio was used for the initial descriptor because of its effectiveness and efficiency [77]. The recursive feature elimination (RFE) scheme was adopted for additional selection, in which the model was repeatedly generated by all but one descriptor. The descriptor, which had the less contribution in predictive performance, was removed after ranking their contributions [78].

#### *2.4. Dataset Selection*

It is not uncommon to identify the outliers and remove them from data collection for model development [79]. As such, outliers were recognized by inspecting molecular distribution in the chemical space [80], which was created by principal components (PCs) using the Diverse Molecules/Principal Component Analysis embedded in Discovery Studio, followed by discovering the outliers.

The remaining molecules were arbitrarily allocated into the training set and test set with an about 4:1 portion as recommended [81] to generate and verify the built model, respectively, using the Diverse Molecules/Library Analysis function within Discovery Studio. Golbraikh et al. have postulated that a sound model can be resulted only when both samples in the training set and test set can show high levels of chemical and biological

similarity [82]. Thus, the data distributions in the training and test set were carefully checked to ensure the high similarity degrees biologically and chemically in both datasets.

2

0.65

0.85

2 2 0.65 and 0.20

**2021**, , x FOR PEER REVIEW 7 of 27

‒

#### *2.5. Hierarchical Support Vector Regression*

Leong et al. originally invented HSVR [54] which was evolved from support vector machine (SVM) proposed by Vapnik et al. [83]. Initially, SVM was designed for classification only and the regression function, termed as support vector regression (SVR), was introduced later [84]. HSVR has a higher level of predictivity and broader applicability domain (AD) as compared with SVR, since it can seamlessly combine the advantages of the local model and global model [56]. More significantly, the superiority of HSVR has been revealed by some studies [1,54–57].

The theory and fulfillment of HSVR have been delineated in detail elsewhere, and the schematic presentation of HSVR can be depicted by Figure 1 of Leong et al. [54]. Basically, an SVR ensemble (SVRE) is used to build an HSVR model, and SVR models in the ensemble are generated from different descriptor combinations and function as local models with their own ADs. Briefly, the svm-train module in *LIBSVM* (software available at http://www.csie.ntu.edu.tw/~cjlin/libsvm/) was employed to build various SVR models using those samples in the training set with different descriptor combinations and SVR run conditions. The module svm-predict in *LIBSVM* was adopted to validate the produced SVR models using the samples in the test set. Radial basis function (RBF) was the designated kernel function due to its simplicity and better functionality [85]. Both *ε*-SVR and *ν*-SVR regression functions were tested. The SVR runtime conditions including *ε*-SVR and *ν*-SVR, their associated *ε* and *ν*, the kernel width *γ*, and cost *C* were tuned by the grid-search technique. ‒

**Figure 1.** The chemical space spanned by three principle components (PCs) displays the distribution of the data samples in the training set (solid circle), test set (gray square), and outlier set (open triangle).

According to the principle of Occam's razor, i.e., the principle of parsimony, the number of descriptors selected to build SVR models should be minimized as much as possible. This principle was also applied to the construction of SVRE, which demanded the minimum number of ensemble members [86]. Initially, the combinations of two SVR models were adopted to generate the HSVR model; this process was repeated until the production of a predictive HSVR. Otherwise, the combinations of three- or even fourmember SVRE were used to develop the HSVR models if the two-SVR ensembles failed to perform well.

#### *2.6. Predictive Evaluation*

The residual yielded by the difference between the observed value (*yi*) and the predicted value (*y*ˆ*i*) for the *i*th molecule was computed based on the following equation:

$$
\Delta\_{\bar{i}} = \mathcal{y}\_{\bar{i}} - \mathcal{y}\_{\bar{i}} \tag{3}
$$

In addition, standard deviation (*s*), maximum residual (∆Max), root mean square error (RMSE), and mean absolute error (MAE) in a dataset with *n* samples were evaluated.

$$\text{RMSE} = \sqrt{\sum\_{i=1}^{n} \triangle\_i^2 / n} \tag{4}$$

$$\text{MAE} = \frac{1}{n} \sum\_{i=1}^{n} \left| \Delta\_i \right| \tag{5}$$

Various statistic metrics were adopted to evaluate the produced models. The squared correlation coefficients including *r* <sup>2</sup> and *q* 2 in the training set and external set, respectively, were computed by the following equation.

$$\langle r^2, q^2 = 1 - \sum\_{i=1}^n (\mathcal{Y}\_i - y\_i)^2 / \sum\_{i=1}^n (y\_i - \langle \mathcal{Y} \rangle)^2 \tag{6}$$

where h*y*ˆ*i*i represents the average predicted value, and *n* is the number of samples in the dataset. The derived models were subjected to the 10-fold cross-validation using the function embedded in *LIBSVM* to give rise to the squared correlation coefficient of 10-fold cross-validation *q* 2 CV. Another internal validation was carried out by the *Y*-scrambling test [87], in which the log *P*app values were randomly permuted and then reapplied to the previous developed model without altering the descriptors. This process was repeated 25 times as suggested [87] to generate the average squared correlation coefficient *r* 2 *s* .

The external dataset was evaluated predictivity by the squared correlation coefficients *q* 2 F1, *q* 2 F2, and *q* 2 F3, and the concordance correlation coefficient (*CCC*) [88–93] using *QSARINS* [94,95].

$$q\_{\rm F1}^2 = 1 - \sum\_{i=1}^{n\_{\rm EXT}} \left( y\_i - \mathfrak{H}\_i \right)^2 / \sum\_{i=1}^{n\_{\rm EXT}} \left( y\_i - \langle y\_{\rm TR} \rangle \right)^2 \tag{7}$$

$$q\_{\rm F2}^2 = 1 - \sum\_{i=1}^{n\_{\rm EXT}} \left( y\_i - \mathfrak{H}\_i \right)^2 / \sum\_{i=1}^{n\_{\rm EXT}} \left( y\_i - \langle y\_{\rm EXT} \rangle \right)^2 \tag{8}$$

$$q\_{\rm F3}^2 = 1 - \left[\sum\_{i=1}^{n\_{\rm EXT}} \left(y\_i - \hat{y}\_i\right)^2 / n\_{\rm EXT}\right] / \left[\sum\_{i=1}^{n\_{\rm EXT}} \left(y\_i - \langle y\_{\rm TR} \rangle\right)^2 / n\_{\rm TR}\right] \tag{9}$$

$$\text{CCC} = \frac{2\sum\_{i=1}^{n\_{\text{EXT}}} \left( y\_i - \langle y\_{\text{EXT}} \rangle \right) \left( \mathcal{g}\_i - \langle \mathcal{g}\_{\text{EXT}} \rangle \right)}{\sum\_{i=1}^{n\_{\text{EXT}}} \left( y\_i - \langle y\_{\text{EXT}} \rangle \right)^2 + \sum\_{i=1}^{n\_{\text{EXT}}} \left( \mathcal{g}\_i - \langle \mathcal{g}\_{\text{EXT}} \rangle \right)^2 + n\_{\text{EXT}} \left( \langle y\_{\text{EXT}} \rangle - \langle \mathcal{g}\_{\text{EXT}} \rangle \right)^2} \tag{10}$$

where h*y*TRi is the averaged observed values in the training set, h*y*EXTi and h*y*ˆEXTi are the averaged observed and predicted values in the external set, respectively; *n*TR and *n*EXT stand for the numbers of samples in the training set and external set, respectively.

In addition, some modified squared correlation coefficients *r* <sup>2</sup> were estimated [96,97]

$$r\_m^2 = r^2 \left(1 - \sqrt{|r^2 - r\_o^2|}\right) \tag{11}$$

$$r'^2\_{\
m} = r^2 \left(1 - \sqrt{\left|r^2 - r'^2\_o\right|}\right) \tag{12}$$

$$
\left\langle r\_m^2 \right\rangle = \left( r\_m^2 + {r'}\_m^2 \right) / 2 \tag{13}
$$

$$
\Delta r\_m^2 = \left| r\_m^2 - r\_m'^2 \right| \tag{14}
$$

$$\left(r^2 - r\_o^2\right) / r^2 < 0.10 \text{ and } 0.85 \le k \le 1.15. \tag{15}$$

To externally evaluate the predictivity of the generated models, the most stringent criteria validation values jointly proposed by Golbraikh et al. [82], Ojha et al. [96], Roy et al. [98], and Chirico and Gramatica [89] were adopted

$$r^2 \text{ } \ q\_{\text{CV}}^2 \text{ } q^2 \text{ } \ q\_{\text{Fn}}^2 \ge 0.70 \tag{16}$$

$$\left|r^2 - q\_{\rm CV}^2\right| < 0.10 \tag{17}$$

$$\left|r\_0^2 - r\_0'^2\right| < 0.30\tag{18}$$

$$r\_m^2 \ge 0.65\tag{19}$$

$$
\left \ge 0.65 \text{ and } \Delta r\_m^2 < 0.20 \tag{20}
$$

$$
\text{CCC} \ge 0.85\tag{21}
$$

where *r* 2 in Equations (15) and (18)-(20) symbolize *r* <sup>2</sup> and *q* 2 in the training set and external set, respectively. The *q* 2 Fn in Equation (16) stands for *q* 2 F1, *q* 2 F2, and *q* 2 F3.

#### **3. Results**

#### *3.1. Dataset Selection*

Of all the molecules enrolled in this study, 104 and 26 molecules were randomly selected as the training set and test set, respectively, giving rise to a ca. 4:1 ratio as suggested [81]. The chemical space with the projection of all molecules is displayed in Figure 1. Three principle components (PCs), which accounted for 97.94% of the variance in the original data, were used to create the chemical space. This figure shows that samples in the training set and test set had similar distribution in the chemical space. The high levels of the biological and chemical similarity between both datasets can be illustrated by the histograms of log *P*app, molecular weight (MW), surface area (SA), polar surface area (PSA), number of hydrogen bond acceptor (HBA), number of hydrogen bond donor (HBD), and *n*-octanol-water partition coefficient (log *P*) in the density form (Figure S1). Thus, it is plausible to assert that the substantial bias did not appear in the data partition.

It is of great significance to characterize the AD of the predictive model and to exclude the outliers from data collection [94]. Various methods to detect outliers have been proposed [99]. The scheme based on the chemical similarity/dissimilarity using principle component analysis (PCA) was adopted in this study [94]. Accordingly, 14 molecules were specified as outliers, which are substantially dissimilar to those ones in both the training and test sets, as shown in the chemical space (Figure 1), from which it can be observed that they are located far from the others. The distinction between the outliers and the others can be actually recognized by the fact that they contain more than nine rings or more than 12 HBAs as compared with the other molecules.

#### *3.2. SVR Models*

Numerous SVR models were generated using different descriptor combinations and runtime conditions. Three SVR models, coined as SVR A, SVR B, and SVR C, were assembled to establish the SVR ensemble, which was successively utilized to generate the HSVR model by another SVR. The optimal runtime conditions of SVR A, SVR B, SVR C, and HSVR are listed in Table S2.

SVR A, SVR B, and SVR C adopted five, five, and seven descriptors, respectively, with different combinations (Table 1). These SVR models in the ensemble were assembled

according to their performances on the molecules and statistical assessments in the training set and test set. Their runtime conditions and their predicted log *P*app values are listed in Tables S1 and S2, respectively. Tables 2 and 3 record their associated statistical evaluations in the training set and test set, respectively.

**Table 1.** The list of ensemble support vector regression (SVR) models and their descriptors, the correlation coefficient (*r*) with *P*app, and their descriptions.


† Selected. ‡ Not applicable.

**Table 2.** Statistic metrics including *r* 2 , ∆Max, mean absolute error (MEA), *s*, root mean square error (RMSE), *q* 2 CV, and *r* 2 s assessed by support vector regression (SVR) A, SVR B, SVR C, and hierarchical support vector regression (HSVR) in the training set.


**Table 3.** Statistic metrics including *q* 2 , *q* 2 F1, *q* 2 F2, *q* 2 F3 *CCC*, ∆Max, MAE, *s*, and RMSE assessed by SVR A, SVR B, SVR C, and HSVE in the test set.


The observed versus the predicted log *P*app values by SVR A, SVR B, SVR C, and HSVR are displayed by the scatter plot in Figure 2, from which it can be observed that SVR A, SVR B, and SVR C predicted the observed values well for the majority of the molecules in the training set, producing small MAE and *s* values consequently (Table 2). Moreover, it can be found from Figure 2 that the points predicted by SVR B are generally closer to the regression line than SVR A and SVR C. SVR B, consequently, gave rise to the lowest

Δ

∆Max (1.19), MAE (0.17), and RMSE (0.32), and the largest *r* 2 (0.77), suggesting that SVR B performed marginally better than SVR A and SVR C in the training set. Δ

**2021**, , x FOR PEER REVIEW 9 of 27

**Figure 2.** Observed log *P*app versus the log *P*app predicted by SVR A (gray circle), SVR B (gray triangle), SVR C (open diamond), and HSVR (solid square) for the training samples. The solid, dashed, and dotted lines represent to the HSVR regression of the data, 95% confidence intervals for the HSVR regression, and 95% confidence intervals for the prediction, respectively.

2 Furthermore, the difference between *r* <sup>2</sup> and *q* 2 CV evaluated by SVR B was 0.58 when subjected to the leave-one-out cross-validation, indicating that SVR B was over-trained which, in turn, can severely limit its application. Over-training was also associated with SVR A and SVR C as manifested by their extremely low *q* 2 CV values. The *r* 2 *s* values produced by SVR A, SVR B, and SVR C were 0.05, 0.03, and 0.03 (Table 2), respectively, when subjected in *Y*-scrambling. These near zero values suggest that there is an almost zero chance correlation associated with those SVR models [87].

The predicted values by SVR A, SVR B, and SVR C are in moderate agreement with the observed values for those test molecules depicted by Figure 3, which shows the scatter plot of observed versus the log *P*app predictions by SVR A, SVR B, SVR C, and HSVR for those samples in the test set. The MAE values generated by SVR A, SVR B, and SVR C increase from 0.28, 0.17, and 0.17 in the training set to 0.42, 0.35, and 0.39 in the test set, respectively (Table 3). RMSE along with the other statistic values also reveal deteriorating performances of these models in SVRE from the training set to the test set (Tables 2 and 3). Moreover, the *q* <sup>2</sup> values produced by SVR A, SVR B, and SVR C were 0.50, 0.58, and 0.60 in the test set, respectively, which are much less than their *r* 2 counterparts in the training set.

**2021**, , x FOR PEER REVIEW 10 of 27

**Figure 3.** Observed log *P*app versus the log *P*app predicted by SVR A (gray circle), SVR B (gray triangle), SVR C (open diamond), and HSVR (solid square) for the test samples. The solid, dashed, and dotted lines represent the HSVR regression of the data, 95% confidence intervals for the HSVR regression, and 95% confidence intervals for the prediction, respectively.

− − The prediction performances of those SVR models in the SVRE were significantly decreased when applied to those samples in the outlier set as suggested by the statistical metrics listed in Table 4. For example, SVR A, SVR B, and SVR C yielded the *q* 2 F2 values of −0.18, −0.41, and 0.16, respectively, which are substantially smaller than the *r* <sup>2</sup> values in the training set and the *q* 2 F2 values in the test set (Tables 2 and 3). Furthermore, the distances between the points and the regression line in the outlier set were much greater than those in the training set shown in Figure 4. As such, it can be asserted that those three models in the SVRE are vulnerable to the outliers that, actually, are not uncommon for most predictive models [100].


**Table 4.** Statistic metrics including *q* 2 , *q* 2 F1, *q* 2 F2, *q* 2 F3, *CCC*, ∆Max, MAE, *s*, and RMSE assessed by SVR A, SVR B, SVR C, and HSVE in the outlier set.

**2021**, , x FOR PEER REVIEW 11 of 27

− −

**Figure 4.** Observed log *P*app versus the log *P*app predicted by SVR A (gray circle), SVR B (gray triangle), SVR C (open diamond), and HSVR (solid square) for the outlier samples. The solid, dashed, and dotted lines represent the HSVR regression of the data, 95% confidence intervals for the HSVR regression, and 95% confidence intervals for the prediction, respectively.

#### *3.3. HSVR Model*

Δ

Δ The HSVR model was generated by the regression of SVRE according to the predictions of all molecules and statistical assessments in the training set (Table S1 and Table 2), and its runtime parameters are recorded in Table S2. HSVR commonly predicted better than SVR A, SVR B, and SVR C for the samples in the training set, as demonstrated by Figure 2, from which it can be noticed that most of predictions by HSVR lie in the range between the largest and the smallest ones predicted by those models in the SVRE. HSVR can improve the predictions in some cases. For instance, the prediction of compound **101** (omeprazole) by HSVR yielded an absolute residual of 0.02, whereas SVR A, SVR B, and SVR C produced the absolute errors of 0.34, 1.10, and 0.18, respectively (Table S1). In addition, HSVR produced the highest *r* 2 (0.91) and *q* 2 CV(0.81) and the lowest ∆Max (0.98), MAE (0.10), *s* (MAE), and RMSE (0.20) values when compared with those models in the SVRE, suggesting that HSVR statistically performed better SVR A, SVR B, and SVR C in the training set. Furthermore, HSVR gave rise to a *r* 2 *s* value of 0.03, indicating that it is least possible that HSVR was created by chance correlation [87].

When applied to the test molecules, marginal performance deteriorations can be found for HSVR. For example, *s* increased from 0.18 in the training set to 0.20 in the test set (Tables 2 and 3). However, ∆Max dropped from 0.98 in the training set to 0.72 in the test set. HSVR still executed better than SVR A, SVR B, and SVR C in the test set as shown in Figure 3. The other statistical parameters listed in Table 3 also assert the performance dominance of HSVR. For instance, the *q* <sup>2</sup> values were 0.50, 0.58, 0.60, and 0.75 generated by SVR A, SVR B, SVR C, and HSVR, respectively. Similarly, HSVR also produced smaller absolute deviations than its counterparts in the SVRE in the test set. For example, the absolute residuals of compound 36 (clozapine) were 0.35, 0.54, 0.35, and 0.03 yielded by SVR A, SVR B, SVR C, and HSVR, respectively (Table S1). HSVR generally produced consistent and small deviations in both training and test sets as asserted by those parameters listed in Tables 2 and 3 in comparison with its counterparts in the SVRE. More importantly, the HSVR model generated the largest *q* 2 (0.75) in the test set and the smallest difference between *r* <sup>2</sup> and *q* 2 CV (0.10), suggesting that it is less likely that HSVR model was over-trained or over-fitted.

HSVR even displayed better performance than the SVR models in the ensemble in the outlier set as depicted by those statistical assessments listed in Table 4. The HSVR model generated the largest *q* <sup>2</sup> value (0.76) and yet SVR A, SVR B, and SVR C yielded 0.45, 0.36, and 0.40, respectively. The superiority of HSVR in the outlier set can also be assured by the other statistical parameters, which is mainly due to the broader application domain of HSVR when compared with its counterparts in the ensemble. That robust HSVR feature makes it more utilizable in practical applications [101].

#### *3.4. Predictive Evaluations*

The scatter plot of residual versus the log *P*app prediction by HSVR for the training, test, and outlier samples is shown in Figure 5, from which it can be found that the residuals are commonly situated on both sides of *x*-axis along with the prediction range in those three datasets, suggesting that it is least likely that systematic error is associated with HSVR. Additionally, the training set, test set, and outlier set had the average residuals of 0.02, −0.13, and 0.06, respectively (Table S1), denoting that there is no biased prediction by HSVR. **2021**, , x FOR PEER REVIEW 13 of 27

**Figure 5.** Residual versus the log *P*app prediction by HSVR in the training set (solid circle), test set (gray square), and outlier set (open triangle).

2

2

Table 5 lists the results when the developed HSVR model was further subjected to the most stringent validation criteria collectively recommended by Golbraikh et al. [82], Ojha et al. [96], Roy et al. [98], and Chirico and Gramatica [89] in the three datasets (Equations (15)–(21)). It can be observed that HSVR completely met those proposed validation requirements in addition to the fact that HSVR exhibited a similar degrees of performance in the training set, test set, and outlier set. As such, it can be asserted that HSVR is an extremely accurate and predictive theoretical model.

**Table 5.** Validation verification of HSVR based on prediction performance of the training, test, and outlier samples.


† Fulfilled; ‡ Not applicable.

#### *3.5. Mock test*

To verify the practical applicability of the generated HSVR model, this model was applied to those drugs measured by Yamashita et al. [30]. There were eight compounds commonly adopted by this study and Yamashita et al., furnishing a sound way to calibrate the challenging system. However, Yamashita et al. assayed the *P*app values at pH 6.0, instead of pH 7.4 used by those compounds collected in this study, suggesting that some *P*app variations can be resulted from both systems (vide supra). These discrepancies make those drugs assayed by Yamashita et al. not appropriate as the second external dataset or the test set because those validation criteria listed in Table 5 cannot be applied to those drugs. The relationship between both different experimental conditions was initially constructed for those eight common compounds, and the resulting scatter plot is exhibited in Figure 6, from which it can be found that both assay systems were reasonably correlated with each other with an *r* value of 0.86), suggesting that this HSVR can be adopted to predict those novel compounds measured by Yamashita et al.

Figure 7 shows the predicted results of seven novel drugs in the mock test. The correlation coefficient *r* value between the predicted log *P*app (pH 7.4) and observed log *P*app (pH 6.0) was 0.86, suggesting that the HSVR model can nearly reproduce the experimental results. In addition, the produced *p*-value was <0.05. This mock test ensured the predictive ability of generated HSVR when applied to the novel compounds with different experimental conditions.

**2021**, , x FOR PEER REVIEW 14 of 27

**2021**, , x FOR PEER REVIEW 14 of 27

**2021**, , x FOR PEER REVIEW 15 of 27 **Figure 6.** Observed log *P*app at pH 7.4 versus observed log *P*app at pH 6.0 for the common drugs in the mock test. The solid, dashed, and dotted lines represent the mock test regression of the observed data, 95% confidence interval for the mock test regression, and 95% confidence interval for the observation, respectively. **2021**, , x FOR PEER REVIEW 15 of 27

s in the mock test. The cor-

‒

‒

−

−

−

−

**Figure 7.** Predicted log *P*app at pH 7.4 versus observed log *P*app at pH 6.0 by the HSVR model for the drugs in the mock test. The solid, dashed, and dotted lines represent the HSVR regression data, 95% confidence interval for the HSVR regression, and 95% confidence interval for prediction, respectively.

#### *3.6. Classification*

It is of interest to verify the qualitative predictivity of HSVR, since a number of qualitative models have been published [25,102]. Accordingly, compounds enlisted in this

‒

‒

*κ* −

*κ* −

study were classified as Caco-2 permeable (Caco-2<sup>+</sup> ) and Caco-2 impermeable (Caco-2- ) based on the threshold value of *P*app (8 × 10 -6 cm/s) as suggested [25,102]. Initially, the confusion matrix was constructed (Table S3), and the Cooper statistics and Kubat's Gmean [103] (Table S4) were employed to qualitatively evaluate the predictivity of HSVR. The results were also compared with predictions made by *admetSAR* [104] (available at website: http://lmmd.ecust.edu.cn/admetsar2/), since *admetSAR* has been adopted by DrugBank (available at: https://go.drugbank.com/) to qualitatively predict Caco-2 permeability. The results are listed in Table 6, from which it can be asserted that HSVR outperformed *admetSAR* in every aspect. For instance, the parameter accuracy was 93.1% produced by HSVR, which is substantially higher than that generated by *admetSAR* (50.7%). The metric MCC is the most distinction between HSVR and *admetSAR* (85.0% vs. −8.0%). Thus, it can be asserted that HSVR is also an accurate and predictive qualitative predictive model.


**Table 6.** Statistical parameters of qualitative predictions by HSVR and *admetSAR*.

#### **4. Discussion**

Caco-2 has been commonly adopted to predict the intestinal permeability in the process of drug discovery because of its morphological and functional similarity with human enterocytes [105]. The mechanism of Caco-2 permeation is rather complex, since it can take place through passive diffusion, which can go through the paracellular and transcellular routes and active transport. The passive diffusion is predominately governed by the concentration gradient, and most hydrophilic drugs prefer to penetrate between cells in a paracellular fashion, whereas hydrophobic drugs are inclined to get across the cells via the transcellular route. Drugs that can permeate the Caco-2 cells by the active transport can interact with the influx and/or efflux transporters expressed on the cell surface [106]. As such, Caco-2 permeability is affected by some physicochemical and physiological properties [106].

Hydrophobicity or lipophilicity plays an important role in passive diffusion through membranes as well as the drug–receptor interactions [17,107,108]. In addition, hydrophobicity, which can represent by the *n*-octanol-water partition coefficient, *viz*. log *P*, is also an important factor affecting the interaction between the molecules and the target protein, since more lipophilic molecules tend to have stronger interactions with both target protein and biological membrane. Therefore, the very lipophilic molecules have poor oral absorption from the stomach [107,109]. Polar and hydrophobic drug must penetrate through the Caco-2 cell membrane [17,110]. In addition, it has been observed that log *P*, hydrogen bond propensity, weight, and volume are closely related with *P*app [43]. As such, log *P* was adopted in this study (Table 1), which is consistent with the fact that numerous published in silico models to predict intestinal absorption, PAMPA permeability [1,111], and Caco-2 permeability also have employed this descriptor [40,112–114]. It can be observed from Figure 8, which displays the average log *P*app for each histogram bin of log *P* for all molecules included in this investigation, that log *P*app increased with log *P* value initially and then decreased afterward, leading to a seemingly bilinear relationship between log *P*app and log *P*. This perplexing dependency can be realized by the fact that the more hydrophobic solutes can easier approach the lipid bilayer to penetrate the membrane. The

opposite relationship between hydrophobicity and permeability will be resulted when the solutes are too hydrophobic due to stronger attractions between solutes and the membrane as well as stronger repulsive forces from the solvent molecules upon the entrance to the solvent environment that can be illustrated by the PAMPA permeability [1,115,116]. Complexity can be even profound when taking into account the fact that P-gp and BCRP, which are efflux transporters in Caco-2 (vide supra), can interact with substrates by hydrophobicity [117], subsequently leading to a low correlation between log *P*app and log *P* (*r* = 0.15). **2021**, , x FOR PEER REVIEW 17 of 27

**Figure 8.** Histogram of average log *P*app versus the distribution of log *P*.

‒ It has been observed that the number of aromatic rings (*n*Ar) has a positive correlation with log *P* with an *r* value of 0.67 [118], suggesting that a predictive model can be overtrained once both log *P* and *n*Ar are adopted simultaneously. However, this issue was not concerned in this study, since only SVR C adopted this descriptor, whereas SVR A and SVR B included log *P* (Table 1). In addition, the aromatic ring is a non-polar group, which can enhance the hydrophobicity [52] and increase the passive diffusion [119,120]. In addition, aromatic ring moieties have been implicated in P-gp substrate recognition and efflux modulation [53], leading to the fact that *n*Ar can be an important factor in Pgp modulation action [121] and BCRP-substrate interactions [122]. As such, *n*Ar plays a complex role in both passive diffusion and active transport in Caco-2 permeability.

‒ ‒ It has been recognized that both PSA and *µ* are associated with passive diffusion [37,123–125]. In addition, these descriptors have been adopted by published in silico Caco-2 permeability models [37,45–49,126–128]. It has been reported in the PAMPA permeability study that larger PSA, *µ*, and polarity can enhance the solute-solute and solute-solvent interactions, which, in turn, require more desolvation energy when the solutes penetrate through the lipophilic membrane to the donor compartment [123,129–132], and conversely decrease the passive diffusion [1], consequently, making permeability less favorable. Therefore, it has been shown that PSA has a negative impact in the permeation rate [133,134]. In addition, Joung et al. have indicated that PSA shows an important role in distinguishing the P-gp substrate from the non-substrates [135]. Accordingly, PSA and *µ* were adopted in this study due to their pivotal roles in Caco-2 permeability.

It is seemingly unusual to include the descriptor |*µ*|max, which is the absolute maximum component of the molecular dipole, in this study, since it has never been employed by any published model before. This inconsistency actually can be manifested by Figure 9, which displays the average |*µ*|max for each histogram bin of *µ*, that the larger *µ*, the larger |*µ*|max, suggesting that they were positively correlated with each other. In addition, *µ* was recruited by SVR A and SVR C, whereas |*µ*|max was enlisted by SVR B only, suggesting that it is less likely to produce an over-trained HSVR, since no single model adopted both two descriptors simultaneously. More importantly, the empirical observation has revealed that HSVR including these selections executed better than the others (data not shown) plausibly because of the descriptor-descriptor interaction [1]. Any other traditional linear or machine learning-based QSAR schemes, conversely, cannot properly render such contradictory descriptor selections. **2021**, , x FOR PEER REVIEW 18 of 27 ‒

**Figure 9.** Histogram of average |*µ*|max versus the distribution of *µ*.

‒ α α α α It has been reported that the molecular size of the solute molecule is of critical importance in the diffusivity of the biological membrane [37,125,136], and the intestinal absorption can decrease with the increase of molecular size [137]. Furthermore, the molecular size also affects passive diffusion through membranes [138,139] and active transport through the P-gp-substrate interactions [121,138]. Molecular size can be represented by a number of descriptors such α, *n*Ring, *V*m, and *n*rot [140–142], which were adopted in the investigation and negatively associated with log *P*app (Table 1). Conversely, Fujiwara et al. adopted the descriptor molecule weight (MW) to develop a theoretical Caco-2 permeability model [37], whereas MW was not included in this study. This discrepancy can be realized by the fact that α was highly correlated with MW with an *r* value of 0.98 for all molecules enlisted in this study, suggesting that it is plausible to replace MW by α in order not to produce an over-trained model. In addition, it has been observed that α is positively correlated to log *P* [143] and is highly associated with absorption [50].

α α The descriptor *n*Ring, which is reportedly related to molecular size [136,141], has never been adopted by any published Caco-2 permeability predictive model and yet was selected by SVR C (Table 1). This disagreement can be recognized by the fact that *n*Ring was greatly correlated with α with an *r* value of 0.78 for all molecules recruited in this study. As such, it is plausible to expect that both *n*Ring and α play similar roles in Caco-2 permeability. The relationship among log *P*app, *n*Ring, and log *P* can be further perplexing as illustrated by Figure 10, which shows the 3D plot of log *P*app, *n*Ring, and log *P*. The relationship between *n*Ring and log *P* has been detailed by Pham-The et al. [125].

**2021**, , x FOR PEER REVIEW 19 of 27

**Figure 10.** The relationship among log *P*app, *n*Ring, and log *P* in 3D presentation.

It has been observed that *V*<sup>m</sup> plays an important role in passive absorption [9,144,145] and it is adopted by a published Caco-2 permeability model [146] as well as in this study. It has been observed in the rat that fewer rotatable bonds, viz. smaller *n*rot, can lead to better oral bioavailability, and *n*rot can also exert a positive effect on the permeation rate [133,143], since more rigid molecules will have smaller *n*rot values that, in turn, can enhance permeability [125]. Furthermore, *n*rot is of importance in intestinal absorption [147], since increased *n*rot can reduce the permeability [133]. Furthermore, a number of published membrane permeability models have also employed the descriptor CSA, which is another feature associated with molecular size and also plays a pivotal role in membrane permeability [70,71]. However, *n*rot was greatly associated with CSA with an *r* value of 0.80 for all molecules enrolled in this investigation, suggesting that using *n*rot in lieu of CSA without producing the over-trained model is plausible. Li et al. also have found that *n*rot is another feature to discriminate P-gp substrates from non-substrates [148]. As such, it is of necessity to recruit *n*rot in model development to properly render Caco-2 permeability as suggested [71,72].

Hydrogen bonding potential, which can be expressed by HBD and HBA, is another important factor in determining the solute–solvent interactions [37], and it is the main contributor for the passive diffusion [143]. It has been observed that Caco-2 permeability is a function of HBD and/or HBA, since more permeable solutes tend to have smaller HBD and/or HBA [130,131,149]. Between HBD and HBA, HBD seemingly shows a more profound effect on Caco-2 permeability as compared with HBA [150] as manifested by the fact that several published in silico models have selected HBD to predict Caco-2 permeability instead of HBA [35,42]. Mechanistically, HBD is one of the features associated

‒

with P-gp-substrate interactions [148,151]. In addition to efflux transport, HDB is one of the features linked to substrate binding with OATP2B1 [7] as well as PepT1 [152]. Thus, it is of necessity to include in Caco-2 predictive models to take into consideration the passive diffusion as well as the active influx/efflux transport.

**2021**, , x FOR PEER REVIEW 20 of 27

The descriptor p*K*a(Max) was selected in this study due to the fact that higher p*K*a(Max) can lead to the lower ionized form of drugs in the donor compartment, which, in turn, can increase the penetration through hydrophobic membrane [153]. Furthermore, it has been recognized that neutral compounds can have higher membrane permeability than the other ion classes [154]. Accordingly, all molecules included in this investigation were categorized into different ion classes based on their p*K*<sup>a</sup> values. In addition, ABC and/or SLC substrates were also identified based on the drug information retrieved from Drug-Bank to understand if the dependence of ion class can be varied by their ion classes. It can be found from Figure 11, which displays the histograms of median log *P*app versus all molecules, ABC substrates, SLC substrates, as well as ABC and SLC substrates for four different ion classes, that the median log *P*app values of neutral compounds are substantially larger than the others, suggesting that neutral compounds exhibit higher Caco-2 permeability regardless of active transporter substrate classes, viz. influx transporter or efflux transporter. This observation actually is very similar to the PAMPA permeability, since the ionized compounds will demand larger desolvation energies, which, in turn, can hinder their penetration [134].

**Figure 11.** The histogram represents the log *P*app versus the molecules belong to ATP-binding cassette (ABC) substrate, solute carrier (SLC) substrate, both ABC and SLC substrate, and total drugs in the acid class, base class, neutral class, and zwitterion class, respectively.

Initially, numerous efforts were made in attempting to build assorted 2-QSAR models by employing the partial least square (PLS) scheme, and yet no productive models were produced (data not shown) [1]. This challenge can be realized by the fact that the correlations between the designated descriptors and log *P*app for all molecules included in this investigation were small, and the largest absolute maximum *r* was only 0.56 between PSA and log *P*app (Table 1), signifying the high non-linearity between them. More significantly, the substantial difference in 2-QSAR development between the passive diffusion, viz. the PAMPA system, and Caco-2 permeability can be greatly attributed to the complex active (influx and efflux) transport. Thus, it is extremely difficult, if not absolutely impossible, to derive a linear Cacao-2 permeability QSAR model. Conversely, the accurate and predictive HSVR model can properly render such non-linear dependence of log *P*app on descriptors.

#### **5. Conclusions**

Intestinal permeability is one of the important ADME/Tox metrics that should be addressed in the process of drug discovery and development. The Caco-2 system has been frequently used as a surrogate to preliminarily investigate the intestinal absorption. An in silico model can be a useful approach to predict Caco-2 permeability in assisting drug discovery and development. However, Caco-2 permeability can occur through passive diffusion and active transport, leading to a complex process. Therefore, it is of necessity to include different descriptor combinations and diverse relationships to address these variations in distinct mechanisms. The innovative machine learning-based HSVR scheme, which possesses the superior features of a local model (greater predictivity) and a global model (larger coverage of the application domain), was employed in this study to construct a theoretical model to predict the Caco-2 permeability. The generated HSVR models unveiled great prediction accuracy for the training, test, and outlier samples. When challenged by a group of drugs assayed at different experimental conditions, the developed HSVR model also executed equivalently well. In addition, HSVR showed excellent qualitative performance in recognizing Caco-2 permeable and impermeable compounds, and the selected descriptors can completely justify the diverse mechanisms related to the passive diffusion and active transport. Thus, it can be assured that this HSVR model can be useful to accurately and swiftly predict the Caco-2 permeability of novel compounds in order to assist drug discovery and development.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/1999-4 923/13/2/174/s1, Table S1. Selected compounds for this study; their names, SMILES strings, CAS numbers, and observed log *P*app values; their predicted values by SVR A, SVR B, SVR C, and HSVR; data partitions; and references; Table S2. Optimal runtime parameters for the SVR models; Table S3. Confusion matrix for the qualitative predictive model; Table S4. The Cooper statistics and Kubat's G-mean calculated from the confusion matrix. Optimal runtime parameters for the SVR models; Figure S1. Histograms of: (A) log *P*app, (B) molecular weight (MW), (C) surface area (SA), (D) polar surface area (PSA), (E) number of hydrogen bond acceptor (HBA), (F) number of hydrogen bond donor (HBD), (G) and the *n*-octanol-water partition coefficient (log *P*) in the training set, test set, and outlier set.

**Author Contributions:** G.H.T., C.-F.W., and M.K.L. conceived and designed the study; G.H.T., C.-S.J., and M.K.L. performed the experiments and analyzed the data; and G.H.T., C.-F.W., and M.K.L. wrote the paper. The final version of manuscript is reviewed and approved by all authors. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Ministry of Science and Technology, Taiwan.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** Parts of calculations were performed at the National Center for High-Performance Computing, Taiwan. The authors are grateful to Paola Gramatica for providing the free license of *QSARINS*.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


#### *Article*
