- Email: [email protected]
Molecular dynamics simulations of spontaneous bile salt aggregation
Colloids and Surfaces A: Physicochem. Eng. Aspects 280 (2006) 182–193
Molecular dynamics simulations of spontaneous bile salt aggregation Dallas B. Warren a , David K. Chalmers b , Keith Hutchison c , Wenbin Dang c , Colin W. Pouton a,∗ a
Department of Pharmaceutical Biology and Pharmacology, Victorian College of Pharmacy, Monash University, Parkville, Vic. 3052, Australia b Department of Medicinal Chemistry, Victorian College of Pharmacy, Monash University, Parkville, Vic. 3052, Australia c Cardinal Health, Pharmaceutical Development, Somerset, NJ, USA Received 15 April 2005; received in revised form 31 January 2006; accepted 6 February 2006 Available online 29 March 2006
Abstract Bile salts are surfactants in bile that facilitate digestion, adsorption and excretion of various compounds. They have planar hydrophobic and hydrophilic faces and therefore exhibit some unusual properties; including the shape and size of the micelles that they form. Molecular dynamics simulations of the spontaneous aggregation of six bile salts (cholate (CHD), glycocholate (GCH), taurocholate (TCH), glycochenodeoxycholate (GCD), glycodeoxycholate (GDX) and glycolithocholate (GLC)) were performed in an aqueous phase to gain insight into their micellar structure. The aggregates that formed spontaneously from a random distribution of molecules ranged in size from 8 to 17 molecules. The structures are highly dynamic in nature and are on average oblate, but can vary from oblate, to spherical or prolate. Intermolecular hydrogen bonding within the micelles was found to be an important factor in determining the micelle size, structure and dynamics. The molecular arrangement within the micelles maximises the hydration of the hydrophilic chains and some favourable orientations for adjacent molecules were acquired. The dynamics of the micelles were investigated using the hydrogen-bond lifetime autocorrelation function correlation time, which exhibited a relationship with the degree of hydroxylation. Comparison of the proposed model to the three literature models showed some features of the disk shaped models of Cary and Small [M.C. Cary, D.M. Small, Arch. Intern. Med. 130 (1972) 506–527] and Kawamura et al. [H. Kawamura, Y. Murata, T. Yamaguchi, H. Igimi, M. Tanaka, G. Sugihara, J.P. Kratohvil, J. Phys. Chem. 93 (1989) 3321–3326], whereas the third, inverted helix model of Giglio et al. [E. Giglio, S. Loreti, N.V. Pavel, J. Phys. Chem. 92 (1988) 2858–2862] can be discounted. The proposed model is better than the existing models, which assumed a rigid and structured molecular arrangement. © 2006 Elsevier B.V. All rights reserved. Keywords: Molecular dynamics; Micelle structure; Spontaneous aggregation; Bile salts; Aggregation number
1. Introduction A detailed knowledge of the surfactant behaviour of bile salts is important for the understanding of drug absorption within the small intestine and determination of the influence of lipid formulations on this process. The current working hypothesis is that the majority of drugs are absorbed from free aqueous solution; however, most drugs in development have low aqueous solubility. For these types of compounds it is thought that vesicles and bile salt/phospholipid mixed micelles act as drug reservoirs. Consequently, the nature of these aggregates plays an important role in drug bioavailability. The current models of the molecular scale aggregation of bile salt surfactants consist
∗
Corresponding author. Tel.: +61 3 9903 9562; fax: +61 3 9903 9638. E-mail address: [email protected] (C.W. Pouton).
0927-7757/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.colsurfa.2006.02.009
of little beyond rudimentary space-filling models. As the foundation for studies of drug transport and distribution within the intestinal lumen we have used molecular dynamics methods to investigate the molecular behaviour of several important bile salt species. Bile salts are curious molecules, being surface active with non-classical surfactant like properties. Classical surfactants (such as sodium dodecyl sulphate) consist of a hydrophobic tail and a hydrophilic head which spontaneously aggregate into micelles in the aqueous phase above the critical micelle concentration. In contrast, bile salts do not have well defined tail and head groups, but exhibit planar polarity with hydroxyl groups generally located on one face and methyl groups on the opposite. A consequence of this planar polarity is that bile salts form smaller micelles than classical surfactants, in the region of 2–9 molecules [1,2]. This observation is perhaps not unexpected, as it is difficult to construct a large aggregate and maintain contact
D.B. Warren et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 280 (2006) 182–193 Table 1 Typical bile salt composition of bile within the gall bladder [6] Bile salt
Mol %
Glycocholate Taurocholate Glycochenodeoxycholate Taurochenodeoxycholate Glycodeoxycholate Taurodeoxycholate Glycolithocholate Taurolithocholate Glycoursodeoxycholate Tauroursodeoxycholate Sulfoglycolithocholate Sulfotaurolithocholate
30.9 12.0 29.3 11.2 9.3 2.2 0.4 0.3 0.4 1.6 1.9 0.5
183
the conjugated amino acid. Conjugation occurs principally with taurine (NH2 CH2 CH2 SO3 H) or glycine (NH2 CH2 COOH) in a ratio of approximately 1:3, with sodium (Na+ ) or potassium (K+ ) present as the counter ions. Of the five major bile acids present within bile, only two are generated by the body in the liver (primary) [5]; cholic (3␣,7␣,12␣-trihydroxy-5-cholan-24-oic acid) and chenodeoxycholic (3␣,7␣-dihydroxy-5-cholan-24oic acid) acids. The remaining three are products of bacterial action within the intestinal tract (secondary) [5]: deoxycholic (3␣,12␣-dihydroxy-5-cholan-24-oic acid), lithocholic (3␣-hydroxy-5-cholan-24-oic acid), and ursodeoxycholic (7hydroxy-5-cholan-24-oic acid) acids. Fig. 1 shows the molecular structures of these compounds as salts. 1.1. Models of the bile salt aggregate structure
between water and all of the hydrophilic faces. In addition, a well defined critical micellar concentration for bile salts does not exist with micelle formation occurring as a continuous association [3]. Bile salts, consisting of a hydroxylated steroid conjugated to an amino acid, are the dominant species present in bile. Bile is secreted from the liver and delivered into the small intestine after storage in the gall bladder. Bile facilitates digestion, adsorption and excretion of compounds such as bilirubin, a by-product of haemoglobin metabolism. It is composed of bile salts (67%, w/w) and phospholipids (22%, w/w), along with small amounts of protein (4.5%, w/w), cholesterol (4%, w/w) and bilirubin (0.33%, w/w) [4]. This ratio remains essentially constant once secreted from the liver until it reaches the intestine, although concentration occurs in the gall bladder [4]. A significant number of bile salt species are found within bile, a typical composition is listed in Table 1. Bile salts differ in the position, stereochemistry and degree of hydroxylation and in the identity of
The structure of bile salt micelles has been the subject of discussion since the 1960s and three models have been proposed within the literature, shown in Fig. 2. These are (i) primary (disk shaped with one or two layers) and secondary micelles (hydrogen bonding between primary micelles) [1]; (ii) disk shaped micelles with a single layer [6]; and (iii) a helical structure with counter ions and water on the axis of the helix [7]. Molecular dynamics (MD) simulations are a powerful technique to investigate aggregates on the molecular scale and have not been performed on bile salts only to date. MD has provided useful insight into the structure and dynamics of micelles, including model [8–10], cationic [11,12], glycoside [13], sodium dodecyl sulphate (SDS) [14,15], phospholipid [16,17] and bile salt/phospholipid [18] systems. To increase the speed of calculations and allow access to large simulation system sizes within a reasonable time, united atom forcefields have been extensively used. United atom forcefields reduce the number of atoms, and therefore the number of calculated interactions, by
Fig. 1. Structures of some bile salts: (CHD) sodium cholate, (TCH) sodium taurocholate, (GCH) sodium glycocholate, (GCD) sodium glycochenodeoxycholate, (GDX) sodium glycodeoxycholate, and (GLC) sodium glycolithocholate.
184
D.B. Warren et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 280 (2006) 182–193
Fig. 2. Representations of the literature proposed structures of bile salt micelles: (i) primary and secondary micelles by Carey and Small [1]; (ii) disk-like model by Kawamura et al. [6], and (iii) helical shaped micelles by Giglio et al. [7].
representing CH, CH2 and CH3 as a single “atom”. Simulations of salt/phospholipid micelles [18], simple surfactant micelles [15,19], and bilayers [19–23] using united atom forcefields successfully reproduce experimental results of interest to this study: the formation and structure of micelles. In this study we report MD simulations of the spontaneous aggregation of the bile salts cholate (CHD), glycocholate (GCH), taurocholate (TCH), glycochenodeoxycholate (GCD), glycodeoxycholate (GDX) and glycolithocholate (GLC) in the aqueous phase, each from a randomly distributed configuration. The structure and dynamics of the resulting aggregates is analysed, compared to the existing literature models and an improved model of bile salt association is proposed. 2. Method 2.1. Software and hardware The GROMACS [24–26] software package (www.gromacs. org), version 3.1.4 was used to perform the MD simulations. Calculations were performed on the following computer systems: a local Linux cluster, the Australian Partnership for Advanced Computing (APAC) super computer (www.apac.edu.au) and the Victorian Partnership for Advanced Computing (VPAC) Linux cluster (www.vpac.org). The local Linux cluster consisted of
AMD 200 MP processors, 1.67 GHz clock speed, with two CPUs and 512 M of memory per box, connected using a 100Base-T switched Ethernet. The APAC SC consisted of a Compaq SC with ES45s, 1 GHz clock speed, with four CPUs and 4 G of memory per box, connected using a Quadrics network. The VPAC LC consisted of Xenon 2.8 GHz clock speed dual CPUs connected using Myrinet. Parallel scaling performance of a sample simulation (25,000 steps) was performed to determine the optimum conditions to operate the simulations under to achieve an efficient use of the CPU time. From these simulations the optimum utilisation was determined to be two CPUs for the local Linux cluster, four for the APAC SC and two for the VPAC LC. The typical CPU time required was 50 CPU hours per ns of simulation, equating to a total of approximately 10,750 CPU hours for the simulations presented. Visualisation of the simulation trajectories was performed using the software package VMD [27] and the solvent accessible surface rendered images were produced by SURF [28] (probe radius = 0.14 nm). The spatial distribution function of atoms relative to the molecules was generated using g sdf [29] (bin width = 0.09 nm), version 1.25. Subsequent analysis scripts were all packaged with GROMACS. The radial distribution function between water atoms was calculated using g rdf (bin width = 0.001 nm). The intermolecular hydrogen bond lifetime was calculated with g hbond (average over all steroid hydroxyl
D.B. Warren et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 280 (2006) 182–193
group hydrogen bonds, donor-hydrogen-acceptor angle = 30◦ , donor-acceptor cut off radius = 0.35 nm, donor = OH, acceptor = O). The rotational autocorrelation functions was calculated using g rotacf (autocorrelation function of cross product of two vectors, n = ij × jk, where i = C3 , j = C12 and k = C15 of steroid backbone, see Fig. 1). The self diffusion coefficient was calculated with g msd. g msd calculates the mean square displacement of atoms/molecules from their initial positions, fits a straight line using the least squares method to the mean square displacement and computes the self diffusion constant using the Einstein relation. 2.2. Bile salt aggregation simulation specifications The united atom molecular simulations utilised the GROMACS united atom forcefield (ffgmx) [26], which is based on the GROMOS-87 force field [30]. A Ryckaert-Bellemans potential [31] was used for the CH2 –CH2 dihedral in the bile acid’s hydrophobic chain and the taurine amino acid group. The Ryckaert-Bellemans potential was required because the energy difference between anti and gauche confirmations using the ffgmx dihedral was too high. The energy difference between the anti- and gauche confirmations of a butane molecule were calculated for the ffgmx and Ryckaert-Bellemans CH2 –CH2 dihedrals as 6.6 and 2.9 kJ mol−1 , respectively. The energy difference between the anti- and gauche confirmations of a butane molecule using the Ryckaert-Bellemans (2.9 kJ mol−1 ) compares much more favourably with the butane ab initio value [32] of (2.6 kJ mol−1 ) than ffgmx (6.6 kJ mol−1 ). Rigid SPC water [33] was utilised as the water model and was constrained using SETTLE [34]. Remaining solute bonds were constrained by the LINCS algorithm [35] and temperature and pressure control was executed by Berensden coupling [36]. Periodic boundary conditions were employed on the cubic simulation box, along with a cut-off distance of 1.5 nm for electrostatic and 1.0 nm for Lennard-Jones non-bonding interactions. These cut-off distances have been found to be suitable in simulations of a bile salt and phospholipid system [18] and were used for the derivation of the force field. To achieve a large simulation time step, the fastest vibrational frequency within the bile salt molecules, exhibited by the polar hydrogens, was reduced by increasing the mass of the hydrogen atom to 4 Da. The attached heavy atom’s mass was decreased by the same amount to conserve mass [37]. The mass redistribution was only applied to the polar hydrogens within the bile salt molecules. This allowed the use of a 4 fs time step, with such large time steps being used increasingly in simulations [16,18,38–41]. All of the bile salt aggregation MD simulations were established and performed using the following procedure. Thirty-one bile salt molecules were randomly distributed within an 8 nm cube, giving a bile salt concentration of approximately 0.10 M. The box was then filled with the appropriate number of water molecules (approximately 1.6 × 104 ). A total of 77 sodium ions and 46 chloride ions were subsequently distributed by random replacement of water molecules, giving 31 counter-ions and a background salt concentration of approximately 0.15 M. The total number of atoms within all simulations was in the region
185
of 50,000. Energy minimisation, using the steepest descent technique [26], was used to remove any bad van der Waals contacts between atoms. This was followed by two short 20 ps simulations. The first short simulation was performed with tight temperature (τ T = 0.01 ps with reference temperature of 310 K) and no pressure coupling, and the second with both tight temperature (τ T = 0.01 ps) and isotropic pressure (τ P = 0.2 ps with reference pressure of 1 bar and 4.5 × 10−5 compressibility) coupling. From this second trajectory and coordinate files the full MD simulation was commenced. A time step of 4 fs for the MD simulation was used, along with Pref = 1 atm, τ P = 0.2 ps, Tref = 310 K and τ T = 0.05 ps. The bulk system properties were stable for the entire length of the simulations, oscillating around the average value. The total energy varied by <1%, while the temperature was maintained at 312 ± 2 K and the density at 990 ± 4 g dm−3 . The density compares favourably with that of pure water, 993 g dm−3 at 37 ◦ C [42]. Simulations continued for 10 to 15 ns after the system had reaching an equilibrium state based on the aggregation and solvent accessible surface area behaviour. 2.3. Forcefield comparison simulation specifications Five smaller scale simulations were performed to validate the utilisation of a united atom forcefield with heavy hydrogens, electrostatic cut-offs and large time step. All forcefield comparison simulations were performed with 6 CHD molecules, 6 sodium counter ions, 1,280 water molecules, started with a 3.5 nm simulation cube and run for 100 ns. Three different united atom forcefield simulations were performed using; 1.5 nm cut-off/4 fs, 1.5 nm cut-off/2 fs, and PME/2 fs. OPLS [43] was utilised as the all atom forcefield with TIP3P [44] the water model using; 2 fs time step, cut-off distance of 0.9 nm for electrostatics and Lennard-Jones non-bonding interactions, and particle-mesh Ewald (PME) [45] method for long range electrostatic interactions. A second OPLS simulation was performed with 1.5 nm cut-off for the electrostatics instead of PME. 3. Results and discussion 3.1. Forcefield comparison The aggregation behaviour simulations performed in this study use the ffgmx united atom forcefield with a 1.5 nm electrostatic cut off and 4 fs time step. To validate this method a series of five simulations were performed, see Table 2. All five simulations rapidly associated within 2 ns, forming a single aggregate of all 6 bile salt molecules that remained intact for the remaining simulation time. Self diffusion coefficients of water and sodium ions were calculated over the entire 100 ns trajectories, with the results summarised in see Table S1, Supplementary Information. Table S1 also contains experimental and molecular dynamics simulation literature values for water and sodium ions. The self diffusion coefficients obtained in this study (SPC 3.87 × 10−5 cm2 s−1 and TIP3P 5.73 × 10−5 cm2 s−1 ) are con-
186
D.B. Warren et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 280 (2006) 182–193
Table 2 Specifications for molecular dynamics simulations performed for method validation Water model
Time step (fs)
Electrostatic treatment
SPC SPC SPC TIP3P TIP3P
4 2 2 2 2
1.5 nm cut-off 1.5 nm cut-off PME 1.5 nm cut-off PME
sistent with literature values under similar MDS conditions (SPC 3.85 × 10−5 cm2 s−1 and TIP3P 5.19 × 10−5 cm2 s−1 ) [46]. Both SPC and TIP3P water models are known to overestimate the diffusivity of water [47], which is 2.30 × 10−5 cm2 s−1 at 25.2 ◦ C and 3.89 × 10−5 cm2 s−1 at 50.2 ◦ C [48]. Increasing the time step from 2 to 4fs for the united atom forcefield only has a minor effect on the self diffusion coefficients, increasing the water and sodium ions values by approximately 5 and 15%, respectively. The spatial distribution functions (SDF) of the water oxygen atoms around the CHD molecules were calculated over the 5 to 100 ns simulation period for the united-atom/1.5 nm cut-off/4 fs and all-atom/PME/2 fs simulations, see Fig. S1, Supplementary Information. The shape and extent of the probability surfaces were similar, with only minor variations. These minor variations are due to differences in the dihedral energy functions of both forcefields and the generally more irregular molecular surface generated by an all atom versus a united atom forcefield. Radial distribution functions (RDF) of oxygen-oxygen, oxygen-hydrogen and hydrogen-hydrogen for the water molecules were calculated for each simulation. Alterations to the electrostatics treatment for both forcefields caused no drift in the location of the RDF peaks, but did increase the peak height or ordering. Increasing the time step to 4 fs and using 1.5 nm cutoff for electrostatics displayed no significant alteration in the RDFs, as displayed by the O-O RDF in Fig. S2, Supplementary Information. Differences in the height and location of the second solvation shell are observed between the two water models, illustrated by the O-O RDF shown in Fig. S3, Supplementary Information. The difference in the RDF around the second solvation shell is due to the water models (SPC versus TIP3P), as has been observed previously [49], rather than the simulation parameters. The average CHD rotational autocorrelation functions were calculated for all five simulations and then fitted with a biexponential function, producing two correlation times. These correlation times are summarised in Table S2, Supplementary Information. The differences between the CHD rotational autocorrelation functions for the united-atom/1.5 nm cut-off/4 fs and all-atom/PME/2 fs are shown in Fig. 3. The two correlation times for these curves are 1.02 and 6.67 × 10−6 ns−1 , and 1.54 and 1.52 × 10−6 ns−1 , respectively. This indicates that micelles simulated using the all atom forcefield are more dynamic on the 500 ps timescale than the united atom forcefield. The effect of increasing the time step on the united atom/1.5 nm cut-off system from 2 to 4 fs is to double both the short and long correlation
Fig. 3. Average rotational autocorrelation function of the 6 cholate molecule simulations: (—) united atom and (— —) all atom forcefield.
times. This increase in correlation time corresponds to a decrease in the rotational freedom of the molecules within the aggregates. The variation in rotational autocorrelation function correlation time simulation conditions is accounted for by considering only the relative values in analysis within this study, not the absolute value. 3.2. Spontaneous aggregation To investigate bile salt aggregation, six bile salt systems were simulated, with the total simulation time performed for each being; CHD (20 ns), GCH (45 ns), TCH (55 ns), GCD (40 ns), GDX (30 ns), and GLC (25 ns). Spontaneous aggregation was observed in all cases. This process is illustrated using GCH in Fig. 4. The starting configuration is presented in Fig. 4(i). After 20 ns (Fig. 4(iii)) all but two of the GCH molecules have aggregated to form a single micelle. At the conclusion of the simulation (Fig. 4(iv)) the free molecules have encountered the micelle and subsequently caused it to split into two separate micelles. The aggregation process was monitored by calculating the number of aggregates and the aggregation numbers over time (shown in Fig. 5). A molecule was defined as belonging to an aggregate when two carbon atoms were within 0.5 nm of a neighbouring molecule’s carbon atoms. This simple rule provided a low noise level and classified only closely associated bile salts as part of an aggregate. The initial aggregation behaviour for all of the bile salts was similar, with spontaneous aggregation complete in approximately 10 ns. For the bile salts CHD, GDX, GCD and GLC, initial aggregation was followed by molecular reorientation to produce a relatively stable, dynamic state. GCH and TCH exhibited more interesting phenomena. Within 16 ns the majority of
D.B. Warren et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 280 (2006) 182–193
187
Fig. 4. Evolution of the glycocholate simulation at time (i) t = 0 ns; (ii) t = 5 ns; (iii) t = 20 ns; and (iv) t = 45 ns. Solid lines indicate the edge of the simulation box. Also displayed are the ions as points and water molecules are omitted for clarity.
the GCH molecules became associated to form a single aggregate, see Fig. 5(ii). After an additional 16 ns, the last two free molecules were absorbed to form a large prolate micelle containing all 31 molecules, which was unstable and short lived. After 2 ns of reorganisation, this micelle split into two approximately equal sized micelles of 15 and 16 molecules, which were then stable for the remaining simulation time. For TCH the aggregation process was far more dynamic, see Fig. 5(vi). Rapid aggregation occurred during the first 10 ns and by 20 ns two micelles containing 10 molecules were present, with the remainder left as free monomers or dimers. Aggregates indicated by the figure to contain around 15 or more molecules were due to transient interactions between the hydrophilic chains of TCH molecules in adjacent aggregates, not between the steroidal backbones. For the remainder of the simulation two micelles containing approximately 10 molecules were maintained, with relatively rapid exchange of monomers and some smaller transient aggregates forming from the remaining molecules. The size and number of aggregates present at the conclusion of the simulations are shown in Table 3, with experimental aggregation numbers provided for comparison. It should be noted that definitive conclusions on aggregate sizes cannot be formed based on the very small data sets presented here. The simula-
tions performed are insufficient in size and time frame to gain an acceptable statistical view of the aggregate size and shape. To be conclusive, significantly larger systems would be required. Moving to a system containing in the order of three hundred bile salt molecules and a simulation time on the microsecond time scale may be sufficient. Nevertheless, provided the limitations of the simulations are considered, these simulations can give initial insight into the behaviour of bile salt micelles on the molecular scale. Table 3 Aggregate size of bile salts from this study and published experimental aggregation numbers Bile salt
Simulation aggregate size
Literature aggregation number
CHD GCH TCH GCD Taurocheno-deoxycholate GDX Taurodeoxy-cholate GLC
31 16, 15 11, 10, 10 × 1 16, 14, 1 – 17, 11, 3 × 1 – 13, 10, 8
15–30 [50], 3.09–3.35 [2] – 3 [51], 2.59 [2] – 2051 – 2451 –
188
D.B. Warren et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 280 (2006) 182–193
Fig. 5. Aggregation behaviour of the molecular dynamics simulation of the bile salts: (i) CHD; (ii) GCH; (iii) GDX; (iv) GCD; (v) GLC, and (vi) TCH. Bubble size is proportional to the number of molecules within aggregates of that aggregation number.
D.B. Warren et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 280 (2006) 182–193
189
Fig. 6. Solvent accessible surface of the various bile salt micelles at the simulation end: (i) CHD; (ii) GCH; (iii) GDX; (iv) GCD; (v) GLC, and (vi) TCH. Surface probe radius = 0.14 nm and each colour indicates a different molecule. (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of the article.)
190
D.B. Warren et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 280 (2006) 182–193
3.3. Aggregate structure The instantaneous arrangements of molecules at the end of each simulation are shown as solvent accessible surfaces in Fig. 6. In general, the bile salts make contact through the steroidal backbone with the hydrophilic chain exposed to the solvent. It is evident that the arrangement of molecules and micellar shapes differ from the three proposed literature models. Furthermore, the micelles are dynamic. For all of the bile salts, the shapes of the micelles were fluid and varied from prolate, through roughly spherical to oblate. It was rare to find an entire hydrophobic face in contact with the solvent for a significant length of time. The hydrophilic face was predominately in contact with the solvent or the hydrophilic face of an adjacent bile salt. The principle radius of gyration for each of the micelles was calculated (atoms explicitly mass weighted) over the final 10 to 15 ns of the trajectories, shown in Table 4. Micelles ranged in size from a long axis radius of 1.31 nm for a 16 molecule GCH micelle, down to 1.00 nm for an 11 molecule GDX micelle. The
Table 4 Aggregation number and principle radii of gyration of MD simulation bile salt micelles Bile salt
Aggregation number
rx (nm)
ry (nm)
rz (nm)
GCH
16 15
1.31 1.08
1.27 1.02
0.86 0.91
GDX
16 11
1.25 1.00
1.19 0.94
0.94 0.83
GCD
16 14
1.26 1.12
1.19 1.06
0.91 0.89
GLC
13 10 8
1.15 1.05 1.11
1.07 0.97 1.02
0.96 0.85 0.81
Balance of molecules (total of 31 per system) are present as un-associated molecules.
principle radii indicate that the average shape of the micelles is oblate, a feature of the model of Kawamura et al., with rx (1.00–1.31 nm) and ry (0.94–1.27 nm) being approximately the same value and the third radius, rz (0.81–0.96 nm), being smaller.
Fig. 7. Four favoured arrangements of GCH molecules within a GCH micelle. Atom colouring is; cyan, carbon; red, oxygen; white, polar hydrogen; and blue, nitrogen. (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of the article.)
D.B. Warren et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 280 (2006) 182–193
An important attribute of the bile salts within the simulated micelles was that the molecules are not arranged parallel to their neighbours, a characteristic of both the disk shape models. Rather, the hydrophilic chain attached to the steroidal backbone maximised the distance from neighbouring chains, allowing maximal hydration and forcing an offset orientation between the molecules. Only a small proportion of parallel alignment occurred. In this case, the molecules were usually aligned with concurrent OH–OH and terminal-terminal chain hydrogen bonding. Parallel alignment involving the hydrophobic face occurred only as transient structures; a result of steric interference from the two methyl groups. A distribution function was calculated for the angle between neighbouring molecules (data not shown), but this showed no substantial ordering, indicating that there is only weak long-range ordering within the micelle. Furthermore, this suggests that the rigid, ordered arrangements of the disk shaped models are overly idealistic. A comparison of bile salt simulations showed that intermolecular H-bonding of the steroid hydroxyl groups was a major factor influencing micelle size, structure and dynamics. It is evident that bile salts with two (GDX and GCD) or three (CHD, GCH and TCH) hydroxyl groups form H-bonds within the micelle in the absence of water. This H-bonding allows the molecules to be buried deeper within the micellar structure, away from the bulk solvent, and to form larger and more spherical micellar structures. In contrast GLC, which has hydrophilic groups located at opposite ends of the molecule, packs differently and adopts cylindrical or disk shaped aggregates to satisfy the interactions required by these hydrophilic groups. This type of stacking reduces the maximum aggregation number for GLC, as shown by the lower aggregation numbers in Table 3. Representative examples of both hydrophilic and hydrophobic intermolecular interactions are shown for GCH in Fig. 7. These examples were obtained by aligning the neighbouring bile salt molecule to high probability regions of the SDF. It was apparent that a large number of alignments are possible, reinforcing the observation that these micelles are dynamic. The hydrophilic (face to face) positioning maximises the number of H-bonds between hydroxyl groups, see Fig. 7(i) and (ii). The hydrophobic alignments have the two methyl groups interlocked to maximise contact, see Fig. 7(iii) and (iv). Similar alignments are seen with the other bile salts but, as the hydroxyl group number decreases, the number of H-bonds reduces, providing more freedom of movement. For the tri-hydroxyl GCH, three anchoring points exist, restricting the molecular movement. In the case of the di-hydroxyls, GDX and GCD, only two points exist, allowing the rest of the molecule to rotate around these more freely. The mono-hydroxyl bile salt, GLC, with only a single anchoring point, has significantly larger conformational freedom while maintaining the single H-bond. The micelle model of Carey and Small [1] proposes that H-bonding between the primary micelles occurs, allowing the formation of secondary micelles that are observed at high bile salt concentrations. Such intermicellar interactions were observed in our simulations and were most evident for GLC and TCH. These interactions between micelles occurred through Hbonding of the hydrophilic tails, not between steroid backbone
191
hydroxyl groups. This was particularly evident in the GLC simulation which separated into three micelles containing 8, 10 and 13 molecules that remained in contact through the hydrophilic tails. The unusual proposed micellar structures of Giglio et al. [7], with the hydrophobic surface exposed to the aqueous phase, are far removed from the observed structures. The proposal that significant hydrophobic area is in contact with the surrounding solvent at low bile salt concentrations is counter-intuitive. The energy reward of having the hydrophilic groups arranged in a helical structure around the counter ions would have to be greater than the energy cost of water structuring from contact with the large hydrophobic surface. Such an arrangement might be possible at high bile salt concentrations, where phase inversion might occur. However, from these simulations it appears to be an unlikely arrangement at micellar phase concentrations. 3.4. Aggregate dynamics The dynamic behaviour of the bile salt aggregates was investigated using a variety of autocorrelation functions. The most insight was obtained from the H-bond lifetime autocorrelation function of the steroid backbone –OH groups, shown in Fig. 8, which exhibited an increasing decay rate as follows: GCH, GCD = CHD, GDX, GLC and TCH. A bi-exponential function was fitted to the autocorrelation functions, with the short correlation time process decaying within 100 ps and showed no correlation to any molecular parameter. Increasing long correlation time is observed across the series of bile salts with the same headgroup (Fig. 9); GCH, GDX, GCD, and GLC. An inverse relationship was established between the long correlation time and the average number of hydroxyl H-bonds per H-bond accep-
Fig. 8. Average hydrogen bond lifetime autocorrelation function of the 31 bile salt molecule simulations: (—) CHD; (···) GCH; (– –) GDX; (– · ·) GCD; (— —) GLC and (– ·) TCH.
192
D.B. Warren et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 280 (2006) 182–193
autocorrelation function correlation time and the degree of hydroxylation. Intermolecular hydrogen-bonding is a major factor influencing size, structure and dynamics of the micelles. This is consistent with the aggregation behaviour of the systems. Intermicellar interactions similar to those proposed by Carey and Small were observed, which could allow formation of larger, secondary micelles. Although the structures observed in these MD simulations possess some components of those previously proposed by Carey and Small1 , and Kawamura et al. [6], they reveal more disordered and dynamic behaviour than that provided by these simple models. Acknowledgements Financial support from Cardinal Health is gratefully acknowledged. We would also like to thank the Australian Partnership for Advanced Computing (APAC) for the Merit Allocation Scheme grant and the Victorian Partnership for Advanced Computing for the provision of computing resources. Appendix A. Supplementary data Fig. 9. Relationship between the correlation time of the H-bond lifetime autocorrelation function and the number of OH H-bonds per number of acceptors: (䊉) CHD, () TCH, () GCH, () GDX, ( ) GCD, and () GLC.
tor, see Fig. 9. Bile salts with lower degrees of intra-micellar Hbonding are able to rotate more freely, and therefore have higher long correlation times. The precise geometry of the hydroxylation is also clearly important to the dynamics, as exhibited by GDX and GCD which have a second hydroxyl group in the 12␣ and 7␣ positions, respectively, and are markedly different in the degree of H-bonding. The effect of the hydrophilic chain on the dynamics is made apparent by the series CHD, GCH and TCH. The addition of the glycine group (CHD to GCH) decreases the correlation time, indicating a restriction in molecular movement. Conversely, the addition of taurine (CHD to TCH) increases the correlation time significantly, indicating a more dynamic system. This observation is consistent with the aggregation process for these two species shown in Fig. 5. The difference between the effect of taurine and glycine might be due changes in the headgroup interactions (steric, hydration and electrostatic); however, the cause requires further investigation. 4. Summary and conclusions MD simulations of six bile salt types were performed to observe their aggregation behaviour in the aqueous phase. The micelles formed were found on average to be oblate with the molecules arranged to maximise hydration of the hydrophilic chain. Contact between the steroidal backbones was mediated through hydrophobic interactions and through hydrogen bonding of hydroxyl groups. Several more common neighbouring molecular orientations were identified but overall there was no favoured angle. The parallel orientations of neighbouring bile salt molecules exhibited by two of the literature models are unlikely. The bile salt micelles were found to be dynamic with an observed relationship between the hydrogen bond lifetime
Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.colsurfa.2006.02.009. References [1] M.C. Carey, D.M. Small, Arch. Intern. Med. 130 (1972) 506. [2] A. Coello, F. Meijide, E.R. Nunez, J.V. Tato, J. Pharm. Sci. 85 (1996) 9. [3] J.P. Kratohvil, W.P. Hsu, M.A. Jacobs, T.M. Aminabhavi, Y. Mukunoki, Colloid Polym. Sci. 261 (1983) 781. [4] C.J. O’Connor, R.G. Wallace, Adv. Colloid Interf. Sci. 22 (1985) 1. [5] A.F. Hofmann, Biliary secretion and excretion: the hepatobiliary component of the enterohepatic circulation of bile acids, in: L.R. Johnson (Ed.), Physiology of the Gastrointestinal Tract, 2, third ed., Raven Press, New York, 1994, p. 1555. [6] H. Kawamura, Y. Murata, T. Yamaguchi, H. Igimi, M. Tanaka, G. Sugihara, J.P. Kratohvil, J. Phys. Chem. 93 (1989) 3321. [7] E. Giglio, S. Loreti, N.V. Pavel, J. Phys. Chem. 92 (1988) 2858. [8] D.W.R. Gruen, J. Phys. Chem. 89 (1985) 146. [9] B.F. Abu-Sharkh, E.Z. Hamad, Langmuir 20 (2004) 254. [10] M.C. Woods, J.M. Haile, J.P. O’Connell, J. Phys. Chem. 90 (1986) 1875. [11] J.-B. Maillet, V. Lachet, P.V. Coveney, Phys. Chem. Chem. Phys. 1 (1999) 5277. [12] H. Shinto, S. Morisada, M. Miyahara, K. Higashitani, Langmuir 20 (2004) 2017. [13] S. Bogusz, R.M. Venable, R.W. Pastor, J. Phys. Chem. B Mater. Surf. Interf. Biophys. 104 (2000) 5462. [14] C.D. Bruce, M.L. Berkowitz, L. Perera, M.D.E. Forbes, J. Phys. Chem. B Mater. 106 (2002) 3788. [15] A.R. Rakitin, G.R. Pack, J. Phys. Chem. B Mater. 108 (2004) 2712. [16] D.P. Tieleman, D. van der Spoel, H.J.C. Berendsen, J. Phys. Chem. B Mater. Surf. Interf. Biophys. 104 (2000) 6380. [17] S.J. Marrink, A.E. Mark, J. Am. Chem. Soc. 125 (2003) 15233. [18] S.J. Marrink, A.E. Mark, Biochemistry 41 (2002) 5375. [19] L.D. Schuler, P. Walde, P.L. Luisi, W.F. van Gunsteren, Eur. J. Biochem. 30 (2001) 330. [20] S.W. Chiu, E. Jakobsson, R.J. Mashl, H.L. Scott, Biophys. J. 83 (2002) 1842. [21] E. Lindahl, O. Edholm, J. Chem. Phys. 115 (2001) 4938. [22] M. Bachar, P. Brunelle, D.P. Tieleman, A. Rauk, J. Phys. Chem. B Mater. 108 (2004) 7170.
D.B. Warren et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 280 (2006) 182–193 [23] A.M. Smondyrev, M.L. Berkowitz, Biophys. J. 80 (2001) 1649. [24] H.J.C. Berendsen, D. van der Spoel, R. van Drunen, Comput. Phys. Commun. 91 (1995) 43. [25] E. Lindahl, B. Hess, D. van der Spoel, J. Mol. Model. 7 (2001) 306. [26] D. van der Spoel, A.R. van Buuren, E. Apol, P.J. Meulenhoff, D.P. Tieleman, A.L.T.M. Sijbers, B. Hess, K.A. Feenstra, E. Lindahl, R. van Drunen, H.J.C. Berendsen, GROMACS User Manual Version 3.1.1, Department of Biophysical Chemistry, Groningen, 2002. [27] W. Humphrey, A. Dalke, K. Schulten, J. Mol. Graph. 14 (1996) 33. [28] A. Varshney, F.P. Brooks, W.V. Wright, IEEE Comput. Graph. Appl. 14 (1994) 19. [29] C. Freudenberger, g sdf; 1.25 ed.; Chemistry, University of Ulm, 2003. [30] W.F. van Gunsteren, H.J.C. Berendsen, Gromos-87 manualGroningen, 1987. [31] J.P. Ryckaert, A. Bellemans, Faraday Discuss. Chem. Soc. 66 (1978) 95. [32] N.L. Allinger, J.T. Fermann, W.D. Allen, H.F. Schaefer III, J. Chem. Phys. 106 (1997) 5143. [33] H.J.C. Berendsen, J.P.M. Postma, W.F. van Gunsteren, J. Hermans, Interaction models for water in relation to protein hydration, in: B. Pullman (Ed.), Intermolecular forces, D. Reidel Publishing Company, Dordrecht, 1981, pp. 331–342. [34] S. Miyamoto, P.A. Kollman, J. Comput. Chem. 13 (1992) 952. [35] B. Hess, H. Bekker, H.J.C. Berendsen, J.G.E.M. Fraaije, J. Comput. Chem. 18 (1997) 1463.
193
[36] H.J.C. Berendsen, J.P.M. Postma, W.F. Gunsteren, A. Dinola, J.R. Haak, J. Chem. Phys. 81 (1984) 3684. [37] K.A. Feenstra, B. Hess, H.J.C. Berendsen, J. Comput. Chem. 20 (1999) 786. [38] G. Colombo, S.-J. Marrink, A.E. Mark, Biophys. J. 84 (2003) 2331. [39] S. Bandyopadhyay, J. Chanda, Langmuir 19 (2003) 10443. [40] S.J. Marrink, E. Lindahl, O. Edholm, A.E. Mark, J. Am. Chem. Soc. 123 (2001) 8638. [41] S.J. Marrink, D.P. Tieleman, Biophys. J. 83 (2003) 2386. [42] G.F.C. Rogers, Y.R. Mayhew, Thermodynamic and Transport Properties of Fluids, fourth ed., Basil Blackwell, Cambridge, 1991. [43] W.L. Jorgensen, D.S. Maxwell, J. Tirado-Rives, J. Am. Chem. Soc. 118 (1996) 11225. [44] W.L. Jorgensen, J. Chandrasekhar, J.D. Madura, R.W. Impey, M.L. Klein, J. Chem. Phys. 79 (1983) 926. [45] U. Essmann, L. Perera, M.L. Berkowitz, D. Tom, H. Lee, L.G. Pedersen, J. Chem. Phys. 103 (1995) 8577. [46] M.W. Mahoney, W.L. Jorgensen, J. Chem. Phys. 114 (2001) 363. [47] B. Guillot, J. Mol. Liquids 101 (2002) 219. [48] K. Krynicki, C.D. Green, D.W. Sawyer, Faraday Discuss. Chem. Soc. 66 (1978) 199. [49] W.L. Jorgensen, C. Jenson, J. Comput. Chem. 19 (1998) 1179. [50] B. Lindman, N. Kamenka, H. Fabre, J. Ulmius, T. Wiellock, J. Colloid Interf. Sci. 73 (1980) 556. [51] N.A. Mazer, M.C. Carey, R.F. Kwasnick, G.B. Benedek, Biochemistry 18 (1979) 3064.
Molecular dynamics simulations of spontaneous bile salt aggregation Dallas B. Warren a , David K. Chalmers b , Keith Hutchison c , Wenbin Dang c , Colin W. Pouton a,∗ a
Department of Pharmaceutical Biology and Pharmacology, Victorian College of Pharmacy, Monash University, Parkville, Vic. 3052, Australia b Department of Medicinal Chemistry, Victorian College of Pharmacy, Monash University, Parkville, Vic. 3052, Australia c Cardinal Health, Pharmaceutical Development, Somerset, NJ, USA Received 15 April 2005; received in revised form 31 January 2006; accepted 6 February 2006 Available online 29 March 2006
Abstract Bile salts are surfactants in bile that facilitate digestion, adsorption and excretion of various compounds. They have planar hydrophobic and hydrophilic faces and therefore exhibit some unusual properties; including the shape and size of the micelles that they form. Molecular dynamics simulations of the spontaneous aggregation of six bile salts (cholate (CHD), glycocholate (GCH), taurocholate (TCH), glycochenodeoxycholate (GCD), glycodeoxycholate (GDX) and glycolithocholate (GLC)) were performed in an aqueous phase to gain insight into their micellar structure. The aggregates that formed spontaneously from a random distribution of molecules ranged in size from 8 to 17 molecules. The structures are highly dynamic in nature and are on average oblate, but can vary from oblate, to spherical or prolate. Intermolecular hydrogen bonding within the micelles was found to be an important factor in determining the micelle size, structure and dynamics. The molecular arrangement within the micelles maximises the hydration of the hydrophilic chains and some favourable orientations for adjacent molecules were acquired. The dynamics of the micelles were investigated using the hydrogen-bond lifetime autocorrelation function correlation time, which exhibited a relationship with the degree of hydroxylation. Comparison of the proposed model to the three literature models showed some features of the disk shaped models of Cary and Small [M.C. Cary, D.M. Small, Arch. Intern. Med. 130 (1972) 506–527] and Kawamura et al. [H. Kawamura, Y. Murata, T. Yamaguchi, H. Igimi, M. Tanaka, G. Sugihara, J.P. Kratohvil, J. Phys. Chem. 93 (1989) 3321–3326], whereas the third, inverted helix model of Giglio et al. [E. Giglio, S. Loreti, N.V. Pavel, J. Phys. Chem. 92 (1988) 2858–2862] can be discounted. The proposed model is better than the existing models, which assumed a rigid and structured molecular arrangement. © 2006 Elsevier B.V. All rights reserved. Keywords: Molecular dynamics; Micelle structure; Spontaneous aggregation; Bile salts; Aggregation number
1. Introduction A detailed knowledge of the surfactant behaviour of bile salts is important for the understanding of drug absorption within the small intestine and determination of the influence of lipid formulations on this process. The current working hypothesis is that the majority of drugs are absorbed from free aqueous solution; however, most drugs in development have low aqueous solubility. For these types of compounds it is thought that vesicles and bile salt/phospholipid mixed micelles act as drug reservoirs. Consequently, the nature of these aggregates plays an important role in drug bioavailability. The current models of the molecular scale aggregation of bile salt surfactants consist
∗
Corresponding author. Tel.: +61 3 9903 9562; fax: +61 3 9903 9638. E-mail address: [email protected] (C.W. Pouton).
0927-7757/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.colsurfa.2006.02.009
of little beyond rudimentary space-filling models. As the foundation for studies of drug transport and distribution within the intestinal lumen we have used molecular dynamics methods to investigate the molecular behaviour of several important bile salt species. Bile salts are curious molecules, being surface active with non-classical surfactant like properties. Classical surfactants (such as sodium dodecyl sulphate) consist of a hydrophobic tail and a hydrophilic head which spontaneously aggregate into micelles in the aqueous phase above the critical micelle concentration. In contrast, bile salts do not have well defined tail and head groups, but exhibit planar polarity with hydroxyl groups generally located on one face and methyl groups on the opposite. A consequence of this planar polarity is that bile salts form smaller micelles than classical surfactants, in the region of 2–9 molecules [1,2]. This observation is perhaps not unexpected, as it is difficult to construct a large aggregate and maintain contact
D.B. Warren et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 280 (2006) 182–193 Table 1 Typical bile salt composition of bile within the gall bladder [6] Bile salt
Mol %
Glycocholate Taurocholate Glycochenodeoxycholate Taurochenodeoxycholate Glycodeoxycholate Taurodeoxycholate Glycolithocholate Taurolithocholate Glycoursodeoxycholate Tauroursodeoxycholate Sulfoglycolithocholate Sulfotaurolithocholate
30.9 12.0 29.3 11.2 9.3 2.2 0.4 0.3 0.4 1.6 1.9 0.5
183
the conjugated amino acid. Conjugation occurs principally with taurine (NH2 CH2 CH2 SO3 H) or glycine (NH2 CH2 COOH) in a ratio of approximately 1:3, with sodium (Na+ ) or potassium (K+ ) present as the counter ions. Of the five major bile acids present within bile, only two are generated by the body in the liver (primary) [5]; cholic (3␣,7␣,12␣-trihydroxy-5-cholan-24-oic acid) and chenodeoxycholic (3␣,7␣-dihydroxy-5-cholan-24oic acid) acids. The remaining three are products of bacterial action within the intestinal tract (secondary) [5]: deoxycholic (3␣,12␣-dihydroxy-5-cholan-24-oic acid), lithocholic (3␣-hydroxy-5-cholan-24-oic acid), and ursodeoxycholic (7hydroxy-5-cholan-24-oic acid) acids. Fig. 1 shows the molecular structures of these compounds as salts. 1.1. Models of the bile salt aggregate structure
between water and all of the hydrophilic faces. In addition, a well defined critical micellar concentration for bile salts does not exist with micelle formation occurring as a continuous association [3]. Bile salts, consisting of a hydroxylated steroid conjugated to an amino acid, are the dominant species present in bile. Bile is secreted from the liver and delivered into the small intestine after storage in the gall bladder. Bile facilitates digestion, adsorption and excretion of compounds such as bilirubin, a by-product of haemoglobin metabolism. It is composed of bile salts (67%, w/w) and phospholipids (22%, w/w), along with small amounts of protein (4.5%, w/w), cholesterol (4%, w/w) and bilirubin (0.33%, w/w) [4]. This ratio remains essentially constant once secreted from the liver until it reaches the intestine, although concentration occurs in the gall bladder [4]. A significant number of bile salt species are found within bile, a typical composition is listed in Table 1. Bile salts differ in the position, stereochemistry and degree of hydroxylation and in the identity of
The structure of bile salt micelles has been the subject of discussion since the 1960s and three models have been proposed within the literature, shown in Fig. 2. These are (i) primary (disk shaped with one or two layers) and secondary micelles (hydrogen bonding between primary micelles) [1]; (ii) disk shaped micelles with a single layer [6]; and (iii) a helical structure with counter ions and water on the axis of the helix [7]. Molecular dynamics (MD) simulations are a powerful technique to investigate aggregates on the molecular scale and have not been performed on bile salts only to date. MD has provided useful insight into the structure and dynamics of micelles, including model [8–10], cationic [11,12], glycoside [13], sodium dodecyl sulphate (SDS) [14,15], phospholipid [16,17] and bile salt/phospholipid [18] systems. To increase the speed of calculations and allow access to large simulation system sizes within a reasonable time, united atom forcefields have been extensively used. United atom forcefields reduce the number of atoms, and therefore the number of calculated interactions, by
Fig. 1. Structures of some bile salts: (CHD) sodium cholate, (TCH) sodium taurocholate, (GCH) sodium glycocholate, (GCD) sodium glycochenodeoxycholate, (GDX) sodium glycodeoxycholate, and (GLC) sodium glycolithocholate.
184
D.B. Warren et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 280 (2006) 182–193
Fig. 2. Representations of the literature proposed structures of bile salt micelles: (i) primary and secondary micelles by Carey and Small [1]; (ii) disk-like model by Kawamura et al. [6], and (iii) helical shaped micelles by Giglio et al. [7].
representing CH, CH2 and CH3 as a single “atom”. Simulations of salt/phospholipid micelles [18], simple surfactant micelles [15,19], and bilayers [19–23] using united atom forcefields successfully reproduce experimental results of interest to this study: the formation and structure of micelles. In this study we report MD simulations of the spontaneous aggregation of the bile salts cholate (CHD), glycocholate (GCH), taurocholate (TCH), glycochenodeoxycholate (GCD), glycodeoxycholate (GDX) and glycolithocholate (GLC) in the aqueous phase, each from a randomly distributed configuration. The structure and dynamics of the resulting aggregates is analysed, compared to the existing literature models and an improved model of bile salt association is proposed. 2. Method 2.1. Software and hardware The GROMACS [24–26] software package (www.gromacs. org), version 3.1.4 was used to perform the MD simulations. Calculations were performed on the following computer systems: a local Linux cluster, the Australian Partnership for Advanced Computing (APAC) super computer (www.apac.edu.au) and the Victorian Partnership for Advanced Computing (VPAC) Linux cluster (www.vpac.org). The local Linux cluster consisted of
AMD 200 MP processors, 1.67 GHz clock speed, with two CPUs and 512 M of memory per box, connected using a 100Base-T switched Ethernet. The APAC SC consisted of a Compaq SC with ES45s, 1 GHz clock speed, with four CPUs and 4 G of memory per box, connected using a Quadrics network. The VPAC LC consisted of Xenon 2.8 GHz clock speed dual CPUs connected using Myrinet. Parallel scaling performance of a sample simulation (25,000 steps) was performed to determine the optimum conditions to operate the simulations under to achieve an efficient use of the CPU time. From these simulations the optimum utilisation was determined to be two CPUs for the local Linux cluster, four for the APAC SC and two for the VPAC LC. The typical CPU time required was 50 CPU hours per ns of simulation, equating to a total of approximately 10,750 CPU hours for the simulations presented. Visualisation of the simulation trajectories was performed using the software package VMD [27] and the solvent accessible surface rendered images were produced by SURF [28] (probe radius = 0.14 nm). The spatial distribution function of atoms relative to the molecules was generated using g sdf [29] (bin width = 0.09 nm), version 1.25. Subsequent analysis scripts were all packaged with GROMACS. The radial distribution function between water atoms was calculated using g rdf (bin width = 0.001 nm). The intermolecular hydrogen bond lifetime was calculated with g hbond (average over all steroid hydroxyl
D.B. Warren et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 280 (2006) 182–193
group hydrogen bonds, donor-hydrogen-acceptor angle = 30◦ , donor-acceptor cut off radius = 0.35 nm, donor = OH, acceptor = O). The rotational autocorrelation functions was calculated using g rotacf (autocorrelation function of cross product of two vectors, n = ij × jk, where i = C3 , j = C12 and k = C15 of steroid backbone, see Fig. 1). The self diffusion coefficient was calculated with g msd. g msd calculates the mean square displacement of atoms/molecules from their initial positions, fits a straight line using the least squares method to the mean square displacement and computes the self diffusion constant using the Einstein relation. 2.2. Bile salt aggregation simulation specifications The united atom molecular simulations utilised the GROMACS united atom forcefield (ffgmx) [26], which is based on the GROMOS-87 force field [30]. A Ryckaert-Bellemans potential [31] was used for the CH2 –CH2 dihedral in the bile acid’s hydrophobic chain and the taurine amino acid group. The Ryckaert-Bellemans potential was required because the energy difference between anti and gauche confirmations using the ffgmx dihedral was too high. The energy difference between the anti- and gauche confirmations of a butane molecule were calculated for the ffgmx and Ryckaert-Bellemans CH2 –CH2 dihedrals as 6.6 and 2.9 kJ mol−1 , respectively. The energy difference between the anti- and gauche confirmations of a butane molecule using the Ryckaert-Bellemans (2.9 kJ mol−1 ) compares much more favourably with the butane ab initio value [32] of (2.6 kJ mol−1 ) than ffgmx (6.6 kJ mol−1 ). Rigid SPC water [33] was utilised as the water model and was constrained using SETTLE [34]. Remaining solute bonds were constrained by the LINCS algorithm [35] and temperature and pressure control was executed by Berensden coupling [36]. Periodic boundary conditions were employed on the cubic simulation box, along with a cut-off distance of 1.5 nm for electrostatic and 1.0 nm for Lennard-Jones non-bonding interactions. These cut-off distances have been found to be suitable in simulations of a bile salt and phospholipid system [18] and were used for the derivation of the force field. To achieve a large simulation time step, the fastest vibrational frequency within the bile salt molecules, exhibited by the polar hydrogens, was reduced by increasing the mass of the hydrogen atom to 4 Da. The attached heavy atom’s mass was decreased by the same amount to conserve mass [37]. The mass redistribution was only applied to the polar hydrogens within the bile salt molecules. This allowed the use of a 4 fs time step, with such large time steps being used increasingly in simulations [16,18,38–41]. All of the bile salt aggregation MD simulations were established and performed using the following procedure. Thirty-one bile salt molecules were randomly distributed within an 8 nm cube, giving a bile salt concentration of approximately 0.10 M. The box was then filled with the appropriate number of water molecules (approximately 1.6 × 104 ). A total of 77 sodium ions and 46 chloride ions were subsequently distributed by random replacement of water molecules, giving 31 counter-ions and a background salt concentration of approximately 0.15 M. The total number of atoms within all simulations was in the region
185
of 50,000. Energy minimisation, using the steepest descent technique [26], was used to remove any bad van der Waals contacts between atoms. This was followed by two short 20 ps simulations. The first short simulation was performed with tight temperature (τ T = 0.01 ps with reference temperature of 310 K) and no pressure coupling, and the second with both tight temperature (τ T = 0.01 ps) and isotropic pressure (τ P = 0.2 ps with reference pressure of 1 bar and 4.5 × 10−5 compressibility) coupling. From this second trajectory and coordinate files the full MD simulation was commenced. A time step of 4 fs for the MD simulation was used, along with Pref = 1 atm, τ P = 0.2 ps, Tref = 310 K and τ T = 0.05 ps. The bulk system properties were stable for the entire length of the simulations, oscillating around the average value. The total energy varied by <1%, while the temperature was maintained at 312 ± 2 K and the density at 990 ± 4 g dm−3 . The density compares favourably with that of pure water, 993 g dm−3 at 37 ◦ C [42]. Simulations continued for 10 to 15 ns after the system had reaching an equilibrium state based on the aggregation and solvent accessible surface area behaviour. 2.3. Forcefield comparison simulation specifications Five smaller scale simulations were performed to validate the utilisation of a united atom forcefield with heavy hydrogens, electrostatic cut-offs and large time step. All forcefield comparison simulations were performed with 6 CHD molecules, 6 sodium counter ions, 1,280 water molecules, started with a 3.5 nm simulation cube and run for 100 ns. Three different united atom forcefield simulations were performed using; 1.5 nm cut-off/4 fs, 1.5 nm cut-off/2 fs, and PME/2 fs. OPLS [43] was utilised as the all atom forcefield with TIP3P [44] the water model using; 2 fs time step, cut-off distance of 0.9 nm for electrostatics and Lennard-Jones non-bonding interactions, and particle-mesh Ewald (PME) [45] method for long range electrostatic interactions. A second OPLS simulation was performed with 1.5 nm cut-off for the electrostatics instead of PME. 3. Results and discussion 3.1. Forcefield comparison The aggregation behaviour simulations performed in this study use the ffgmx united atom forcefield with a 1.5 nm electrostatic cut off and 4 fs time step. To validate this method a series of five simulations were performed, see Table 2. All five simulations rapidly associated within 2 ns, forming a single aggregate of all 6 bile salt molecules that remained intact for the remaining simulation time. Self diffusion coefficients of water and sodium ions were calculated over the entire 100 ns trajectories, with the results summarised in see Table S1, Supplementary Information. Table S1 also contains experimental and molecular dynamics simulation literature values for water and sodium ions. The self diffusion coefficients obtained in this study (SPC 3.87 × 10−5 cm2 s−1 and TIP3P 5.73 × 10−5 cm2 s−1 ) are con-
186
D.B. Warren et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 280 (2006) 182–193
Table 2 Specifications for molecular dynamics simulations performed for method validation Water model
Time step (fs)
Electrostatic treatment
SPC SPC SPC TIP3P TIP3P
4 2 2 2 2
1.5 nm cut-off 1.5 nm cut-off PME 1.5 nm cut-off PME
sistent with literature values under similar MDS conditions (SPC 3.85 × 10−5 cm2 s−1 and TIP3P 5.19 × 10−5 cm2 s−1 ) [46]. Both SPC and TIP3P water models are known to overestimate the diffusivity of water [47], which is 2.30 × 10−5 cm2 s−1 at 25.2 ◦ C and 3.89 × 10−5 cm2 s−1 at 50.2 ◦ C [48]. Increasing the time step from 2 to 4fs for the united atom forcefield only has a minor effect on the self diffusion coefficients, increasing the water and sodium ions values by approximately 5 and 15%, respectively. The spatial distribution functions (SDF) of the water oxygen atoms around the CHD molecules were calculated over the 5 to 100 ns simulation period for the united-atom/1.5 nm cut-off/4 fs and all-atom/PME/2 fs simulations, see Fig. S1, Supplementary Information. The shape and extent of the probability surfaces were similar, with only minor variations. These minor variations are due to differences in the dihedral energy functions of both forcefields and the generally more irregular molecular surface generated by an all atom versus a united atom forcefield. Radial distribution functions (RDF) of oxygen-oxygen, oxygen-hydrogen and hydrogen-hydrogen for the water molecules were calculated for each simulation. Alterations to the electrostatics treatment for both forcefields caused no drift in the location of the RDF peaks, but did increase the peak height or ordering. Increasing the time step to 4 fs and using 1.5 nm cutoff for electrostatics displayed no significant alteration in the RDFs, as displayed by the O-O RDF in Fig. S2, Supplementary Information. Differences in the height and location of the second solvation shell are observed between the two water models, illustrated by the O-O RDF shown in Fig. S3, Supplementary Information. The difference in the RDF around the second solvation shell is due to the water models (SPC versus TIP3P), as has been observed previously [49], rather than the simulation parameters. The average CHD rotational autocorrelation functions were calculated for all five simulations and then fitted with a biexponential function, producing two correlation times. These correlation times are summarised in Table S2, Supplementary Information. The differences between the CHD rotational autocorrelation functions for the united-atom/1.5 nm cut-off/4 fs and all-atom/PME/2 fs are shown in Fig. 3. The two correlation times for these curves are 1.02 and 6.67 × 10−6 ns−1 , and 1.54 and 1.52 × 10−6 ns−1 , respectively. This indicates that micelles simulated using the all atom forcefield are more dynamic on the 500 ps timescale than the united atom forcefield. The effect of increasing the time step on the united atom/1.5 nm cut-off system from 2 to 4 fs is to double both the short and long correlation
Fig. 3. Average rotational autocorrelation function of the 6 cholate molecule simulations: (—) united atom and (— —) all atom forcefield.
times. This increase in correlation time corresponds to a decrease in the rotational freedom of the molecules within the aggregates. The variation in rotational autocorrelation function correlation time simulation conditions is accounted for by considering only the relative values in analysis within this study, not the absolute value. 3.2. Spontaneous aggregation To investigate bile salt aggregation, six bile salt systems were simulated, with the total simulation time performed for each being; CHD (20 ns), GCH (45 ns), TCH (55 ns), GCD (40 ns), GDX (30 ns), and GLC (25 ns). Spontaneous aggregation was observed in all cases. This process is illustrated using GCH in Fig. 4. The starting configuration is presented in Fig. 4(i). After 20 ns (Fig. 4(iii)) all but two of the GCH molecules have aggregated to form a single micelle. At the conclusion of the simulation (Fig. 4(iv)) the free molecules have encountered the micelle and subsequently caused it to split into two separate micelles. The aggregation process was monitored by calculating the number of aggregates and the aggregation numbers over time (shown in Fig. 5). A molecule was defined as belonging to an aggregate when two carbon atoms were within 0.5 nm of a neighbouring molecule’s carbon atoms. This simple rule provided a low noise level and classified only closely associated bile salts as part of an aggregate. The initial aggregation behaviour for all of the bile salts was similar, with spontaneous aggregation complete in approximately 10 ns. For the bile salts CHD, GDX, GCD and GLC, initial aggregation was followed by molecular reorientation to produce a relatively stable, dynamic state. GCH and TCH exhibited more interesting phenomena. Within 16 ns the majority of
D.B. Warren et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 280 (2006) 182–193
187
Fig. 4. Evolution of the glycocholate simulation at time (i) t = 0 ns; (ii) t = 5 ns; (iii) t = 20 ns; and (iv) t = 45 ns. Solid lines indicate the edge of the simulation box. Also displayed are the ions as points and water molecules are omitted for clarity.
the GCH molecules became associated to form a single aggregate, see Fig. 5(ii). After an additional 16 ns, the last two free molecules were absorbed to form a large prolate micelle containing all 31 molecules, which was unstable and short lived. After 2 ns of reorganisation, this micelle split into two approximately equal sized micelles of 15 and 16 molecules, which were then stable for the remaining simulation time. For TCH the aggregation process was far more dynamic, see Fig. 5(vi). Rapid aggregation occurred during the first 10 ns and by 20 ns two micelles containing 10 molecules were present, with the remainder left as free monomers or dimers. Aggregates indicated by the figure to contain around 15 or more molecules were due to transient interactions between the hydrophilic chains of TCH molecules in adjacent aggregates, not between the steroidal backbones. For the remainder of the simulation two micelles containing approximately 10 molecules were maintained, with relatively rapid exchange of monomers and some smaller transient aggregates forming from the remaining molecules. The size and number of aggregates present at the conclusion of the simulations are shown in Table 3, with experimental aggregation numbers provided for comparison. It should be noted that definitive conclusions on aggregate sizes cannot be formed based on the very small data sets presented here. The simula-
tions performed are insufficient in size and time frame to gain an acceptable statistical view of the aggregate size and shape. To be conclusive, significantly larger systems would be required. Moving to a system containing in the order of three hundred bile salt molecules and a simulation time on the microsecond time scale may be sufficient. Nevertheless, provided the limitations of the simulations are considered, these simulations can give initial insight into the behaviour of bile salt micelles on the molecular scale. Table 3 Aggregate size of bile salts from this study and published experimental aggregation numbers Bile salt
Simulation aggregate size
Literature aggregation number
CHD GCH TCH GCD Taurocheno-deoxycholate GDX Taurodeoxy-cholate GLC
31 16, 15 11, 10, 10 × 1 16, 14, 1 – 17, 11, 3 × 1 – 13, 10, 8
15–30 [50], 3.09–3.35 [2] – 3 [51], 2.59 [2] – 2051 – 2451 –
188
D.B. Warren et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 280 (2006) 182–193
Fig. 5. Aggregation behaviour of the molecular dynamics simulation of the bile salts: (i) CHD; (ii) GCH; (iii) GDX; (iv) GCD; (v) GLC, and (vi) TCH. Bubble size is proportional to the number of molecules within aggregates of that aggregation number.
D.B. Warren et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 280 (2006) 182–193
189
Fig. 6. Solvent accessible surface of the various bile salt micelles at the simulation end: (i) CHD; (ii) GCH; (iii) GDX; (iv) GCD; (v) GLC, and (vi) TCH. Surface probe radius = 0.14 nm and each colour indicates a different molecule. (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of the article.)
190
D.B. Warren et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 280 (2006) 182–193
3.3. Aggregate structure The instantaneous arrangements of molecules at the end of each simulation are shown as solvent accessible surfaces in Fig. 6. In general, the bile salts make contact through the steroidal backbone with the hydrophilic chain exposed to the solvent. It is evident that the arrangement of molecules and micellar shapes differ from the three proposed literature models. Furthermore, the micelles are dynamic. For all of the bile salts, the shapes of the micelles were fluid and varied from prolate, through roughly spherical to oblate. It was rare to find an entire hydrophobic face in contact with the solvent for a significant length of time. The hydrophilic face was predominately in contact with the solvent or the hydrophilic face of an adjacent bile salt. The principle radius of gyration for each of the micelles was calculated (atoms explicitly mass weighted) over the final 10 to 15 ns of the trajectories, shown in Table 4. Micelles ranged in size from a long axis radius of 1.31 nm for a 16 molecule GCH micelle, down to 1.00 nm for an 11 molecule GDX micelle. The
Table 4 Aggregation number and principle radii of gyration of MD simulation bile salt micelles Bile salt
Aggregation number
rx (nm)
ry (nm)
rz (nm)
GCH
16 15
1.31 1.08
1.27 1.02
0.86 0.91
GDX
16 11
1.25 1.00
1.19 0.94
0.94 0.83
GCD
16 14
1.26 1.12
1.19 1.06
0.91 0.89
GLC
13 10 8
1.15 1.05 1.11
1.07 0.97 1.02
0.96 0.85 0.81
Balance of molecules (total of 31 per system) are present as un-associated molecules.
principle radii indicate that the average shape of the micelles is oblate, a feature of the model of Kawamura et al., with rx (1.00–1.31 nm) and ry (0.94–1.27 nm) being approximately the same value and the third radius, rz (0.81–0.96 nm), being smaller.
Fig. 7. Four favoured arrangements of GCH molecules within a GCH micelle. Atom colouring is; cyan, carbon; red, oxygen; white, polar hydrogen; and blue, nitrogen. (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of the article.)
D.B. Warren et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 280 (2006) 182–193
An important attribute of the bile salts within the simulated micelles was that the molecules are not arranged parallel to their neighbours, a characteristic of both the disk shape models. Rather, the hydrophilic chain attached to the steroidal backbone maximised the distance from neighbouring chains, allowing maximal hydration and forcing an offset orientation between the molecules. Only a small proportion of parallel alignment occurred. In this case, the molecules were usually aligned with concurrent OH–OH and terminal-terminal chain hydrogen bonding. Parallel alignment involving the hydrophobic face occurred only as transient structures; a result of steric interference from the two methyl groups. A distribution function was calculated for the angle between neighbouring molecules (data not shown), but this showed no substantial ordering, indicating that there is only weak long-range ordering within the micelle. Furthermore, this suggests that the rigid, ordered arrangements of the disk shaped models are overly idealistic. A comparison of bile salt simulations showed that intermolecular H-bonding of the steroid hydroxyl groups was a major factor influencing micelle size, structure and dynamics. It is evident that bile salts with two (GDX and GCD) or three (CHD, GCH and TCH) hydroxyl groups form H-bonds within the micelle in the absence of water. This H-bonding allows the molecules to be buried deeper within the micellar structure, away from the bulk solvent, and to form larger and more spherical micellar structures. In contrast GLC, which has hydrophilic groups located at opposite ends of the molecule, packs differently and adopts cylindrical or disk shaped aggregates to satisfy the interactions required by these hydrophilic groups. This type of stacking reduces the maximum aggregation number for GLC, as shown by the lower aggregation numbers in Table 3. Representative examples of both hydrophilic and hydrophobic intermolecular interactions are shown for GCH in Fig. 7. These examples were obtained by aligning the neighbouring bile salt molecule to high probability regions of the SDF. It was apparent that a large number of alignments are possible, reinforcing the observation that these micelles are dynamic. The hydrophilic (face to face) positioning maximises the number of H-bonds between hydroxyl groups, see Fig. 7(i) and (ii). The hydrophobic alignments have the two methyl groups interlocked to maximise contact, see Fig. 7(iii) and (iv). Similar alignments are seen with the other bile salts but, as the hydroxyl group number decreases, the number of H-bonds reduces, providing more freedom of movement. For the tri-hydroxyl GCH, three anchoring points exist, restricting the molecular movement. In the case of the di-hydroxyls, GDX and GCD, only two points exist, allowing the rest of the molecule to rotate around these more freely. The mono-hydroxyl bile salt, GLC, with only a single anchoring point, has significantly larger conformational freedom while maintaining the single H-bond. The micelle model of Carey and Small [1] proposes that H-bonding between the primary micelles occurs, allowing the formation of secondary micelles that are observed at high bile salt concentrations. Such intermicellar interactions were observed in our simulations and were most evident for GLC and TCH. These interactions between micelles occurred through Hbonding of the hydrophilic tails, not between steroid backbone
191
hydroxyl groups. This was particularly evident in the GLC simulation which separated into three micelles containing 8, 10 and 13 molecules that remained in contact through the hydrophilic tails. The unusual proposed micellar structures of Giglio et al. [7], with the hydrophobic surface exposed to the aqueous phase, are far removed from the observed structures. The proposal that significant hydrophobic area is in contact with the surrounding solvent at low bile salt concentrations is counter-intuitive. The energy reward of having the hydrophilic groups arranged in a helical structure around the counter ions would have to be greater than the energy cost of water structuring from contact with the large hydrophobic surface. Such an arrangement might be possible at high bile salt concentrations, where phase inversion might occur. However, from these simulations it appears to be an unlikely arrangement at micellar phase concentrations. 3.4. Aggregate dynamics The dynamic behaviour of the bile salt aggregates was investigated using a variety of autocorrelation functions. The most insight was obtained from the H-bond lifetime autocorrelation function of the steroid backbone –OH groups, shown in Fig. 8, which exhibited an increasing decay rate as follows: GCH, GCD = CHD, GDX, GLC and TCH. A bi-exponential function was fitted to the autocorrelation functions, with the short correlation time process decaying within 100 ps and showed no correlation to any molecular parameter. Increasing long correlation time is observed across the series of bile salts with the same headgroup (Fig. 9); GCH, GDX, GCD, and GLC. An inverse relationship was established between the long correlation time and the average number of hydroxyl H-bonds per H-bond accep-
Fig. 8. Average hydrogen bond lifetime autocorrelation function of the 31 bile salt molecule simulations: (—) CHD; (···) GCH; (– –) GDX; (– · ·) GCD; (— —) GLC and (– ·) TCH.
192
D.B. Warren et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 280 (2006) 182–193
autocorrelation function correlation time and the degree of hydroxylation. Intermolecular hydrogen-bonding is a major factor influencing size, structure and dynamics of the micelles. This is consistent with the aggregation behaviour of the systems. Intermicellar interactions similar to those proposed by Carey and Small were observed, which could allow formation of larger, secondary micelles. Although the structures observed in these MD simulations possess some components of those previously proposed by Carey and Small1 , and Kawamura et al. [6], they reveal more disordered and dynamic behaviour than that provided by these simple models. Acknowledgements Financial support from Cardinal Health is gratefully acknowledged. We would also like to thank the Australian Partnership for Advanced Computing (APAC) for the Merit Allocation Scheme grant and the Victorian Partnership for Advanced Computing for the provision of computing resources. Appendix A. Supplementary data Fig. 9. Relationship between the correlation time of the H-bond lifetime autocorrelation function and the number of OH H-bonds per number of acceptors: (䊉) CHD, () TCH, () GCH, () GDX, ( ) GCD, and () GLC.
tor, see Fig. 9. Bile salts with lower degrees of intra-micellar Hbonding are able to rotate more freely, and therefore have higher long correlation times. The precise geometry of the hydroxylation is also clearly important to the dynamics, as exhibited by GDX and GCD which have a second hydroxyl group in the 12␣ and 7␣ positions, respectively, and are markedly different in the degree of H-bonding. The effect of the hydrophilic chain on the dynamics is made apparent by the series CHD, GCH and TCH. The addition of the glycine group (CHD to GCH) decreases the correlation time, indicating a restriction in molecular movement. Conversely, the addition of taurine (CHD to TCH) increases the correlation time significantly, indicating a more dynamic system. This observation is consistent with the aggregation process for these two species shown in Fig. 5. The difference between the effect of taurine and glycine might be due changes in the headgroup interactions (steric, hydration and electrostatic); however, the cause requires further investigation. 4. Summary and conclusions MD simulations of six bile salt types were performed to observe their aggregation behaviour in the aqueous phase. The micelles formed were found on average to be oblate with the molecules arranged to maximise hydration of the hydrophilic chain. Contact between the steroidal backbones was mediated through hydrophobic interactions and through hydrogen bonding of hydroxyl groups. Several more common neighbouring molecular orientations were identified but overall there was no favoured angle. The parallel orientations of neighbouring bile salt molecules exhibited by two of the literature models are unlikely. The bile salt micelles were found to be dynamic with an observed relationship between the hydrogen bond lifetime
Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.colsurfa.2006.02.009. References [1] M.C. Carey, D.M. Small, Arch. Intern. Med. 130 (1972) 506. [2] A. Coello, F. Meijide, E.R. Nunez, J.V. Tato, J. Pharm. Sci. 85 (1996) 9. [3] J.P. Kratohvil, W.P. Hsu, M.A. Jacobs, T.M. Aminabhavi, Y. Mukunoki, Colloid Polym. Sci. 261 (1983) 781. [4] C.J. O’Connor, R.G. Wallace, Adv. Colloid Interf. Sci. 22 (1985) 1. [5] A.F. Hofmann, Biliary secretion and excretion: the hepatobiliary component of the enterohepatic circulation of bile acids, in: L.R. Johnson (Ed.), Physiology of the Gastrointestinal Tract, 2, third ed., Raven Press, New York, 1994, p. 1555. [6] H. Kawamura, Y. Murata, T. Yamaguchi, H. Igimi, M. Tanaka, G. Sugihara, J.P. Kratohvil, J. Phys. Chem. 93 (1989) 3321. [7] E. Giglio, S. Loreti, N.V. Pavel, J. Phys. Chem. 92 (1988) 2858. [8] D.W.R. Gruen, J. Phys. Chem. 89 (1985) 146. [9] B.F. Abu-Sharkh, E.Z. Hamad, Langmuir 20 (2004) 254. [10] M.C. Woods, J.M. Haile, J.P. O’Connell, J. Phys. Chem. 90 (1986) 1875. [11] J.-B. Maillet, V. Lachet, P.V. Coveney, Phys. Chem. Chem. Phys. 1 (1999) 5277. [12] H. Shinto, S. Morisada, M. Miyahara, K. Higashitani, Langmuir 20 (2004) 2017. [13] S. Bogusz, R.M. Venable, R.W. Pastor, J. Phys. Chem. B Mater. Surf. Interf. Biophys. 104 (2000) 5462. [14] C.D. Bruce, M.L. Berkowitz, L. Perera, M.D.E. Forbes, J. Phys. Chem. B Mater. 106 (2002) 3788. [15] A.R. Rakitin, G.R. Pack, J. Phys. Chem. B Mater. 108 (2004) 2712. [16] D.P. Tieleman, D. van der Spoel, H.J.C. Berendsen, J. Phys. Chem. B Mater. Surf. Interf. Biophys. 104 (2000) 6380. [17] S.J. Marrink, A.E. Mark, J. Am. Chem. Soc. 125 (2003) 15233. [18] S.J. Marrink, A.E. Mark, Biochemistry 41 (2002) 5375. [19] L.D. Schuler, P. Walde, P.L. Luisi, W.F. van Gunsteren, Eur. J. Biochem. 30 (2001) 330. [20] S.W. Chiu, E. Jakobsson, R.J. Mashl, H.L. Scott, Biophys. J. 83 (2002) 1842. [21] E. Lindahl, O. Edholm, J. Chem. Phys. 115 (2001) 4938. [22] M. Bachar, P. Brunelle, D.P. Tieleman, A. Rauk, J. Phys. Chem. B Mater. 108 (2004) 7170.
D.B. Warren et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 280 (2006) 182–193 [23] A.M. Smondyrev, M.L. Berkowitz, Biophys. J. 80 (2001) 1649. [24] H.J.C. Berendsen, D. van der Spoel, R. van Drunen, Comput. Phys. Commun. 91 (1995) 43. [25] E. Lindahl, B. Hess, D. van der Spoel, J. Mol. Model. 7 (2001) 306. [26] D. van der Spoel, A.R. van Buuren, E. Apol, P.J. Meulenhoff, D.P. Tieleman, A.L.T.M. Sijbers, B. Hess, K.A. Feenstra, E. Lindahl, R. van Drunen, H.J.C. Berendsen, GROMACS User Manual Version 3.1.1, Department of Biophysical Chemistry, Groningen, 2002. [27] W. Humphrey, A. Dalke, K. Schulten, J. Mol. Graph. 14 (1996) 33. [28] A. Varshney, F.P. Brooks, W.V. Wright, IEEE Comput. Graph. Appl. 14 (1994) 19. [29] C. Freudenberger, g sdf; 1.25 ed.; Chemistry, University of Ulm, 2003. [30] W.F. van Gunsteren, H.J.C. Berendsen, Gromos-87 manualGroningen, 1987. [31] J.P. Ryckaert, A. Bellemans, Faraday Discuss. Chem. Soc. 66 (1978) 95. [32] N.L. Allinger, J.T. Fermann, W.D. Allen, H.F. Schaefer III, J. Chem. Phys. 106 (1997) 5143. [33] H.J.C. Berendsen, J.P.M. Postma, W.F. van Gunsteren, J. Hermans, Interaction models for water in relation to protein hydration, in: B. Pullman (Ed.), Intermolecular forces, D. Reidel Publishing Company, Dordrecht, 1981, pp. 331–342. [34] S. Miyamoto, P.A. Kollman, J. Comput. Chem. 13 (1992) 952. [35] B. Hess, H. Bekker, H.J.C. Berendsen, J.G.E.M. Fraaije, J. Comput. Chem. 18 (1997) 1463.
193
[36] H.J.C. Berendsen, J.P.M. Postma, W.F. Gunsteren, A. Dinola, J.R. Haak, J. Chem. Phys. 81 (1984) 3684. [37] K.A. Feenstra, B. Hess, H.J.C. Berendsen, J. Comput. Chem. 20 (1999) 786. [38] G. Colombo, S.-J. Marrink, A.E. Mark, Biophys. J. 84 (2003) 2331. [39] S. Bandyopadhyay, J. Chanda, Langmuir 19 (2003) 10443. [40] S.J. Marrink, E. Lindahl, O. Edholm, A.E. Mark, J. Am. Chem. Soc. 123 (2001) 8638. [41] S.J. Marrink, D.P. Tieleman, Biophys. J. 83 (2003) 2386. [42] G.F.C. Rogers, Y.R. Mayhew, Thermodynamic and Transport Properties of Fluids, fourth ed., Basil Blackwell, Cambridge, 1991. [43] W.L. Jorgensen, D.S. Maxwell, J. Tirado-Rives, J. Am. Chem. Soc. 118 (1996) 11225. [44] W.L. Jorgensen, J. Chandrasekhar, J.D. Madura, R.W. Impey, M.L. Klein, J. Chem. Phys. 79 (1983) 926. [45] U. Essmann, L. Perera, M.L. Berkowitz, D. Tom, H. Lee, L.G. Pedersen, J. Chem. Phys. 103 (1995) 8577. [46] M.W. Mahoney, W.L. Jorgensen, J. Chem. Phys. 114 (2001) 363. [47] B. Guillot, J. Mol. Liquids 101 (2002) 219. [48] K. Krynicki, C.D. Green, D.W. Sawyer, Faraday Discuss. Chem. Soc. 66 (1978) 199. [49] W.L. Jorgensen, C. Jenson, J. Comput. Chem. 19 (1998) 1179. [50] B. Lindman, N. Kamenka, H. Fabre, J. Ulmius, T. Wiellock, J. Colloid Interf. Sci. 73 (1980) 556. [51] N.A. Mazer, M.C. Carey, R.F. Kwasnick, G.B. Benedek, Biochemistry 18 (1979) 3064.