Molecular Dynamics Simulations on Interindividual Variability of Intestinal Fluids: Impact on Drug Solubilization.

Efficient delivery of oral drugs is dependent on their solubility in human intestinal fluid, a complex and dynamic fluid that contains colloidal structures composed of small molecules. These structures solubilize poorly water-soluble compounds, increasing their apparent solubility, and possibly their bioavailability. In this study, we conducted coarse-grained molecular dynamics simulations with data from duodenal fluid samples previously acquired from five healthy volunteers. In these simulations, we observed the self-assembly of mixed micelles of bile salts, phospholipids, and free fatty acids. The micelles were ellipsoids with a size range of 4-7 nm. Next, we investigated micelle affinities of three model drugs. The affinities in our simulation showed the same trend as literature values for the solubility enhancement of drugs in human intestinal fluids. This type of simulations is useful for studies of events and interactions taking place in the small intestinal fluid.

Introduction

Oral formulations such as tablets and capsules are preferred for delivery of small-molecule drugs because of their cost efficacy, stability, and patient compliance. For the drug to reach systemic circulation, the oral formulation needs to disintegrate and dissolve in the gastrointestinal tract, transfer through the mucus, and permeate the gut wall. Any of these steps can be rate limiting for the drug absorption process, but it is critical that the compound (molecularly) dissolves as only dissolved drug molecules will be absorbed through the intestinal barrier and reach the systemic circulation. Ideally, the water solubility of candidate drug molecules should be evaluated early in the formulation pipeline.

Typically, solubility is tested by in vitro dissolution to give an idea of what type of administration routes are possible and how much formulation effort is required to design the final product. However, since orally administered drugs have their main absorption in the small intestine, the solvent of greatest relevance for the solubility assessment is the human intestinal fluid (HIF). HIF is composed of a mixture of water and bile secreted by the bile ducts. There can be large differences in solubility in water compared to that observed in HIF because the latter has a far more complex composition. Understanding the solubilizing ability of HIF is therefore of great interest for researchers in drug development and formulation.

The small molecular components have been thoroughly quantified in HIF [1-3], and commercial simulated intestinal fluids (SIFs) are available for in vitro experiments. HIF changes dynamically over time as lipid digestion, water and bile secretion, and absorption of components occur simultaneously. The composition of HIF is also heavily dependent on the prandial state of the individual. While bile salts, phospholipids, and free fatty acids are available in the fasted state, higher concentrations of them are present in the fed state, in addition to glycerides of different forms, depending on the type of food consumed.[1] The higher apparent solubility of many compounds in HIF compared to water is linked to the colloidal structures—micelles and vesicles of different shapes and sizes—formed by the HIF components.[2−4] Techniques to characterize HIFs and SIFs include microscopy (cryogenic transmission electron microscopy and atomic force microscopy) and scattering (dynamic light scattering, small-angle X-ray scattering, small-angle neutron scattering).[5−7]

How these colloidal structures interact with drug molecules and excipients of drug formulations is key to accurately estimate drug solubility in a patient population. This information identifies formulation strategies that are likely to be successful for specific drug candidates. Molecular dynamics (MD) simulations are an alternative to experimental assessment of these molecular interactions. MD simulations are based upon physics equations that describe the movement of molecules within a defined system. MD simulations have long been used in studies of amphiphilic aggregating systems to gain insights into the morphology and assembling of micelles[8,9] and bilayers.[10] They have also been used to investigate drug partitioning in colloidal systems.[11,12]

Bile components are one group of molecules that have been simulated because of their self-assembling properties.[13,14] For example, Mark and Marrink simulated the self-assembly of cholate and palmitoyloleoyl phosphatidylcholine,[15] which gave insights into the morphology of bile salt micelles and mixed micelles in the gastrointestinal tract. Birru et al. explored aggregation behavior as a function of bile salt and phospholipids for different degrees of digestion of fatty acids. They observed micelles, wormlike micelles, and phase separations.[16] All of the simulations above were conducted at atom- or untied atomic resolution, simulations that limit accessible length- and timescales due to the high cost of computational power needed. Therefore, the system size is often small. It is difficult to simulate physiologically relevant systems with a sufficient number of molecules to assemble the multitude of colloidal structures present in, e.g., intestinal fluids.

Depending on the MD method, different resolutions and accuracies are available.[17] To simulate larger systems of bile components, there are lower-resolution techniques such as coarse-grained (CG) molecular dynamics and dissipative particle dynamics (DPD). These reduce the cost of computing resources, but of course give lower resolution than all-atom and united-atom simulations. Both CG and DPD have been used to simulate micellar systems of bile salts, mixed bile salts phospholipid micelles,[18−21] and small phospholipid vesicles.[22]

In this study, our aim was to investigate human duodenal fluid in the fasted state, with a focus on the assembly of the colloidal structures. We used CG MD simulations of previously collected experimental data from the aspirated duodenal fluids of five healthy volunteers. CG MD successfully identified the impact of interindividual variability of HIF components on the colloidal structures, in terms of size, shape, concentration, composition, and solubilization capacity.

Method

Composition of Human Intestinal Fluid in the Fasted State

Simulations were made using data collected by Riethorst et al.[1] for aspirated duodenal fasted-state samples in fasted and fed states of 20 healthy volunteers. Samples contained eight different bile salts and phospholipids, free fatty acids (8–18 carbon atoms), cholesterol, monoglycerides (20–28 carbon atoms), diglycerides (30–42 carbon atoms), and triglycerides (>42 carbon atoms). In the fasted state, all of the mentioned components except glycerides were present. The simulations in this study were based on the concentration data from 5 of the 20 healthy volunteers (HVs). These five were selected to represent a wide range in components; concentrations for each type of component are listed in Table .

Table 1

Concentrations of Bile Components in Aspirated Samples from Five Healthy Volunteers (HVs)a

	HV3	HV6	HV9	HV16	HV20
bile salts (mM)	6.5	1.2	6.3	1.8	3.2
phospholipids (mM)	0.4	1	1.1	0.8	0.5
free fatty acids (mM)	1.5	1.5	3.2	1.1	1.7
cholesterol (mM)	0.1	0.1	0.04	0.1	0.1
bile salts–phospholipids ratio	16.3	1.2	5.7	2.3	6.4

Concentrations were used for building five HIF simulations, representing one HV each.

Setup of the Molecular Dynamics Simulations

Simulations were performed using the Gromacs software version 2016[23] with the CG Martini force field,[24] in which each CG bead typically represents three to four heavy atoms. Topologies for all components in our simulations were taken from the Martini website, except for the bile salts for which the topology input is available in the Supporting Information. Phospholipids were represented by two different structures: 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) and palmitoyl phospatadylcholine (PPC). The free fatty acids were represented with a carbon chain length corresponding to oleic or palmolic acid. To represent bile salts, the four following types were chosen: glycocholate (GC), glycodeoxycholate (GDC), taurocholate (TC), and taurodeoxycholate (TDC) (Figure ). TC and TDC were parameterized in a previous study.[21] In that study, the topologies were based on the Martini cholesterol, with additional modifications that make the bile salts more polar. For this study, we did a similar procedure to produce topologies for GC and the slightly less polar GDC; specific bead types for all bile salts can be seen in Figure S1. We decided to use only two deoxycholated versions of the bile salts in this study: TDC and GDC. The reason for this was that the CG MD methodology does not capture the small chemical difference very well between deoxycholated bile salts (here, glycochenodeoxycholate, glycoursodeoxycholate, taurochenodoxycholate, and tauroursodeoxycholate). Except for the conjugated amino acid, the difference is only the position of a single hydroxyl group. To represent water, the standard Martini water model was used, in which one bead represents four water molecules. Positively charged sodium ions were added to neutralize the system from the negatively charged bile salts and free fatty acids, which were all represented in deprotonated forms.

Figure 1

Components used in the MD simulations of aspirated human intestinal fluid (HIF): two-dimensional (2D) structure and the corresponding CG representation. Bile salt is colored gray; cholesterol, pink; free fatty acids, green; phospholipid head groups, red and yellow; and their tails, green.

The first set of simulations mimicked the five aspirated samples. Molecules, representing the molar concentrations from the five aspirated samples, were randomly distributed in a cubic simulation box. The box size was selected after an early estimate of how the size would influence the colloidal structures in the system. Several cubic boxes were assessed at different box lengths from 25 to 70 nm, with an average HIF concentration. The number of micelles appeared linear with the volume of the boxes between box lengths 40 and 70 nm. To save computational resources, the box length was therefore set to 45 nm.

Simulations were performed with the recommended electrostatic and van der Waals interactions settings for Gromacs, in combination with standard Martini water beads. The v-rescale thermostat was set to 310 K, and the Berendsen isotropic barostat to 1 bar. For all simulations, an energy minimization was applied using the steepest step algorithm for at least 1000 steps. This was followed by equilibration with an increasing time step from 1 to 30 fs for a duration of 12 ns, followed by an additional equilibration at 30 fs for a duration of 12 ns. Production runs were performed at a time step of 30 fs for a total duration of 3 μs. The total time spent per HV system during a production run was roughly 21 000 CPU hours.

Micelle Determination and Shape Factor Analysis

For the assignment of molecules to individual micelles, an in-house Python script was used, which performed a closest-neighbor search. Each molecule was assigned to a micelle (cluster) if the distance between it and any micelle molecule was below 5 nm. The minimum aggregation number (Nagg) was set to four molecules for the cluster to count as a micelle. Micelles were characterized in terms of shape and size. The shape factors of the micelles were determined by dividing the largest and the smallest eigenvalues from the radius of gyration tensor. According to this calculation, a perfect sphere would have a shape factor equal to 1, and ellipsoidal micelles >1. Micelle diameter was measured as the maximum distance between atoms in one micelle (dmax).

Coarse-Grained Representation of Model Drugs

For the micelle–drug interaction simulations, we used three nonionizable model drugs of different solubilities and hydrophobicities (prednisolone, fenofibrate, and probucol) (Figure ). For prednisolone, we used the topology published by Estrada-López et al.[25] The CG topologies for fenofibrate and probucol were built with a combination of the automated tools auto_martini.py[26] and PyCGTOOL.[27] Additional manual editing of force constants, angles, and bonds was done as needed to make the molecules stable for simulations at 30 fs time steps, while keeping angle and bond conformations in agreement with matching all-atom simulations. To evaluate the topology files, fenofibrate was simulated in water, both as all-atom and as CG. Based on the trajectories, the CG fenofibrate rapidly formed “sticky” aggregates in a manner not observed in the all-atom simulations. In addition, preliminary simulations of fenofibrate with HIF micelles resulted in strong self-aggregation of fenofibrate in water and micelles. This indicates that the hydrophobic forces between the fenofibrate beads were too strong.

Figure 2

Molecular structure and CG representation of (a) prednisolone (published by Estrada-López et al.[25]), (b) fenofibrate, and (c) probucol. To inhibit strong self-aggregation, the Lennard-Jones potentials were reduced between beads of the same drug.

This type of behavior was also observed for probucol. Similar observations of exaggerated self-aggregation are reported for saccharides and proteins.[28] One technique to reduce this self-aggregation is to reduce the depth (ε) of the Lennard-Jones potential for specific bead interactions.[29] Using this technique, we reduced self-aggregation in the fenofibrate and probucol by different amounts. Thereafter, concentrations below and above the solubility limit of fenofibrate[30] were simulated. Because the solubility of fenofibrate is very low in water, it was not feasible to simulate such low concentrations, and therefore propanol was substituted as solvent for these studies. This enabled appropriate reduction in hydrophobicity so that strong self-aggregation did not occur below the solubility limit. Probucol was evaluated in the same manner, but with ethanol instead of propanol. (Ethanol was used because the vendor data were available from Cayman Chemical. Snapshots from simulations at different concentrations can be found in Supporting Information Figure S2.)

Micelle Simulations with Model Drugs

The second set of simulations used the micelle data from the first set of simulations to examine the solubilization of poorly water-soluble drugs. For each HV, the largest micelle in the last frame from the first simulation was isolated and placed in a cubic box with a length of 20 nm. The box was filled with water beads, and sodium ions were added to keep the system neutral in charge. Energy minimization and equilibration were applied as described above. To ensure that the micelles were kept intact during the first equilibration steps (when the box pressure fluctuates), absolute position restraints were applied to the micelle molecules. The simulation box containing a single micelle was replicated with the addition of six concentrations (0.8, 1, 1.5, 2, 3, and 4 mM) for each of the three model compounds, resulting in 18 simulation boxes per HV. The new simulation boxes were equilibrated again and simulated for a duration of 3 μs.

The last 190 ns of each simulation were analyzed by calculating the number of contacts between drug–water and drug–micelle entities. Contacts were calculated with gmx mindist using a cutoff of 0.6 nm. Here, we used the number of contacts as a measure of how many interactions a drug molecule would have with either micelle molecules or water, on the assumption that the number of contacts between drug and micelle reflects the extent of drug solubilization. To do this, we divided the number of drug–micelle contacts (contactsd–m) with the number of drug–water ones (contactsd–w); see eq . This value gives an estimate of the affinity of the drug molecules for either micelles or water. Thus, a high number indicates a preference for micelles, hereafter referred to as micellar affinity (AM). Next, we normalized the micellar affinity with the ratio of the volume fraction of micelle and water (ϕm,w); see eq . This provides an estimate of the drug–micelle affinity, unbiased by the volume fraction in the system (KAM). KAM was used to estimate the differences between the larger (45 nm box length) HV systems. To achieve this, we extrapolated to the affinity for the whole simulated HIF systems (AM,45 nm) by dividing KAM with the volume fraction of the larger systems (ϕm,w,45 nm); see eq . This approach provided relative values for the solubilization of the three drugs.

Results and Discussion

Self-Assembly to Micelles

Five systems, representing five different healthy individuals in the fasted state, were simulated for a duration of 3 μs. The simulations resulted in the self-assembly of colloidal structures in the form of ellipsoidal micelles for each HV (Figure a). Each system followed the same trend in the development of micelles during the duration of the simulations. As seen in Figure b, the number of micelles peaked within the first 100 ns. At this point, there were several undeveloped, small micelles with a low aggregation number. The number of micelles then declined over time, as micelles fused and single bile salts transferred from smaller micelles to larger ones. The number of micelles in the system declined during the remaining duration of the simulation, with none or very few dynamic changes in the last 0.5–1 μs. This indicates that the simulations reached equilibrium, with no further (drastic) changes within a reasonable time frame for the simulation. It should be noted that the time frames discussed here were only from the production run, post-equilibration. The micelles started to form already during the equilibration; thus, the starting time points varied in their number of micelles. In the initial configuration before equilibration, there were no micelles present.

Figure 3

(a) Snapshots from the final frame of MD simulations of human intestinal fluid representing concentration from five healthy volunteers. Colloidal structures in terms of micelles can be seen. (b) The number of micelles (N) changes during the simulations. At each time point, (N) has been normalized to the number of micelles in the last frame for each system. (c) Shape factors for each micelle show their tendency to form ellipsoids rather than spheres. (d) Micelle size, described as the maximal diameter.

The total number of monomers (mostly bile salts) decreases more quickly than the number of micelles in the systems, reaching a base value already at 400 ns. However, the monomers vibrate more because the bile salts dynamically shift between being part of micelles or being free monomers during the simulations (Figure S3). The number of monomers that were not part of the micelles seemed dependent on the number of glycocholate (GC) and taurocholate (TC) molecules in the systems since GC and TC were a clear majority of the total monomers present; see Figure S4.

The difference in behavior of bile salts in the simulations is due to the higher polarity of TC/GC. The TC/GC bile salts have an additional SP1 bead, representing the additional hydroxyl group in their molecular structure. Their behavior in the simulation is also reasonable, as their critical micelle concentration (CMC) values are quite different from each other. TC and GC have an up to 5 times higher CMC than their respective deoxycholated forms, although the values change drastically when components such as phospholipids are available in the system.[31] Since there are a lot of the deoxycholated forms in HIF, the difference in aggregation behavior could be relevant to consider, since SIFs often represent bile salts with only TC. The micelle composition and number of free bile salt monomers may be considerably different for SIF and HIF compositions; however, previous experimental studies have shown that the type of bile salt used in SIFs has a minor effect on drug solubility.[32]

Micelle Size, Shape, and Conformation

Even though there was a difference in the HIF compositions of the HVs, there was no clear interindividual variability in micelle shape. The shape factors for micelles in the simulations varied between 1.2 and 1.9 (Figure c), indicating ellipsoidal micelles, where the lower shape factors describe a “more spherical” ellipsoidal shape. The slightly more spherical micelles were in HV6 and HV16; these individuals had intestinal fluids with a lower ratio of bile salts to phospholipids. The micelle sizes, seen in Figure d, ranged between 2.3 and 7.3 nm, with average sizes between 4.4 and 5.7 nm for the five simulated HV systems. The size of the micelles—calculated as the maximum distance between two beads within a micelle—did not significantly differ among the five HVs, despite the great variety in bile salt, phospholipid, and free fatty acid concentrations.

The micelles in the simulations were slightly smaller than what is reported in experimental research, where the colloidal structures are 10–50 nm in diameter in fasted-state HIF.[33−35] In comparison, the size of the larger colloidal structures in our simulations (5.7–7.3 nm) is close to the lower fraction of micelles observed experimentally. This suggests that the CG simulation methodology performs well in describing the smaller-sized fraction of the colloidal structure landscape in fasted HIF. Even though we are using a CG approach, the size of the simulation box size and the time scale of the simulations may still hinder formation of larger structures.

Differences between HVs of largest and smallest micelles are best exemplified by HV3 and HV6. Not only did these systems have the highest and lowest micelle concentrations but also the greatest range of bile salt-to-phospholipid ratios. HV6 had the lowest micellar Nagg of all five HVs, and the simulated intestinal fluid also contained the fewest number of micelles. This could be due to the low bile salt-to-phospholipid ratio and the low concentration of molecules. The lowest aggregation numbers were between 4 and 16 molecules per micelle (Figure S4). Still, the largest micelle in HV6 was in the same range as for micelles in the other HV systems. HV3, in contrast to HV6, had the highest ratio of bile salts to phospholipids, with 10-fold more monomeric GC and TC and almost twice the number of micelles in the last frames. The largest micelle of HV3 was also slightly more elongated. This shape is probably due to the rough stacking of bile salt molecules that occurs at high concentrations of bile salt.

A closer look at the micelles in the simulations showed that the positioning of bile components within a micelle followed a clear trend. The bile salts formed a hollow shell around a core of phospholipid tails, with their head groups pointing out from the micelle (Figure ). This is in agreement with results from similar techniques used for simulation of phopholipid and bile salt systems.[20,21] The free fatty acids were positioned similarly to the phospholipids, with their tails deeply inserted into the micelle and their head groups pointing outward, accessible for interactions with polar water beads. The ratio of bile salts to phospholipids in the HIF systems varied from 0.7 to 16 (Table ). For micelles with a high bile salt-to-phospholipid ratio, the formation of a lipid core was not as definite as for those with a low ratio (see, e.g., HV3; Figure b). Since the free fatty acids are part of the more hydrophobic core, their carbon chain lengths strongly influence their positioning in the micelles. In our simulations, the free fatty acids were represented by five beads, corresponding to the structure of oleic or palmoleic acid with the force field used. A fatty acid with a shorter length than the phospholipids could possibly alter the packing and the conformation of a micelle. It could be of interest to investigate the effect of the carbon chain lengths of the free fatty acids on micelle formation.

Figure 4

(a) Largest micelles for each human intestinal fluid simulation. Molecular compositions are displayed as pie charts. (b) A single micelle displaying how the phospholipids (green tails, red and yellow head groups) and free fatty acids (green) form a denser core than the bile salts (gray). It is also clear that cholesterol, if present, resides in the core rather than the shell.

Drug–Micelle Interactions

To evaluate differences in drug solubilization of the five HV systems, three neutral, model drugs were added to the largest micelle from each HV. The micelle affinity for each simulated drug relative to each micelle was calculated using eq (Figure ). The estimated micelle affinity was greater than 1 for all model drugs, meaning the affinity of the drug was stronger for micelles than for water, and therefore solubilization would be likely. The model drugs showed a clear rank order of micelle affinity. Probucol had by far the highest micelle affinity, followed by fenofibrate and finally prednisolone. This rank order follows the hydrophobicity as described by the calculated and pH-adjusted partition coefficients between octanol and water (logP values). These are 1.4 for prednisolone, 5.3 for fenofibrate, and 10.0 for probucol.[36] The micelle affinity did deviate between simulations, but did not show any clear concentration-related correlation. In Figure a, the micelle affinities are displayed as average values from all six concentrations for each drug and HV. Since probucol showed a higher micelle affinity in all five micelles, relative to the other drugs, the simulations suggest that probucol would have a higher degree of partitioning into micelles and presumably a significantly increased solubility in the HIF compared to water.

Figure 5

Micelle affinities from simulations calculated as ratios of contacts between drug and micelle and between drug and water. (a) Micelle affinity from the largest micelle in each simulation, extracted to, and simulated in, 20 nm boxes. (b) Micelle affinities adjusted for the volume fraction of micelle and water in the larger, previously simulated, human intestinal fluid systems for five HVs (Methodsection, eq ). Values for each HV were sampled from six simulations with drug concentrations ranging from 0.8 to 4 mM, displayed here as average values with standard deviations (bar and whiskers).

The literature reports that solubility enhancement in HIF (compared to water) is greatest for probucol, followed by fenofibrate and finally prednisolone.[36,37] This reported order of increased solubility is in accordance with the trend for micelle affinities seen in our simulations, supporting that simulation methodology can differentiate the degree of micellar solubilization of neutral drugs with different lipophilicities. However, it should be noted that we are only simulating the smaller fraction of HIF colloids. Large colloids, such as vesicles and droplets, are not taken into account in the simulations; these can greatly affect drug solubilization, but on the other hand, these are not nearly as abundant in fasted-state intestinal fluids as in fed-state ones.[38]

The difference in micelle affinity of the model drugs was clear; however, we could only see small tendencies of interindividual variability between the HVs. For fenofibrate, the micelle affinity was similar in all HV micelles. The differences were amplified by extrapolating to the overall concentrations in HIF of HVs (eq ). This gave higher values for HV3 relative to HV9, due to the larger micellar volume in HV3. Interestingly, probucol seemed to have a stronger affinity for micelles in HV3, HV9, and HV20, and weaker affinity for HV6 and HV16 micelles. This difference was also amplified upon extrapolation. This could be explained by a drug-specific solubilization mechanism in the simulations.

In micelles from HV3, HV9, and HV20, probucol positioned itself deeper within the micelle core. In most cases, this would lead to an inner core of pure probucol, with very few to none contacts between probucol and water, and hence, a greatly increased ratio of contacts. Snapshots from the last frames in all of the drug–micelle simulations are presented in Figure S6, and contact ratio for drug–drug and drug–water for probucol can be seen in Figure S7. Similar solubilization behavior, but for a less hydrophobic drug (celecoxib), has been investigated by Elvang et al. in different SIF preparations.[39] They observe that micelles seem to swell with the addition of celecoxib; therefore, they propose that celecoxib behaves as a lipid in the bile-salt-rich micelle. Our simulations with probucol in micelles from HV6 and HV16 showed that the micelles did not manage to entrap probucol in their cores to the same extent as the other three HVs. Instead, probucol resided closer to the surface of the micelles, giving more contacts with water and therefore a lower contact ratio with the micelle. This poor drug entrapment could be linked to the composition of the micelles. HV6 and HV16 had a much lower bile salt–phospholipid ratio (2.1 and 1.6, respectively) than the other three micelles, which had 6–10 (Figure a). This indicates that the micelle probably requires an initial dense bile salt shell to incorporate probucol into the core. This was possible for HV3, HV9, and HV20 because their micelles had a high bile salt-to-phospholipid ratio. In summary, these simulations suggest that the micelle composition and drug lipophilicity affect the drug solubilization mechanism of the model drugs.

Conclusions

We have attempted to mimic fasted-state human intestinal fluids—focusing on the colloidal structures—by simulating concentrations of small-molecule components. These concentrations were based on previously acquired in vivo data from five healthy volunteers. The simulation methodology used the popular and available CG Martini force field. The self-associated colloidal structures were ellipsoidal micelles ranging from 2 to 7 nm. Some of these were then selected for further study of how well they could solubilize three model drugs (of neutral charge, but different hydrophobicities). The estimated order of relative micelle affinity for the model drugs was in line with solubility data from the literature. The solubilization mechanism differed, depending on the type of drug and the composition of the micelles. This type of simulation could improve understanding of drug solubility, solubilization, and other events taking place in HIF, in particular, substance-specific solubilization behavior.