Partitioning into Colloidal Structures of Fasted State Intestinal Fluid Studied by Molecular Dynamics Simulations.

We performed molecular dynamics (MD) simulations to obtain insights into the structure and molecular interactions of colloidal structures present in fasted state intestinal fluid. Drug partitioning and interaction were studied with a mixed system of the bile salt taurocholate (TCH) and 1,2-dilinoleoyl-sn-glycero-3-phosphocholine (DLiPC). Spontaneous aggregation of TCH and DLiPC from unconstrained MD simulations at the united-atom level using the Berger/Gromos54A7 force fields demonstrated that intermolecular hydrogen bonding between TCH molecules was an important factor in determining the overall TCH and DLiPC configuration. In bilayered systems, these intermolecular hydrogen bonds resulted in embedded transmembrane TCH clusters. Free energy simulations using the umbrella sampling technique revealed that the stability of these transmembrane TCH clusters was superior when they consisted of 3 or 4 TCH per bilayer leaflet. All-atom simulations using the Slipids/GAFF force fields showed that the TCH embedded in the bilayer decreased the energy barrier to penetrate the bilayer (ΔGpen) for water, ethanol, and carbamazepine, but not for the more lipophilic felodipine and danazol. This suggests that diffusion of hydrophilic to moderately lipophilic molecules through the bilayer is facilitated by the embedded TCH molecules. However, the effect of embedded TCH on the overall lipid/water partitioning was significant for danazol, indicating that the incorporation of TCH plays a crucial role for the partitioning of lipophilic solutes into e.g. lipidic vesicles existing in fasted state intestinal fluids. To conclude, the MD simulations revealed important intermolecular interactions in lipidic bilayers, both between the bile components themselves and with the drug molecules.

Introduction

The hunt for highly potent new drug molecules, often targeting receptors with lipophilic molecular requirements, has resulted in a general trend of increasing number of highly lipophilic, poorly water-soluble drugs.[1−3] Between 75 and 90% of all new drug molecules have a solubility which is too low to allow complete dissolution of the dose in the intestinal fluid after oral administration.[4,5] Therefore, the bioavailability of the drug is compromised due to low and erratic absorption. In the gastrointestinal tract lipophilic drug molecules are solubilized in mixed lipid aggregates in the intestinal fluids. In the fasted state, these are composed of bile salts and phospholipids. For lipophilic drug molecules, the interaction of the drug with the bile salts and phospholipids can have a significant effect on the partitioning into these structures and, hence, impact the solubilization and, consequently, the bioavailability. The composition of the intestinal fluid shows high interindividual variability in the fasted state,[6,7] and the large variation in bile secretion and composition adds another unpredictable element to drug delivery. Unpredictable bioavailability can be critical when the drug has a limited therapeutic window; i.e., a low dose may be ineffective while a high dose results in toxicity.

The two primary components in simulated human intestinal fluid are the bile salt taurocholate (TCH) and the phospholipid phosphatidylcholine 1,2-dilinoleoyl-sn-glycero-3-phosphocholine (DLiPC).[8,9] Phospholipids have polar headgroups and long (nonpolar) acyl chains (Figure ). They therefore tend to self-assemble into ordered and extended bilayered membranes and micelles in water to minimize the exposure of their nonpolar groups to water. Bile salts are amphiphilic surfactants that easily self-assemble into small micelles in water solutions at millimolar concentrations,[10−16] with the capacity to solubilize both polar and nonpolar compounds.[17,18] Bile salts originate from cholesterol in the liver and contain a hydrophilic side chain attached to a large, rigid, and boat-shaped tetracyclic ring system. The ring system is characterized by a hydrophilic and a hydrophobic face, where the overall amphiphilic properties stem from the hydrophilic hydroxyl and hydrophobic methyl groups located on different sides of the ring system (see Figure ).[19]

Figure 1

United-atom molecular dynamics snapshots and Lewis structures illustrating the deprotonated bile salt taurocholate (TCH, top) and the phospholipid 1,2-dilinoleoyl-sn-glycero-3-phosphocholine (DLiPC, bottom). Note the stereochemistry of the TCH hydroxyl groups.

Experimental characterization and drug dissolution studies show that mixtures of DLiPC and TCH produce different mixed lipid aggregates that may increase the solubility of certain drugs by several tenfold.[20] The conditions of formation of these mixed lipid aggregates and their colloidal structure (micellar, uni-multilamellar vesicles) are hence of fundamental importance for the solubilization of drug molecules in the intestine. The aggregate size and the type of the structures formed in the three-component systems of water/bile salt/phospholipid differ significantly compared to two-component systems (water/bile or water/phospholipid) that typically form micelles or vesicles.[11,21−25] In this study, we used MD simulations to study drug partitioning into colloidal structures typically found in fasted intestinal fluid. As the first step, we investigated spontaneous aggregation of mixed TCH and DLiPC systems and simulated a range of systems with different water/TCH/DLiPC ratios on the united-atom (UA) level, starting from either random or ordered bilayer structures. Then, to quantify and further investigate the mixed lipid aggregates, we performed simulations of embedded TCH clusters in DLiPC bilayers using free energy calculations via umbrella sampling (US) simulations. Finally, we investigated the free energy barriers of the interaction of selected model drug compounds with representative DLiPC bilayer membranes, either without or in the presence of embedded TCH at a TCH:DLiPC ratio of 1:2, using US simulations on the all-atom level.

Results and Discussion

Interaction between Taurocholate and Phospholipids

Initially, several different attempts to probe the tendency for TCH and DLiPC molecules to spontaneously aggregate into ordered bilayers were performed on the united-atom level using the Berger/Gromos54A7 force fields by simulating systems having random initial configurations and different TCH and DLiPC compositions and sizes. These initial simulations were performed in view of the findings by Marrink et al., who showed that several phospholipids (1,2-dipalmitoyl-sn-glycero-3-phosphocholine (DPPC), 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC), 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC), and 1,2-dioleoyl-sn-glycero-3-phosphoethanolamine (DOPE)) spontaneously aggregate into an equilibrated bilayer phase in water within a few tens of nanoseconds, using random starting configurations.[26] In contrast to the reported systems and simulations, only a handful of the simulated DLiPC and mixed DLiPC/TCH systems performed in this study from random starting configurations yielded well-defined lipid aggregates after more than 100 ns of simulations. As for DLiPC, this can be explained by the polyunsaturated acyl chains yielding a low gel to liquid phase transition temperature (Tm) of −57 °C.[27] For comparison, this is approximately 100 °C lower than for the saturated 16- and 18-carbon analogues, DPPC and DSPC (Tm = 41 and 55 °C, respectively[28]). With this in mind, it appeared that the full equilibration of the self-assembly between mixed DLiPC/TCH systems may require much longer time scales than those used in this study and that by Marrink et al. In line with previous studies, however,[29] the obtained lipid structures were found to be highly dependent on the properties of the periodic boundary conditions (PBC) and pressure control. For instance, 64 DLiPC molecules (with or without 16 TCH molecules saturated with 55 water molecules per lipid) did form a bilayered phase when semi-isotropic pressure control was used. However, cylindrical micelles were formed with isotropic pressure control. It is interesting to note that both types of structures can be found experimentally and are known to coexist under certain conditions.[30−33]

Given that fully equilibrated self-assemblies of the DLiPC/TCH lipid systems typically appeared to be outside the time scale of our simulations and likely subjected to PBC artifacts, the stability of preformed bilayer structures of DLiPC with TCH molecules were also investigated for six systems having 128 DLiPC molecules and TCH:DLiPC ratios of 0:1, 1:4, 1:2, 1:1, 2:1, and 4:1, respectively, containing 55 water molecules per lipid (i.e., fully water saturated conditions). The resulting bilayer structures displayed remarkably high stability on the time scale of the simulations. Figure shows the density profiles of mixed systems during 100 ns. Interestingly, even at a 1:1 ratio, the DLiPC sustained a seemingly normal bilayer structure with the polar headgroups extending toward the aqueous phase, well beyond the TCH polar side chain. This was further confirmed by the similar positions of the DLiPC main functional groups (choline, phosphate, glycerol) in the system with a 1:2 TCH:DLiPC ratio (Figure , center) as compared with the pure DLiPC system (not shown). As for TCH (Figure , bottom), most of the sterol groups are buried deep into the DLiPC acyl chains in the bilayer interior, and the polar side chains face the aqueous phase, although not beyond the DLiPC polar headgroups. Interestingly, with increasing TCH concentration in the simulated systems, increasing amounts of water was also found within the bilayer interior, increasing from 0 wt % to 0.2, 1, 1.7, 3.3, and 12 wt % with increasing TCH:DLiPC ratio.

Figure 2

Top: density profiles of mixed TCH/DLiPC united-atom systems containing 0–512 TCH (green) and 128 DLiPC (black) molecules at ratios of 0:1, 1:4; 1:2, 1:1, 2:1, and 4:1 TCH to DLiPC hydrated with 55 water molecules per lipid. Arrows indicate the effect of increasing TCH concentration. Dashed lines show the pure DLiPC system. Middle and bottom: density profiles of the 32:128 TCH:DLiPC system, displaying the relative positions of the molecular moieties of DLiPC (middle) and TCH (bottom). For clarity, the dotted lines show the total density profiles for DLiPC and TCH. Add the following at the end: (Abbreviations in the middle panel refer to the nitrogen group (N grp), the phosphate group (P grp) and the glycerol linker (Glyc grp).

Stability of Bilayers

The apparent stability of the mixed TCH:DLiPC bilayers under the given conditions might be due to the ability of the TCH molecules to organize into transmembrane clusters in the normal direction within the DLiPC bilayers (Figure ). At low TCH:DLiPC ratios, most TCH molecules were found slightly below the DLiPC headgroups; the nonpolar groups faced the bilayer interior, and the polar groups interacted with the DLiPC headgroups and water molecules. However, with increasing TCH concentrations within the bilayer, the number of transmembrane clusters increased, extending throughout the bilayer in the normal direction. This behavior can be explained by the amphiphilic nature of TCH, which maximizes the hydrophilic and hydrophobic interactions of the polar and nonpolar groups simultaneously.

Figure 3

Left: snapshot of a mixed TCH:DLiPC united-atom bilayer system. The equilibrated structure from a bilayered system contains 32 deprotonated TCH and 128 DLiPC molecules. Water molecules are omitted for clarity; charge-balancing Na+ ions are shown in blue. Right: an isolated 3/4 TCH cluster having three and four TCH molecules in the top and bottom bilayer leaflet, respectively. The molecules are stabilized by hydrogen bonds between the TCH molecules within the same bilayer leaflet and between TCH molecules in the opposing leaflet. Note the single water molecule in the upper TCH cluster; it penetrates deep into the bilayer interior due to the interaction with the buried hydrogen donors and acceptors of the TCH. PDB files with these structures are available as Supporting Information.

To study the stability of the transmembrane TCH clusters as a function of number of TCH molecules per cluster, a set of simulations were performed with systems of 64 DLiPC lipids containing a single isolated TCH cluster. One to four TCH molecules were kept in the upper leaflet and four in the lower one. The TCH transmembrane clusters in systems with 1 and 2 TCH in the upper leaflet were not stable during the timespan of the simulation, and the cluster disaggregated. These systems had a lower probability of transmembrane hydrogen bonding, in contrast to the systems with 3 or 4 TCH upper leaflet TCH molecules; hence, the TCH clusters were not stabilized. With only 1 or 2 TCH molecules in the upper leaflet, the initial hydrogen bonds with the TCH in the opposing bilayer leaflet were disrupted, and the TCH molecules positioned themselves adjacent to and below the DLiPC headgroups after a few tens of nanoseconds. However, the systems with 3 or 4 TCH in the upper leaflet remained virtually unchanged even after 200 ns, indicating superior stability over the smaller clusters. Hydrogen bond analysis (with the maximum bond length and angle between the hydrogen bond donor and acceptor being 0.35 nm and 30°) further revealed a slightly larger number of H-bonds per TCH molecule between both TCH/TCH and TCH/DLiPC for clusters with 3 TCH compared to 4 TCH molecules, indicating greater stability of the former type (Figure ).

Figure 4

Hydrogen bond participation per united-atom TCH in transmembrane clusters containing 3 (blue) and 4 TCH (red) molecules per bilayer leaflet. The number of inter-TCH hydrogen bonds (approximately 3.5) is highest in clusters with 3 TCH molecules; in general, this number of TCH molecules also forms more hydrogen bonds with its surroundings.

The TCH clusters with 3 TCH molecules in the same bilayer leaflet were often slightly tilted relative the normal direction of the bilayer, with the polar −OH groups in the sterol moiety facing each other. The clusters with 4 TCH in the same bilayer leaflet typically had one peripheral and partly excluded TCH, or one to two TCHs slightly twisted around their principal axis, exposing −OH groups toward the DLiPC acyl chains. Analysis of the external surface area of the entire TCH transmembrane clusters with a 0.14 nm probe (i.e., the equivalent of the water molecule radii) was also performed on the hydrophobic and hydrophilic parts of those clusters with atom charges less or greater than ±0.2. This analysis revealed that the 3 TCH clusters had an overall 15% larger hydrophobic area and 4% larger hydrophilic area per molecule than the clusters with 4 TCH. However, analysis of the external surface area of the sterol moieties only (without the polar side chain) resulted in a 12% increase and 10% decrease in the hydrophobic and hydrophilic areas, respectively. This shows that the sterol groups of the 3 TCH cluster display a more optimal packing geometry when embedded in the DLiPC bilayer than the corresponding 4 TCH cluster, by exposing more hydrophobic and less hydrophilic atoms toward the lipid bilayer interior.

Spontaneous formation of the transmembrane TCH clusters was also investigated, i.e., TCH cluster formation independent of the possible artifacts introduced by using ordered starting configurations. For this, we performed simulations using a pre-equilibrated micelle composed of 8 TCH embedded in the center of 64 DLiPC bilayer. Upon initial equilibration, predominately pairs of TCH molecules interacting with their hydrophilic groups formed within 1–2 ns, with some orienting the headgroup toward the center of the bilayer. However, during the production run (500 ns), a 3 + 3 transmembrane TCH cluster formed within a few nanoseconds with the two residual TCH molecules positioned adjacent to the DLiPC headgroups on each side of the bilayer. This result further indicated that TCH transmembrane clusters could form naturally in bilayers of DLiPC.

To further investigate the TCH and DLiPC interactions in mixed lipid bilayers, the potential of mean force (PMF) for the TCH molecule perpendicular the bilayer (z-direction) was determined using US simulations followed by the weighted histogram analysis.[34] A PMF indicates the affinity of a TCH molecule for a particular position. The overall PMF profile also reveals the free energy barriers along the chosen reaction coordinate. Figure shows the PMF profiles for TCH molecules in different transmembrane clusters over half the bilayer. For a TCH molecule not interacting with any other TCH molecules (denoted 1/0 in Figure ), the PMF profile is even positive (+8 kJ/mol) in the center of the bilayer. This clearly shows that the interior of the bilayer is an unfavorable environment for the TCH alone. However, with 4 TCH molecules in the opposing bilayer leaflet (denoted as 1/4 in the figure), the corresponding value decreased to −68 kJ/mol, indicating a gain in stability by formation of a transmembrane cluster. It is also notable that the minimum of the PMF profile (the optimal position of the molecule COM) of the 1/4 system was more shifted toward the bilayer center and away from the DLiPC headgroups than the 1/0 system. With additional TCH molecules present in the same bilayer leaflet (i.e., systems with 2/4 and 3/4 TCH molecules in the upper and lower bilayer leaflets, respectively), there were even deeper PMF profiles, demonstrating a greater overall stability of the TCH clusters. For clusters with up to 3 TCH molecules in the upper bilayer leaflet, the minimum of each PMF curve shifted away from the bilayer center due to the increasing interactions with adjacent TCH. For the same systems, this trend was accompanied by an increasing free energy penetration barrier, Gpen (defined analogously to eq ), as indicated in Figure . However, the same trend was broken upon additional TCH, suggesting that clusters with >3 TCH are less rigid due to the small energy gain achieved by adding an additional TCH to the otherwise stable and rigid 3 TCH-molecule clusters.

Figure 5

Density and PMF profiles for TCH/DLiPC united-atom systems. Top: simulation snapshot showing the system setup for the US simulations with black arrows indicating the approximate position of the on average 40 umbrella windows used. Middle: density profiles of DLiPC and the TCH cluster, with five separate PMF curves for TCH along the Z-direction, normal to the DLiPC bilayer. Note the shift in the position of the PMF minimum with a changing number of TCH molecules. Black arrows indicate ΔGlipid/water and ΔGpen for the 3/4 cluster. Bottom: density of DLiPC and TCH and the PMF profiles for a 3/4 TCH cluster without and with additional (xTCH) of surface TCH.

To better represent physiological conditions where free TCH molecules are also present, an additional system was created by adding 64 TCH to the aqueous bulk phase. After equilibration, the aqueous phase TCH mostly accumulated at the DLiPC/water interface, which drastically changed the PMF profile of the probe TCH molecule (Figure , bottom). Nevertheless, the minimum value and position of the PMF profile strongly resembled systems without bulk phase TCH. Hence, TCH was energetically favored to interact with the transmembrane TCH clusters in the lipid bilayer rather than interacting with the TCH-rich water phase in contact with the DLiPC bilayer.

Drug Interactions with Bile Salt and Phospholipids

Bile salts and phospholipids such as TCH and DLiPC solubilize drugs[35] and alter lipid membrane properties.[19] Therefore, additional US simulations were performed to investigate the interactions of drug molecules with the mixed lipid bilayers of DLiPC and TCH (Figure ). The simulations were performed at the all-atom level using the Slipids and GAFF force fields, as recommended by a recent benchmarking study,[36] and to better account for specific small molecule–lipid interactions. Because of the computational expense of all-atom US simulations for three-component systems,[37−39] the simulations were limited to five different molecules: water, ethanol, carbamazepine, felodipine, and danazol (Figure ). These molecules were chosen to cover a wide range of lipophilicity, and the resulting log P values, as obtained from eq , are summarized in Table .

Figure 6

Density profiles of a mixed lipid bilayer containing 48 DLiPC and 24 TCH molecules.

Figure 7

Left: molecular structures of ethanol (top left), carbamazepine (top right), felodipine (middle), and danazol. Right: symmetrized PMF profiles of water (red), ethanol (green), carbamazepine (blue, denoted Car), felodipine (purple, denoted Fel), and danazol (black, denoted Dan) over half an all-atom bilayer (centered at Z = 0) containing only DLiPC or (dotted line) mixed TCH/DLiPC (solid line).

Table 1

The log P Values Calculated from All-Atom Simulations with Slipids/GAFF Using Eq

	TCH/DLiPCa	DLiPCb	expc
water	–3.1 ± 0.4	–4.6 ± 0.1
ethanol	1.1 ± 0.1	0.75 ± 0.05	–0.3d
carbamazepine	1.9 ± 0.2	2.3 ± 0.2	2.32e
felodipine	4.2 ± 0.3	5.3 ± 0.2	5.58f
danazol	8.1 ± 0.5	6.8 ± 0.3	4.53f

log Plipid/water in mixed TCH/DLiPC bilayer.

log Plipid/water in pure DLiPC bilayer.

Experimental log Poct.

Reference (40).

Reference (41).

Reference (20).

For water, there was a clear decrease in the otherwise positive PMF profile in the mixed TCH/DLiPC system compared to the pure DLiPC bilayer when investigating deep partitioning into the membrane (>1 nm). Inspection of the trajectories revealed—in line with the results in Figure and the united-atom simulations—occasional penetration of water molecules into the transmembrane TCH clusters, facilitated by hydrogen bonding with the TCH sterol hydroxyl groups. The overall free energy difference (ΔGlipid/water) of water in the mixed TCH/DLiPC and the DLiPC reference was 18.6 ± 2.3 and 27.6 ± 0.6 kJ/mol, respectively, well in line with previous studies.[42] This demonstrates that the water partitioning of a DLiPC bilayer increases by a factor of 30 when TCH is present (calculated from ΔΔGlipid/water and eq ).

The PMF profiles over the mixed TCH/DLiPC and pure DLiPC bilayer revealed only small differences in partitioning of the drug molecules and ethanol (Table ).

However, intricate details governing the molecular interactions between the different compounds and TCH/DLiPC were found by interaction energy analysis (separating out electrostatic and van der Waals interactions). For the lipophilic compound danazol (log Poct of 4.53),[20] the corresponding log Plipid/water value was higher in the mixed TCH/DLiPC and the pure DLiPC systems, respectively, than that observed for octanol. The increased partitioning of danazol in the mixed TCH/DLiPC bilayer and the shift of the minimum value in the PMF profile toward the center of the bilayer are explainable by the van der Waals interactions with the embedded TCH molecules; these were much stronger than the electrostatic interactions. This was further indicated by the higher total density of atoms in the mixed bilayer than in the pure bilayer (Figure ). Despite these findings, the simulated trajectories did not point to greater extent of binding between danazol and TCH as compared to danazol and DLiPC.

For felodipine (log Poct of 5.58),[20] the corresponding log Plipid/water values were lower in the mixed TCH/DLiPC than the pure DLiPC system. The decrease in partitioning of felodipine in the mixed TCH/DLiPC bilayer resulted from both weaker van der Waals and electrostatic interactions with DLiPC. Unlike in the danazol system, these were not compensated for by the interactions with the embedded TCH. Figure S1 of the Supporting Information shows the overlap of the umbrella sampling windows, as well as the effect of the symmetrization procedure on the resulting PMF.

The differences in the PMF profiles for the TCH/DLiPC and pure DLiPC bilayers were smaller for the more hydrophilic carbamazepine and ethanol (log Poct values of 2.32 and −0.30, respectively[40,41]). An inspection of the respective trajectories identified that these molecules interacted with the TCH −OH groups in the center region, which was highly accessible due to low particle density in the bilayer center. Although the differences in the PMF profiles were small, the embedded TCH may accelerate drug diffusion over the lipid bilayer by lowering the ΔGpen. The resulting log Plipid/water for carbamazepine was slightly lower in the mixed TCH/DLiPC compared to the pure DLiPC system, which was equivalent to the reported log Poct value. The corresponding log Plipid/water values for ethanol were, on the other hand, slightly higher in the mixed TCH/DLiPC compared to the pure DLiPC system; i.e., there was a slightly increased partitioning of ethanol when TCH was embedded.

Conclusions

The MD simulations strongly suggest that in mixtures of TCH and DLiPC, reflecting the composition of fasted state simulated intestinal fluid, embedded transmembrane TCH clusters constitute an integral component in phospholipid bilayers composed of DLiPC. This new finding complements existing models of mixed bile–phospholipid micelle and bilayer systems[21,23] by revealing the preferred number and configuration of bile salt molecules embedded in the phospholipid bilayers. Both unconstrained and US free energy united-atom simulations showed that the transmembrane TCH clusters were (i) stabilized by mutual hydrogen bonds (OH3α···OH3α) across the bilayer center and (ii) preferably consisted of 3 or 4 (3 + 1) adjacent and clustered TCH molecules per bilayer leaflet. In the latter type, the additional TCH constituted a more peripheral member to the overall cluster.

From simulations on the all-atom level, the effect of embedded TCH in DLiPC bilayers on the partitioning of five model substances (water, ethanol, carbamazepine, felodipine, and danazol) showed that only water displayed a positive free energy profile over the entire bilayer regardless of embedded TCH. For water this positive free energy barrier was however reduced in the presence of TCH due to penetration into the transmembrane clusters. For the other hydrophilic solutes—ethanol and carbamazepine—the overall lipid/water partitioning into the lipid bilayer was similar to and without embedded TCH. The effect of embedded TCH on the overall lipid/water partitioning (i.e., from ΔGlipid/water) was most significant for danazol, demonstrating that the embedded TCH may facilitate partitioning into, and hence, solubilization of lipophilic solutes in the mixed bilayer. However, this was not observed for the other highly lipophilic drug felodipine, which displayed a slightly reduced log Plipid/water in the presence of embedded TCH.

Methods

Simulation Setup and Parameters

The GROMACS[43,44] suite of programs was used for all simulations. Unless otherwise stated, all lipid simulations were conducted at a physiological NaCl concentration of approximately 0.15 M. For the united-atom simulations (SPC/Gromos54A7/Berger[45,46]) and the all-atom simulations (Tip3P/GAFF/Slipids[47]) we used a universal short-range cutoff at 1.2 nm for both electrostatic and Lennard-Jones interactions. Long-range electrostatic interactions were treated by particle mesh Ewald (PME[48]). A long-range dispersion correction was applied to both energy and pressure.[49] All covalent bond lengths were constrained using LINCS,[50] and a time step of 2 fs was used in all production simulations.

Simulation systems of DLiPC with and without TCH were constructed to have either random starting configurations or ordered bilayers. The latter type was obtained by initially distributing the lipid monomers on a square grid and pre-equilibrating bilayered systems for >50 ns while weakly restraining the lipid movements in the direction normal to the bilayer.

DLiPC and TCH Parametrization

To model the mixed lipid systems at the united-atom level, Berger force field parameters for DLiPC were adapted from a DOPC topology, by enforcing a cis-conformation of the double bonds in the acyl chains.[45] For TCH, the topology and all interaction parameters were obtained from the Automated Topology Builder server using the Gromos54A7 force field.[51] However, the partial charges were calculated from restricted electrostatic potential fitting using the pyRED server as OPLS compatible charges (RESP-O1), i.e., by Gaussian09 calculations at the HF/6-31G* level followed by RESP.[52,53]

Similarly, a Slipids compatible topology for DLiPC was adapted from a DOPC topology.[54] This was done by removal of two hydrogen atoms and by enforcing a cis-conformation of the two neighboring double bonds and by increasing the partial charges of the CH2 atoms in the adjacent methylene group by +0.03 (to +0.06) to maintain charge neutrality. For the all-atom TCH and drug molecules, the general Amber force field (GAFF) was used, having molecular topologies constructed by Antechamber software through the acpype.py script (available online: http://code.google.com/p/acpype/) applied to molecular coordinates obtained from the ZINC database.[55]

All regular and unconstrained MD productions runs were typically preceded by (i) energy minimization using the steepest descent algorithm with a tolerance of 1000 kJ mol–1 nm–1, followed by (ii) a 50 ps simulation in the canonical (NVT) ensemble using position restraints on all solute particles, followed by (iii) a 5 ns volume equilibration simulation in the isothermal–isobaric (NPT) ensemble. The latter simulations used the Berendsen[56] barostat for semi-isotropic pressure control, whereas all production runs used the Parrinello–Rahman[57] barostat with semi-isotropic pressure coupling. For temperature control, the modified Berendsen and velocity-rescaled thermostats were used for equilibration and production runs, coupling the non-water and water molecules to separate thermostats. By analogy, removal of center of mass (COM) translation was applied to the lipids and water molecules independently; for bilayer simulations, the COM of the upper and lower bilayer leaflets were controlled separately.

Pure DLiPC and mixed TCH:DLiPC US systems—with either a random or an ordered bilayered initial configuration—were simulated for a minimum of 100 ns, with either 64, 128, or 256 DLiPC molecules and with a TCH:DLiPC ratio from 0 to 4:1. To ensure full water saturation, the systems contained a minimum of 30 water molecules per lipid, in contrast to 21.5, the reported experimental value for DLiPC.[58]

In the following studies we focused on free energy calculations making use of the bilayers since it is well-known that unilamellar vesicles are naturally present in fasted intestinal fluids.[11] Free energy profiles for TCH molecules over the bilayers, i.e., the potential of mean force (PMF), were computed from US simulations. To generate the starting configurations for the sampling, a single TCH molecule was pulled to the aqueous bulk phase from the center of a representative TCH cluster embedded in a bilayer consisting of 64 DLiPC molecules. The systems contained 0–4 adjacent TCH molecules within the same bilayer leaflet and 4 TCH in the opposing bilayer leaflet. The total number of lipids in each leaflet was set constant by removal of DLiPC molecules from the opposing bilayer leaflet. Each starting structure was equilibrated with 90 water molecules per lipid, after which a TCH molecule was pulled by a harmonic force at a rate of 0.00025 nm/ps with the Gromacs pull code. For the US simulations, more than 40 configurations were chosen from an initial set of 2500. In these 40 configurations, the probed TCH molecule was displaced approximately 0.15 nm along the reaction coordinate, i.e., in the normal direction (z-axis) to the bilayer. Each configuration was then simulated for 5 ns with a harmonic force constant of 800 kJ mol–1 nm–2. To ensure proper overlap of the sampled umbrella windows, each umbrella simulation was performed twice with different starting configurations. From the resulting potential of mean force (PMF) profiles an equivalent bilayer/water partitioning (log Plipid/water) was calculated aswhere ΔGlipid/water is the free energy of transfer from the aqueous phase to the lipid bilayer, expressed as the difference in the free energy of solubilization, ΔGwater – ΔGlipid, i.e., analogous to calculations of the traditional octanol and water partition coefficient log Poct/water.

The free energy profiles over the bilayers for the small molecules (danazol, felodipine, carbamazepine, water, and ethanol) were similarly computed from US simulations, this time using all-atom systems. To evaluate the effect of TCH in the mixed lipid bilayers, umbrella simulations were also performed using pure DLiPC bilayers without TCH. The mixed lipid bilayer was constructed with 24 TCH and 48 DLiPC molecules, whereas 64 DLiPC molecules were used in the systems without TCH. Similarly to the TCH-only US simulations, representative configurations were chosen with a spacing of 0.15 nm along the pull vector. Each chosen configuration was simulated for 20 ns using a harmonic force constant of 700 kJ mol–1 nm–2. Unlike the TCH US simulations, however, 2 (or 4 for ethanol) equally spaced probed molecules were pulled simultaneously through the bilayer, i.e., sampling both sides of the bilayer leaflet at the same time. This was possible because of the small size and neutral charge of these compounds. To avoid overlap of the inserted molecules, a slow growth protocol over 50 consecutive simulation steps was applied, after which the complete drug/solvent/TCH/DLiPC structure was energy minimized and equilibrated. Because of this, a lower pull rate of 0.000 05 nm/ps was used to avoid structural artifacts in the bilayer caused by the pulled drug molecules.

To obtain the final PMF along the reaction coordinate, i.e., the free energy profiles, all US simulations were analyzed with the weighted histogram analysis method as implemented in the GROMACS g_wham utility. Standard deviations were computed from bootstrapping over 200 runs using a tolerance of 1 × 10–5.[34] In the mixed TCH/DLiPC system, the PMF with calculated standard deviations for water was obtained from the averaged Boltzmann-weighted density profiles of water in the drug molecule simulations, i.e., from unconstrained conditions and a total simulation time of approximately 2 μs.