Prediction of Passive Membrane Permeability by Semi-Empirical Method Considering Viscous and Inertial Resistances and Different Rates of Conformational Change and Diffusion.

Membrane permeability is an important property of drugs in adsorption. Many prediction methods work well for small molecules, but the prediction of middle-molecule permeability is still difficult. In the present study, we modified a classical permeability model based on Fick's law to study passive membrane permeability. The model consisted of the distribution of solute from water to membrane and the diffusion of solute in each solvent. The diffusion coefficient is the inverse of the resistance, and we examined the inertial resistance in addition to the viscous resistance, the latter of which has been widely used in permeability prediction. Also, we examined three models changing the balance between the diffusion of solute in membrane and the conformational change of solute. The inertial resistance improved the prediction results in addition to the viscous resistance. The models worked well not only for small molecules but also for middle molecules, whose structures have more conformational freedom.

Introduction

Permeability is one of the most important factors in a drug's adsorption and target‐binding properties in cells. The understanding and predicting membrane permeability of molecules have been studied for last few decades. It is still one of the hot topics, especially under circumstances where the molecular weights of drug molecules have been increasing and larger molecules often face the lower permeability than smaller drug molecules do. There have been a number of reports on permeability.1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26 The main permeability problems are adsorption in human intestine, extraction from kidney, penetration of the blood‐brain barrier, skin permeability, and the permeability of the cell membrane to approach target proteins in cells. Caco‐2 cells and MDCK cell systems are two of the model systems that mimic human intestine adsorption and extraction from kidney, respectively. Parallel artificial membrane permeability assay (PAMPA) systems have been popular in vitro assay systems for the past 20 years.26, 27, 28, 29, 30, 31, 32 PAMPA systems have been improved to mimic in vivo permeability by trying various membrane materials, pH levels of donor and acceptor liquids and the other conditions. Certain mechanisms underlie permeability.2, 3, 26 Namely, solute molecules penetrate the cell membrane by diffusion (transcellular), the solute molecules go through the tight junction (paracellular), and transporters and channel proteins work in the influx and efflux processes. Among these mechanisms, PAMPA permeability represents transcellular passive permeability.

Recent advances in molecular dynamics (MD) simulations enable the understanding and evaluation of transcellular passive permeability.4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15 In these calculations, MD simulations have been applied to explicit atomic models of membrane molecules with solvent water molecules. Since permeation is a very slow process, biased MD simulations have been popular in this analysis. MD simulations have shown that the distribution of the existence probability of solute and the diffusion of solute in the membrane gave the permeability constant.

On the other hand, many approaches have adopted quantitative structure‐property relationship (QSPR) models for passive permeability based on Fick's law.2, 3, 26 Previous works have shown the efficiency of this basic theory, and some extensions from this theory have improved prediction accuracy.16, 17 Fick's law explains that permeability is a combination of the transfer of solutes from the donor water into the membrane and the diffusion of the solutes from the donor to the acceptor sides in the membrane. Only the neutral molecule moves into the lipid layer so that the pK a and the partition coefficient (LogP) or the distribution coefficient (LogD) determine the distribution of solute between the donor water and the membrane. The permeability P app is given as follows.

where P app, D, M, and h represent the apparent membrane permeability, distribution coefficient, diffusion constant of the solute, and thickness of the membrane, respectively. The above MD simulations support this assumption. Namely, LogD corresponds to the probability distribution of existence and M corresponds to the diffusion of solute obtained by the MD simulation, respectively.

The diffusion constant M is estimated by the Einstein‐Stokes relation

where k B, T, μ, and R are the Boltzmann constant, temperature, viscosity and radius of the solute, respectively. Since LogD corresponds to the free energy of the transfer from water to membrane, this value could be approximated by the accessible surface area (ASA) method, the generalized Born (GB) method, polar surface area (PSA), and so on.33 Thus, the following linear regression model can estimate LogP app.

where x, c, and N descriptor are the molecular descriptors, the fitting parameters, and the total number of descriptors, respectively, and the fitting parameters represent the characteristic properties of the membrane.

Recent advances in pharmaceutical research have increased the molecular size of drugs, and middle‐molecule drugs, with molecular masses >500 Da, have become popular. In the last few decades, the major drug targets have been receptors and enzymes. Pharmaceutical companies have released several thousands of drug molecules; however, they are now facing a severe depletion in the druggable targets with conventional strategies. Research activities focusing on protein‐protein interaction (PPI) inhibitors have been started instead of the projects on finding ligands to receptors and enzymes.34, 35 To inhibit PPIs, larger molecules are often preferable than the small molecules. However, the hydrophobic and often insoluble physical properties of middle molecules that are distinct from small ligands cause serious problems in their development stages.

While one of the major problems of middle molecules is the low membrane permeability, the synthetic accessibility in chemical modifications is also limited. The synthetic processes of middle molecules are more complicated and time‐consuming than that of small molecules. Thus, we need the mechanistic understandings to unveil the dominant factors of membrane permeability, rather than black‐box permeability predictions, to guide a rational modification of the middle molecules.

A number of reports suggested that conformational change is essential in permeability and cyclization as well as in the methylation of main‐chain amide groups, and some chemical modifications have improved permeability.18, 19, 20, 21, 22 Many studies have challenged this problem by using QSPR models and have successfully predicted membrane permeability, including that of middle molecules.23 However, some studies have suggested that it is still difficult to predict the membrane permeabilities of middle molecules.

Diffusion constant M in eq. 1 consists of viscous resistance and inertial resistance. Equation 2, which is the inverse of viscous resistance, assumes the slow migration of solute in equilibrium. Permeability is non‐equilibrium process in general, and the density of the solute keeps changing in the experiment. When solute molecules push solvent molecules, the inertial resistance is proportional to the cross section of the solute. In general, we have ignored inertial resistance without examining the actual experimental data. In a non‐equilibrium state, the density of solute in the membrane cannot reach equilibrium. The solute can pass through the membrane before the conformational ensemble of the solute reaches equilibrium distribution (see Figure 1A).4 In this case, the distribution of solute is different from the distribution coefficient D observed in equilibrium.

Figure 1

Permeability mechanism Models A and B. The rings represent the solute molecules and the different shapes represent the different conformers. Thick solid and thin dotted rings represent the large and small populations of molecules, respectively. The buffer solvent region is separated by the lipid bilayer in the middle (gray). Orange and green arrows represent the diffusion in the membrane and conformer change in the membrane depicted in gray. The open forms (3 and 4) are major rather than the closed form (1 and 2) in the buffer solvent. In Model A, there is no conformation change in the membrane. In Model B, the solute molecules in membrane show conformation change and the closed form becomes major.

In the present study, we examined the diffusion process by using a classical QSPR model. We considered inertial resistance the same as viscous resistance. Also, we developed formulas considering the above two cases in which the conformational change occurs faster or slower than the diffusion in the membrane.5

Methods

We classify the membrane permeability process into two models according to the speed of permeability.5 If the membrane permeation is faster than the conformational change of the solute in the membrane, each conformer penetrates the membrane while keeping the 3D structure of each conformer, and the dominant conformer in water could be the main contributor to the permeability (permeability model A in Figure 1A). If the membrane permeation is slower than the conformational change of the solute in the membrane, the most stable conformer in the membrane could contribute to the permeability (permeability model B in Figure 1B). In this article, we apply Fick's law to membrane permeability, since the physical meaning of Fick's law is clear and simple.

Permeability Model A

Since each solute penetrates the membrane without conformational change, the distribution of conformers in the membrane is equal to that in water solvent. The total permeability ratio (P app all) is the summation of the permeability of all conformers.

Here, P app(a), d wat(a), and N conformer are the apparent permeability of the a‐th conformer, the distribution of the a‐th conformer in water, and the total number of conformers, respectively. The following relation is obvious.

where the brackets < >wat stand for the average over the distribution in water. The fraction of the a‐th conformer in water (=d wat(a)) is given by

where E wat(a) and n(a) are the energy and the degeneracy of the a‐th conformer, respectively.

When we apply Kubo's cumulant expansion to eq. 5, the first two terms of the expansion are as follows.36

The first term is the average LogP app and the second is the deviation of LogP app.

A linear regression model that is a weighted summation of the descriptor values approximates the LogP app of the a‐th conformer as follows,

where c and x(a) are the fitting parameters and the descriptor of the a‐th conformer, respectively.

If the descriptors x A and x B are independent from each other, =+ and also σ(x A+x B)=σ (x A)+σ (x B), where σ stands for the deviation. Thus, if all the descriptors are independent of each other, eq. 7 becomes as follows.

In permeability model A, the deviations of all the descriptors, including the ASA and the radius of the solute, contribute to LogP app all in addition to the average values of these descriptors.

Permeability Model B

If the membrane penetration is much slower than the conformation change, the distribution of conformers reaches equilibrium in the membrane and in water. The distribution constant D is given by the partition function of the molecule as follows.

where Z°ct, Z wat, and n(b) are the partition functions in octanol and in water and the degradation number of the b‐th conformer. E°ct and E wat are the energy of the conformer in octanol and water, respectively.

The D value gives the density of molecules on the donor−‐membrane interface of the membrane. The summation of diffusions of all the conformers gives the total diffusion.

The fraction of the a‐th conformer in the membrane is given by

As in eq. 7, here we apply Kubo's cumulant expansion to eq. 11, giving

In permeability model B, the deviations of the descriptors relating to the diffusion contribute to LogP app all besides the average values of the descriptors. Namely, the deviation of the radius of the solute should contribute to the prediction of LogP app.

Regression and Prediction

Our physical‐property prediction method was a principal component regression (PCR) with an L2 regularization term based on the molecular descriptors.37, 38 The regression model was the same as that used in our previous study. Namely, the principal component (PC) analysis projected each compound into each point in a chemical space of the PC, and a multiple linear regression was applied to the molecular coordinates in the chemical space. The principal component axes were selected to minimize the root‐mean‐square error (RMSE) in regression. The addition of the L2 term reduced the RMSE in the prediction.

The physical property of the i‐th molecule V(i) is estimated based on the molecular descriptors {s i } where b represents the b‐th descriptor as follows.

where c, p i , and d b are the intercept parameter of the linear function, the PC vector, and the loading vector, respectively. The PC analysis of the descriptor s gives the loading vector d b and the PC vector (axis) p i . The values of {c} are determined by multiple linear regression (MLR). N axis (N axis

Here, the parameter set {c} is determined to minimize the O value. The weighting parameter λ is fixed to a constant (λ=0.000001) in the present study.

We apply this regression model to P app, LogP, LogD, and pK a. Since the permeability depends on the LogP, LogD, and pK a, the P app prediction model should include the descriptors for LogD, LogP, and pK a prediction. All the prediction models shared the same molecular descriptors. We applied 4‐fold CV tests to evaluate the accuracy, and all the values in the tables were predicted in the 4‐fold CVs. The CV tests showed RMSE and Q2 values. The definitions of Q2 and RMSE are determined as follows.

Here, Y pred and Y exp represent the predicted value in validation and the experimental Y value, respectively. In the present study, we do not compare the Q2 values of different properties, since the variances of the experimental data (denominator in eq. 17) differed from each other among the different data sets.39

Models A and B use the degeneracy and energy of the conformers in eqs. 6 and 9 so that we must sample the conformers. Although an exhaustive search of conformers is very difficult and time‐consuming, many heuristic conformer generation methods have been proposed. The major approaches are random search and systematic search.40, 41, 42, 43, 44, 45, 46, 47 Especially, the conformer generation of ring systems is more difficult than that of nonring systems. For ring systems, some specialized methods, such as corner and edge flips, have been used.40, 41 We applied a Monte Carlo random‐sampling method.42 In conformer generation, a difficult problem is the estimation of degeneracy (n(a)) in eqs. 6 and 9. Most conformer generation programs focus on the reproduction of the most stable conformers precisely, but the degeneracy is unclear. Long‐time MD simulations could measure the degeneracy but are not realistic in this work, since there are too many molecules.

One of the simplest approaches should be a uniform sampling of conformers without solvent bias. If the conformer search was a random sampling of dihedral angles of rotatable bonds of a molecule, the generated structural ensemble should mimic the uniform sampling (see some discussions in APPENDIX A in the supporting information). Unfortunately, this approach could not work in the sampling of large or cyclic compounds, since the ring systems could be opened by the random rotations of the dihedral angles and an atomic collision could occur. Thus, we applied a force field to close the ring systems and to avoid atomic collisions. Our conformer generation program transforms the ring systems of a molecule to nonring systems by removing one of the bonds in the ring. The following energy optimization closes the rings. The force field used was an AMBER‐like one including 1–2, 1–3, 1–4, and 1–5 interactions. The 1–4 and 1–5 interaction terms consisted of only van der Waals interactions without electrostatic interactions, since the dielectric constant of water is very different from that of membrane.

A clustering analysis of the conformers generated above estimates the degeneracy number (n(a) in eqs. 6 and 9). The clustering threshold is the heavy atom RMSD<1.5 Å.

Descriptors

The molecular descriptors consisted of physical and chemical ones. The physical descriptors represent mainly the size of a molecule that is related to the diffusion. The accessible surface area (ASA) is a useful idea with which to represent the solute−solvent interaction. The chemical properties represent mainly the hydrophobicity/hydrophilicity of a molecule that is related to the distribution between water solvent and membrane. In general, charged molecules cannot penetrate a membrane whose dielectric constant is small, but neutral molecules can. The total charge of a molecule is determined by the pK a and the MACCS key,48 which is a set of substructures that recognizes the functional groups closely related to pK a.48 The atomic charges were AM1‐BCC charges obtained by MOPAC7.49, 50, 51 The hydrogen bonds are determined for the hydrogen atoms of OH and NH groups and for the acceptor atoms with lone pairs (O, N, S). The GBSA method was used to calculate the LogD by eq. 10.33 Since the descriptors include the ASA, the ASA and the LogD term double‐counted the solute‐solvent interaction. The PCR could combine these dependent terms and thus reduce the number of independent fitting parameters. The atomic solvation parameters in the GBSA method were set to 10 cal/mol/Å2 and −5 cal/mol/Å2 for water and membrane, respectively.

Models A (eq. 9) and B (eq. 13) include both the average and deviation terms. The average and deviation σ(A) of property A are determined as follows. The conformer generation program in section 2.4 generates the conformers.

where sol represents the solvent (water or membrane). E s°l(i) is the relative energy of the i‐th conformer in the solvent compared to that in vacuum.

Since membrane permeability is related to LogD, or to LogP and pK a values as described in eqs. 1, 3, 8, and 13, the descriptors that could approximate LogD, LogP, and pK a values should contribute to P app values. Dissociations of hydrogen atoms depend on the chemical bonds of the functional groups and the electrostatic field generated by the charge distribution of the solute. A previous work proposed a linear correlation between the proton chemical shift and the pK a value, and the chemical shift depends on the electron density on the nucleus.52 Thus, we adopted the descriptors used in the LogD and pK a predictions. Namely, we used the numbers of hydrogen donors and acceptors, the atomic charges of the hydrogen atoms of the NH and OH groups, and the MACCS key to represent the chemical structures of the solutes.

Data Preparation

The PAMPA permeability data and molecular structures were extracted from ChEMBL database version 24.53 The ChEMBL assay and compound IDs used are summarized in Table S1. Most of the experimental conditions included pH 7.4 and the observation times were several hours, but the details of the conditions varied. In addition, a variety of membrane materials have been used.15, 26, 27, 28, 29, 30, 31, 32 A desirable assay data set consists of the P app values obtained under the unique experimental conditions.18 The careful data selection limited the number and diversity of examined compounds. While there have been differences in the P app values observed through the PAMPA assay, the systems were designed to reproduce P app values of the Caco‐2 assay system. To evaluate the wide variety of compounds, the present data set consisted of 737 compound data obtained under the heterogeneous experimental conditions.; There were 70 pairs of the same solute with different P app values; for these 70 pairs, the average difference was 0.401 and the RMSD was 0.500. Thus, we discussed the prediction accuracy of LogP app >0.5. As with the P app, we also prepared data on the octanol‐water distribution coefficient (LogD, 4215 compounds),54 and the dissociation coefficient (pK a, 240 compounds).55, 56

Astellas Pharmaceutical kindly provided additional experimental data, since there were fewer data on macrocyclic compounds than on acyclic compounds. The observed experimental property data were P app, LogD, and LogS data. The test compounds were 58 selected commercially available macrocyclic compounds. In addition, we asked Enamine Ltd. to evaluate the P app values of two cyclic peptides using the PAMPA assay. The molecular structures and experimental values are summarized in Tables S1 and S2 in the supporting information.57

We prepared three‐dimensional molecular structures of the examined compounds in the present study. The computational procedure was summarized as APPENDIX B in the supporting information. The initial molecular structures were 2D planar structures without hydrogen atoms in the SD file format. A file transfer program, Hgene/myPresto, was used to prepare the three‐dimensional molecular structures with hydrogen atoms in electrically neutral forms including acids and amines. The atomic charges were the AM1‐BCC charges provided by the MOPAC program. A molecular dynamics simulation program, cosgene/myPresto, gave the energy‐optimized structures with the general AMBER force field (GAFF) assigned by the tplgeneL/myPresto topology generator.58, 59 The energy optimization gave one stable or meta‐stable molecular structure for each molecule in vacuum. All the molecular structures were prepared in the Sybyl mol2 file format.

The total data set consisted of 795 data points from the ChEMBL database and the unpublished assay data. Table 1 shows the distributions of molecular weight and the total number of atoms as well as the numbers of heavy atoms, ring structures, rotatable bonds, and atoms included in the biggest ring system of each molecule.

Table 1

Statistical features of the data set for permeability prediction.

	Average	Deviation
Molecular weight (Da)	424.1	160.4
No. of all atoms	54.5	23.1
No. of heavy atoms	29.6	11.4
No. of rings	3.2	1.3

Statistical features of the data set for permeability prediction.

Average

Deviation

Molecular weight (Da)

424.1

160.4

No. of all atoms

54.5

23.1

No. of heavy atoms

29.6

11.4

No. of rings

3.2

1.3

Results and Discussion

Prediction Models for Models A and B

Our LogP app prediction model represents the dissociation of solute in water, the distribution of solute from donor water (donor) to membrane, and the diffusion of solute in the membrane towards the water (acceptor) phase. Diffusion depends on the solute molecular radius R, and many previous works adopted R as the descriptor in LogP app prediction.2, 16, 23, 24 The ES term (eq. 2) represents the viscous resistance, while the real solute molecule feels both the viscous resistance and the inertial resistance. The inertial resistance is proportional to the product of the cross section of the solute and the square of the velocity of the molecule. The Taylor series of the total diffusion becomes as follows.

where R, f v, and f are the average radius of the solute and the respective coefficients of the diffusions related to the viscous and inertial resistance. In the actual calculation, we used the average R (denoted as ) instead of R. Thus, we described the diffusion term of LogP app as a weighted summation of , 2, 3, and Log. Also, the 3 term represented the permeant's volume, since the transfer of the molecule from the donor phase to the membrane needed a volume change of the membrane, and this change generated additional energy corresponding to pressure in the membrane. In addition to these terms, we adopted the σ(R) term of the cumulant expansion (eqs. 7 and 13).

We discussed only a passive permeability. The Fick’ law (eq. 1) suggested that the permeability process consisted of several processes in a sequential order (the diffusions in water, the partitioning between water and membrane, the diffusion in membrane) and the LogP app could be given by a summation of the effects of these processes. In addition, our modifications (eqs. 7, 13 and 21) could be represented by linear regression models. Thus, we started from simple linear regression models as the QSPR equations based on Models A (eq. 22) and B (eq. (22).

and

where A, B, q(OH), q(NH), N atom, N rot, N HA, N HD, and MACCS represent the ASA of the i‐th atom type, the number of intramolecular hydrogen bonds, the maximum atomic charge of the H atom in the OH groups, that in the NH groups, the number of atoms, the number of rotational bonds, that of hydrogen donors, that of acceptors, and the MACCS key. The coefficient c’s (c D, c AW, c AM, c B, c R, c S, c V, c r, c l, c°H, c NH, c at°m, c r°t, c HA, c HD, c MACCS, and c 0) are the fitting parameters determined by the regression. The atom type was Sybyl mol2.

Prediction Results by Models A and B: Conformer Dependence

To examine the adequacy of Models A and B, we calculated the conformer dependence of each. Before the validation of the prediction models, we examined the conformer generation. Table 2 shows the average values of wat and mem, σ(A)wat, and σ(A)mem over all of the 795 solute molecules in the dataset, when the number of conformers was set to 100. The conformer generation represented the solute‐size change in water and in membrane. The results supported previous findings.18, 19, 20, 21, 22, 25 Namely, the solutes in membrane were folded into smaller compact structures than those in water (mem<wat), and the solutes in water fluctuated more than those in membrane (σ(R)mem<σ(R)water). These results showed a consensus with the previously reported phenomena.18, 19, 20, 21, 22, 25 Thus, we applied this conformer generation method in the present study.

Table 2

Statistics on diffusion‐related values of 795 compounds at N structure=100.

Term	Average	Deviation	Min	Max
<R>wat a	4.27	0.80	1.73	7.07
<R>mem a	3.85	0.75	1.73	6.44
σ(R)wat a	0.28	0.17	0.00	0.64
σ(R)mem a	0.03	0.05	0.00	0.51
<B>wat b	0.02	0.13	0.00	1.41
<B>mem b	0.02	0.16	0.00	2.00

a: units in Å. b: number of intramolecular hydrogen bonds.

Statistics on diffusion‐related values of 795 compounds at N structure=100.

Term

Average

Deviation

Min

Max

wat a

4.27

0.80

1.73

7.07

mem a

3.85

0.75

1.73

6.44

σ(R)wat a

0.28

0.17

0.00

0.64

σ(R)mem a

0.03

0.05

0.00

0.51

wat b

0.02

0.13

0.00

1.41

mem b

0.02

0.16

0.00

2.00

a: units in Å. b: number of intramolecular hydrogen bonds.

We applied Models A and B to the ensemble of conformers by restricting the maximum number of generated conformers up to 300. Then the 4‐fold CVs of LogP app predictions were used to estimate the Q2 and RMSE values. Figure 2 and Table S3 show the conformer‐number dependence of the LogP app prediction. Both Models A and B worked well, and the increase in conformers improved the prediction accuracy. Some previous works showed that the conformer change increased the membrane permeability.3, 6, 7, 24, 25 Namely, if the molecule has some conformers that show hydrophilic or hydrophobic surfaces, it can select stable conformers depending on the solvent (hydrophobic conformer in lipid and hydrophilic conformer in water).

Figure 2

Conformer‐number dependence of predicted properties. X axis is in Log scale. Red, green, and black circles represent the RMSE values obtained by Models A, B, and AB, respectively. (A) LogP, (B) LogD, (C) pK, (D) LogP without σ(R) terms, and (E) LogP values with n(i)=1, respectively.

Conformer‐number dependence of predicted properties. X axis is in Log scale. Red, green, and black circles represent the RMSE values obtained by Models A, B, and AB, respectively. (A) LogP, (B) LogD, (C) pK, (D) LogP without σ(R) terms, and (E) LogP values with n(i)=1, respectively.

There are two problems that make the validation of the models difficult. One of the problems is the lack of sufficient experimental observations of the conformational dynamics of solute in permeability process. Solution NMR experiments could determine the conformations in different environments including those mimic membrane interiors. However, available NMR experimental data that directly compare aqueous and lipid bilayer environments are rather limited. Often, the limited solubilities of middle molecules prevent the observation of NMR signals in aqueous solution.

In addition, most of the permeability data are for small molecules, which are less flexible than middle molecules and have smaller number of possible conformers. There are only small number of experimentally determined permeability data on middle molecules to establish the dynamics‐activity relationships. The limited number of flexible molecules with membrane permeability data also causes the problem in validating of the models based on the prediction accuracy. Recent progress of the solution NMR and molecular dynamics simulations of membrane systems should elucidate and validate the permeability mechanism in near future.

Table S4 also shows the average and maximum CPU times elapsed for one molecule. The longest CPU time was less than 62 minutes. The calculation speed should be faster than the precise MD simulation.

The results obtained by Model A were close to those by Model B. Thus, we combined the models into one by merging their descriptors. Model AB represents the united model, and the fitting parameters were determined in the same way as for Models A and B. The results obtained by Model AB were slightly better than those by Models A and B (see Figure 2 and Table S3). The RMSE values were about 0.5 and close to the deviation of the experimental data (about 0.4). Since these results showed that both Models A and B were possible permeability processes, we could not clarify the ratio of the speeds of conformer change and diffusion. The ratio should be different for each molecule in the data set.

The conformer generation represented two effects. One was the molecular‐size change depending on the solvent. The radius of the molecule in water was different from that in membrane. One conformer rich of intra‐molecule hydrophobic contacts could be dominant in water and the other conformer rich of intra‐molecule hydrogen bonds could be dominant in water. The other effect was the structural fluctuation of the molecule. The dynamics of the molecule in water were different from those in membrane. In this case, multiple conformers existed in water and in membrane, and the population of conformers in water could be different from that in membrane. The deviation terms (σ(R)) in eqs. 22 and 23 correspond to this structural fluctuation of the solute molecule. To evaluate the effect of this deviation term, we removed the deviation terms from eqs. 22 and 23. Figure 2D shows the RMSE values obtained by Models A, B, and AB without the σ(R) terms. The σ(R) terms did not improve the accuracy drastically as the number of conformers increased. These results suggested that the structural change of the solute molecule was important in the LogP app prediction, but the dynamic structural fluctuation of the solute molecule in each condition was not so important in the LogP app prediction.

We examined which terms in eqs. 22 and 23 reflect the conformer number dependence. The permeability process includes the dissociation of solutes in water and the distribution of electrically neutral solute molecules from water to membrane. These processes relate to the LogD, LogP, and pK a values. Thus, the permeability‐prediction models (eqs. 22 and 23) should predict the LogD, LogP, and pK a values by adjusting the coefficients {} of eqs. 22 and 23. Figures 2B and 2 C show the conformer dependences of the predicted LogD and pK a values, respectively. The prediction models worked and the RMSE values of these predictions were similar to those previously reported. These predicted values did not show clear dependence on the number of conformers. Using the same model, only the predicted LogP app values showed the conformer‐number dependence among these properties. Considering eq. 1, these results suggested that the conformer‐number dependence of LogP app originated from the diffusion process (M in eq. 1) in the present models and that the terms in eqs. 22 and 23 should contribute to the conformer‐number dependence of LogP app.

We also examined the RMSE and Q2 values when all the degeneracies of conformers were identical (n(i)=1 in eqs. 19 and 20). Figure 2E shows the conformer dependence of the predicted LogP app values with n(i)=1. The results were worse than those with estimated degeneracy, and increasing the number of generated conformers did not improve the RMSE. These results suggested that the conformer generation worked properly and that the estimation of conformer populations was important in the LogP app prediction.

Contribution of Diffusion Terms to LogP Prediction

We examined the diffusion process in terms of viscous (ES term, eq. 2) and inertial resistances. Table 3 shows the Q2 and RMSE values for various combinations of diffusion terms. The diffusion was described by the Log, , 2, 3, and σ(R) terms in eqs. 22 and 23. Since the Taylor series in eq. 21 includes higher‐order terms than 3, we examined the effect of 4 too. 3 was proportional to the volume of solute molecule, but it did not originally refer to volume. 2 and 3 correspond to the cross section of the solute molecule and the cross section of it that causes inertial resistance. The and Log terms represented viscous resistance. The 2 and 3 terms improved the accuracy, as did the ES term that corresponded to and Log. The higher‐order term (4) and the deviation σ(R) did not improve the results so much.

Table 3

LogP prediction results obtained by various diffusion terms at Nstructure=100.

Diffusion term including R	Model A		Model B		Model AB
	RMSE	Q2	RMSE	Q2	RMSE	Q2
None	0.54	0.77	0.54	0.77	0.54	0.77
<R>	0.54	0.77	0.54	0.78	0.54	0.77
Log<R>	0.54	0.78	0.53	0.78	0.54	0.78
σ(R)	0.54	0.77	0.54	0.77	0.54	0.77
R, σ(R)	0.54	0.77	0.53	0.78	0.54	0.78
Log<R>, σ(R)	0.53	0.78	0.53	0.78	0.53	0.79
<R>, Log<R>	0.53	0.78	0.53	0.78	0.53	0.78
<R>, Log<R>, σ(R)	0.53	0.78	0.53	0.78	0.53	0.79
{<R>n; n=1–2}, Log<R>, σ(R)	0.53	0.78	0.53	0.78	0.53	0.79
{<R>n; n=1–3}, Log<R>, σ(R)	0.53	0.78	0.52	0.79	0.52	0.79
{<R>n; n=1–2}, Log<R>	0.53	0.78	0.54	0.78	0.53	0.78
{<R>n; n=1–3}, Log<R>	0.53	0.78	0.53	0.78	0.52	0.79

LogP prediction results obtained by various diffusion terms at Nstructure=100.

Diffusion term including R

Model A

Model B

Model AB

RMSE

Q2

RMSE

Q2

RMSE

Q2

None

0.54

0.77

0.54

0.77

0.54

0.77

0.54

0.77

0.54

0.78

0.54

0.77

Log

0.54

0.78

0.53

0.78

0.54

0.78

σ(R)

0.54

0.77

0.54

0.77

0.54

0.77

R, σ(R)

0.54

0.77

0.53

0.78

0.54

0.78

Log, σ(R)

0.53

0.78

0.53

0.78

0.53

0.79

, Log

0.53

0.78

0.53

0.78

0.53

0.78

, Log, σ(R)

0.53

0.78

0.53

0.78

0.53

0.79

{n; n=1–2}, Log, σ(R)

0.53

0.78

0.53

0.78

0.53

0.79

{n; n=1–3}, Log, σ(R)

0.53

0.78

0.52

0.79

0.52

0.79

{n; n=1–2}, Log

0.53

0.78

0.54

0.78

0.53

0.78

{n; n=1–3}, Log

0.53

0.78

0.53

0.78

0.52

0.79

Molecular Weight and Ring‐Size Dependences

We examined how Models A and B worked in the prediction of the LogP app of middle molecules. If the prediction model was adequately constructed, the predicted results should not depend on the molecular size. Since molecules with MW>500 Da and those with N ring (the number of ring‐member atoms of the biggest ring) >12 are so‐called middle molecules and macrocyclic molecules, respectively, we examined the MW and N ring dependences of the predicted data distributions. Figures 3 and 4 show the prediction results obtained by the 4‐fold CV of Model B at N structure=100. Figure 3A and 4 A show the correlations between the predicted and experimental LogP app values, and Figure 3B and 4B show the principal component analysis results as the chemical space of the data points. Since Models A and AB showed the same trends as Model B, these results were omitted.

Figure 3

Predicted and experimental LogP app (A), chemical spaces (B) based on the descriptors in eq. 23 and (C) based on the Mordred descriptors in term of MW, respectively. The model used is Model B at N struct=100. Red, orange, gray, green, and blue spheres represent the molecules with 0

Figure 4

Predicted and experimental LogP app (A) and chemical spaces (B) based on the descriptors in eq. 23 and (C) based on the Mordred descriptors, respectively. The model used is Model B at N struct=100. Gray, green, and blue spheres represent the molecules with N C<12, 12

Predicted and experimental LogP app (A), chemical spaces (B) based on the descriptors in eq. 23 and (C) based on the Mordred descriptors in term of MW, respectively. The model used is Model B at N struct=100. Red, orange, gray, green, and blue spheres represent the molecules with 0

Predicted and experimental LogP app (A) and chemical spaces (B) based on the descriptors in eq. 23 and (C) based on the Mordred descriptors, respectively. The model used is Model B at N struct=100. Gray, green, and blue spheres represent the molecules with N C<12, 12

The colors in Figure 3 represent molecular weight (MW). The overall trend was that the set of small molecules showed high permeability and that of large molecules showed low permeability, but the overlaps among different color points were spread across a wide area. The overlap among the smaller and bigger molecules showed wide distribution in the chemical space. The RMSE did not depend on the MW so much. These results suggested that the present prediction model should be robust against the MW difference.

The colors in Figure 4 represent the number of ring‐member atoms of the biggest ring (N ring). The distributions of the same color data points were not localized in the LogP app prediction and chemical space, or the data points of similar size compounds spread across a wide area. The distributions of middle molecules and macrocycles overlapped that of small molecules. The RMSE did not depend on the N ring so much. In addition, Table 4 shows the dependence of the RMSE on the number of rotational bonds (N rot). N rot represents the flexibility of a molecule. As with the MW and N ring, the N rot dependence of LogP app was weak. These results suggested that the prediction model worked well in the wide molecular‐size range and should be robust against various molecular sizes and shapes.

Table 4

RMSE and Q2 values in terms of molecular size at N structure=100.

Molecular features	No. of mols	Model A		Model Bb		Model AB
		RMSE	Q2	RMSE	Q2	RMSE	Q2
MW<150	37	0.47	0.36	0.45	0.43	0.46	0.40
150<MW<500	571	0.52	0.78	0.50	0.80	0.49	0.80
500<MW	187	0.40	0.77	0.39	0.78	0.38	0.78
Nrot<10	242	0.48	0.73	0.45	0.76	0.45	0.77
10<Nrot<20	448	0.50	0.81	0.49	0.82	0.48	0.82
20<Nrot	99	0.46	0.84	0.43	0.86	0.44	0.85
Nring<6	511	0.52	0.81	0.50	0.82	0.50	0.82
6<Nring<11	171	0.41	0.71	0.39	0.73	0.37	0.76
11<Nring	111	0.47	0.80	0.42	0.84	0.44	0.82

RMSE and Q2 values in terms of molecular size at N structure=100.

Molecular features

No. of mols

Model A

Model Bb

Model AB

RMSE

Q2

RMSE

Q2

RMSE

Q2

MW<150

37

0.47

0.36

0.45

0.43

0.46

0.40

150

571

0.52

0.78

0.50

0.80

0.49

0.80

500

187

0.40

0.77

0.39

0.78

0.38

0.78

N<10

242

0.48

0.73

0.45

0.76

0.45

0.77

10

448

0.50

0.81

0.49

0.82

0.48

0.82

20

99

0.46

0.84

0.43

0.86

0.44

0.85

N<6

511

0.52

0.81

0.50

0.82

0.50

0.82

6

171

0.41

0.71

0.39

0.73

0.37

0.76

11

111

0.47

0.80

0.42

0.84

0.44

0.82

We depicted the PCA plots based on Mordred descriptors to examine the chemical space of the collected molecules in the present study.61 Mordred consists of 1826 molecular descriptors and it has been widely used in the chemoinformatics. Figures 3C and 4 C show the results. As same as Figures 3B and 3 C, the distributions of molecules were widely spread, and the groups of molecules colored according to their size and shape distributed contiguously except very small molecules <150 Da (colored in red). Thus, most of the data should be suitable for the present analysis..

Verifications of Models A and B

To verify Models A, B and AB, we also examined the L1 regularization (eq. 24) instead of the L2 regularization in eq. 16.

The L1 regularization could show the descriptors that show the major contributions to the results. The accuracy was almost equivalent to that obtained by the L2 regularization. Namely, Models A, B, and AB showed the RMSE and Q2 values of 0.73 and 0.60, respectively, in the 4‐fold cross validation tests (Table 5). The important descriptors suggested by the L1 regularization were the , Log, 2, some ASA, mainly sub‐structures including O and N atoms, and q(OH)/q(NH) among 324 descriptors including 166 MACCS keys. The results were summarized in Table S5 and Figure S1 in the supporting information. The results supported the results obtained by the present regression models.

Table 5

RMSE and Q2 values obtained by eq. 24 (L1 regularization).

	Model AB	Model A	Model B
RMSE	0.59	0.59	0.59
Q2	0.72	0.73	0.73

RMSE and Q2 values obtained by eq. 24 (L1 regularization).

Model AB

Model A

Model B

RMSE

0.59

0.59

0.59

Q2

0.72

0.73

0.73

We examined the robustness of the present method by using hold‐out tests and compared the results obtained by the Mordred‐descriptor set as an alternative method. We examined five hold‐out tests and compared the coefficient sets {} of generated regression models, in addition to the comparisons of predicted LogP app values. In each hold‐out test, the molecules were sorted by the molecular features and the top 25 % bigger molecules form the hold‐out set and the other smaller 75 % of molecules were used for the prediction model construction. The considered molecular features were the molecular weight (MW), number of atoms (N), number of ring structures (N), number of rotational bonds (N), and the number of member atoms of the maximum ring system of molecule (N). In addition to these biased sets, we prepared a non‐biased set (None) by a random selection of molecules to make a reference prediction model.

In each hold‐out test, we applied the 3‐fold cross validations to the teaching sets (the smaller molecules) and generated the prediction models. The prediction models estimated the LogP app values of the molecules in the hold‐out set (the bigger molecules). Table 6 summarizes the correlation coefficients (R) and the RMSE between the predicted and experimental LogP app values of the hold‐out set. These results show that both the present and Mordred descriptors worked well (see Figure S2 in the supporting information).

Table 6

RMSE and R values obtained by the hold‐out tests.

Molecular feature	The present descriptors (eq. 23)
	L1 (eq. 24)a		L2 (eq. 16)a
	R	RMSE	R	RMSE
Noneb	0.84	0.64	0.89	0.54
MW	0.43	0.73	0.78	0.54
Natom	0.81	0.78	0.90	0.58
Ncycle	0.45	0.62	0.59	0.53
Nrot	0.82	0.80	0.89	0.63
Nring	0.68	0.71	0.84	0.51

a: Model B was used. b: The reference models used in Table 7.

RMSE and R values obtained by the hold‐out tests.

Molecular feature

The present descriptors (eq. 23)

L1 (eq. 24)a

L2 (eq. 16)a

R

RMSE

R

RMSE

Noneb

0.84

0.64

0.89

0.54

MW

0.43

0.73

0.78

0.54

N

0.81

0.78

0.90

0.58

N

0.45

0.62

0.59

0.53

N

0.82

0.80

0.89

0.63

N

0.68

0.71

0.84

0.51

Molecular feature

Mordred descriptors

L1

L2

R

RMSE

R

RMSE

Noneb

0.81

0.68

0.89

0.53

MW

0.58

0.72

0.83

0.47

N

0.84

0.74

0.91

0.55

N

0.48

0.59

0.66

0.49

N

0.84

0.76

0.90

0.61

N

0.71

0.68

0.81

0.55

a: Model B was used. b: The reference models used in Table 7.

Then, we examined the robustness of the model construction. Table 7 shows the correlation coefficients between the coefficients {} of the prediction model obtained by the hold‐out test and that of the reference model. The {} based on the present descriptors did not depend on the difference of the teaching sets so much. On the other hand, the {} based on the Mordred descriptors depended on the difference of the teaching sets strongly. The PAMPA systems represent a passive permeability only. It means that the {} does not depend on the choice of the teaching sets. Thus, the present method should be realistic and useful rather than the collection of many descriptors.

Table 7

Correlation coefficients (R) among {} of the defferent prediction models in the hold‐out tests.

Molecular feature	The present descriptors (eq. 23)		Mordred descriptors
	L1 (eq. 24)	L2 (eq. 23)	L1	L2
NW	1.00	0.94	0.70	0.90
Natom	1.00	0.95	0.71	0.84
Ncycle	1.00	0.62	0.33	0.60
Nrot	1.00	0.90	0.75	0.86
Nring	1.00	0.88	0.81	0.63

Correlation coefficients (R) among {} of the defferent prediction models in the hold‐out tests.

Molecular feature

The present descriptors (eq. 23)

Mordred descriptors

L1 (eq. 24)

L2 (eq. 23)

L1

L2

NW

1.00

0.94

0.70

0.90

N

1.00

0.95

0.71

0.84

N

1.00

0.62

0.33

0.60

N

1.00

0.90

0.75

0.86

N

1.00

0.88

0.81

0.63

Conclusion

We proposed a QSPR method for for evaluating the apparent membrane permeability (P app) based on an analysis of the diffusion process and the partition function calculation with conformer sampling. This method generated conformers of the solute by a random structural sampling and the following structure optimization. The molecular descriptors were calculated based on the structural ensemble of the solute. Namely, the descriptors were the average and deviation values of the calculated ASA, the LogD and number of hydrogen bonds in water and in membrane, and the MACCS key. In addition, the descriptor set included a diffusion coefficient that is the inverse of the resistance in diffusion. To estimate the diffusion of solute, we examined the inertial resistance in the diffusion of solute in addition to the viscous resistance. We assumed two types of permeability models of solute with multiple conformers. One was the fast diffusion process (Model A), in which the solute diffused in the membrane with a fixed conformer, whose fractions were the same as those in water. The other was the slow diffusion process (Model B), in which the solute changed the conformer in the diffusion process and the fractions of conformers followed the most stable distribution in the membrane. The present QSPR models represented Models A and B based on the molecular descriptors mentioned above.

This prediction method worked in the P app prediction of the middle molecules and macrocycles, the same as with that of the small molecules. The results suggested that the inertial resistance should be important in the diffusion, as is the viscous resistance known as the Einstein‐Stokes equation. The ensemble of conformers improved the prediction accuracy. This study supported both Models A and B, and the permeability process could be a combination of Models A and B.

Conflict of Interest

None declared.

Supporting Information

The appendices, Tables S1–S5, Figures S1 and S2 were supplied as described in the Supporting Information.

As a service to our authors and readers, this journal provides supporting information supplied by the authors. Such materials are peer reviewed and may be re‐organized for online delivery, but are not copy‐edited or typeset. Technical support issues arising from supporting information (other than missing files) should be addressed to the authors.

Supplementary

Click here for additional data file.