Kinetics of formation of bile salt micelles from coarse-grained Langevin Dynamics simulations

As mentioned in the main text, investigating the mechanism(s) of micelle formation requires simulations where the number of fission and fusion events is sufficiently high to allow proper sampling of the configuration space of the system. In the main text we show that we achieve sufficient sampling by comparing the micelle size distribution obtained in our NVT simulations to that obtained using Grand Canonical Langevin Dynamics. Here we show further evidence indicating that our simulations achieve sufficient sampling: we demonstrate that individual molecules sample a wide range of micellar environments and that the reaction rates calculated from the micelle size distribution are identical to those directly measured in the simulation.


Introduction
Bile salts and bile acids are surfactants that play a key role in the digestion of fats by humans.2][3][4][5][6][7][8] After solubilization takes place, DMMs carry nutrients to the intestinal wall where they are released and absorbed by intestinal cells. 2,3Bile salts are also used in the formulation of both hydrophilic and hydrophobic drugs, because of their ability to enhance drug transport across tissues. 7,9][15][16] The molecular structure of bile salts differs from those of well-studied surfactants like SDS or the DMM component phosphatidylcholine.Instead of the hydrophilic head-hydrophobic tail (head-tail) structure that leads to a prolate shape typically associated with surfactants, bile salts have a rigid steroid group with four rings attached to a short and flexible tail, as shown in Fig. 1(a) for the bile salt taurochenodeoxycholate.The hydrophilic character of one of the steroid faces results from the presence of two or three hydroxyl groups, which are absent from the hydrophobic face.In the upper intestinal tract, the bile salt tail predominantly ends in a hydrophilic glycine or taurine group which, depending on pH, may be charged.Further down in the intestinal tract the taurine or glycine conjugation is removed and the bile salts become more hydrophobic.8][19][20][21][22][23][24][25][26][27] The low average aggregation numbers of pure bile micelles compared to those of canonical surfactants like SDS stem from the bile salts' unusual molecular structure.This structure enables them to form oligomers (with aggregation numbers, n = 2, 3, 4. .][34] While the static properties of pure bile micelles are reasonably well characterized as summarized above, considerably less is known about their mechanism of formation.Understanding this mechanism is important because the rate and extent to which lipophilic nutrients are solubilized in DMMs are directly related to the rate and extent of uptake by the body. 35Empiric knowledge about the role of DMMs in nutrient uptake has been used for several years to influence the solubilization of particular nutrients or drugs by DMMs for nutritional or health benefits, with only partial success. 2,3,5Further progress has been hampered in part by our lack of detailed knowledge of the molecular scale mechanism of formation of DMMs.To gain such knowledge, we must first understand the solution behavior of the individual components of DMMs.Bile salts are the least understood of all DMM components, so here we focus on them. 7][38][39][40][41] For ionic surfactants with the typical head-tail geometry like SDS, theory and simulation 38,[40][41][42] suggest that micelle size changes occur primarily via the addition and removal of monomers, and that addition and removal of dimers, trimers and other small fragments should play a very small role, consistent with the expected very low concentration of these fragments in solution.Only at high counterion concentrations should fission and fusion of fragments become frequent.In contrast with this canonical picture, Ekwall and Small propose a two-step mechanism of micelle formation for bile salts, made possible by the abundance of bile salt oligomers in solution even at relatively low counterion concentrations: above the CMC, bile salts initially form smaller micelles (containing between 2 and 9 molecules) by stepwise addition; 17,43 larger (called secondary) micelles are then assumed to form by fusion of the primary ones via intramicellar hydrogen bonds.5][46] At present, direct molecular scale information on the mechanism of formation of pure bile micelles is unavailable, so it remains controversial whether bile micelles form primarily via addition and removal of monomers or of oligomers.
In this paper we address this on-going controversy by performing Langevin dynamics simulations of conjugated dihydroxy bile salts using a coarse-grained model developed by us. 31The model satisfactorily reproduces the average aggregation number of bile salts determined from experiment as a function of concentration and the molecular scale structure of bile micelles determined from simulations using all-atom models.The modest computational cost of this model allows for the indispensable long simulations at concentrations close to the CMC, where fission and fusion of both monomers and oligomers can be wellsampled.We identify the dominant mechanisms involved in size change of pure bile micelles, calculate their rate constants, relate this kinetic information to free energy profiles along an appropriate reaction coordinate and relate our observations to the biological role of bile salts.

Approach
Using Langevin dynamics simulations, we investigate the formation of micelles of conjugated dihydroxy bile salts at physiological temperature and NaCl concentration for bile salt concentrations 1.8, 8.9 and 18 times larger than the CMC (0.000442 molecules/s 2 ); these concentrations are denoted as 1.8 CMC, 8.9 CMC and 18 CMC.These concentrations are reported as a multiple of the CMC because different experimental methods lead to CMCs differing by a factor as large as 3 for the same bile system; 31 expressing concentrations as a multiple of the CMC thus facilitates comparisons with experiment.The concentration range used in this study is sufficiently broad to examine the dependence of micelle formation mechanisms on this variable 47 and is physiologically relevant.We define the CMC as the concentration at which 50% of all molecules exist as monomers, a definition which is commonly employed in simulations. 42,48,49onomers are distinguished from aggregates using a distance criterion: 31 two molecules are considered part of the same aggregate if the distance between their centers of mass is o3s, where s is the unit of length in our simulations and is defined in more detail below.Full details of the model used in this study are given in a prior publication 31 so below we offer only a brief description.

Coarse-grained model of bile salts
The coarse-grained model of dihydroxy taurine-or glycinesubstituted bile salts used here is shown in Fig. 1.In the same figure we show the correspondence between the coarse-grained model and the molecular structure of taurochenodeoxycholate; This journal is © The Royal Society of Chemistry 2016 we emphasize, however, that the model is also representative of other conjugated dihydroxy bile salts.
Significant overlap between any two beads is prevented through a Weeks-Chandler-Andersen potential.The water and the salt ions present in aqueous solutions of bile salts are modeled through effective interactions.To mimic the effective attraction experienced by hydrophobic groups in water we use a cosine square potential.This potential has a finite range and has been shown to capture the characteristic short-range attraction of hydrophobic interactions necessary for self-assembly of micelles and bilayers. 50,51We specify no net attraction or repulsion between any two hydrophilic beads or between a hydrophilic and a hydrophobic bead beyond volume exclusion because these interactions are comparable to the interactions of these groups with water, and are much weaker than hydrophobic interactions. 31he only exceptions are the charged beads, which interact with other charged beads through a Debye-Hu ¨ckel potential.This potential accounts for the screening effect brought by the dielectric constant of water and the presence of NaCl on the electrostatic interactions.The model is parameterized to reproduce the critical micellar concentration and average aggregation number of conjugated dihydroxy bile salts near the CMC at physiological temperature and concentration of NaCl but it can easily be adapted to represent bile salts under other conditions, or even molecules similar to bile salts like the triterpenoids asiatic acid and madecassic acid. 52,53

Simulation details
The model is implemented in the ESPResSo molecular dynamics simulation package. 54For computational convenience the implementation is done using reduced units, i.e. we chose suitable units of length, mass and energy and express all quantities as multiples of these units.The unit of mass is m E 8 Â 10 À26 kg, of energy is e = 8.6 Â 10 À21 J = 1.2 kcal mol À1 = 2k B T and of length is s E 4.0 Å as described elsewhere. 31It follows that our unit of temperature is e/k B , where k B is the Boltzmann constant.The unit of time t is estimated by equating the diffusion coefficient of cholate ions (5 Â 10 À10 m 2 s À1 ) 55 to the diffusion coefficient of monomers (0.04s 2 /t) from our simulation at the lowest concentration here tested.We find t E 12 ps.We note that this value represents only a coarse estimate of the unit of time, so in this work we focus on understanding the relative dynamics of the various processes involved in changes of micelle size and the sequence of events involved in such processes.The trends seen in the simulations are thus directly comparable to experiment despite the uncertainty in our estimate of t.
We investigate micellization using Langevin dynamics simulations at constant volume, temperature and number of particles.The initial configuration for each simulation is obtained from a prior simulation (described in ref. 31) performed with Grand Canonical Langevin Dynamics (GCLD).The initial configuration thus already has a micelle size distribution characteristic of equilibrium, and further equilibration is not necessary.Cubic boxes with periodic boundary conditions are used for all simulations.The dimensions of the simulation boxes and the number of molecules used at each concentration are given in Table 1.
All simulations are performed at a temperature of 0.5e/k B , enforced through a Langevin thermostat.The equations of motion are integrated using a velocity Verlet algorithm with a time step of 0.01t.Production runs last 30 000 000 steps, with configurations saved every 100 steps.Saved configurations are visualized using VMD-Visual Molecular Dynamics. 56Results and discussion

Simulations properly sample the system configuration space
Investigating the mechanism(s) of micelle formation requires simulations where the number of fission and fusion events is sufficiently high to allow proper sampling of the system configuration space.We demonstrate that our NVT simulations meet this condition by: (i) verifying that the NVT simulations sample the correct micelle size distribution at each concentration by comparison to simulations using grand canonical Langevin dynamics; (ii) confirming that individual molecules sample a wide range of micellar environments by calculating the mean square displacement in micelle size space and comparing its saturation value with predictions from diffusion theory 57 (see ESI †); and (iii) for the most frequent events -the addition and removal of monomers, dimers or trimers to/from micelles, as shown belowby showing that the rates calculated from the micelle size distribution are identical to those directly calculated in the simulation (see ESI †).
The micelle size distribution is quantified here through the probability, P(n), of finding micelles with n molecules (i.e., with aggregation number n): where N n (or N i ) indicate the number of micelles composed of n (or i) molecules, N mol is the total number of molecules in the simulation, i.e., the maximum micelle size, and the averages are calculated over all saved configurations for each concentration.This probability mass function is shown in Fig. 2 for the three different concentrations investigated here.For comparison, we also show the equivalent distributions obtained with GCLD simulations performed during the course of prior work. 31The distributions obtained with GCLD do not suffer from sampling deficiencies and are thus our reference.The two sets of P(n) distributions are similar, indicating that the current simulations are sufficiently long to sample the full distribution of micelle sizes.Fig. 2 also shows that at the lowest concentration the system is principally composed of oligomers whereas at the highest concentration mostly larger micelles are present.
At the intermediate concentration both oligomers and larger micelles abound.This set of three concentrations is thus adequate to investigate whether oligomers make a large contribution to the mechanism of formation of micelles.

Micelle size changes occur primarily via monomer removal and addition at all concentrations
To evaluate which processes contribute most to changes in micellar size we calculate the reaction rates, l, for addition and removal of monomers and larger fragments.Fission and fusion events are described by quasichemical reactions where A n refers to a cluster with aggregation number n, or a monomer if n = 1.For convenience, in what follows n is the size (i.e.aggregation number) of the initial micelle (or monomer) whereas m is the size of the fissioned or fusioned fragment.The rate constants for fusion (k + ) and fission (k À ) are related to the reaction rates according to the well-known expressions with C n the number density of micelles with aggregation number n.
We calculate the reaction rates as the number of fission or fusion events occurring per unit time and per unit volume.Fission and fusion events are identified by detecting a change in aggregate size between two consecutive saved configurations.The reaction rates calculated this way necessarily depend on the time interval between saved configurations but we expect this dependence to be small because the rate constants are sufficiently low (e.g., the fission rate constant k À , is O(10 À2 t À1 )), that there is only a negligible chance that the same micelle has been involved in more than one fission or fusion event between two saved configurations (separated by 1t).
The fission reaction rates, l À , as a function of the aggregation number of the initial micelle, n, are shown in Fig. 3. Almost identical curves (not shown) are obtained for the fusion reaction rates, l + , as expected for a system in equilibrium.Inspection of Fig. 3 indicates that the fission reaction rates depend strongly on concentration.This dependence follows that expected from eqn (5) assuming a rate constant, k À , that is independent of the concentration.In the ESI † we show that this assumption indeed holds for our systems.Fig. 3 clearly shows that fusion and fission of monomers are by far the most frequent events even at concentrations far from the CMC: the reaction rates for fission and fusion of dimers, trimers etc. are at least one order of magnitude lower than for monomers.Events involving fission or fusion of larger aggregates are nevertheless still frequent.At the lowest concentration investigated (1.8 CMC) these events  represent 10% of all fission and fusion events; at the highest concentration (18 CMC) they represent 25%.Most of the events with m 4 1 are restricted to fission and fusion of dimers or trimers (m = 2, 3): even at the highest concentration, fission and fusion events with m 4 3 represent only 13% of the total number of events.
These results indicate that the dominant mechanism for changes of micellar size is monomer fission and fusion for both large and small micelles, even at concentrations much higher than the CMC.This observation is supported by recent ultrasonic relaxation measurements of micellar kinetics in bile salt solutions. 58In that work, the experimental results were wellfitted by analytical models of micelle kinetics relying solely on monomer addition and removal, indicating that fission and fusion of larger micelles should play a secondary role.The observed importance of monomer addition and removal is in line with simulations of model coarse-grained ionic head-tail surfactants 42 and, more generally, with the expected behavior of ionic surfactants at concentrations near the CMC and for low to moderate counterion concentration. 40Our results thus do not support the mechanism initially proposed by Small and Ekwall, 17,43 that larger micelles form primarily by fusion of smaller ones.In apparent contradiction with our observations and the ultrasonic relaxation measurements 58 we mentioned above, recent conductivity measurements 34 have been interpreted to support the twostage mechanism proposed by Small and Ekwall. 17,43It is our opinion, however, that such a conclusion may not be reached from those measurements.In conductivity measurements, the CMC is estimated from the position of a break point in a graph of conductivity vs. concentration.In the study by Kumar et al., that break point is not sharp, which was interpreted to suggest the presence of oligomers as well as micelles near the CMC.We agree with that interpretation; however, we point out that the mere presence of oligomers does not constitute evidence that larger micelles form via the fusion of small oligomers, as suggested by Kumar et al.

Models of micellar kinetics developed for head-tail surfactants apply for bile salts
7][38] Our results confirm that relaxation times in bile salt solutions at physiologically relevant concentrations may be interpreted using those models, as was recently successfully done by Ravichandran et al., 58 as discussed above.For improved quantitative agreement, however, fission and fusion of dimers and trimers should be included in the models at high concentrations, where these events become more frequent.7][38] In the ESI † we show that rate constants depend weakly on the total bile salt concentration, that the fusion rates are indeed independent of micelle size for all but the smallest aggregates, but that the fission rates depend on micelle size, and we discuss the implications of these findings in the context of the theoretical models used to interpret chemical relaxation times in micellar systems.
3.4 Rate constants and free energy landscapes 3.4.1 Addition and removal of monomers from bile micelles are well-described by a single reaction coordinate: the distance between two fragments.Comparing micellar kinetics between bile salt solutions and solutions of other surfactants is important to understand the biological role of bile salts, but is not easily done between different simulation studies because of the large uncertainty associated with t, as mentioned above in the Introduction.We avoid this difficulty here by obtaining the free energy landscape associated with fission and fusion of micelles and comparing the free energy barriers associated with fission and fusion in bile salt solutions with the equivalent quantities found in the literature for head-tail surfactants.
For head-tail surfactants, simulations 47,59 suggest that addition and removal of a monomer from a micelle is well-described via a single reaction coordinate, the distance, r, between monomer and micelle.The removal process is activated, being associated with large free energy barriers (8k B T, including the contribution of the 2 ln(r) term 47 ).In contrast, monomer addition is essentially barrierless: the rate constants for monomer fusion directly measured in simulations are large compared to the fusion rate for barrierless attachment estimated from the Smoluchowsky model (k + = 4p(R n + R 1 )(D n + D 1 ), where R are the radii and D are translational diffusion coefficients of a monomer or of a micelle of aggregation number n). 59 To assess whether r is a reasonable reaction coordinate for bile salts as well, we compute the free energy, A(r), as a function of the distance between two aggregates, and relate this landscape with the fission and fusion rate constants, shown in the ESI.† We first calculate the probability, P n,m (r), of finding monomers, dimers or trimers (m = 1, 2, 3) around aggregates of size n and then obtain A(r) = Àk B T ln(P n,m (r)).We make the assumption that the free energy profile is not significantly affected by small changes in n, and calculate P n,m (r) averaged over a narrow range of micelle sizes (n = 9-13).To obtain a free energy landscape that gives insight into fission and fusion, the function P n,m (r) must include not only pairs of distinct aggregates with sizes m and n but also pairs of subaggregates of sizes m and n that are a part of a micelle of size n + m; therefore, P n,m (r) is defined as Here N n,m (r) is the number of pairs of aggregates with sizes m and n at separation r between their centers of mass, and N n,m (r) is the number of pairs of subaggregates of size m and n separated by distance r and belonging the same micelle (with size n + m); N n,m and N n,m are the analogous quantities summed over all r.Each micelle with size n + m is arbitrarily divided into a single pair of subaggregates of sizes m and n, so N n,m is equal to the number of micelles of size n + m in the simulation box.We note that, because P n,m (r) does not include a 1/(4pr 2 ) term, A(r) includes the ideal gas contribution to the free energy, i.e., the entropic term 2 ln(r).Fig. 4 shows the free energy profiles as a function of the distance r between micelles and monomers, for different bile salt concentrations.Both fusion and fission of monomers to/from micelles are associated with low free energy barriers.For monomer fusion, the free energy difference between r = 10s and the maximum in A(r) (at r = 5-6s, depending on the concentration) is of order 1-2k B T in all cases.This low barrier is consistent with the fact that the rates we measure are large compared to the fusion rate for barrierless attachment estimated from the Smoluchowsky model.The fission of monomers from micelles is also associated with low free energy barriers: e.g., at 8.9 CMC, A(r = 5.5s) À A(r = 2.5s) E 2k B T. For both fission and fusion events, the free energy barriers have large contributions from the ideal gas term (2 ln(r)).The free energy barrier associated with fission of monomers from micelles increases slightly with increasing concentration -less than 0.5k B T between 1.8 and 18 CMC -an increase which is quantitatively consistent with the observed decrease in monomer fission rates shown in the ESI; † the opposite trend is seen for the free energy barrier of fusion of monomers to micelles, consistent with the larger monomer fusion rates observed at higher concentrations.The reasonable agreement between the trends in fission and fusion rate constants and those in the one dimensional free energy profiles shown in Fig. 4 thus suggests that to understand addition or removal of monomers to/from micelles, it is sufficient to consider only the reaction coordinate r.In the ESI, † we show that this coordinate is also sufficient to describe addition and removal of oligomers to/from micelles.
3.4.2Rate constants of bile salt removal from bile micelles are much larger than for typical head-tail surfactants.Having established that r is an appropriate reaction coordinate to describe fusion and fission of monomers to/from micelles, both for bile salts and for head-tail surfactants, we can compare the kinetics of monomer addition and removal between these two types of systems by comparing their respective free energy barriers associated with fission and fusion.Bile salts behave similarly to head-tail surfactants when it comes to addition of monomers: for both types of molecules, addition of monomers to micelles is associated with negligible free energy barriers (see Fig. 4). 47In contrast, for all concentrations investigated here, the free energy barrier for fission of bile monomers from micelles is only 2k B T, much lower than the equivalent barrier for micelles of head-tail surfactants (8k B T, as mentioned above 47 ).These results suggest that bile salts may leave bile micelles at much faster rates than those typically found for head-tail surfactants.These differences in free energy barriers are qualitatively consistent with differences in the monomer removal rate estimated from experiment: for sodium deoxycholate, 58 k À = O(10 7 s À1 ); for head-tail surfactants 38 with CMC similar to that of sodium deoxycholate, k À = O(10 6 s À1 ); for head-tail surfactants with lower CMC, k À becomes even lower.

Concluding remarks
The coarse-grained model of bile salts used in this work has modest computational cost, enabling the investigation of the mechanisms by which pure bile micelles change aggregation number at concentrations near and far the critical micellar concentration via straightforward Langevin dynamics simulations.1][62] These models may also be used to understand the general trends leading to improved solubilization of drugs and nutraceuticals in artificial DMM-like micelles for oral drug delivery, and the subsequent interaction of these carrier vehicles with DMMs. 7,63][68] Our results indicate that the dominant mechanism of micelle size change in bile salt solutions is monomer addition or removal, even at concentrations 18 times higher than the critical micellar concentration, in agreement with recent ultrasonic studies of the kinetics of bile salt micelles. 58The simulations do not support the popular assumption that bile micelles predominantly form in a two-stage process, by formation of oligomers and subsequent oligomer fusion to produce larger micelles.Addition of monomers to micelles is nearly barrierless, similarly to head-tail surfactants, and the free energy barrier to monomer removal is of order 2k B T, i.e., much lower than the 8k B T found for head-tail surfactants. 47his journal is © The Royal Society of Chemistry 2016 We propose that the low free energy barrier associated with removal of monomers from bile micelles plays a role during digestion.Solubilization of fats by surfactant micelles occurs typically via two mechanisms: 69,70 (i) for fats with sizable solubility in water, micelles can capture fat molecules dissolved in the solution; (ii) when fats are poorly soluble, the surfactant must first adsorb at the surface of the fat droplet, and subsequently a micelle containing both surfactant and fat is formed and desorbs from the fat droplet.This second mechanism, however, is only effective for non-ionic surfactants or ionic surfactants at high ionic strength; when ionic surfactants are present in solutions with low ionic strength, the electrostatic repulsion between surfactant micelles and partially covered fat droplets prevents further adsorption of the micelles.We speculate that the low free energy barrier associated with fission of bile monomers combined with the high CMC of bile salts may increase the effectiveness of solubilization of fats via mechanism (ii): these two properties of bile salts enable them to more rapidly cover the surface of fat droplets, because the electrostatic repulsion between single bile monomers and a fat droplet partially covered with bile salts is much smaller than that between the fat droplet and a micelle.A recent simulation study of the solubilization of a phospholipid vesicle by sodium cholate using a coarse-grained model similar to ours suggests that adsorption of monomeric bile salts to a phospholipid vesicle is indeed fast. 66

Fig. 1
Fig. 1 Correspondence between (a) taurochenodeoxycholate and (b) the coarse-grained model for dihydroxy bile salts.Beads 0 to 6 are hydrophobic, beads 7, 9 and 10 are hydrophilic and bead 8 is charged.

Fig. 2
Fig.2Probability mass function, P(n), of finding micelles composed of n molecules at concentrations 1.8 CMC (black), 8.9 CMC (red) and 18 CMC (green), from Langevin dynamics simulations in the canonical ensemble.For comparison, distributions obtained from grand canonical Langevin dynamics simulations at the same concentrations are also shown (gray).For each concentration, the P(n) distributions obtained with both methods approximately coincide.Notice the change of scale of the y-axis, to better display the probability distribution in the micellar region.

Fig. 3
Fig. 3 Fission reaction rates, l À , as a function of the size of the initial micelle, n, and of the departing fragment, m, for the three different concentrations here investigated.Fission events with n 4 30 or m 4 7 are undersampled so the corresponding reaction rates are not shown.

Table 1
Number of molecules (N mol ) and length (L) of the edge of the cubic simulation boxes used in each simulation