From microemulsion phase diagrams to hydrophilicity and hydration controlled adsorption: a dissipative particle dynamics modelling study of phospholipid assembly in bio oils

We explore here the assembly and adsorption response of a ternary bio oil–phospholipid–water system via dissipative particle dynamics (DPD) simulations. The mesoscale, particle-based modelling approach allows the examination of large-scale self-assembly response of dipalmitoylphosphatidylcholine (DPPC) phospholipids in a model bio oil solvent (modelled by triglycerides) in the presence of varying amounts of water. We report the reverse micellar and microemulsion assembly phase diagrams of the ternary mixture, verifying the model against literature data. The results show water content vs phospholipid concentration dependent transitions in reverse micellar to network-like and various lamellar phases in bulk assembly. Examination of the DPPC adsorption to smooth, homogeneous adsorbate surfaces of differing polarity reveals that phospholipid adsorption response transitions between discrete assemblies on polyethylene-like hydrophobic to continuous coating on mica-like hydrophilic substrates as the function of phospholipid and water concentrations. The significance of the work is that the presented model for phospholipid assembly in apolar solvents is capable of predicting accurately large scale assembly response and morphology changes including adsorption response in terms of system variables. The presented parametrization and verification information of the model enable readily extending the approach to other systems. The work provides computational access for tuning lipid-based microemulsion systems and their adsorption.


Introduction
Biobased microemulsions exist in key roles in biological systems but also have high technological relevance in pharmaceutic materials, cosmetics, food additives, and drug delivery systems.In biology, for example lipid droplets that have a central role in cellular lipid storage and homeostasis, correspond to the emulsion phase. 1,2The technological benefits of biobased emulsions rise from biocompatibility, aqueous green chemistry processing, and high tunability of the materials response by additives and conditions. 3,46][7] In such phospholipid based microemulsions, the hydrophobic tails of the phospholipid molecules reside in the low polarity solvent phase and the hydrophilic head groups face the water phase. 6][14][15] The high water sensitivity means that the water/surfactant molar ratio w is a key structural parameter to control the characteristics of microemulsions and reverse micellar systems. 16Presence of water typically promotes surfactant aggregation, causes swelling of the micelles, and, at high enough w, results in phase separation. 179][20][21] In biobased systems, trace amounts of water readily persist due to the hygroscopic nature of many bio-based surfactants and solvents. 22Water binds to the polar surfactant head group and, at a threshold water content, reverse micelles with a water pool at the micelle core form.Further water addition leads to anisotropic swelling and elongation of the reverse micelles, and eventually a transformation into cylindrical or worm-like assemblies. 23,24t elevated w, typically large lamellar precipitates or hexagonal phases emerge. 25,26The sensitivity of the surfactant assembly to water content, and the eventually induced phase separation, have technological uses in, e.g.bio oil purification processes, but also tuning rheological and chemical properties of microemulsions. 27he sensitivity of phospholipid reverse micellar assemblies to apolar solvent species is also well demonstrated with the assembly response extensively studied in various organic solvents, 8 including cyclohexane, [28][29][30] benzene 31 and toluene. 324][35] The existing works on the effect of solvent leave a gap over assembly response in bio oils, or triglycerides, as the nonpolar media.Works on phospholipid reverse micelles or emulsification in such solvents remain sparse, 17,27,36,37 and focus on specific conditions.Interesting for the assembly response, bio oils can solubilize large amounts of water compared to most common organic solvents, and this is reflected as curious variations in their ternary system phase response. 13,38To address this gap, we focus here on the assembly response of phospholipids in bio oils at varying water contents, i.e. the microemulsion assembly response dependencies on the system composition.
These dependencies, assembly phase diagrams, are important in designing systems with desirable colligative, rheological, and spectroscopic properties, but also in engineering microenvironments, e.g., for chemical reactions control and catalysis. 14However, detailed characterization of microemulsions phase behaviour and their microscale structural features is challenging. 39,40Particularly challenging to capture are the smooth transitions between different types of microemulsions, characteristic to reverse micellar systems, see e.g.ref. 3 and 6.Often several different experimental techniques, including NMR, 41 conductometry, electron microscopy, 42 or scattering experiments, 41,43,44 are used in tandem.Also labelling techniques provide information on assembly. 12][46] Besides bulk assembly response, also the adsorption response of ternary bio oil -phospholipid -water mixtures is broadly interesting.Phospholipids at air-water and oil-water interfaces have been widely studied, see ref. 6 for a comprehensive review.The adsorption to solid substrates is relevant, for example, in bio oils processing, where phospholipids and other impurities can be removed from the oil through adsorption onto substrates such as silica, [47][48][49] or adsorbent clays. 50,51Surface adsorption is also key to interfacial lubrication properties, 52 but also to charge transfer at electrified emulsion interfaces. 53,54][57] Molecular modelling methods offer a versatile tool for accessing surfactant aggregation and dynamics on the molecular scale in welldefined molecular conditions.Particularly, atomistic molecular dynamics approaches have been used to examine water dynamics and phospholipid hydration in reverse micelles 11,58,59 and the effect of cosurfactants, such as fatty acids, on reverse micelle morphology and dynamics. 37Due to time and length scales involved, most molecular modelling studies of reverse micelles and microemulsions focus on pre-built aggregates 11,37,58 or planar oil-surfactantwater interfaces. 60,613][64] Additional challenges to molecular level modelling of surfactants in apolar solvents are posed by the low dielectric solvent medium and typically high viscosity of the oils as these lead to heterogeneous screening of electrostatics and slow structural relaxation (in comparison to water solutions), as well as, the poor transferability of parameters from aqueous systems to nonpolar oils. 10,63,65,66At equilibrium state level, thermodynamic modelling and theory approaches for the aggregation and adsorption response can be formulated, see e.g.ref. 67, but such approaches loose molecular level structural and assembly information.9][70] While the approach allows modelling assembly and dynamics at time scales in the microsecond to millisecond range and describing system sizes in the micrometer range which brings access to large scale assembly morphologies and their transitions [71][72][73][74][75] and allows also mapping system dependencies on adsorption phases, [76][77][78][79] the trade-off is the abstraction of chemical features and molecular scale interactions, e.g.hydrogen bonding and solvent shell formation.However, for surfactant systems, well argumented DPD parametrization schemes result in micellization response with good match to experimental observables [80][81][82][83] and established scaling laws. 84ncouraged by the success of the mesoscale modelling approach to colloidal materials modelling and extracting predictions for large scale assembly in them, we develop here a DPD model for a ternary bio oil-phospholipid-water system to examine the assembly response of a phospholipid (dipalmitoylphosphatidylcholine, DPPC) in bio oil (triglyceride solvent) in the presence of varying amounts of water.The modelled system parameters correspond to a range in which the system is expected to form (w/o) microemulsions and reverse micellar systems.The model is verified against expected literature response of the assembly, obtaining a good match.Finally, adsorption response based on the developed model on smooth, homogeneous surfaces of differing polarity is mapped.The examined polarity changes correspond to changes in adsorbent hydrophilicity, or alternatively, chemical modification of the surface.Importantly, the results highlight a difference in phospholipid adsorption regime with the polarity of the surface This journal is © The Royal Society of Chemistry 2023 as well as the role of water in determining the packing of surfactant on the adsorbent and modifying aggregation behaviour.

Computational methods
DPD is a particle-based, mesocale coarse-grained (CG) simulations method that relies on describing the molecular system by beads that capture the effective characteristics of regions of the system. 85The time evolution of the beads are obtained via integration of Newton's equations of motion.Additionally, the setup conserves momentum and preserves hydrodynamics.DPD is computationally efficient because of dissipative forces, a short cut-off distance, and pair-wise additive calculation between the DPD beads that interact via a soft-core potential.
The DPD beads i interact via a total force 68,86 In this, F C ij is the conservative force, F D ij the dissipative force, and F R ij the random force.The forces are pair-wise additive, resulting in the summation j a i running over all neighbouring beads j within a cutoff r c .At distances r ij Z r c , each force component is zero.In this, r ij = |r ij |, and r ij = r i À r j .
The DPD conservative force where a ij is the conservative force repulsion coefficient between beads i an j.The unit vector r ˆij = r ij /r ij determines the force direction.The random F R ij and dissipative F D ij forces are and respectively.Here, v ij is the relative velocity between beads i and j.The magnitude of the force components is determined by the friction coefficient g and noise amplitude s, together with the distance dependent weight functions for the dissipative force and random force o D (r ij ) and o R (r ij ), as well as the random number x ij , with zero average and unit variance (normalized Gaussian white noise).The coupling of the dissipative and random forces of the DPD approach via and where T is the absolute temperature and k B the Boltzmann constant, forms the DPD thermostat via a pair-wise Brownian dashpot.It also captures the effect of Brownian motion into the larger length scale.The coupling is also responsible for momentum conservation in the system.
To construct a DPD model for the phospholipid/oil/water ternary system, we focus on tuning the DPD conservative force F C ij repulsion parameters a ij , which relate to the chemical nature of the modelled system.We consider DPPC as the model phospholipid to be parametrized, the oil solvent to be bio oil, modelled by triglycerides, and water as the polar additive.For the DPD parametrization and bead coarse-graining, these components of the phospholipid/oil/water ternary system are considered to be composed of molecule regions corresponding chemically to four different DPD bead types, CHOL, GLYS, TAIL, and Water, see Fig. 1.The chosen degree of coarsegraining N CG = 8 is to allow construction of a model that can predict equilibrium assembly responses at length and time scales corresponding to large scale assembly morphology changes.For the adsorbent surface modelling, a DPD bead SURF corresponding to varying hydrophilicity to model chemically specific adsorbent surfaces and another DPD bead species BARRIER are introduced.The bead BARRIER is for the inner parts of the adsorbate and has a high repulsion to all other bead components in the system (a BARRIER,j = 160, for all bead species j other than the BARRIER beads).This provides a penetration barrier, needed to prevent species diffusion into and through the adsorbent, otherwise allowed by the soft core repulsion of DPD particles.
The self-repulsion parameter a ii for GLYS, TAIL, Water and the SURF beads is set to 25k B T, derived from water isothermal compressibility at DPD particle density r = 3, and a common choice also in DPD modelling of other solvent species. 80,87,88or CHOL beads, a ii = 30 is used to better reflect the charged nature and increased self-attraction resulting from the zwitterionic DPPC head group.Notably, no electrostatic interactions are explicitly present in the model and the modification to CHOL self-repulsion is purely an effective modification.DPD electrostatics approaches based on a smeared charge do exist, 84,[89][90][91] and would be an option here.However, the molecular self-assembly response is dependent on the smearing. 84dditionally, charge smearing implies a uniform dielectric screening environment and here the solvent oil phase, the (hydrated) reverse micelle cores, and the adsorbent oil interface provide heterogenous, directionally varying screening.In considering the results of the current model, the lack of direct charge-charge interactions should be kept in mind, as e.g., electric double layer and electric dipole interactions do not rise from the model setup.
The DPD conservative force repulsion parameters a ij between different bead types i and j connect with the Flory-Huggins mixing parameter w ij as In this, a = 0.101 AE 0.001 for DPD particle density r Z 3. The parameter a is a constant related to the radial distribution function between the beads in a single component system.The Flory-Huggins mixing parameter was obtained from Hildebrand solubility parameters d i and d j 92 for the components i and j as Here, V bead is the DPD bead volume.The Hildebrand solubility parameter d for each of the components corresponding to the DPD beads is defined based on cohesive energy E coh determined here via all-atom molecular dynamics (MD) simulations via the energy difference between the condensed phase and the gas phase of each bead component.Details of the parametrization MD simulations including the non-bonded energies and the resulting solubility parameters (Table S1) are provided in the ESI.† The Hildebrand solution theory functions well for non-polar molecules, but fails to accurately describe molecules with high polarity, hydrogen bonding capability, or charged moieties.Other solubility theories, e.g.Hansen solubility theory, address this by identifying dispersion, dipole-dipole, and hydrogen bonding contributions in the cohesion energy.However, dividing the cohesion energy into such separate contributions lacks a physical basis.The parameters rise from an empirical estimate based on experimental data, group contribution methods, or from molecular modelling where explicit energy terms are typically associated for dispersion, dipole-dipole, and H-bonding interactions. 80,93,94In this work, the DPD interaction parameters derived from Hildebrand solubility theory were fine-tuned to match structural data of atomistic MD simulations of binary mixtures corresponding to molecules modelled by the different bead types.The protocol including the radial distribution function matching data are provided in ESI.† The obtained parametrization, i.e. the resulting set of DPD conservative force repulsion coefficients a ij is presented in Table 1.Notably, also Hansen derived DPD interaction parameters require scaling to accurately replicate interfacial dynamics in real chemical systems. 80ditionally, to constrain the DPD beads into molecule-like compounds, see Fig. 1, harmonic potentials to bind two consecutive beads together by a spring and constraint angles between three consecutive beads are included in the DPD models corresponding to Fig. 1.These contribute to additional spring force contributions to the force terms of eqn (1) and Here, K b and K a are the bond and angle spring coefficients, and r 0 and y 0 the corresponding equilibrium bond length and angle value.i, j, and k refer to consecutive beads corresponding to the bond or the angle.Table 2 presents a summary of the angle and bead-bead interactions in the parametrization.
The DPD simulations employ a DPD particle density r = 3.0.The dissipative coefficient g = 4.5.Following standard DPD formulation, all interacting beads, regardless of bead type or model chemical species, have the same mass m i and diameter r c .Reduced units, such that r c = m i = k B T = 1, are employed throughout the formulation and simulations.For the presented parametrization, the DPD cut-off distance r c = 0.891 nm when converted to real units based on the volume of 8 water molecules.

Simulated systems
Two types of simulation set-ups were examined: phospholipid/ triglyceride/water bulk mixtures (no adsorbent surface), and bulk mixtures on a planar adsorbent surface.For bulk simulations, the DPD equivalent of 50 mM, 80 mM, 100 mM, 200 mM, or 500 mM DPPC was solvated in triglyceride solvent in a 50r c Â Table 1 DPD conservative force repulsion parameters a ij for the different bead types i and j in the phospholipid/oil/water ternary system.The barrier beads in substrate have a BARRIER,j = 160 and the adsorbent surface bead interactions a SURF,j are detailed by eqn (8) for all bead species j other than the BARRIER or the SURF beads For the adsorption study, to prevent unwanted solvent species diffusion into the surface as result of the soft interactions, the centermost slab of the relaxed surface slab was designated as highly repulsive DPD particles BARRIER that have a cross-repulsion coefficient of a ij = 160.0 with all other DPD beads in the system.Sandwiching this, the outermost 1.2r c shell of the surface comprised of SURF beads that vary in their cross-repulsion coefficients a SURF,j according to the hydrophilicity variation captured by a hydrophilicity index l.SURF beads have repulsive interactions Here a Water, j and a TAIL, j are the cross-repulsion coefficients between bead type j and the Water or the TAIL beads (hydrocarbon), respectively.Hence, l provides a scaling for the hydrophilicity of the SURF beads between that of the DPD water bead and the TAIL bead.Values of 0.2, 0.4, 0.6, and 0.8 for l were used.6][97] After relaxation of the substrate slab, in the actual simulation containing also the bulk mixture, all surface beads SURF and BARRIER were frozen to ensure an intact, non-morphing surface throughout the simulation.The structure of the adsorbent surface is visualized in Fig. 2.
All described DPD systems are simulated using a time step of Dt = 0.01 for 5 Â 10 6 time steps.Simulation trajectories were analyzed for the last 2 Â 10 5 time steps.Simulations were performed by the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) molecular dynamics program. 98,99Results and discussion

Phase behaviour of DPPC/water/triglyceride ternary systems
The phase diagram of the DPPC/triglyceride/water ternary system as a function of DPPC and water concentration, as predicted by the DPD modelling, is provided in Fig. 3.In the absence of water, the model predicts a CMC between 50 mM and 100 mM DPPC concentration, with concentrations exceeding the CMC resulting in small aggregates limited in size by the lack of water core.Notably, at 50 mM concentration, DPPC aggregation remains negligible and the surfactants form a nearly isotropic solution in the triglyceride.Literature reports the CMC of DPPC in soy-bean oil being B0.085 mM. 100 The difference can be explained by the presence of trace water, but also impurities acting as nucleation centers in the experiments, as well as potential underestimation of the localized DPPC head grouphead group interactions in the DPD mesoscale coarse-grained model.The introduction of water into the system even in the coarse-grained DPD model results in stronger assembly and more diverse phase behaviour.At low water concentrations, the model predicts DPPC to form elongated, ellipsoidal aggregates around a water core.Increase in hydration or DPPC concentration both promote the aggregates to elongate into worm-like reverse micelles that may form a network-like, connected structure.Finally, at sufficiently high water content, here w = 16, large lamellar aggregates and bicontinuous lamellar assemblies form.
The simulations indicate that DPPC forms aggregates in triglyceride even in the absence of water.The size distribution of the aggregates at different DPPC concentrations is presented in Fig. 4. The exponential shape of distribution indicates that aggregation follows step-wise, open association, meaning that the aggregates do not have a clearly preferential size.Such response is characteristic to step-wise surfactant aggregation in apolar solvents.Opposed to this would be the closed association leading to bias in preferred aggregate sizes and corresponding to the micellization peak characteristic to surfactant micellization aggregate size distributions in aqueous solutions.Notably, bias in aggregate size, indicating preference in aggregate size, may emerge even for surfactants in apolar solvents, assuming that sufficient attraction between the head groups exists: the bias has been observed even for non-ionic surfactants with sufficient hydrogen bonding capability at high concentrations. 12,63However, the effect is most relevant for strongly aggregating species, such as ionic surfactants, or in the presence of a gelating agent such as water, or other small polar molecule.
In prior works, similar observations have been made for reverse micelles of other charged surfactants.For example, for AOT in isooctane, Urano et al. 101 observed a rapid decrease in the free energy of surfactant insertion into a reverse micellar aggregate with increasing aggregation number for small aggregates.The free energy of surfactant insertion was reported to have a minimum near aggregate size of 20 AOT molecules and increase for aggregates larger than 30 AOTs. 101These findings support the existence of a biased step-wise regime (and consequently a preferential aggregate size, with likely also an aggregate geometry corresponding to the preferred size) for the strongly aggregating ionic surfactant, such as AOT, in an apolar solvent.Consistent with this, the data of Fig. 4 shows the size distribution deviates from the exponential distribution.Notably, at c(DPPC) Z 100 mM, instead of the exponential decay indicated by a straight line on the logarithmic scale, the distributions show a clear bend and favoring of the larger aggregates.
In the simulations, at c(DPPC) = 50 mM DPPC in the absence of water forms small, tightly packed, spherical aggregates.
The structure of the formed aggregates was analysed by considering the radial distribution function (RDF) calculated between the CHOL beads of DPPC in the triglyceride solvent at different DPPC concentrations c(DPPC), see Fig. 5. Coinciding with the concentration induced bend in the distribution, associated with the bias in preferred size, the RDFs of Fig. 5 show the emergence of a second structural length scale.This is a telltale sign of the aggregates becoming elongated, i.e. instead of the RDF showing just the length scale corresponding to spherical aggregates, also the elongated dimension of the aggregate contributes to the RDF signal.The bend in the CHOL-CHOL RDF at r E 1.5r c becomes more prominent with concentration.Additionally, the RDF data shows that the first  This journal is © The Royal Society of Chemistry 2023 peak in Fig. 5 shifts to larger r with increased phospholipid concentration, suggesting less dense packing in the reverse micelle core with increasing concentration.This connects directly with the increase in mean size indicated by the size distributions in Fig. 4.
Based on literature, an increase in phospholipid concentration in apolar solvent such as the triglyceride here should lead to formation of an isotropic, viscoelastic, dense micellar phase. 38,102While this phase does not exhibit dominant lamellar, hexagonal, or cubic phase ordering, the presence of closedpacked structures for 60 wt% phospholipid in 50 : 50 hexane : oil mixture has been observed. 38n our simulations, at 500 mM (B57 wt%) the DPPC aggregates start exhibiting some ordered packing, as indicated by the subtle long range ordering peak at r E 5.8r c in the CHOL-CHOL RDF curve following the dip corresponding to the individual aggregates, see the inset of Fig. 5.The long-range ordering peak indicates not only correlations in the small aggregate positions but also depletion of solvent triglyceride as replaced by the aggregates.Notably, such long range order is absent for c(DPPC) o 500 mM, indicating the aggregates remain freely solvated in the triglyceride solvent.The location of the observed long range ordering peak r E 5.8r c E 5.1 nm is in excellent agreement with d-spacing of 5.45 nm from SAXS measurements of 60 wt% phospholipid in 50 : 50 hexane : oil mixture for the dense reverse micelle phase. 38Opposed, at low phospholipid concentrations, such long range order is not present. 38,102,103This means that the DPD model of this work captures the transition from a dilute micellar phase to a dense micellar phase with increased DPPC concentration.Notably, in the presented modelling, the extent of long-range ordering is limited by the heterogeneous size of the aggregates and transition from spherical to ellipsoidal micelles, as well as, the finite simulations box size.Despite the simulation box being extensive 50r c Â 50r c Â 50r c , the periodicity of the simulation box sets a constraint for the long-range assembly.This affects both capturing long range order and the periodic structures formed.
The transitions in reverse micellar shape and aggregation response in the presence of water in the system can be considered via changes in the surfactant packing parameter, that is, effective shape of the molecule determined by the head group size vs. size of the hydrophobic tails.Hydration increases the proportional size of the head group, which in this case promotes first elongation of the micelles and eventually leads to lamellar packing structures.
In oils, phospholipids can be expected to have a strongly bound water shell of 6-10 water molecules, hydrating the phospholipid head group and greatly limiting the growth of reverse micellar aggregates at small water concentrations. 38owever, even trace amounts of water may trigger the elongation of phospholipid reverse micelles. 8,28,104Introduction of water in the phospholipid-oil system first leads to micelles transitioning from spherical to rod-like, and further addition of water promotes the growth of cylindrical micelles.Cylindrical micelles transition to worm-like micelles once their size exceeds their persistence length.Also branching may occur.The elongation and entanglement of the aggregates can cause an abrupt change in solution viscosity, potentially by several orders of magnitude. 8,105The amount of water needed for the formation of worm-like micelles depends strongly on the solvent; for phosphatidylcholine lipids in hydrocarbon solvents, transition of the solution into a gel occurs at w E 7. 8 Specifically for DPPC in triglyceride, a slightly higher transition water content is expected as the hydrophilic moiety of the triglyceride molecules also binds water, meaning that some of the water partitions into the solvent as well.
In the DPD simulations, at water-to-DPPC ratio w Z 16 the assembly morphologies correspond to large disk-like aggregates and lamellar phases.The formation of lamellar phases is also linked to changes in water dynamics, particularly greater mobility of the water DPD beads.Fig. 6 presents the scaled diffusion coefficients of the DPD water beads at different waterto-DPPC ratios.Notably, water diffusivity at DPPC concentrations r100 mM remains systematically low and largely independent of the degree of hydration (water-to-DPPC ratio w).This is expected as the assembly phases correspond to formation of either small rod-like or elongated and worm-like reverse micelles, in which water remains mostly bound tightly around the DPPC polar head groups and confined to the aggregate cores.The confinement limits the diffusion based movement.Conversely, for DPPC concentrations Z200 mM, water mobility increases clearly as a function of hydration, with the water-to-DPPC ratio w controlling the diffusion.This matches expectations for networklike reverse micellar aggregate phases and lamellar structures as water has pathways for long-range propagation in such assembly morphologies.At elevated concentrations, in such assembly morphologies, water can be expected to have even free movement at the scale of the system.Previously, for similar microemulsion forming ternary systems, an increase in water mobility was linked to formation of bicontinuous inverse lamellar phases. 106ater molecules in microemulsion droplets or reverse micelles can be classified into either two or three different types based on how tightly the water is bound.1][112] Similar bound and unbound water populations have been observed in e.g., polyelectrolyte multilayers 112 but also in water adsorbed in doped carbon nanopores. 113The RDFs of water beads around the DPPC heads in Fig. 6 indicate that the increased water mobility results mainly from the amount of bulklike water increasing, opposed to the bound water around the DPPC head group remaining at constant level.Analogous distinction between bound and bulk water has been observed in reverse micelles with sizeable water pools 114 and also mapped in atomistic simulations. 11Molecular level interactions of water and the lipid head groups determine largely the binding strength and size of the head group water shell.For the DPD model here, and coarse-grained models in general, modelling the hydration number is mainly by the bead size and its polarity (conservative interaction parameter).Notably, the microemulsion droplet-like reverse micelles that have relatively large water cores at c(DPPC) = 50 mM and w = 32 have overall low diffusivity of the water.This is a direct consequence of the model being a coarse-grained mesoscale model and the diffusion coefficient capturing mobility only inside the water pool, now limited in size by the coarse-graining.The fine grained structure and detailed dynamics of water is lost.For comparison with atomistic detail simulations results, see ref. 11.
We conclude that the DPD model reproduces at qualitative level the DPPC assembly response in triglyceride in the absence of water, but also when water is added.Notably, the presence of CMC and small aggregate formation is predicted accurately, as is the emergence of dense micellar phase at higher DPPD concentrations.Furthermore, the transitions to elongated and wormlike aggregates, as expected based on literature, emerge upon addition of water, and further water addition transitions the assembly phases to various lamellar phases, as expected.

Adsorption of DPPC on adsorbents of differing polarity
As the DPD model reproduces bulk assembly response to good degree, we move to employing the model to map much less characterized phospholipid assembly response, that of surface adsorption from ternary DPPC/triglyceride/water system.First the adsorption of DPPC and water on planar surfaces of varying polarity was examined.
Adsorption acts as a competing process for DPPC and water aggregation.Importantly, the results reveal a difference in adsorption regime on hydrophobic and hydrophilic surfaces, see Fig. 7 and 8. Fig. 7 presents the density of DPPC head groups on the adsorbent surface as a function of surface hydrophilicity parameter l and the water-to-DPPC ratio w in the oil, modelling the oil-surfactant system hydration.Fig. 8 presents the corresponding adsorbed phospholipid and water quantifications for different DPPC concentrations.In the absence of water at w = 0, DPPC adsorption at the surface increases with DPPC concentration and surface hydrophilicity.On hydrophilic surfaces (l = 0.8, comparable to a mica-like surface) a DPPC monolayer forms with coating saturation dependent on DPPC concentration, see Fig.  3 nm 2 per DPPC lipid have been reported for adsorption from oil onto hydrophilic surfaces, including sepiolite, silica, silica-alumina and bleaching earths. 50,57The tighter packing in the simulations can be explained by the idealized setup, for example the flat, ideal surface.Notably, the resulting packing density predicted by the model at the hydrophilic substrate is much closer to what would be expected of phosphatidylcholine monolayers (0.6 nm 2 ). 61Notably, the packing density is likely to be overestimated here due to insufficient repulsion between phosphatidylcholine head group DPD beads.
In the ESI, † Fig. S1 summarizes the densities of the system constituents near the adsorbent.On hydrophilic surfaces (l = 0.8, comparable to a mica-like surface), both DPPC head groups and water adsorb directly onto the surface.Here, water acts as a competing adsorbate, as shown by the decrease in DPPC adsorption in systems containing water.On hydrophobic surfaces (l = 0.2, comparable to a polyethylene surface) both the water and DPPC head group density peaks shift away from the surface, while the DPPC tail peak shifts closer to the surface.This indicates either aggregate or lamellar like structures on the surface.On hydrophobic surfaces also the triglyceride solvent exhibits structuring at the interface, likely due to the Interestingly, the density mappings of Fig. 7 reveal the presence of clustered DPPC aggregates on hydrophobic surfaces (l = 0.2, comparable to a polyethylene surface) at 100 mM DPPC concentration.These aggregates interact only weakly with the hydrophobic surface.The assembly at surface corresponds to the hydrophobic hydrocarbon tails facing both adsorbent surface and the solvent phase.Such reverse micellar aggregates on hydrophobic surfaces support previous experimental findings of minimal phospholipid adsorption on smooth hydrophobic polymer surfaces and increased adsorption on porous surfaces of similar chemistry. 56The findings could be explained by surface reverse micellar aggregates being trapped in the pore geometry of the adsorbate.Even step edges align and bind surfactants 115 and any pore edges may have chemically different character, promoting aggregate pinning: on patterned and stepped surfaces, the length scale of the surface defects and the surfactant molecular size, as well, as chemical heterogeneity of the adsorbent surface leads to deviation from the planar adsorption response, see e.g., ref. 116-118.Notably, adsorbent surfaces decorated by adsorbed micelles have been previously reported for adsorption of both non-ionic and ionic surfactant onto silica in water and such assemblies may be described in literature as admicelles, bilayers, patchy bilayers, or interdigitated bilayers. 119,120In aqueous systems, an S-shaped adsorption isotherm response with initial hemimicelle adsorption at dilute surfactant concentration, followed by admicelle adsorption above CMC, when bulk aggregation occurs, can be expected. 119n the oil system, a more simple adsorption response is expected due to the lack of CMC.
The data of Fig. 9 shows that with increasing adsorbent hydrophilicity, surface coverage shifts from island like aggregate coating to the formation of an evenly distributed monolayer.This difference in adsorption regimes is indicated by the RDF curves for DPPC heads and tail with respect to the surface particles: the DPPC heads gradually shift away from direct contact with the adsorbent surface with increased adsorbent hydrophilicity parameter l.An opposite trend is observed for DPPC tails, with direct adsorption of the tail groups with adsorbent on hydrophobic surfaces (l = 0.2).
The adsorption regimes observed in the absence of water persist also in the presence of water, see Fig. morphs the shape of the spherical surface aggregates into elongated ellipsoids.The ellipsoids undergo further growth with increasing water concentration.This results in increased DPPC adsorption on hydrophobic surfaces with increased water content, see Fig. 10.The aggregates present on the surface largely reflect those present in bulk solution.For example, at elevated water concentration, water-to-DPPC ratio w = 16, large inverse lamellar structures form on the surface, consistent with bulk phase behaviour.Here, due to adsorption being driven by diffusion, the time scale of the modelling may not accurately capture equilibrium state adsorption of larger aggregates, particularly slowly diffusing lamellar precipitates.
The competitive adsorbate character of water to DPPD on hydrophilic surfaces results in depletion of DPPC from the surface with increasing system water content, see Fig. 8.As most of the adsorbed water is bound to a phospholipid head group, the depletion of DPPC rises from the adsorbed hydrated DPPC head groups occupying a larger surface area than head groups with less water.The increase in head group radius can be expected to scale with w 1/3 as it is a volumetric quantity.Experimentally, the effect of water on the adsorption of both hydratable and non-hydratable phospholipids has been shown to be minimal. 37,57For example, Lehtinen et al. showed negligible effect of up to 1.0 wt% water on DPPC adsorption from rapeseed oil onto acid activated sepiolite.However, the denser monolayer packing predicted by the model here is more sensitive to changes in the DPPC head group size with hydration.We hypothesise that the sparser packing of DPPC in experiments 37 results from electrostatic repulsion between the lipid heads, which is largely unaffected by the formation of the first, strongly bound hydration layer around the head group.The employed DPD model considers the effect here as coarsegrained effective size change.

Conclusions
Here, we presented a DPD model for a DPPC/triglyceride/water ternary system and examined its bulk phase behaviour and adsorption properties.The bulk assembly response predicted by the model agreed very well with experimental assembly morphologies in large scale.We demonstrated that the system undergoes initial aggregation, with aggregate growth greatly limited in water-less systems.However, even for dry binary mixtures of DPPC and triglyceride solvent, onset of cooperative   aggregation with increased phospholipid concentration was observed.Additionally, the model captures structural transition from a dilute micellar phase to a dense micellar phase.This coincides with elongation of the aggregate shape profile.When water is present, aggregates undergo growth and elongation from ellipsoidal to worm-like.Continued aggregate growth results in the formation of network like organogel phase and lamellar phase depending on oil water content.For the relatively high degree of coarse-graining (N CG = 8) used here to enable modelling large scale assembly and morphology changes in equilibrium, the agreement with experimental verification data is very good.
For surface adsorption, we reported two different adsorption regimes depending on surface hydrophilicity: monolayer formation (hydrophilic adsorbent) and adsorption as aggregates (hydrophobic adsorbent).Water acted as a competing adsorbate on hydrophilic adsorbents but also promoted aggregate growth and adsorption of admicelles on hydrophobic surfaces.2][123] Here, we generalized this adsorption response for surfactants in oil.It is also interesting to speculate what response would other additives besides water cause.Small polar additives can be expected to behave like water, whereas e.g.hydrogen bonding ability will compete with both the solvent, surface and DPPC, see e.g.ref. 37.In prior work, we examined nitrogen compound adsorption from bio oil, resolving molecular mechanisms for related biospecies. 65e also reported that increasing hydration pushes some of the DPPC lipids from the surface.This depletion of phospholipids could largely be explained by the increased area of hydrated head groups.The adsorption of phospholipids is largely insensitive to water content in experiments mainly due to sparse packing of the lipid monolayer on the surface by electrostatic repulsion of the head groups. 50Here, the DPD model disregards electrostatic interactions, resulting in a denser packing structure of phospholipid monolayers despite the increase self-repulsion between DPPC head groups.
Overall, the work shows that DPD offers a versatile mesoscale modelling method for examining equilibrium structures and assembly phases for even complex, colloidal multicomponent systems.The obtained significant extension of the temporal and spatial resolution in comparison to, e.g. even coarsegrained classical atomistic MD approaches makes the approach powerful in terms of modelling reach.Additionally, the presented work provides guidelines for designing both microemulsion and reverse micellar bulk systems but also for predicting phospholipid adsorption response from bio oil solutions at interfaces of varying polarity.

Fig. 1
Fig. 1 Coarse-graining scheme for the DPD model of DPPC.An additional TAIL bead was added to the phospholipid hydrocarbon tails to capture better the molecular asymmetry.Model for triglyceride solvent was constructed from the hydrocarbon TAIL beads and a central GLYS bead.The water bead corresponds to B8 H 2 O molecules.

Table 2
Parameters for harmonic bonds and angles between DPD beads in the model.The bead-bead bonds have equilibrium length r 0 and spring constant K b and the equilibrium angle y 0 has a spring constant K a .Water, SURF, and BARRIER beads are not connected to other beads via the harmonic interactions in the model Bond r 0 [r c ] This journal is © The Royal Society of Chemistry 2023 50r c Â 50r c simulation box with 0, 0.5, 1, 2, or 4 water beads per DPPC molecule.Due to coarse-graining at 8 water molecules per DPD bead, this corresponds to 0-32 water molecules per DPPC molecule.All molecules were set initially randomly and uniformly distributed into the simulation box.The planar DPD surface was constructed by equilibrating a 40r c Â 40r c Â 5r c slab of identical DPD particles with density r = 3 and self-repulsion coefficient a ii = 25.0.This slab was then centered perpendicular to the z-axis into a simulation box of size 40r c Â 40r c Â 100r c , resulting in the surface being periodic along the x,y-plane, and the z-axis being perpendicular to the plane.The remainder of the box was filled with DPD equivalent of 50 mM, 100 mM, or 200 mM DPPC solvated in triglyceride, with w = 0, 0.5, 1, or 2 water beads per DPPC molecule.Example visualizations of the simulated systems are presented in Fig. 2.

Fig. 2
Fig.2Visualizations of example systems representing the simulated bulk self-assembly and surface adsorption systems.Coloring of DPD beads follows Fig.1for the beads corresponding to DPPC and water.Triglyceride solvent is omitted in visualization for clarity.The adsorbent surface consists of two types of beads with the chemically specific surface beads SURF in orange and the highly repulsive barrier beads BARRIER in white.

Fig. 3
Fig. 3 Assembly phase diagram for phospholipid (DPPC) in triglycerides solvent as a function of phospholipid concentration and water-to-DPPC ratio.The snapshots show representative visualizations from each assembly phase.The water-to-DPPC ratio w = N(water)/N(DPPC) has been converted from the DPD beads to number of molecules, i.e. one DPD water bead corresponds to 8 water molecules.

Fig. 4
Fig. 4 DPPC aggregate size distributions for different DPPC concentrations in triglyceride in the DPD simulations.The data corresponds to water-to-DPPC ratio w = 0.

Fig. 5
Fig. 5 Radial distribution function RDF calculated between the CHOL beads of DPPC in triglyceride solvent in the absence of water (w = 0).The inset presents the RDF for c DPPC = 500 mM with the long range delocalized 2nd ordering peak marked with arrow.

8
. The formed monolayer coverages correspond with lipid surface areas of 0.85r c 2 E 0.68 nm 2 for c(DPPC) = 50 mM, 0.57r c 2 E 0.45 nm 2 for c(DPPC) = 100 mM, and 0.50r c 2 E 0.40 nm 2 for c(DPPC) = 200 mM.The values compare well with prior experimental characterizations in which areas of 1.2-1.

Fig. 6
Fig. 6 Top panel presents the scaled diffusion constant vs. water-to-DPPC ratio w = N(water)/N(DPPC) for different DPPC concentrations.Bottom panel presents the radial distribution function RDF calculated for DPD water beads around the DPPC head (CHOL + GLYS beads) for 200 mM DPPC in triglyceride solvent at varying water-to-DPPC ratios w = N(water)/N(DPPC).

Fig. 8
Fig. 8 Top panel presents the surface density of adsorbed DPPC and the bottom panel the corresponding data for adsorbed water molecules on hydrophilic (l = 0.8, mica-like substrate) adsorbent surface as a function of water-to-DPPC ratio in the system.

Fig. 9
Fig. 9 Radial distribution function RDF calculated between DPPC heads and surface beads (top panel) and tails and surface beads (bottom panel) for adsorbent surfaces of differing hydrophobicity parameter l.The presented data corresponds to DPPC concentration 100 mM and w = 0. Inset cartoons point out the transition from aggregate adsorption to monolayer coverage with increasing l.

Fig. 10
Fig. 10 Surface density of adsorbed DPPC on hydrophobic (l = 0.2) adsorbent surface as a function of system water content.