Structure-based prediction of drug distribution across the headgroup and core strata of a phospholipid bilayer using surrogate phases.

Solvation of drugs in the core (C) and headgroup (H) strata of phospholipid bilayers affects their physiological transport rates and accumulation. These characteristics, especially a complete drug distribution profile across the bilayer strata, are tedious to obtain experimentally, to the point that even simplified preferred locations are only available for a few dozen compounds. Recently, we showed that the partition coefficient (P) values in the system of hydrated diacetyl phosphatidylcholine (DAcPC) and n-hexadecane (C16), as surrogates of the H- and C-strata of the bilayer composed of the most abundant mammalian phospholipid, PC, agree well with the preferred bilayer location of compounds. High P values are typical for lipophiles accumulating in the core, and low P values are characteristic of cephalophiles preferring the headgroups. This simple pattern does not hold for most compounds, which usually have more even distribution and may also accumulate at the H/C interface. To model complete distribution, the correlates of solvation energies are needed for each drug state in the bilayer: (1) for the H-stratum it is the DAcPC/W P value, calculated as the ratio of the C16/W and C16/DAcPC (W for water) P values; (2) for the C-stratum, the C16/W P value; (3) for the H/C interface, the P values for all plausible molecular poses are characterized using the fragment DAcPC/W and C16/W solvation parameters for the parts of the molecule embedded in the H- and C-strata, respectively. The correlates, each scaled by two Collander coefficients, were used in a nonlinear, mass-balance based model of intrabilayer distribution, which was applied to the easily measurable overall P values of compounds in the DMPC (M = myristoyl) bilayers and monolayers as the dependent variables. The calibrated model for 107 neutral compounds explains 94% of experimental variance, achieves similar cross-validation levels, and agrees well with the nontrivial, experimentally determined bilayer locations for 27 compounds. The resulting structure-based prediction system for intrabilayer distribution will facilitate more realistic modeling of passive transport and drug interactions with those integral membrane proteins, which have the binding sites located in the bilayer, such as some enzymes, influx and efflux transporters, and receptors. If only overall bilayer accumulation is of interest, the 1-octanol/W P values suffice to model the studied set.

Introduction

The interactions with the headgroup (H) and core (C) strata of phospholipid bilayers play key roles in membrane protein folding[1] and channel function,[2] as well as in absorption and distribution of peptides, nucleic acids, drugs, and other small organic molecules in organisms.[3]

The bilayer H- and C-strata differ in the ability to shield electrostatic charges and in the types of interactions, in which the dissolved molecules participate. Despite intense thermal motion,[4] the strata are clearly distinguishable in neutron[5] and small-angle X-ray[6] diffraction experiments, as well as in molecular dynamics simulations[7] of the phosphatidylcholine (PC) bilayer. The main H- and C-strata can be further subdivided based on density,[7] voids,[8] and the average positions of PC fragments.[9] While most studies refer to individual layers as regions, we prefer to use the term “strata” so that the term “regions” can be reserved for distinct lateral areas (also called lipid rafts) caused by immiscibility of some lipids in bilayer[10] or interactions of lipids with proteins,[11] cholesterol,[12] and other compounds.[13]

Among almost 200 known human lipid species, PC represents 42–49% of phospholipids and 0.7–2.5% w/w of fresh tissue (containing 70–80% of water), for skeletal muscle,[14,15] heart, and liver.[15] Linoleic acid (18:2) with 37.7% is the most abundant fatty acid in PC of human skeletal muscle, followed by palmitic acid (16:0) with 27.3% and oleic acid (18:1) with 14.0%.[16] In the adipose tissue, the contents are 10%, 23%, and 48%, respectively.[15] Considering a usual pattern, a saturated acid and an unsaturated acid esterifying the PC glycerol’s O-1 and O-2 atoms, respectively, it can be inferred that 1-palmitoyl-2-linoleoyl-sn-glycero-3-phosphocholine and 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine are the most abundant phospholipid molecules in skeletal muscle and adipose tissue, respectively. While the lipid diversity plays significant roles in physiological processes, most partition data have been obtained with PC containing saturated fatty acids because of easier synthesis and resistance to oxidation. Among them, dimyristoylphosphatidylcholine (DMPC) was used most frequently thanks to a lower transition temperature, oxidation stability, and availability, although a recent study pointed out that DMPC and other PC bilayers may differ in accumulation of some compounds.[17]

Several complementary experimental methods, mainly wide-angle[18,19] (WAXS) and small-angle X-ray scattering (SAXS)[20,21] and neutron diffraction,[17,22−24] various nuclear magnetic resonance (NMR)[25−28] and electron paramagnetic resonance (EPR)[29] techniques, surface plasmon resonance,[30] as well as fluorescence[31] and fluorescence quenching,[32] have been used for characterization of the preferred location of compounds in the phospholipid bilayer.[33] Except WAXS, SAXS, and neutron diffraction, these methods usually do not provide a complete intrabilayer distribution profile. All methods are tedious and not amenable to a higher throughput. On a routine scale, typically only the overall bilayer partitioning is measured, without obtaining the details pertinent to individual bilayer strata.

To extrapolate the valuable bilayer location data to other compounds and quickly obtain estimates of solvation free energies of chemicals in the bilayers or their strata, organic solvents have been examined as surrogates. The partition coefficients P in two solvent systems, which are capable of similar interactions with the studied compounds, are related according to the Collander equation:[34]The parameters α and β are obtained by the fit to experimental data. The approach is widely used in the design of bioactive compounds,[35] computational chemistry,[36] protein folding,[37] and other areas.[38] The Collander equation (eq 1) is one of the numerous extrathermodynamic relationships, which include linear free energy relationships (LFERs),[39] and cover a broad area of physical organic chemistry. The LFERs describe correlations between free or activation energy changes with chemical structure, in two or more processes. Such changes in reference processes can be used to define LFER parameters characterizing individual substructures, such as the Hammett constants[40] and fragment solvation characteristics,[35] or entire molecules such as solvatochromic parameters.[39,41] The requirement of similarity of the two processes holds for simple linear correlations, such as the logarithmized Collander equation (eq 1), but may be waived if multiple LFER parameters are used.

Additive-constitutive nature of the partition coefficient P provides for its deconvolution into the solvation characteristics (f) of structural fragments. This trait allows prediction of the P values from structure and the modeling of partitioning of amphiphilic molecules, which may interact with the H/C interface and protrude in both H- and C-strata. The f values in the surrogate phases for individual strata can be summed up to estimate the overall solvation free energy. The f values represent enthalpy and entropy changes upon transfer of individual fragments between the strata. The additional entropy changes associated with the entire molecule are not treated explicitly but are, to some extent, accommodated by the enthalpy–entropy compensation, which is a common phenomenon in partitioning.[42−44] The two attributes of the P values, Collander scaling (eq 1) and additive-constitutive nature, create the basis for the development of the reference systems, such as the widely used 1-octanol/water (O/W) system, and also lie at the heart of the present approach.

The key issue for the selection of a proper surrogate is the similarity of interactions with the imitated system. In this aspect, diacetyl phosphatidylcholine (DAcPC) contains all structural fragments and represents the closest match to the PC headgroup. High aqueous solubility of DAcPC allows its use in the hydration state close to that of a fluid PC bilayer at room and body temperatures, represented by 6–16 water molecules per a headgroup,[45−48] and making hydrated DAcPC a potential surrogate phase of the H-stratum of the PC bilayer. Whether the lacking anisotropy and the racemic nature of hydrated DAcPC affect its ability to emulate solvation behavior of a PC bilayer remains to be seen.

Hydrocarbons are appropriate core surrogates as evidenced by experimental observation of similar molecular packing[49] and dynamics[50] of the fatty acyl chains in the core of the bilayer and in bulk liquid alkanes. The studies with black lipid membranes show that the selectivity of the bilayer in the partitioning of the H-bonding groups lies somewhere between those of 1-octanol and alkanes[51,52] or alkenes[53] but is significantly closer to the latter solvents for the studied compounds.

DAcPC, when dissolved in water in molar ratio 1:14 used in our experiments, forms a homogeneous, isotropic, and only slightly viscous solution.[54] The DAcPC molecules do not measurably partition from this solution into n-hexadecane (C16) in the two-phase system. We showed that the C16/DAcPC partition coefficients perfectly discriminate between the lipophiles accumulating in the core and cephalophiles preferring the H-stratum, in contrast to the O/W and C16/W partition coefficients alone or their ratio.[55] This observation indicates the suitability of hydrated DAcPC and C16 as the surrogate phases of the H- and C-strata.

As criteria for selecting a proper surrogate, early studies used partitioning into cells and transport through multicellular biological films, which are both affected by protein binding, multiplicity of bilayers, and other factors. These problems are avoided when using unilamellar phospholipid liposomes, which are more suitable for evaluation of the surrogate phases. Liposome partitioning depends on the temperature, which should be well above the gel–fluid phase transition temperature,[56] fluidity,[57] phospholipid surface density,[58] and other related factors.[59]

Computational approaches for studying macromolecules and drugs in the bilayer have used different levels of approximation for the bilayer structure, ranging from implicit continuum models, to coarse-grained mesoscale models,[60] to fully atomistic representations.[61] Our approach formally represents a simple model belonging to the first category.

Continuum approaches started with conformational optimization using a classical force field in a pure hydrocarbon slab immersed in water, with different dielectric constants for the aqueous and nonpolar phases.[62] Initial insertion of chemicals and peptides was guided by the hydrophobic and hydrophilic centers of the molecules, which were determined using the atomic O/W transfer energies.[63] A hydrophobic potential energy function was devised to represent the bilayer environment[64] and combined with a cavity formation term.[65] The depth-dependent dielectric constant and atomic hydrophobicity parameters were used to estimate the solute solvation energies in the bilayer.[66] A restrain function was developed to generate peptide conformations typical for bilayers.[67] All these studies relied on the O/W and alkane/W partitioning data, which, as we have shown,[55] do not correctly emulate solvation in bilayer strata as they fail to identify the cephalophiles and lipophiles accumulating in the H- and C-strata of the PC bilayer. Solvation energies can be estimated using two components: (1) intramolecular electrostatic energy, calculated by classical or quantum mechanical techniques and (2) the “nonpolar” energy spent for cavity formation and released by interactions with the solvent, approximated by techniques based on the atomic solvent-accessible surface area (SASA) combined with empirical atomic surface tensions, derived from partition experiments.[68] A well-known application to different solvents is represented by the series of the SMx models.[69] The two-component concept was applied to bilayer solvation in several studies differing in sophistication of the model setup and used computational techniques. A simple hexane slab was used to study the preferred bilayer location of steroids,[70] with the finite-difference solutions to the Poisson equation providing the electrostatic part,[71] and a linear relation to the SASA defining the nonpolar part. A multistrata dielectric model setup was used to study membrane proteins in combination with generalized Born electrostatics and the SASA multiplied by a depth-dependent function, which was obtained from the free-energy profile of oxygen in the bilayer, generated by molecular dynamics (MD) simulation.[72] Similar techniques have been implemented in the CHARMM package.[73] A multistrata, continuum COSMOmic (mic for micelle) approach[74] is built around the density functional theory (DFT) calculations combined with conductor-like screening model for realistic solvation (COSMO-RS).[75] The bilayer details are obtained from atomistic MD simulations, allowing for study of mixed bilayers.[76,77] A detailed PPM (Positioning of Proteins in Membranes) approach combines the trans-bilayer profiles of the dielectric constant and some solvatochromic parameters[39,41] with the electrostatics calculated using the dipole moments rather than atomic charges for nonpolar solutes and the Born equation for ions.[78,79]

Mesoscale models were mainly used in the MD studies of structural aspects of bilayer dynamics and were deemed unsuitable for partitioning of small molecules.[80] This situation is gradually changing with the development of more sophisticated coarse-grained force fields, which can be smoothly transformed to atomistic models in the areas of interest.[81−83]

Atomistic MD simulations were used to analyze bilayer interactions with small molecules, with the first studies focusing on benzene,[84,85] dimethyl 1,4-dihydro-2,6-dimethyl-4-(2-methylphenyl)-3,5-pyridinedicarboxylate (a nifedipine analogue),[86] trichloroethylene,[87] and a 29 amino acid fragment of corticotropin-releasing factor.[88] Numerous other MD studies revealed several important aspects of trans-bilayer diffusion, such as varying diffusion rates in individual strata, a faster speed of lateral than transversal diffusion,[84] and a preferred longitudinal orientation of transported molecules along the bilayer normal in the most dense interfacial stratum, unless strong H-bonds are formed.[89] Several studies emphasized preferred molecular positions and orientations during the diffusion process.[90,91] Free-energy estimates were obtained using various biasing MD techniques, to reduce the computing time.[92] These data are getting more precise with the development of specialized force fields for lipids.[93−95] The force fields were compared for their ability to reproduce the bilayer characteristics and the overall bilayer/water partition coefficients of several compounds.[96]

Here we present a simple continuum approach that uses a coarser representation of the PC bilayer than the most detailed continuum models: the homogeneous H- and C-strata, connected through a planar H/C-interface passing right below the fatty acid ester groups of PC. This simplification allows us to directly use the solvation energies of molecules and their fragments, which were obtained from partitioning experiments with relevant surrogate phases. We hypothesize that quantitative estimates of the intrabilayer drug distribution can be obtained by modeling measured overall partitioning of compounds into the DMPC bilayers in unilamellar liposomes and the DMPC monolayers adsorbed on nonporous, octadecylated silica microparticles. All reasonable interfacial poses of each molecule are examined in the model. To reduce complexity in generating the poses, simple molecules with reduced conformational flexibility and avoiding ionization under physiological conditions were selected. The correlation also provides, for each compound, estimated distribution in bilayer strata, including the preferred pose(s) at the H/C interface. The predictions of preferred location in the PC bilayer are compared with experimental data for the compounds, where this information exists.

Experimental Section

Chemicals

DAcPC was obtained from Euticals (Prime European Therapeuticals S.p.A, Lodi, Italy). All studied solutes (Table 1) and C16 were purchased from Sigma-Aldrich (St. Louis, MO).

Table 1

Studied Compounds, Partition Coefficients (P) in the Solvent and PC Systems, Number of Possible Interfacial Poses (NP), and Predicted Fractions in Bilayer Strata

		logP								relative amount (%)h
					monolayer/W		bilayer/W
no.	compound name	C16/Wa	C16/DAcPCb	O/Wc	exptld	calcde	exptlf	calcde	NPg	H	I	C
1	1,2-dichlorobenzene	3.52i	2.24	3.43	3.28o	3.28	3.33	3.28	13	10	20	70
2	1,3,5-trichlorobenzene	4.48i	2.89	4.19	3.71o	4.15	3.95q	4.12	10	2	17	81
3	1,3-dichlorobenzene	3.69i	2.46	3.53	3.22o	3.50	3.03	3.51	15	5	35	60
4	1,3-dimethyl-2-nitrobenzene	2.48	1.80	2.95	2.52	2.69	2.99	2.78	12	22	55	23
5	1,4-dichlorobenzene	3.70i	2.41	3.44	3.22o	3.43	3.56	3.40	7	7	16	77
6	1,4-dimethylbenzene	3.25	2.59	3.15	3.49o	3.00	2.98q	2.98	7	14	9	77
7	1-bromo-2-phenylethane	2.35i	2.07	3.09	2.89	2.81	2.71	2.93	16	13	75	12
8	1-bromo-3-phenylpropane	3.66	2.18	3.72	3.95	3.60	3.04	3.65	18	4	55	40
9	1-bromo-4-chlorobenzene	3.78	3.34	3.54	3.24	3.50		3.47	14	4	17	79
10	1-bromonaphthalene	4.17	3.43	4.06	3.96	3.91		3.90	18	2	30	68
11	1-hexanol	0.44	–0.01l	2.03		2.55	1.91q	2.69	12	24	76	0
12	1-naphthaldehyde	2.21	1.36	2.67n	2.62o	2.41	2.80o	2.49	16	47	29	25
13	1-naphthol	0.44i	–1.18	2.85	2.82	2.43	2.80	2.58	17	56	43	0
14	2,2′,3,3′,4,4′-hexachlorobiphenyl	5.98	8.36	7.32		5.53	6.09q	5.47	9	0	5	95
15	2,2′,3,3′,4,5,5′,6-octachlorobiphenyl	7.25	10.37l	8.16n		6.71	6.18q	6.66	16	0	1	99
16	2,2′,3,3′,4-pentachlorobiphenyl	5.98	7.35	6.30	5.94	5.53		5.47	16	0	5	95
17	2,2′,3,3′,6,6′-hexachlorobiphenyl	5.98	8.36	7.12		5.53	5.90q	5.47	8	0	5	95
18	2,2′,3,3′-tetrachlorobiphenyl	5.35	6.35	6.18	-	4.93	5.40q	4.88	8	0	5	95
19	2,2′,3,4,4′,5,5′-heptachlorobiphenyl	7.25	8.90	8.16n	6.91p	6.72	6.33q	6.66	16	0	2	98
20	2,2′,3,5,6′-pentachlorobiphenyl	6.14	6.89	6.72n	6.21	5.72		5.69	16	0	19	81
21	2,2′,3,5′-tetrachlorobiphenyl	5.83	5.88	5.81	5.47o	5.38		5.33	8	0	5	95
22	2,2′,3,6-tetrachlorobiphenyl	5.35	6.35	6.01n	4.99o	4.99		4.96	15	0	22	78
23	2,2′,3-trichlorobiphenyl	5.04	5.34	5.60	5.26o	4.68		4.66	14	0	20	80
24	2,2′,4,4′,6,6′-hexachlorobiphenyl	6.49j	7.43	7.29		6.06	6.00p	6.04	9	0	22	78
25	2,2′,4,5,5′-pentachlorobiphenyl	6.62	6.42	6.85	5.68o	6.14		6.10	16	0	10	90
26	2,2′,4,5-tetrachlorobiphenyl	5.83	5.88	6.27	5.96p	5.43	5.55q	5.40	16	0	20	80
27	2,2′,4,5′-tetrachlorobiphenyl	6.30	5.42	6.36	6.26	5.84		5.80	9	0	9	91
28	2,2′,4-trichlorobiphenyl	5.51	4.88	5.76	5.06o	5.15	5.03	5.14	16	0	25	74
29	2,2′,5,5′-tetrachlorobiphenyl	6.27j	5.42	6.26	5.69o	5.80		5.75	8	0	6	94
30	2,2′,5-trichlorobiphenyl	5.51	4.88	5.60	4.92p	5.14		5.12	15	0	22	78
31	2,2′,6-trichlorobiphenyl	5.04	5.34	5.47	4.70o	4.68		4.65	8	0	19	81
32	2,2-diphenylethanol	1.05i	–0.08	2.77n	2.99	2.32	2.34	2.46	8	58	40	2
33	2,2′-dichlorobiphenyl	4.72	4.33	4.97	4.37o	4.38		4.35	7	1	19	81
34	2,3,3′,6-tetrachlorobiphenyl	5.83	5.88	6.26n	5.47o	5.38		5.33	13	0	6	94
35	2,3,4,5,6-pentachlorobiphenyl	6.62	7.82l	6.74		6.12	5.90q	6.06	16	0	1	99
36	2,3,4,5-tetrachlorobiphenyl	5.83	6.35	6.41		5.39	5.44q	5.35	16	0	9	91
37	2,3,4′,6-tetrachlorobiphenyl	5.83	5.88	6.26n	6.07	5.40		5.35	13	0	10	90
38	2,3′,4-trichlorobiphenyl	5.99	4.41	5.92n	5.10o	5.55		5.50	10	0	10	90
39	2,3′,5-trichlorobiphenyl	5.99	4.41	5.76	4.74o	5.54		5.48	10	0	6	94
40	2,3,6-trichlorobiphenyl	5.04	5.34	5.67	5.00o	4.69		4.67	15	0	23	77
41	2,3′,6-trichlorobiphenyl	5.51	4.88	5.67	5.00o	5.14		5.12	15	0	23	77
42	2,3-diaminonaphthalene	–0.22	–1.70m	1.54	2.39	2.30		2.44	21	72	28	0
43	2,3-dichlorobiphenyl	4.72	4.33	4.99	4.85o	4.40		4.37	14	0	23	77
44	2,3′-dichlorobiphenyl	5.20	3.87	5.02	4.53p	4.84		4.82	10	0	22	77
45	2,4,5-trichloroaniline	2.14	1.04	3.69	3.31o	3.36	3.58p	3.50	19	5	93	2
46	2,4,5-trichlorobiphenyl	5.74i	4.88	5.90		5.33	5.16q	5.29	16	0	15	85
47	2,4,6-trichlorobiphenyl	5.51	4.88	5.71		5.15	4.71q	5.14	16	0	26	74
48	2,4′,6-trichlorobiphenyl	5.51	4.88	5.75	5.26o	5.15		5.14	16	0	26	74
49	2,4-dichlorobiphenyl	5.20	3.87	5.30	4.91o	4.86		4.84	16	0	26	74
50	2,4′-dichlorobiphenyl	5.20	3.87	5.10	4.85o	4.86		4.84	15	0	26	74
51	2,5-dichlorobiphenyl	5.20	3.87	5.16	4.91o	4.84		4.82	14	0	23	77
52	2,5-dimethylphenol	0.37i	0.36	2.33	2.32	2.03		2.17	13	65	34	1
53	2,6-dichlorobiphenyl	4.72	4.33	4.98	4.37o	4.38		4.36	11	1	20	80
54	2,6-dimethoxyphenol	–2.28	–1.14	1.15	1.99	1.64	1.65	1.78	9	89	11	0
55	2-bromonaphthalene	4.17	3.43	4.18n	4.05	3.92		3.92	20	2	32	66
56	2-bromostyrene	3.73	3.32	3.73n	3.20	3.49		3.48	12	4	27	69
57	2-bromotoluene	3.50	3.46	3.50	3.40	3.29		3.29	16	5	30	65
58	2-chlorobiphenyl	4.41	3.33	4.53	4.24p	4.09	4.09	4.07	7	1	20	79
59	2-hydroxybiphenyl	0.75	–0.15	3.09	3.15o	2.27	3.21p	2.42	15	57	42	1
60	2-methylanthracene	5.05	3.41	5.00	4.15	4.76		4.76	23	0	36	64
61	2-nitroaniline	0.22i	–0.89	1.85	1.99	2.49	2.45	2.64	15	38	62	0
62	2-nitrotoluene	2.25i	1.22	2.30	2.34	2.48	2.41q	2.57	15	43	34	23
63	3-bromoaniline	1.04k	–0.15	2.10	2.10	2.44	2.36	2.58	17	45	53	2
64	3-nitroaniline	–0.68i	–1.41	1.37	2.18o	2.02	2.17q	2.16	16	94	6	0
65	3-nitrotoluene	1.81	1.29	2.45	2.60	2.41	2.56	2.53	18	37	54	9
66	4,4′-dichlorobiphenyl	5.67	3.41	5.58		5.21	4.78q	5.25	9	0	11	89
67	4-amino-3-methylphenol	–2.63	–1.19m	0.70n	1.80o	1.53	1.41p	1.68	15	99	1	0
68	4-aminoacetophenone	–0.70i	–1.95	0.86	1.67	2.11		2.26	12	97	3	0
69	4-aminophenol	–3.34	–3.04	0.04	1.76	1.77	1.82	1.91	14	100	0	0
70	4-biphenylcarboxaldehyde	3.01	2.18	3.38n	3.08o	2.99	3.14r	3.04	23	13	47	40
71	4-bromoaniline	0.65	–0.26	2.11	2.92	2.31	2.47	2.45	14	53	46	1
72	4-bromobenzophenone	2.46k	2.43	4.12n	2.58	2.53	2.66	2.59	14	25	42	34
73	4-bromophenol	–0.01i	–1.20	2.59	2.56	2.40	3.01	2.55	14	48	51	0
74	4-chloro-3-methylphenol	0.22k	–1.08	3.10	2.73o	2.97	3.19p	3.12	16	14	86	0
75	4-chlorobenzophenone	2.86k	2.49	3.79	3.45	2.78	2.57p	2.81	14	18	33	49
76	4-chlorotoluene	3.48	2.47	3.33	3.06o	3.25		3.24	14	9	21	70
77	4-dimethylamino-benzaldehyde	2.50	0.05	2.01	1.83	2.62		2.69	14	65	6	29
78	4-nitroaniline	–0.72k	–2.23	1.39	2.16	2.17	2.01	2.32	14	96	4	0
79	4-nitrotoluene	1.81	1.20	2.37	2.10	2.37	2.40	2.49	14	42	47	10
80	4-n-pentylphenol	2.11	0.89l	4.24		4.17	4.31q	4.02	22	1	98	0
81	4-phenoxybutyl bromide	3.45	2.18	3.50n	4.18	3.42	3.82	3.47	22	6	55	39
82c	9-anthracenemethanol	2.34k	–0.42	3.04	3.21o	3.18	3.50r	3.32	8	18	77	5
83	acetophenone	1.14i	0.62	1.58	1.73	2.08	2.01	2.21	14	75	20	5
84	aniline	–0.10k	–0.52	0.90	2.32o	1.99	1.63q	2.13	12	87	13	0
85	anisole	2.09i	1.44	2.11	2.30	2.26	2.31	2.35	14	59	14	27
86	anthracene	4.20i	2.53	4.45	4.08o	4.03		4.06	13	2	47	51
87	benzaldehyde	1.06i	0.79	1.48	1.80	1.96		2.09	12	88	7	5
88	benzene	2.15i	1.27	2.13	2.43o	2.31	2.10	2.39	5	60	12	27
89	benzylalcohol	–0.43i	–1.18	1.10	1.85	2.24	2.03	2.38	14	57	43	0
90	biphenyl	4.08i	2.32	4.01	3.61o	3.80		3.78	7	4	21	75
91	bisphenol A	–0.91	–2.70	3.32	2.87	2.24	2.95	2.39	10	94	6	0
92	chlorobenzene	2.84i	2.74	2.89	3.15o	2.77	2.81q	2.81	12	16	37	47
93	dibutyl phthalate	4.32	5.24	4.72	3.54	4.00		3.97	13	1	19	81
94c	ethanol	–2.19i	–2.08	–0.31	1.81	1.81	1.96s	1.96	4	99	1	0
95	ethylnicotinate	0.55k	0.27	1.32	1.79	1.97	1.91	2.11	16	85	14	2
96	fluoranthene	5.00	2.75	5.16	4.93	4.66		4.64	12	1	24	75
97	methyl-4-chloro-2-nitrobenzoate	0.95	0.98	2.38n	2.13	1.88	2.27	2.01	12	91	4	5
98	methyl-4-nitrobenzoate	0.63	0.90	1.85n	2.26	1.82		1.96	12	92	5	3
99	methylbenzoate	1.56i	1.81	2.12	1.88	1.96		2.07	14	71	13	16
100	N,N-dimethylaniline	2.17	1.20	2.31	2.03	2.48	2.78	2.57	14	41	40	19
101	n-propylbenzene	3.84	2.98	3.72	3.72	3.64	3.29	3.65	16	3	37	59
102	naphthalene	3.41i	1.90	3.30	3.19o	3.22		3.22	9	12	26	62
103	nitrobenzene	1.47k	0.77	1.85	1.80o	2.18	2.01q	2.30	12	67	26	8
104	phenanthrene	4.74	2.53	4.46	4.28p	4.41		4.38	11	1	20	79
105	phenol	–1.08i	–1.68	1.47	1.71	2.04	1.87	2.18	12	84	16	0
106	pyridine	–0.45k	–0.89	0.65	1.89	1.98	1.86s	2.12	11	89	11	0
107	toluene	2.76i	1.93	2.73	2.43o	2.68	2.56	2.71	12	28	23	49

Estimated by the fragment method,[97] if not marked otherwise.

Measured by us and published,[55] if not marked otherwise.

Experimental values from the ClogP database,[119] if not marked otherwise.

Measured in this study, if not marked otherwise.

Predicted values from the present model-logarithmized eq 10 combined with eq 11, with the coefficient values given in the last line of Table 3.

Measured in this study by isothermal titration calorimetry (ITC) incorporation method,[102] if not marked otherwise.

Number of possible poses at the interface.

Fraction of the bilayer amount in headgroups (H), at interface (I), and in core (C). More details in Table S1 in the Supporting Information.

Experimental values published.[120]

Published.[111]

Measured by us and published.[97]

Estimated by fragment method.[55]

Measured by us and published.[54]

Estimated by the ClogP approach.[119]

Measured by us and published.[98] Estimated values from this reference were not used.

Measured in this study and supersedes the published measured or estimated value.[98]

Published values for polychlorinated biphenyls (PCB)[57] and for other compounds.[107]

Measured by us using the ITC incorporation method and published.[98]

Measured in this study using the ITC solvent null method.[104]

C16/DAcPC and C16/W Partitioning

The experiments in the pre-equilibrated, stirred two-phase systems were described in detail previously.[55,97] The kinetics of partitioning was monitored, in both phases when feasible, to make sure that the equilibrium concentrations were determined. The compounds were analyzed by UV−Vis, HPLC, or GC/MS after SPME extraction as appropriate. Evaporation was accounted for and the compounds were checked for instability and self-aggregation under experimental conditions.

Partitioning to Phospholipid Monolayers

Preparation of DMPC monoloayers on octadecylated, nonporous silica particles of uniform size was described previously, along with the partition experiments.[98]

Partitioning to Liposomes

Liposome Preparation

Unilamellar DMPC liposomes with diameters of about 80 nm were prepared by the extrusion method[99] using a mini-extruder (Avanti Polar Lipids, Alabaster, AL) with 0.1 μm filter (Nucleopore polycarbonate Track-Etch membrane, Whatman, Florham Park, NJ) at 37 °C. Before extrusion, powdered DMPC (Avanti Polar Lipids, Alabaster, AL) was hydrated in borate buffer (5 mM Na2B4O7, 180 mM H3BO3, and 18 mM NaCl, pH = 7.4) for 48 h at 37 °C with occasional vortexing. The suspension was extruded for 29 times to achieve even distribution of vesicle sizes. After extrusion, the exact DMPC concentration was determined using phosphate assay.[100]

Isothermal Titration Calorimetry (ITC) Measurement of Partitioning of Chemicals to Liposomes

The DMPC bilayer/water partition coefficients in liposome suspensions were measured using the VP calorimeter (MicroCal, Northampton, MA). The calorimeter cell has a volume of 1.4347 mL, and the syringe has a volume of 294.26 μL. Before loading in the cell, all solutions were subjected to 10 min degassing in order to get rid of the air bubbles, which may cause errors in the solution volume in the cell. The temperature of all used solutions was kept at 37 °C, i.e., above the main transition temperature of DMPC, in order to keep the bilayer in the fluid state. A stirring speed of 310 rpm was set for all the experiments. MicroCal Origin software was used to perform all the fitting procedures.[101]

The setup of the experiments was chosen to maximize precision. On the basis of the affinities of the tested chemicals for the bilayer, two different setups were used: the incorporation protocol and the solvent null approach.

Incorporation Protocol

This method,[102] using the titration of the solution of the chemicals by the liposome suspension, was used to determine the bilayer/water partition coefficients for the majority of the studied chemicals. The suspension of DMPC liposomes was loaded in the syringe, and the cell was filled with the solution of the chemical in borate buffer (pH 7.4). Liposome suspension was injected in the aliquots of 5 or 10 μL. For each compound, three background titrations were performed: (1) the buffer into the solution of the chemical to account for the heat of dilution and to confirm the absence of self-association of the chemical, (2) liposome suspension into the buffer to account for heat of dilution of the lipid, and (3) the buffer into the buffer to account for apparent heat caused by possible small differences between the temperatures of the solutions in cell and in syringe.

The corresponding heat generated after the ith injection, δhi, can be expressed aswhere cD is the concentration of the chemical in the cell after the ith injection, ΔHD is the transfer enthalpy of the chemical from the aqueous phase to the bilayer, K is the mole-ratio partition coefficient, δcL is the increase of the lipid concentration in the cell due to each injection of lipid suspension, and Vcell is the volume of the cell. Two parameters, ΔHD and K, were optimized by the fit of eq 2 to experimental data. The mole-ratio partition coefficient, K, was then converted to molarity-based partition coefficient, P, usingassuming that the concentrations of the chemical in the bilayer (cb) and aqueous phase (cw) are both sufficiently low to approximate activities. The subscripts D and L indicate total concentrations (c) or amounts (n) of the chemical and lipid, respectively. The density of the DMPC in liposomes[103] of ∼1 g/mL was used to calculate the bilayer volume (Vb).

Solvent Null Protocol

This method[104] uses the titration of the solutions of the chemical with different concentrations by the equilibrated mixture of the liposomes with the chemical. If the chemical concentration in the calorimeter cell does not match the free concentration of the chemical in the suspension loaded in the syringe, heat will be either absorbed or released upon mixing. No heat will be generated only when the concentration in the cell is equal to the free concentration in the suspension. This so-called solvent null concentration is used to calculate the actual amount of the chemical bound to lipid.

The solvent null method was applied to those chemicals that have very low partition coefficients so that the curvature of the δhi plots in the incorporation protocol was insufficient to allow for a good fit by eq 2. Aliquots of 10 μL DMPC liposome chemical-buffer suspension were injected to various concentrations of chemical solutions in borate buffer (pH 7.4) loaded in the cell. Before loading, the liposome suspension was combined with the solution of the chemical at room temperature, and the mixture was incubated at 37 °C for 1 h. For each experiment, the corresponding controls were performed: (1) the chemical in the null concentration titrated into all concentrations of the chemical in the cell; and (2) the neat lipid titrated into buffer. The molarity-based partition coefficient, P, was calculated aswhere n0 is the amount of the chemical mixed with the liposome suspension; nf is the free amount determined through the null concentration, and Vb and Vcell represent the volumes of the bilayer and the cell, as defined previously.

Modeling Drug Distribution in Bilayer Strata

Experimental information about overall partitioning of drugs into PC monolayers or bilayers can be conveniently obtained using common analytical techniques, some of which are described above. The measurements can be made directly on the suspension of monolayer particles or liposomes in equilibrated drug solution (by UV–Vis difference spectroscopy,[105] isothermal titration calorimetry,[102,104] and pH-metric titration for ionizable drugs[106]) or, alternatively, on the equilibrated drug solution after separation of liposomes or monolayer particles. The molarity-based partition coefficients (P) between the PC phase and water are calculated using the estimated volumes of monolayer or bilayer PC phases.

The partitioning of a compound into the PC monolayer or bilayer comprises its accumulation in the H- and C-strata, and the interface between them. The accumulation at the interface can be described using the surface-based interface/water equilibrium constant Si/w (dimension 1/length), representing the ratio of the surface concentration si (per unit area) and the aqueous concentration. The surface of the interface is practically equal to the surface of the H-stratum (Ah) because of the small thickness (hh) of the H-stratum as compared to the curvature radius of the liposome or the particle. The drug amount at the interface is ni = siAh = ciVh. The second equality comes from the conversion of the surface-based concentration to the volume-based concentration, with the reference to the H-stratum volume. This step needs to be done because the fragment solvation energies, which will be used to describe the interface partitioning, are given for the volume-based concentrations, and the volumes are used to express the drug amounts in the H- and C-strata. The relationship between the surface-based equilibrium constant Si/w and the volume-based partition coefficients Pi/w = ci/cw is

The partition coefficients (P) values for the entire monolayer or bilayer can be expressed as the ratio of drug concentration in the PC phase and water:where c are the concentrations, and V are the volumes. The subscripts x, w, h, c, and i denote the entire monolayer or bilayer (x = m or b, respectively), the aqueous phase including the internal aqueous phase of liposomes, the H-stratum, the hydrocarbon C-stratum, and the H/C interface, respectively. The interface/water partitioning is characterized by a concentration-based partition coefficient, with the concentration of the interface-bound species expressed with regard to the headgroups’ volume, to maintain consistency with the descriptions of the other two partitioning processes.

The surrogate phases for the H- and C-strata are introduced into the model using the Collander eq 1. The core is imitated by the C16 phase. The surrogate partition coefficient for the headgroups, PDacPC/W, cannot be measured directly but can be calculated as the ratio of two other partition coefficients:The regression coefficients α and β have different values for the H-stratum, the C-stratum, and the H/C interface, the latter expressed using eq 8 (subscripts h, c, and i, respectively).

Various poses (orientations) of a compound at the H/C interface may contribute to the observed Pi/w value:Here, the subscripts ij mark the volume-based interface concentration (cij) and the partition coefficient Pij of the jth pose at the interface. The Pij is defined using two sums of the logP-related fragment solvation characteristics: one using the DAcPC/W characteristics for the fragments located in the H-stratum and the other using the C16/W characteristics for the fragments located in the C-stratum.

Results

The data set assembled for testing the hypothesis about estimating drug distribution in bilayer strata from modeling overall bilayer partitioning data is presented in Table 1. For 107 compounds, there are 95 Pm/w values for the monolayer partitioning and 63 Pb/w values for the bilayer partitioning. The Pm/w and Pb/w values for the same compound, where available, do not differ much because the headgroup and core phases of monolayers and bilayers behave in similar way, and the difference is only caused by the volumes of the two strata. The published logPb/w were measured in unilamellar DMPC liposomes[57,107] at 35–37 °C. Other published DMPC partitioning data were not used because either the experiment temperature[108−113] was too close to the DMPC gel–fluid transition temperature (24 °C)[114] or the compound characteristics for eqs 10 and 11 below were missing.[42,115−117] The used compounds either cannot ionize in water or have negligible ionization in the used medium with pH = 7.4, i.e., the bases have pKa ≤ 5.4 and the acids have pKa ≥ 9.4 (estimated by ACD Percepta,[118] see Table S2 in Supporting Information). Most C16/W and all C16/DAcPC partition coefficients were either measured in our lab or estimated using the solvation characteristics for fragments obtained using the ClogP fragmentation scheme,[35] which we published previously.[55,97] Bilayer partitioning data were measured in DMPC liposomes by us or others, and in DMPC monolayers adsorbed on octadecylated nonporous silica microspheres by us. All available data are included in Table 1. Some of the data from our laboratory were previously published. Some predictions from the calibrated model, the partition coefficients and the intrabilayer distribution, are also listed in Table 1.

Table 3

Regression Coefficients and Statistical Indices of Calibrated Model (eqs 10 and 11)

αc	βc	αh	βh	αi	βi	r2	SD	F
1	1	1	1	1	1	0.792	0.607	149.38
4.229 ± 2.560	0.783 ± 0.059	1	1	1	1	0.821	0.563	173.60
1	1	209.3 ± 48.90	0.123 ± 0.108	1	1	0.893	0.436	289.96
1	1	1	1	21.37 ± 4.061	0.353 ± 0.048	0.913	0.392	357.94
1	1	265.7 ± 48.37	0.231 ± 0.078	0.768 ± 0.221	0.832 ± 0.170	0.924	0.366	410.13
0.741 ± 0.407	0.973 ± 0.045	1	1	20.71 ± 3.589	0.381 ± 0.043	0.928	0.356	434.52
0.222 ± 0.218	1.023 ± 0.073	225.308 ± 43.302	0.129 ± 0.088	1	1	0.924	0.368	407.23
0.954 ± 0.465	0.946 ± 0.039	272.0 ± 42.97	0.215 ± 0.069	0.551 ± 0.129	0.918 ± 0.129	0.941	0.324	523.22

For individual variables in Table 1, the minimum and maximum values, the average, and the standard deviation are −3.34, 7.25, 3.16, and 2.49 for logPC16/W; −3.04, 10.37, 2.62, and 2.85 for logPC16/DAcPC; −3.12, 2.76, 0.54, and 1.02 for logPDAcPC/W (calculated as logPC16/W – logPC16/DAcPC, used in the model but not shown in Table 1); −0.31, 8.16, 3.88, and 1.93 for logPO/W; 1.67, 6.91, 3.51, and 1.35 for logPm/w; 1.41, 6.33, 3.25, and 1.36 for logPb/w.

Distributions for all variables are shown in Figure 1. Only the logPC16/DAcPC values exhibit normal distribution at the 0.05 probability according to the Shapiro–Wilk test,[101] although the deviations from normality for other variables are not dramatic. The distributions for the C16/W and O/W partition coefficients are slightly skewed because of a number of nonpolar aromatic molecules, which have the logP values of ∼6. These compounds have lower partition coefficients in the C16/DAcPC system thanks to the interactions of aromatic rings with DAcPC, and this fact keeps the C16/DAcPC distribution normal. The calculated values for the DAcPC/W system continuously rise to a sharp maximum at logPDAcPC/W = 1 but then exhibit a fast decline leading to a slightly asymmetrical distribution. The maximum logPDAcPC/W = 1 ± 0.5 comprises about a half of the studied set. In contrast to the other two solvent/water systems, the logPO/W values lack the negative values because of a higher partitioning of hydrophilic compounds in 1-octanol than in nonpolar solvents thanks to the high water content and H-bonding in 1-octanol. The monolayer and bilayer data show incomplete normal distributions; they lack the lower logP values, which are difficult to measure. For the solvent/water systems, the low logP measurements can be performed using varying volume ratio of the used phases.

Figure 1

Distribution of independent variables (logPC16/W, green; logPC16/DAcPC, magenta; logPDAcPC/W, black; logPO/W, blue) and dependent variables (logPm/w, red; logPb/w, brown) (Table 1). The bins are one unit wide and centered at the integer logP magnitudes.

Model of Drug Distribution in Bilayer

The concentration-based, overall partition coefficient (P) between the entire monolayer or bilayer (subscript x = m or b, respectively) and water (subscript w) is described as a weighted sum of the partition coefficients in individual strata (H-stratum − headgroups, C-stratum − core, and H/C interface denoted by subscripts h, c, and i, respectively), where the weights are relative volumes of the H- and C-strata (see eq 6 for details):Here, V are the volumes of the phases indicated by the subscripts. The interface/water partitioning is characterized by a concentration-based partition coefficient (Pi/w), with the concentration of the interface-bound species expressed with regard to the headgroups’ volume, to maintain consistency with the descriptions of the other two partitioning processes, and to account for the double interface area of the bilayer as compared to the monolayer. This description is permissible as long as the surface density of the interface-bound molecules is sufficiently low to ensure independent binding. More details are provided in the text accompanying eq 6.

Using the Collander equation (eq 1), the surrogate phases for the H- and C-strata are introduced into eq 9. The headgroups are imitated by the hydrated DAcPC phase and the core by the C16 phase:The surrogate partition coefficient for the headgroups is calculated as PDAcPC/W = PC16/W/PC16/DAcPC. Compounds may interact with the H/C interface in various poses (orientations, conformations, and positions). Therefore, the observed Pi/w value in eq 10 is expressed as the sum of the partition coefficients for each of the poses (see eq 8 for more details):The subscript ij indicates the relationship to the jth pose at the interface. The partition coefficient (Pij) for each pose is defined as the partition coefficient obtained from fragments (Pf), scaled by the Collander coefficients as in eq 1. Each logPf value is generated using two sums of the logP-related fragment solvation characteristics (f, summarized in Table 2), which were determined previously for the fragmentation used in the ClogP system[35] by regression analysis of the DAcPC/W and C16/W partition coefficients.[55,97] In the first sum (k-summation), the DAcPC/W characteristics are used for the fragments interacting with the headgroups (H-location). The second sum (l-summation) contains the C16/W characteristics for the fragments located in the core (C-location). The interfacial poses exhibit the Boltzmann distribution, with the fraction of the jth pose being equal to Pij/ΣPij.

Table 2

Fragment Solvation Characteristics for C16/DAcPC, C16/W, and DAcPC/W Systems

fragment		fragment solvation parameters
symbol	full name	C16/DAcPCa	C16/Wb	DAcPC/Wc
Ca	aromatic [a] isolating carbon	0.084	0.077	–0.007
CA	aliphatic [A] isolating carbon	0.336	0.108	–0.228
H	hydrogen on isolating carbon	0.134	0.287	0.153
OHa	hydroxyl [a]	–2.651	–2.788	–0.137
OHA	hydroxyl [A]	–2.788	–3.656	–0.868
OHZ	hydroxyl [Z]	–2.976	–2.967	0.009
OAa	ether [Aa]	–0.106	–0.955	–0.849
EsAa	ester [Aa]	–0.144	–1.299	–1.155
Ala	aldehyde [a]	–0.283	–1.052	–0.769
COAa	carbonyl [Aa]	–1.278	–1.774	–0.496
COaa	carbonyl [aa]	–0.096	–1.562	–1.466
NH2a	primary amine [a]	–1.866	–2.166	–0.300
NAAa	tertiary amine [AAa]	–0.646	–1.120	–0.474
NO2a	nitro [a]	–0.768	–0.592	0.176
Naraa	aromatic nitrogen [aa]	–1.340	–1.611	–0.271
Cla	chloride [a]	0.698	1.078	0.380
Bra	bromide [a]	1.401	1.096	–0.305
BrA	bromide [A]	–0.247	0.069	0.316
		Correction Factors
BB	benzyl bond to simple aromates	0.139	–0.177	–0.316
Chain	chain	–0.197	–0.057	0.140
CHBr	chain and cluster branch	–0.290	–0.125	0.165
FAC	extended aromatic iso-C	0.001	0.149	0.148
HB	pair-H bond	0.627	1.323	0.696
NOrtho1d	ortho substitution 1	0.038	0.147	0.109
NOrtho2d	ortho substitution 2	0.436	–0.476	–0.912
NOrtho3d	ortho substitution 3	0.3103	–0.27	–0.580
PCCY	phenyl-fragment pair	–0.065	0.244	0.309
PIWR1e	potential interaction within ring 1	–0.52	0.092	0.612
PIWR2e	potential interaction within ring 2	–0.632	0.114	0.746
PIWR3	potential interaction within ring 3	–0.211	0.370	0.581

Published.[97]

Published.[55]

Calculated as the difference between the C16/W and C16/DAcPC values.

For the interacting fragments with the f > −0.12 (Northo1), −0.25 ≤ f ≤ −0.12 (Northo2), and f < −0.25 (Northo3). The average is used if the f values are from different categories.[97]

For the interacting fragments with the f < 0.15 (PIWR1), 0.15 ≤ f ≤ 0.5 (PIWR2), and f > 0.5 (PIWR3). The average is used if the f values are from different intervals.[97]

Individual poses were created for each compound using the following algorithm, illustrated for ethanol (94, Table 1) in Figure 2: (1) generate the most extended conformation; (2) orient the molecule in this conformation so that its long axis is parallel to the bilayer normal, to minimize bilayer disturbance as observed in MD simulations;[89,121] (3) create individual poses by moving the oriented molecule from one phase to the opposite phase across the H/C-interface in small, preferably atom-by-atom, steps leading to different numbers of atoms in the H- and C-strata; and (4) repeat step 3 with inverted molecule, so that both possible orientations along the long axis are examined. If a molecule cannot adopt a clearly prolonged shape, and the long axis is not significantly longer than other axes, repeat steps 2–4 also for the axis perpendicular to the long axis.

Figure 2

Set of generated interface poses for ethanol (94). Individual logPf values, used in eq 11, are obtained by summing up the fragment solvation characteristics (Table 2) for the phase, where the fragment is located. The chain factor, correcting for the chain length, was scaled based upon the location of the two carbons forming the chain. This principle was applied to other correction factors, which are applied to a larger substructure: they were composed of the contributions from each phase, sized by the proportion of the substructure present in the given phase.

The model represented by logarithmized eq 10 combined with eq 11 was calibrated by nonlinear regression analysis[122] based on experimental data (Pm/w and Pb/w in Table 1) using the independent variables PDAcPC/W = PC16/W/PC16/DAcPC and PC16/W (Table 1) and f (Table 2). The regression coefficients α and β may have different values for the H- and C-strata and the H/C interface (subscripts h, c, and i, respectively). On the basis of the experimental data[48] and molecular dynamics simulations,[123] the volume ratio of both headgroup strata to the core in the DMPC bilayer was estimated as Vhb/Vc = 7:13 and the H/C volume ratio in monolayers was estimated as Vhm/Vc = 1:3. These volume ratios were kept fixed in the optimization process.

The optimized values of regression coefficients α and β for the fit of logarithmized eq 10 combined with eq 11 to experimental data given in Tables 1 and 2 are summarized in Table 3. The number of experimental data points was n = 158 in all cases.

Interestingly, even without calibration, the model explains 79% of experimental variance (line 1 in Table 3). This is a typical model quality for calibrated models based on frequently used topological descriptors, which often bear no physical relationship to the solvation process. Optimization of any pair of the Collander coefficients α and β (lines 2–4) enhances the fit, with most improvement seen for the optimized i/w coefficients leading to the squared correlation coefficient value of r2 = 0.913 and a 35% reduction in the value of the standard deviation (SD = 0.392). When two pairs of the Collander coefficients are optimized (lines 5–7), the agreement with the experiment increases, up to r2 = 0.928 and SD = 0.356 for the best combinations of the c/w and i/w coefficients. The best result, 94% of explained variance and SD = 0.324, is achieved when all six values of the α and β coefficients are optimized (line 8). Apparently, no reduction in the number of coefficients leads to a better fit, as indicated by the increasing F-test values (F). The standard errors of coefficients are acceptable, especially considering that all of them are nonlinear coefficients in the logarithmized eq 10 combined with eq 11.

Given the number of data points (158) and the number of optimized coefficients (6), there is low probability of a chance correlation. The stability of the correlation is confirmed by the cross-validation based on the omission of 21 (13.3%) data points at a time, which was repeated ten times. The omission of compounds was performed by a random selection of three compounds from each of the seven bins of about equal size for the data points ordered by the magnitude of the Px/w values (x = m or b), to ensure that the entire Px/w range is represented. This binning does not cluster similar compounds because the overall partitioning depends on the affinity for the H-, C-, and H/C strata, and a small structural change may cause a large change in the Px/w value. The average predictive squared correlation coefficient, q2 = 0.924, corresponds well with r2 = 0.941 for the full fit. The average predictive sum of squares of deviations PRESS = 2.315 leads to a value of the corresponding SD equivalent (∼square root of 2.315/21) equal to 0.332, which is commensurate with the SD = 0.324 for the full training set.

The calibrated model, logarithmized eq 10 combined with eq 11, with the optimized coefficient values in the last line of Table 3, predicts the drug distribution in individual strata of the DMPC bilayer (Table 1, the last three columns, and Table 4). For the interfacial poses, more detailed characterization is provided in Figure 3 and Table S1 in Supporting Information, listing relative drug amounts in individual interfacial poses as well as in both H- and C-strata.

Table 4

Predicted Bilayer Locations for Compounds with Known Bilayer Locations

		log P
					bilayer/water			relative amount (%)f,g
no.	compound	C16/Wa	C16/DAcPCb	O/Wc	exptld	calcd	NPe	H	I	C	exp. locationg
108	1-propanol	–1.53	–1.23	0.25		1.94	6	94	6	0	H
109	3-methylindole	0.81	0.88	2.17		2.11	14	71	26	3	H
110	9H-carbazole	2.21	0.74	3.06		2.64	10	45	37	18	Hi
88h	benzene	2.15a	1.27b	2.13c	2.10	2.39	5	60	12	27	H
89h	benzylalcohol	–0.43a	–1.18b	1.10c	2.03	2.38	14	57	43	0	H
91h	bisphenol A	–0.91a	–2.70b	3.32c	2.95	2.39	10	94	6	0	H
94h	ethanol	–2.19a	–2.08b	–0.31c	1.96	1.96	4	99	1	0	H
111	indole	0.75	0.14	1.67		2.26	10	70	28	2	H
112	methanol	–2.77	–2.05	–0.77		1.83	2	100	0	0	H
106h	pyridine	–0.45a	–0.89b	0.65c	1.86	2.12	11	89	11	0	H
113	1-butanol	–0.81	–0.83	0.88		2.00	9	96	4	0	H/Ij
114	1-pentanol	–0.81	–0.42	1.35		2.25	10	44	56	0	H/Ij
82h	9-anthracenemethanol	2.34a	–0.42b	3.04c	3.50	3.32	8	18	77	5	Ik
11h	1-hexanol	0.44a	–0.01b	2.03c	1.91	2.69	10	24	76	0	I
115	1-heptanol	1.06	0.39	2.41		3.19	14	9	91	0	I
116	1-octanol	1.69	0.80	2.94		3.74	16	3	97	0	I
117	4-tert-octylphenol	3.49	0.96	5.16		4.65	15	1	96	3	I
6h	1,4-dimethylbenzene	3.25a	2.59b	3.15c	2.98	2.98	7	14	9	77	C
118	1,6-diphenyl-1,3,5-hexatriene	5.83	4.25	5.64		5.44	17	0	27	73	C
119	1-methyl-4-(6-phenyl-1,3,5-hexatrien-1-yl)-benzene	6.33	4.93	6.14		5.95	30	0	32	68	C
120	9-methylanthracene	5.05	3.26	5.07		4.66	6	1	18	82	C
121	9-ethylanthracene	5.68	3.67	5.52		5.22	6	0	13	87	C
122	ethylbenzene	3.20	2.46	3.15		3.26	14	11	36	53	C
123	n-decane	7.01	4.74	5.98		6.45	9	0	7	93	C
124	n-hexane	4.49	3.11	3.90		4.06	5	2	8	90	C
101h	n-propylbenzene	3.84a	2.98b	3.72c	3.29	3.65	16	3	37	59	C
107h	toluene	2.76a	1.93b	2.73c	2.56	2.71	12	28	23	49	C

Predicted values according to ref (97) except the marked values, which are listed in Table 1.

Predicted values according to ref (55) except the marked values, which are listed in Table 1.

Estimated by the ClogP approach[119] except the marked values, which are listed in Table 1.

Listed in Table 1.

Number of possible poses at the interface.

Fraction of the amount in the bilayer strata (total amount in the bilayer = 100%).

Headgroups, H; interface, I; core, C.

Listed in Table 1 and used in model calibration, which did not utilize the location information.

See Discussion.

Published as I only; see Discussion for explanation.

Published as H, but I is more plausible; see Discussion for explanation.

Figure 3

Distribution of the prevalences of interfacial poses of individual compounds predicted by the calibrated model. The poses were generated using the algorithm described in Figure 2 and the accompanying text. The compounds are arranged in the order of decreasing interface accumulation (left to right). The poses are listed in the order of decreasing prevalence. White color indicates the poses, which were not generated (PNG).

Inside the calibrated model, logarithmized eq 10 combined with eq 11, the sum of the Pij values for individual poses is dominated by prevalent poses, and contributions of weak poses are marginalized. After calibration, individual Pij values provide information about the prevalence of the given pose. Individual poses for the used compounds are summarized in Table S1 in Supporting Information, along with the fraction of the compound present in each pose. The number of poses for each compound is also listed in Table 1. Inclusion of configurational entropy using the number of individual poses was attempted but did not lead to an improvement of the fit.

The number of considered poses for a compound increases with the length of the molecule in the extended conformation or with the absence of such conformation, when the poses along two axes need to be created. For individual compounds, 4–23 possible interfacial poses were generated, and all were simultaneously fed into the model. The distribution of predicted prevalences of poses for individual compounds is depicted in Figure 3. The model effectively reduces the number of considered poses. The number of significant poses (>2%) is four or less for 56 (52%) compounds. Note that the first three shades mark low, insignificant prevalences and are only shown to illustrate the numerical aspect of the approach. Figure 3 also shows that the number of considered poses does not vary substantially between compounds exhibiting a high preference for the interface and the rest of the compounds. Expected amphiphilicity of a compound, or the lack of thereof, is not considered in generating possible interfacial poses.

Model Validation

The presented model was calibrated using the overall partition coefficients and predicts, in addition to these characteristics, intrabilayer distribution. These latter predictions were chosen for independent validation of the model because they are more challenging than the overall partition coefficient (see Discussion). For several dozen compounds, the preferred locations in the fluid PC bilayer were experimentally determined, as summarized previously, along with the used techniques.[33] For all compounds, for which experimental or predicted independent variables (PC16/W and PC16/DAcPC) and the fragment solvation characteristics (f) for both C16/W and C16/DAcPC systems were available, preferred locations are compared with predictions from our model (Table 3, the last line) in Table 4. Compounds with available Pb/w values in Table 4 were included in the set used for model calibration (Table 1), but the information about the preferred location was not used in that process.

The validation data set is not trivial and the locations are in many cases hard to guess using standard chemical intuition. The preferred location depends on the character of fragments forming a molecule. Small aromates prefer the headgroups (benzene, 88; pyridine, 106), as do larger aromates containing heteroatoms (9H-carbazole, 110; indole, 111), contrary to the common perception that they should accumulate in the core because of high hydrophobicity as indicated by their O/W partition coefficients (Table 4). Aromates lose the preference for the headgroups by addition of alkyl substituents, as seen in toluene (107), ethylbenzene (122), n-propylbenzene (101), and 1,4-dimethylbenzene (6), or an alkenyl chain as in 1,6-diphenyl-1,3,5-hexatriene (118) and its methyl derivative (119). Naturally, alkanes such as n-hexane (124) and n-decane (123) then reside almost exclusively in the core. Halogenation of aromates increases their affinity for the core thanks to the positive values of the Cla and Bra fragment parameters listed in Table 2. A single hydroxyl group can pull a molecule from the core into the interface, as seen for the 9-methylantracene (120)/9-anthrancenemethanol (82) pair, or even into headgroups, as illustrated by the toluene (107)/benzylalcohol (89) pair. However, there are no simple qualitative rules for estimation of preferred bilayer location, which is determined by a delicate balance of fragment solvation characteristics in the molecule. For instance, the addition of a single methyl group is sufficient to change the preference from the headgroups to the core for the benzene (88)/toluene (107) pair but not for the indole (111)/3-methylindole (109) pair.

The calibrated model faithfully captures this complex behavior (Table 4). The agreement between the model predictions and experimental data is very good taking into account that the experimentally determined preferred location means prevalent but not exclusive accumulation in that location. Only for two (1-butanol, 113; 1-pentanol, 114) out of 27 compounds are the experimentally determined preferred locations different from the prevalent location indicated by the model, and the experimental data are not definite in these two cases (see Discussion). The experimentally determined preferred strata have more than 50% and often significantly higher prevalence predicted in all cases except 9H-carbazole (110), which is still predicted correctly as having the highest accumulation in the headgroups.

Discussion

The presented approach to prediction of the bilayer distribution of compounds is unique thanks to the use of hydrated DAcPC as the H-stratum surrogate. Solvatochromic analysis[55] showed that the DAcPC phase, in comparison to bulk water, maintains higher H-bond acceptor ability but has diminished H-bond donor ability. These observations are contradictory to the often used assumption about similar solvation properties of the headgroup stratum and bulk water,[70,124−127] stemming from the high headgroup hydration.

The calibrated model (Table 3, the last line) describes well the overall partition coefficients of 107 compounds, which do not ionize under experimental conditions. Extensive cross-validation shows that the model is stable and predictive. The Collander coefficients have the optimized values close to unity for the core (both αc and βc) and the H/C interface (βi only), indicating high similarity and collinearity, respectively, between the natural and surrogate partition coefficients. For headgroups, there is a significantly nonlinear relationship between the natural and surrogate partition coefficients, described by the low βh value, which is compensated by the high αh value. This fact could indicate that the Collander coefficients for headgroups account for some phenomena, which were not included in the model, e.g., electrostatic interactions, entropy changes, and/or cavity creation (although the latter aspect is partially covered by the fragment solvation characteristics).

Interpretation of Some Model Validation Data

Experimentally determined preferred locations of 27 compounds were used for independent validation of the model. The preferred location was predicted correctly for 25 compounds, and the experimental data were not conclusive for the remaining two compounds, 113 and 114 (Table 4). Two more compounds, 82 and 110 (Tables 1 and 4), also deserve some discussion.

The preferred location of 9-anthracenemethanol (82) was determined by the parallax analysis, which uses fluorescence quenching fragments covalently bound at different positions along the PC molecule.[32,128] The method does not provide much structural detail and places the center of the molecule at the extrapolated peak quenching depth. In the published pose, with the long axis perpendicular to the bilayer normal, the molecule of 82 is located at and above the ester groups of fatty acids, a position qualifying for the headgroup location. The published pose, however, does not minimize the disturbance of the interface and is therefore less probable than a pose with the long axis parallel to the bilayer normal, as observed in MD simulations for other compounds.[89,121] In the more probable pose, the atoms of one of the outer rings would reach into the core, and the pose would be classified as amphiphilic, in contrast to our previous assessment,[33] which used the conclusion about the headgroup location made in the original studies.[32,128] As seen in Table 4, our model prefers the interface location (74%), along with some headgroup accumulation (20%).

The 9H-carbazole (110) is the only correctly predicted compound where the predicted prevalence did not exceed 50%. Location of 110 in the bilayer was also determined by the parallax analysis[128,129] as being about 3 Å closer to the headgroups than that of 82. By applying the same considerations as for 82 in the preceding paragraph, 110 with the long axis parallel to the bilayer normal would be classified, in accordance with the model prediction (Table 4), as a cephalophile (predicted prevalence 44%), being very close to the interface (predicted prevalence 36%).

For n-alkanols, the preferred accumulation in headgroups was determined independently for ethanol (94).[130] It can be expected that with the chain elongation, n-alkanols become more amphiphilic, and at some point, the interactions with the interface prevail. There are two reports pinpointing the chain length when this change happens in 1-palmitoyl-2-oleoyl-PC vesicles. The lifetime and intensity ratios of fluorescence probes in the presence of n-alkanols indicated that the first three congeners (112, 94, and 108) prefer headgroups and that longer congeners (113, 114, 11, 115, 116) accumulate at the interface.[131] Photolysis of rhodopsin, an integral membrane GPCR protein, in the presence of n-alkanols placed the break point in the series on 1-hexanol (11), for which the effect was negligible, while a promoting effect was observed for shorter n-alkanols (112, 94, 108, 113, and 114), and an inhibitory effect was seen for longer n-alkanols (115 and 116).[132] Our model captures the transition from the headgroup preference to the interface preference with the alkyl chain extension but puts the transition point on 1-pentanol (114), i.e., between the outcomes of the two studies.

Predicted Distribution in Bilayer Strata

The calibrated model, logarithmized eq 10 combined with eq 11, with the optimized coefficient values given in the last line of Table 3, provides a prediction of drug distribution in individual strata of the DMPC bilayer for the compounds, for which the input data, the C16/W and the C16/DAcPC partition coefficients, and fragment solvation characteristics (Table 2) are available. The interfacial poses are generated using the algorithm described in Figure 2 and the accompanying text. The fragment-based partition coefficients of individual poses, Pij, are calculated (eq 11) using the fragment solvation characteristics f in the DAcPC/W and C16/W systems for the fragments embedded in the headgroups and the core, respectively. The interfacial poses exhibit the Boltzmann distribution, with the fraction of the jth orientation being equal to Pij/ΣPij.

The predicted values of relative drug amount in the headgroups, core, and at the interface (Tables 1 and 4) are plotted against the C16/DAcPC partition coefficients in Figure 4. The sigmoidal dependences marking, for the majority of compounds, the maximum fraction in the H- and C-strata can be discerned. The sigmoidal function for the H-stratum (FH, in %) as dependent on the respective logP values (the red line in Figure 4) is based on the definition of the C/H partition coefficient:[55]The exponent β was added to account for a different composition of the bilayer strata and the surrogate systems using the Collander equation,[133] but it was not necessary because its value could be kept fixed as β = 1. The fraction in the C-stratum is, obviously, equal to 100 – FH (the black line in Figure 4). The functions adopt reasonable shapes if the coefficient is set to α = 0.05.

Figure 4

Relative amounts accumulated in headgroups, core, and the headgroups/core interface (red, black, and blue points, respectively) plotted against the C16/DAcPC partition coefficients (Tables 1 and 4). Points represent the predicted values, and stars mark the compounds with known bilayer location. The lines (see text for the equations) indicate the maximum fraction in the headgroup and core strata (red and black, respectively) attainable for the given P value.

As seen in Figure 4, the PC16/DAcPC value alone has a limited use as a predictor of bilayer location; it seems to only recognize cephalophiles having almost 100% accumulation in the H-stratum as compounds with the logPC16/DAcPC < −3. Lipophiles with logPC16/DAcPC > 3 prefer to accumulate in the C-stratum but also exhibit a minor but significant (up to 30%) adsorption at the interface, probably because of the aromatic character of the used lipophilic compounds. The compounds with −3 < logPC16/DAcPC < 3 show widely varying distributions between individual bilayer strata. The asymmetric peak of blue points marking the interface accumulation starts rising at logPC16/DAcPC = −3, reaches the expected maximum for 0 < logPC16/DAcPC < 1, but then decreases slowly, forming a long shoulder with the prevalence around 20% for up to logPC16/DAcPC ≈ 8. This shoulder contains a significant number of halogenated aromatics.

For several compounds, the prediction of the overall bilayer distribution, including the significant interfacial poses, is depicted in Figure 5. The shown molecular geometries and orientations approximate real interfacial poses, while absolutely no geometry or orientation is specified for the molecules accumulated in the H- and C-strata. The relative amounts in Figure 5 represent fractions of the total amount present in the bilayer. These bilayer amounts differ for individual compounds, if they are added to a bilayer/water system in equal doses. Drug distribution between the bilayer and aqueous phases can be calculated using the Pb/w values in Tables 1 and 4 and assuming no protein binding. Using a realistic estimate for muscles, liver, and heart,[14,15] the volume of bilayers represents 2.5% of the total tissue volume. Then, the fractions of the total amount in the tissue that are present in the bilayer are 81% for indole (111, Figure 5A); 68% for 1-hexanol (11, B); 99% for 9-anthracenemethanol (82, C); 97% for 4-chloro-3-methylphenol (74, D); and 93% for 1-bromo-2-phenylethane (7, E).

Figure 5

Predicted fraction (%) of the overall bilayer amount accumulated in the headgroup and core strata (pink and gray areas, respectively) and at the interface of the fluid PC bilayer for: A, indole (111); B, 1-hexanol (11); C, 9-anthracenemethanol (82); D, 4-chloro-3-methylphenol (74, one interface pose with prevalence of 16% is omitted for clarity); E, 1-bromo-2-phenylethane (7). The shown interfacial poses coarsely correspond to expected experimental poses.

For indole (111, Table 4 and Figure 5A), peculiar experimental data are available showing bimodal distribution, with predominant binding to the headgroup phosphates but also some adsorption to the H/C-interface.[134] Interestingly, the presented approach is able to capture both headgroup accumulation and interfacial poses, which are very close to observed positions, and even approximate their prevalences (70% and 28%, respectively, Table 4), if all interfacial poses are lumped together. Figure 5B and 5C show well-separated prevalences for the H-stratum accumulation, three interfacial poses, and low core partitioning, which are predicted for 1-hexanol and 9-anthracenemethanol (11 and 82, Table 4), respectively. The predictions are in good agreement with experimental data, classifying these compounds as amphiphilic (Table 4). A similar pattern holds for 4-chloro-3-methylphenol (74), as shown in Figure 5D, although this compound exhibits smaller differences in the preferences for the H-stratum accumulation and individual interfacial poses (one pose with 16% prevalence, Table S1 in Supporting Information, is omitted for clarity). For 1-bromo-2-phenylethane (7, Figure 5E), the bilayer prevalences are much more comparable than for the other shown compounds. Interestingly, the interfacial poses are not similar as in the previous cases. The pose with the lowest shown prevalence resides in the headgroups except the benzene ring end, in contrast to the higher ranked poses, which only have the bromine end interacting with the headgroups.

The importance of intramolecular fragment positions for amphiphilicity is nicely illustrated using 2-, 3-, and 4-nitroaniline (61, 64, and 78, respectively, Table 1), having the overall logPb/w values in the same range (2.45, 2.11, and 2.01, respectively). However, their intrabilayer distribution is rather different: 64 and 78 partition almost completely to the headgroups, while 61 shows higher preference for the interface. This fact cannot be explained by the DAcPC/W partition coefficient (=P16/W/P16/DAcPC), which actually has a higher value for 61 than for 64. All three compounds contain cephalophilic amino and nitro groups, with the C16/DAcPC fragment solvation parameters f(NH2) = −2.166 and f(NO2) = −0.592 (Table 2, the suffix a denotes an aromatic skeleton). While the ortho-derivative 61 can have both groups located in the headgroups and the rest of the molecule crossing the interface into the core resulting in the interface prevalence of 62%, this option is greatly diminished for the meta-derivative 64 (6%) and the para-derivative 78 (4%).

Intrabilayer Distribution: Comparison with Other Approaches

Numerous computational studies mostly focused on qualitative prediction of bilayer distribution for one or two compounds at a time. For a comparison with our very good model predictions (Table 4), we selected studies predicting the preferred location of several compounds.

A multistrata, continuum approach, PPM (Positioning of Proteins in Membranes), was also applied to bilayer partitioning of small molecules.[79] None of the six compounds with experimentally determined bilayer locations had the location reproduced correctly. These results are a bit surprising given the level of sophistication of the approach and can probably be improved by refocusing the approach on partitioning of small molecules.

Atomistic depth-constrained MD simulations with several all-atom and united-atom lipid force fields were compared for their ability to predict partition coefficients of 11 compounds (see below).[96] The preferred locations of individual molecules (Table S6, Supporting Information)[96] varied between force fields, indicating that this characteristic is more difficult to reproduce than the overall partition coefficient. Experimentally determined location data are only available for a few compounds. For 1-butanol (113, Table 4) preferring the headgroups[132] or the interface,[131] the free energy minimum was placed at 8.9–11.1 Å from the bilayer center by four of the five tested force fields, and only one force field provided the more realistic distance of 19.8 Å. Benzylalcohol (89, Tables 1 and 4) prefers the headgroups,[26] but all force fields placed the free energy minimum in the more ordered core portion at 8.9–11.1 Å from the bilayer center. For 1,4-dimethylbenzene (6, Tables 1 and 4) accumulating in the bilayer center,[135] two force fields made the correct prediction, and three force fields placed the free energy minimum in the more ordered portion of the core, 6.1–10.0 Å from the center.

The same study[96] also used a multistrata, continuum COSMOmic approach[74] to examine bilayer partitioning of 11 compounds. Predictions of the free energy minima for the three compounds above (6, 89, and 113) did not agree with experimental data.

Overall Bilayer Partitioning: Comparison with Other Approaches

The monolayer and bilayer logP values, as predicted by the present model (Table 3, the last line), are summarized in Table 1. The used approach takes into account the volume of individual strata, so the monolayer and bilayer data are described by the same model. Overall accumulation in both systems is predicted satisfactorily: the fit has the r2 = 0.941 and the SD = 0.324 for 158 data points. We also examined other approaches that are used for modeling of overall accumulation in biological systems: the partition coefficients for the 1-octanol/water system using the logarithmized Collander equation (eq 1)[34] and the solvatochromic approach.[136] The correlations of the monolayer and bilayer accumulation for 95 and 63 compounds, respectively (Table 1), needed to be performed separately because the approaches do not include the volumes of the strata.

The correlations of the logPm/w and logPb/w values with the logPO/W values are very good, as indicated by the statistical indices, summarized in the first two lines of Table 5. The simplest approach using the PO/W values works well for neutral compounds, in accord with previous results.[33,137] In the overall linear trend, more precise linear correlations for polar and nonpolar chemicals, as well as for congeneric series were distinguishable.[33,107] The PO/W values are very well suited for QSAR modeling involving overall accumulation, especially if the PO/W parameter is a part of a more complex disposition function.[138] The situation may change for compounds with PO/W < 0.01 (not present in our data set), which may be affected by the high water content (4% vol) in the water-saturated 1-octanol,[139] depending upon how much of this water can be used for hydration of drug molecules. The PO/W parameter did not perform well in correlations for liposome partitioning of ionizable compounds, forming mixtures of ionized and neutral molecules at physiological conditions because of different interactions of neutral and charged molecules.[140,141]

Table 5

Overall Accumulation (logP) in the PC Systems Correlated with logPO/Wa and the Solvatochromic Parameters

	coefficients
PC system	a	b	s	e	v	const.	n	r2	SD	F
monolayera	0.729 ± 0.020					0.882 ± 0.082	95	0.965	0.355	1339
bilayera	0.654 ± 0.022					1.057 ± 0.086	63	0.966	0.351	895.2
monolayer	0.0651 ± 0.197	–2.665 ± 0.235	–0.642 ± 0.210	0.693 ± 0.149	2.443 ± 0.181	0.934 ± 0.158	95	0.934	0.346	262.5
bilayer	0.382 ± 0.220	–2.672 ± 0.277	–0.875 ± 0.230	0.754 ± 0.242	2.093 ± 0.244	1.444 ± 0.192	63	0.934	0.344	167.9

For the correlation with logPO/W using logarithmized eq 1, a is the slope and const is the intercept.

The solvatochromic approach was applied to the logPm/w and logPb/w data using the overall H-bond acidity, A; overall H-bond basicity, B; dipolarity/polarizability, S; excess molar refraction, E; and characteristic volume V, which were obtained using the Absolv software.[118] The estimated and experimental values of the solvatochromic properties A–V are summarized in Table S2 in the Supporting Information. The logPm/w and logPb/w values are expressed as the linear combination of solvatochromic parameters, with the coefficients scaling individual parameters being denoted by the same letters, just using lowercase, as the parameters. The results, summarized in the last two lines of Table 5 demonstrate very good correlations, in accordance with a previous report.[137] Standard deviation of the regression coefficient a for the monolayer (Table 5, line 3) is larger than the coefficient itself. Removal of a would provide similar statistical indices, but it was not performed to allow a comparison with other data sets. The solvatochromic approach is well suited for a direct modeling of overall accumulation in the DMPC liposomes. The higher number of optimized parameters makes its use in the disposition function[138] for QSAR models of biological effects more cumbersome.

The PPM approach was mainly developed for prediction of bilayer interactions with proteins and peptides but was also applied to partitioning of small organic molecules in the bilayer.[79] For 23 neutral and charged compounds, the model explained 65% of variance in the experimental binding free energies. The compounds were collected from several studies,[98,107,142−144] based on the criteria of conformational rigidity and composition limited to the C, H, N, O, and S atoms, for which the parameters were available. The performance may improve with more focus on partitioning of small molecules, distinguishing the interactions, which dominate in the bilayer partitioning of neutral and charged molecules, and expanding the set of parametrized atoms.

The COSMOmic approach[74] was used to predict the logPb/w values for more than 200 neutral compounds from two published sets.[145,146] The model explains 60% and 80% of experimental variance, respectively, and provides leave-one-out cross-validation results at similar levels. The use of MD simulations to obtain representative lipid geometries improved the performance to 81% to 92% of explained variance for the sets of 42–66 neutral compounds.[76,77]

Atomistic depth-constrained MD simulations with several force fields were compared for their ability to predict partition coefficients of 11 compounds.[96] The explained variance varies from 72% to 93% for individual force fields.

Applications of the Intrabilayer Distribution Model

The comparisons with other approaches show that our model predicts both the overall partition coefficient and the intrabilayer distribution. Other approaches focus mostly on the former aspect, and if they handle the distribution, they do not predict it well. The intrabilayer distribution is an important phenomenon, which plays key roles in drug action. It strongly affects at least two phenomena playing significant roles in drug development: interaction of compounds with the membrane proteins, which have binding sites only accessible from inside the bilayer and the trans-bilayer transport.

The connection intrabilayer distribution and membrane interactions are straightforward: only the fraction of the compound present at a proper bilayer depth is able to readily interact with the intrabilayer binding site. We expect that quantitative predictions of intrabilayer distribution will have a profound impact on modeling of drug interactions with the membrane proteins having the binding sites accessible from the bilayer, such as P-glycoprotein,[147] cytochromes P450,[148] cyclooxygenases,[149] and other targets.

The connection with the transport rates is less obvious, albeit logical. We hypothesized that a high transport rate requires a comparatively even intrabilayer distribution, in addition to sufficient overall bilayer accumulation. Too strong or too weak interactions with either H- or C-stratum or the H/C-interface create the bottlenecks, which slow down the transport.[3,33,150] Similar considerations for headgroups and core were used to justify the parabolic dependence of transport rates on the reference partition coefficients.[151] The hypothesis is also supported by the description of the steady-state permeability coefficient: its inverse value is equal to the sum of resistances in individual strata.[152,153] Each resistance is expressed as the ratio of the thickness and the product of the partition coefficient and the diffusion coefficient, all for the given stratum. As the drug–stratum interactions strengthen, the partition and diffusion coefficients increase and decrease, respectively. Consequently, optimum interaction strength could be expected as a condition for fast trans-bilayer diffusion. A slow-down may arise from a combination of a high partition coefficient and low diffusion coefficient or vice versa. Although preliminary data support our hypothesis,[3] a rigorous validation remains to be performed.

Trans-bilayer diffusion is a key factor for the rates of drug absorption and distribution in the body. While the importance of the absorption rate for oral bioavailability is well recognized, the role of distribution rate is unduly underappreciated. Most drugs are developed for oral administration, which often goes hand-in-hand with fast, whole-body distribution. However, for drugs with significant side effects, a toxicity reduction could be achieved by distribution that is restricted to the target body region, e.g., tumors for anticancer drugs. A limited distribution can be achieved by drugs with both targeted administration and slow trans-bilayer transport. While targeted administration is an intensively studied issue, the need for limited distribution is not generally recognized.

Conclusions

The calibrated model, logarithmized eq 10 combined with eq 11, with the optimized coefficient values given in the last line of Table 3, provides a prediction of both overall partitioning and drug distribution in individual strata of the DMPC bilayer for the compounds, for which the input data, the C16/W and the C16/DAcPC partition coefficients, and fragment solvation characteristics (Table 2) are available. The model was created for neutral compounds or compounds having negligible ionization under physiological conditions. The interfacial poses are generated using the algorithm described in Figure 2 and the accompanying text. The fragment-based partition coefficients of individual interfacial poses, Pij, are calculated (eq 11) using the fragment solvation characteristics f in the DAcPC/W and C16/W systems for the fragments embedded in the headgroups and the core, respectively. The interfacial poses exhibit the Boltzmann distribution, with the fraction of the jth orientation being equal to Pij/ΣPij.

The presented model is robust and stable, and has good predictive ability. The use of the established ClogP fragmentation scheme[35] allows the extension of these predictions for compounds, for which the ClogP fragment solvation parameters in the C16/DAcPC[55] and C16/W[97] are available. To the best of our knowledge, this is currently the only approach that allows reliable predictions of both overall partitioning and intrabilayer distribution for new compounds with available input data.

Overall, accumulation alone was also described well by other approaches: the O/W partition coefficients, the solvatochromic approach, COSMOmic combined with MD simulations, and atomistic MD simulations with some specialized force fields.

Outlook

For the studied set of neutral molecules, the calibrated model showed good descriptivity, predictivity, and stability. In the future, if needed, the model can be expanded by including, e.g., electrostatics, molecular cross-section determining the disturbance of the most dense interface region, and entropy changes. The last aspect may be of lower significance thanks to the enthalpy–entropy compensation, which is a common phenomenon in partitioning.[42−44] Inclusion of electrostatic interactions, which may significantly contribute to bilayer partitioning of charged molecules,[43,143,154] is probably necessary because most drug molecules ionize under physiologic conditions.[155]

The model has been calibrated for the DMPC bilayer and needs to be expanded to other bilayer compositions. Extension to the bilayers consisting of PC with other saturated fatty acids should be accomplished by the use of adequate volume fractions for the headgroups and the core (eq 10). We assume that for bilayers composed of PC containing unsaturated fatty acids, the model could be directly recalibrated, with the main change expected in the Collander α and β coefficients for the core. This expectation is based on good Collander-type correlations for the partition coefficients of compounds measured in saturated and unsaturated hydrocarbons[156] and similar coefficients of the solvatochromic equation for the two-phase systems of n-hexadecane,[55] hexadec-1-ene,[157] and deca-1,9-diene[157,158] with water. The situation may be more complex, though, as indicated by experimental data for α-tocopherol showing the preferred location varying with PC fatty acid chains.[17] The effects of admixtures of other phospholipids and cholesterol to PC as the main phospholipid in mammalian membranes on bilayer distribution of compounds will require more complex model modifications.