Aggregation Behavior of Medium Chain Fatty Acids Studied by Coarse-Grained Molecular Dynamics Simulation

Medium chain fatty acids (MCFA) are digestion products of lipid-rich food and lipid-based formulations, and they are used as transient permeability enhancers in formulation of poorly permeable compounds. These molecules may promote drug absorption by several different processes, including solubilization, increased membrane fluidity, and increased paracellular transport through opening of the tight junctions. Therefore, understanding the aggregation behavior of MCFAs is important. A number of studies have measured the critical micelle concentration (CMC) of MCFAs experimentally. However, CMC is highly dependent on system conditions like pH, temperature, and the ionic strength of the buffer used in different experimental techniques. In this study, we investigated the aggregation behavior of four different MCFAs using the coarse-grained molecular dynamics (CG-MD) simulations with the purpose to explore if CG-MD can be used to study MCFA interactions occurring in water. The ratio of deprotonated and non-charged MCFA molecules were manipulated to assess aggregation behavior under different pH conditions and within the box sizes of 22 × 22 × 44 nm3 and 44 nm3 for 1 μs. CMCs were calculated by performing CG-MD simulations with an increasing number of MCFAs. The resulting aggregate size distribution and number of free MCFA molecules were used to determine the CMC. The CMCs from simulations for C8, C10, and C12 were 1.8–3.5-fold lower than the respective CMCs determined experimentally by the Wilhelmy method. However, the variation of MCFA aggregate sizes and morphologies at different pH conditions is consistent with previous experimental observation. Overall, this study suggests that CG-MD is suitable for studying colloidal systems including various MCFAs. Electronic supplementary material The online version of this article (10.1208/s12249-018-1289-4) contains supplementary material, which is available to authorized users.


INTRODUCTION
One of the most successful approaches in promoting the absorption of poorly soluble drugs through the intestinal epithelium has been the use of lipid-based enhancers primarily based on medium chain fatty acids (MCFA) [1]. MCFAs can increase the epithelial flux through transient transcellular perturbation as well as paracellular pathways [1][2][3][4][5]. The effect of MCFAs on the intestinal epithelium is transient and rapidly reversible [6,7]. Therefore, MCFAs are ideal components for transient permeability enhancers for peptides, proteins, macromolecules, and oral formulations [1,6]. Transient permeability enhancers have become an essential component of a number of oral formulations including oral peptide and oligonucleotide dosage forms currently undergoing clinical trials [8][9][10][11][12]. MCFA and MCFA derivatives are used in solid dosage forms, and is commonly known as a gastrointestinal permeation enhancement technology (GIPET) [11]. It has been tested with various doses with poorly absorbed drugs in 16 Phase I studies for which oral bioavailability of 5-13% was observed when GIPET dosage forms were used in comparison to the bioavailability of less than 1% that is typically achieved [11]. Furthermore, MCFAs are common in food products and may be generated from glycerides during lipid digestion. No significant toxicity effect Guest Editor: Sanyog Jain Electronic supplementary material The online version of this article (https://doi.org/10.1208/s12249-018-1289-4) contains supplementary material, which is available to authorized users. has been attributed to the use of these lipid-based transient permeability enhancers, even with repeated dosing [8,11,13].
In aqueous solution, fatty acid molecules are able to selfassemble and micelles of various sizes and structures can form depending on fatty acid concentration, which may affect the efficiency of MCFAs as transient permeability enhancers, because free fatty acid monomers present near the intestinal absorption site can directly impact the extent of transient permeability [1]. Therefore, it is important to estimate the number of fatty acid monomers in both free and aggregate forms at a given concentration in order to identify the threshold concentration required to enhance absorption. Typically, the critical micelle concentration (CMC) is considered a good criterion for threshold concentration [1,14]. However, the aggregation behavior and CMC value of different MCFAs are highly dependent on various system properties such as pH, temperature, and ionic strength [1,14]. A large number of studies have investigated aggregation behavior and determined the CMC values of MCFAs using different experimental techniques such as surface tension measurement [14], the dye micellization method [15], the relative viscosity-based method [16], high resolution ultrasonic spectroscopy [6], ion activity measurement, and conductimetry techniques [17]. One key observation regarding the experimentally measured values is that CMC values also depend on the methods of measurement: the difference between CMC values measured using different techniques at similar or specific system conditions is relatively large.
Molecular dynamics (MD) simulations have shown significant promise for investigating the aggregation behavior of different surfactant molecules [18][19][20]. However, the atomistic all atom (AA) model [21,22] in which each atom is represented individually has limited capabilities for measuring equilibrium properties like CMC and aggregate size distribution, since these methods are generally confined to nanosecond time scales and nanometer length scales [20]. Sammalkorpi et al. [23] estimated the CMC of sodium hexyl sulfate (S6S) using an atomistic MD simulations with a united atom model. In this system, 200 ns was sufficient for the surfactant model to reach an equilibrium state. For MCFAs and surfactants like sodium dodecyl sulfate (SDS), such a short time scale is insufficient for fully equilibrating the system [23]. One approach for extending the box size and time scale of experiments is using coarse-grained (CG) MD simulations [24]. In the CG-MD approach, a number of atoms are represented by a single CG bead, which significantly reduces the number of interaction sites in a simulation [24][25][26]. Therefore, these simulations can be performed for relatively longer times and greater length scales, but at the cost of sacrificing the atomic resolution allowed by the AA model [27]. CMC values determined using CG-MD for different surfactants like dodecylphosphocholine (DPC) [28], and S6S [29], sodium nonyl sulfate [29], and SDS [29,30] were found to have reasonable agreement with experimental values. Results using the dissipative particle dynamics (DPD) also agreed quantitatively with experimental CMC and aggregation number for several surfactants such as octaethylene glycol monooctyl ether, dodecyldimethylamineoxide (DDAO), and N-decanoyl-N-methyl-D-glucamine [31,32]. However, an important limitation of DPD is the limited availability of verified model parameters. In CG-MD, the parameterization issue is greatly reduced by, e.g., the Martini model developed by Marrink et al. [25,26]. In the Martini model framework, a specific molecule can be constructed relatively easily by using the different particle types and interactions described in the model [25]. The model was used to determine the CMC value of H 1 T 4 and tetra(ethylene) glycol monoalkyl ether (C 8 E 4 ) and was found to overpredict the experimental values by a factor of 5 and 4 respectively [20]. However, the model was able to capture the aggregate size distribution and average aggregate size for H 1 T 4 and C 8 E 4 [20], as well as complex morphological changes of the oleic acid aggregates based on the protonation state of the oleic acid head groups [33] and complex phase behavior of oleic acid including the pKa shift from 4.8 in water to 6.5 in small micelle and 6.6 in bilayer [34].
The aim of this study was to investigate CMC and aggregation behavior of four MCFAs (C 8 , C 10 , C 12 , and C 14 ) as a function of protonation state of those respective MCFAs and the ionic strength of the buffer using CG-MD simulations with the Martini model. To verify how closely the CG-MD reproduced the experimental CMC values, we also measured the CMCs of C 8 , C 10 , and C 12 experimentally at different representative system conditions by mimicking our simulation conditions as closely as possible. The aggregation behavior obtained from the simulation was also compared with different experimental observations available in the literature.

CG Model and Simulation
The simulations of MCFA aggregates were run using the CG Martini v2.0 model [26]. In this model, typically two to four heavy atoms are represented by one CG bead. To model C 14 (which has 14 carbon chain with a carboxyl headgroup), we used the available MCFA topology from the cgmartini.nl website, with three C1 beads used to represent the carbon chain and one bead (P4 or Qa) was used to represent the carboxylic head group, either as neutral or with a negative charge. The model of C 10 was obtained from C 14 by excluding one CG bead from the carbon chain, since one CG bead typically represented four heavy atoms. Therefore, the CG model of C 10 included three beads with P4/Qa for the headgroup and two beads for representing the carbon chain. The four-to-one mapping of heavy atoms into beads is both a strength and a weakness of the Martini model. One consequence is that it is not possible to separate, e.g., fatty acids with similarly long carbon chains from each other. In this context, we, however, sought to investigate whether it was possible to circumvent and obtain relevant CMC values. To model C 12 which has 12 carbon chains with a carboxyl headgroup, we used four CG beads similar to C 14 . However, to represent the reduced chain length of C 12 compared to C 14 , the effective inter-particle distance of the beads was reduced from 0.47 to 0.37 nm, and a similar change was made going from C 10 to C 8 . No other changes were made to the force field. The ratio of non-charged and deprotonated molecules was varied to simulate the effect of pH on the MCFA aggregation (Table I).
MD simulations were performed in Gromacs 2016.4 [36] with a 30-fs time step. Isotropic pressure coupling with a reference pressure of 1 bar (1 bar = 100 kPa) was maintained with the Berendsen coupling method. Temperature was set at 37°C for most of the simulation. One set of simulation for C 12 was performed at 50°C for the purpose of mimicking the experimental condition. Non-bonded interactions were set according to the recommendations for the Martini force fields and Gromacs 2016.4 version. For the simulations of C 8 and C 10 aggregates, the length, width, and height of the box was 22 nm, 22 nm, and 44 nm respectively. To represent the concentration of 1 mM of MCFA in the simulation box, 13 fatty acid molecules were needed. Approximately 192,000 water beads were used to fill up the box with water. When MCFA were used in the deprotonated state, an equal amount of positively charged sodium ions were also included to make the system neutral. To represent an ionic strength of 140 mM NaCl, 1820 Na + and Cl − ions were added in the simulation box. For a number of simulations with C 12 and C 14 , a larger box of 44 nm 3 was used, because the aggregation occurs at much lower concentrations and the larger box size was required for reliable calculation of the number of free and aggregated molecules. To determine the CMC of a particular system, a number of simulations with increasing number of MCFA molecules were performed initiated from random dispersion of molecules in the box. Prior to the production run, energy minimization was performed for 10,000 steps using the steepest descent algorithm for each system followed by four short equilibration runs (50,000 steps) with time steps of 1, 2, 5, and 20 fs, respectively.

MCFA Aggregates Determination
To determine the MCFA aggregates, two MCFA molecules were considered to be in the same aggregate if their constituent beads were within a specific distance from each other. This cutoff distance was determined to be 0.6 nm for the Martini model [20]. The cutoff size for the free MCFA molecules was selected from the aggregate size distribution and considered to lie between the first two peaks of the distribution [20]. Typical aggregate size distribution for C 8 and C 10 is shown in Fig. 1a, indicating a cutoff aggregate size of 5 to be reasonable. The aggregate size distribution for C 12 and C 14 is shown in Supplementary Fig. S1. It is possible that the boxes did not include enough molecules for obtaining a complete aggregate size distribution, but Supplementary Fig.   S1 shows a trend similar to that of Fig. 1a, so cutoff size 5 was used for C 12 and C 14 as well.

CMC Definition and Determination
Different definitions of CMC are currently in use, including the concentration at which micelles start to appear, or the concentration at which half of the molecules are in the aggregate form and half are free [20]. This study is based on the second definition. The number of free MCFA molecules  were calculated over the simulation trajectory as shown in Fig. 1b. The fraction of free molecules (averaged over 20 snapshots within 0.9 to 1 μs of the simulation time) was used to determine CMC.

CMC Surface Tension Measurements
The surface tension measurement utilizing the Wilhelmy plate is the most common CMC measurement technique [14,15,37]. This method is also known as the Wilhelmy method. In this method, the surface tension of aqueous solution with surfactants is measured by using a platinum plate placed inside the solution. The surface tension value decreases with increasing surfactant concentration in the solution and becomes constant above the CMC [15].
In our experiment, the surface tension of phosphate buffer containing MCFA was measured on a Sigma 70 instrument to determine CMC. Measurement temperatures were 37°C for C 8 , C 10 , and 50°C for C 12 in order to stay above the Krafft point. Buffer solutions were designed to have ionic strengths of 140 mM. Two different pH values of 7 and 12 were used when measuring the surface tension to mimic the desired system conditions. All measurements were performed in triplicate.

RESULTS
To determine CMC using CG-MD, we first calculated the fraction of free molecules at different concentrations for each MCFA in specific system conditions (Fig. 2). The free fraction of molecules was generated based on 20 snapshots from the simulations, and the average value obtained from these snapshots was used for CMC calculation. Each set of values was fitted with a nonlinear power-law equation which was then used to determine the concentration at which 50% of the molecules are free in the simulation, and this concentration value was taken as the CMC value (see Methods; Table I). The R-squared values for fitting equations were 0.941-0.995.
For each MCFA, the largest CMC value was obtained when all molecules were in the deprotonated state and no other ions were present in the system (Table I and Fig. 2). Including 140 mM NaCl with the deprotonated MCFA molecules decreases the CMC values. The lowest CMC value was observed for all four MCFAs when all the molecules were at the non-charged state in the presence of 140 mM NaCl (note that using all the molecules in the deprotonated and non-charged state mimics high and low pH conditions respectively). CMC values for the systems with 1:1 ratios of deprotonated and non-charged molecules were higher than CMC values for the system with 100% non-charged molecules, but lower than the system with 100% deprotonated molecules, in line with what has been experimentally observed [18]. The CMC values of the different MCFAs decreased with increasing carbon chain length. This relationship was consistent for all four different systems. The CMC values obtained for C 8 and C 12 using the reduced interparticle distance showed good consistency with respect to the values obtained for C 10 and C 14 which used the actual bead inter-particle distance of the Martini model.
We also measured the CMCs of C 8 , C 10 , and C 12 experimentally to validate the exact model configuration used in this study, with particular interest in capturing the effect of decreased bead interaction distance used in the C 8 and C 12 model. These experimentally determined CMC values are in Table I. To determine these CMCs using the Wilhelmy method, the surface tensions as a function of concentration were measured (Fig. 3). The linear parts of the data sets for each MCFA were fitted using linear regression by least square. The corresponding concentration of the intersection point between the two fitted lines was considered to be the CMC.
The experiments performed at pH = 12 with 140 mM NaCl corresponds to the simulated system with 100% deprotonated MCFA molecules with ionic strength. The pKa values for C 8 and C 10 are slightly above 7 [35]. Therefore, the experiments performed at pH = 7 and 140 mM NaCl corresponds to the simulated system with a 1:1 ratio of deprotonated and non-charged fatty acid molecules and NaCl. To compare experimentally determined CMC for C 12 , we performed a new set of simulations at 50°C with 100% deprotonated molecules and 140 mM NaCl. Based on the fraction of free molecules obtained in the new simulations ( Supplementary Fig. S2), the CMC value was calculated to be 1.71 mM (Table I). Figure 4 shows snapshots of simulations of C 10 at different concentrations for the system with 100% deprotonated molecules and 140 mM NaCl. For this system, calculated CMC was 14.89 mM ( Table I). The number of free C 10 molecules is clearly very high in the simulation snapshot (Fig. 4a) of 10 mM C 10 . As the concentration of C 10 increases, the decreasing number of free molecules and increasing number of C 10 aggregates are clearly indicated by the simulation snapshots.
Simulation snapshots of C 10 at different concentration levels for the system with a 1:1 ratio of deprotonated and non-charged molecules and 140 mM NaCl are shown in Fig. 5. The CMC of C 10 for this system was calculated to be 8.27 mM ( Table I). The snapshots demonstrate that the fraction of free molecules decreases as the concentration increases (similar to the 100% deprotonated system). However, this system has fewer aggregates and a larger aggregate size compared to the 100% deprotonated system at the same concentration level. The snapshots with the C 10 concentration higher than CMC (Fig. 5c, d) also indicate that a majority of the free molecules are deprotonated.

DISCUSSION
Previous studies have used the free monomer concentrations from simulations performed with high surfactant concentrations to determine CMCs, on the assumption that the monomer concentration does not vary much above the CMC [28,31], but this approach has been shown to not be suitable for surfactants with relatively low CMC values (less than 60 mM) [29]. A theoretical model was derived from the free energy density of a system containing free monomers, micelles, and counterions to calculate CMC; however, this approach was also shown to underestimate CMC values of various surfactants [30]. CMC has also been predicted from its relationship with Gibbs free energy and temperature [20]. Using this method, CMC values for nonionic surfactants were closer to the measured values, but for zwitterionic surfactants, CMCs were again significantly underestimated [20].
In this work, we used a more direct approach based on [20], using CG-MD simulations and a free fraction of molecules of 50% as the CMC value. This method resulted in CMC values consistent with previously reported results such as the observed decrease in CMC values with increasing carbon chain length has been observed in a number of experimental measurements [6,37,38]. The decrease in CMC values with increasing ionic strength has also been experimentally determined for various surfactants, including C 12 [39], SDS [39], DDAO [40], and including bile salts [41]. We observed that the protonation state of the MCFA concentrations affects CMC values. At the same ionic strength, when the system changed from 100% molecules with deprotonated carboxylic head groups to a 1:1 ratio of deprotonated to protonated carboxylic head groups, CMC values decreased 56% for C 8 , 44% for C 10 , 14% for C 12 , and 11% for C 14 . The systems with all molecules with neutral and deprotonated head groups mimic low and high pH, while the systems with 1:1 charged and non-charged head groups mimics the pH equal to the pKa of the molecule [33]. The effect of pH has previously been experimentally investigated for SDS [42] and decanoic acid/decanoate [18], where it was found that CMC  [43]. A number of studies have investigated the transition between vesicles and micelles at various system pH values and ratios of deprotonated and neutral molecules, and have found that MCFAs form bilayers at low and moderate pH values and micelles at high pH values [44,45]. We noticed similar pHdependent morphologies of the aggregates for C 8 (bilayers at low/intermediate pH, micelles at high pH; snapshots in Supplementary Fig. S3).
In addition to the simulations described above, we also measured CMC of C 8 , C 10 , and C 12 at two different system conditions (Table I, section 3) using the Wilhelmy method. The measured CMC values also showed decrease with increasing chain length of the MCFA molecules and decreasing pH of the system. Our simulations underestimated the experimentally measured CMCs by a factor of 1.8 to 3.5. This discrepancy was probably due to limitations of the CG force field used in this study, which is not capable of truly representing a specific pH condition of the buffer solution (water in this case). Note that the CMC is highly dependent on the buffer solution of the system. Different experimental methods used to measure the CMC also provided large   5. Snapshots of the simulations using a 7.5, b 10, c 15, and d 25 mM C 10 in the system with a 1:1 ratio of deprotonated and non-charged molecules and 140 mM NaCl. The red and green color represent the headgroup beads as negatively charged and neutral, respectively. The cyan color represents the carbon chain beads variation in results mainly due to the variation of the measurement methods and the use of different buffer solutions [6,14,16,45].
The key finding of this study was to demonstrate the viability of the CG-MD simulation technique in determining aggregation and CMCs of MCFAs at different system conditions. The close correlation between the simulation results and experimental data suggests that the CG-MD technique may be a suitable candidate for studying fatty acids as formulation candidates and permeability enhancers for the delivery of oral macromolecular drugs. This technique captured the small size difference among the various MCFAs by the very simple change of introducing a reduced CG beadbead distance. This modification of the CG Martini model seems to provide results consistent with the experimental values, but this approach requires careful validation for other systems and research questions. The major advantage of the CG-MD technique is that it requires fewer computational resources than the AA-MD simulations. Therefore, relatively large systems can be simulated over longer time scales. Our simulations needed around a microsecond to reach equilibrium so that the number of molecules in the free and aggregated forms were stable. We also needed a larger than usual box size (44 nm 3 ) for the simulations of C 12 and C 14 in order to maintain an adequate number of MCFA molecules so that aggregation could occur. We performed a few simulations by varying the box sizes but with the same ratio of MCFA molecules and box size to determine the fraction of free molecules in each case. We were able to verify that five times more molecules than the cutoff value of the aggregate size is sufficient to obtain reliable fraction of free molecules.

CONCLUSIONS
In this study, we investigated the aggregation behavior of four MCFAs commonly used as transient permeability enhancers for oral delivery of macromolecules. Our CMC measurements from CG-MD simulations were 1.8-to 3.5-fold lower than the experimentally measured values. The aggregate size, numbers, and morphologies as a function of carbon chain length, MCFA concentration, and ratio of deprotonated and non-charged MCFA molecules was consistent with experimental data. The calculated CMC values also ranked the compounds similarly as the experimental values. However, in order to reproduce the experimental CMC values more closely, improved force field parameterization to better represent the exact pH value of the buffer system might be needed. Among the current available models, polarizable water model can be used to better represent the electrostatics of the system. Finally, this study suggests that the CG-MD Martini model is promising for evaluating MCFAs as formulation candidates and transient permeability enhancers.