The significance of lipid membrane curvature and lipid order for nano-C60 uptake

Abstract

Interactions that invasive mesoscopic nanoparticles (NPs) experience with the variety of lipid membranes present in our body have the potential to trigger diverse biological responses and endocytic pathways. From a safety perspective, it is therefore crucial to comprehend how these NPs are absorbed by complex biological environments and how they subsequently disperse and disassemble within them. Previous computational research has focussed on evaluating the compositional dependence of NP binding and translocation to ‘model’ lipid membranes, i.e. membranes that are in equilibrium, flat, fluid and composed of a few lipid species. The observation of strong local variation of curvature and phase along many membranes in biology, however, suggests that these features are important for membrane-related processes. In particular, the thermodynamic factors underlying curvature formation and phase separation appear to be essential for regulating intricate membrane functions such as insertion and dispersion in a precise and controllable manner. In this study, we have employed coarse-grained simulation to shed light on the role of curvature and phase on the adsorption and dispersion of a typical hydrophobic NPs like nano-C60 (or n-C60), for a model setup of a lung membrane at biologically relevant scales. The effect of curvature on n-C60 uptake, which is particularly striking when the membrane is in the gel phase, is explained in terms of lipid defects, and we find that uptake modulates the fluidity and defect characteristics of the resulting membrane in a state-dependent fashion, and will affect protein binding. Our discovery offers key functional insight in the role of the particular curvature and phase, i.e. a combination of a flat gel state and curved fluid state, that the lung membranes exhibits as a primary barrier.

---

1 Introduction

We inhabit a world characterised by hierarchical structures that engage in diverse interactions. Nanoparticles (NPs) stand out as a class of structures that is unique in their elevated surface-to-volume ratio as well as in their ability to spontaneously or actively internalise and disperse within larger and more complex structures such as the human body. While these unique characteristics provide avenues for intentional interventions in animals and humans, for instance via targeted drug delivery or local manipulation through external fields, NPs also have the potential to disrupt the intricate network of structures and interactions that are crucial for maintaining key biological functions. Consequently, gaining a deeper insight in the insertion of NPs, via knowledge of how NPs and biomembranes interact, is essential from the point of view of materials design.

Within this class of materials, carbon-based NPs like fullerenes and nanotubes stand out in terms of produced quantity, controllability of structural properties, and experimental characterisation effort. Fullerenes or C60 are special in that they are utterly reproducible, generated both by industrial fabrication methods and natural processes, and small enough (diameter ≈1 nm) in their monomeric form to spontaneously enter a lipid membrane. Insertion, taking place via a process that is anticipated to involve short-lived micropores,¹ makes them settle in regions of high density without significantly perturbing chain conformations of the surrounding lipid matrix. Whereas monomeric C60s interact with protein via binding pockets, and may thereby alter protein function, albeit in a fairly a-selective fashion,² small aggregates of C60 called n-C60 that form in poor solvent like water are usually stable at a size of 50–150 nm,³ and have been identified to carry a protein corona in blood serum⁴ just like other hydrophobic NPs of comparable sizes. Inside the membrane, high C60 loading is known to change mechanistic membrane properties, which can be beneficial for C60-filled carrier vesicles. Membrane destabilisation, which has been proposed as one of the mechanisms for C60 toxicity, has only been observed experimentally at very high concentrations of C60 and was never reproduced in silico.^5,6 The effective hydrophobicity of fullerenes, which in combination with aggregate size regulates their membrane binding affinity, can be tuned via surface functionalization. For pristine C60, the thermodynamic driving force for aggregation in water is substantial.² As a consequence, both in experimental and computational setups special conditions are needed to control the aggregate size prior to membrane binding and insertion.

Their reproducibility and broad scope for design of binding interactions via surface functionalisation, both in terms of hydrophobicity and specificity, their sensitivity to external factors and to their local environment, and their natural propensity to enter hydrophobic domains like membranes spontaneously and perturb structural properties in a semi-controlled fashion enables them to perform very diverse and unique roles.⁷ For instance, C60 were previously functionally applied as photosensitizer⁸ and as a scavenger of free radicals. Yet, the oxidative stress that nano-C60 is also known to induce, in combination with cell membrane damage at high concentrations, poses a severe problem in practice, particularly since C60 is known to disperse throughout the body including the brain. Precise detail of what causes toxic effects at a molecular level are absent, and thereby the means for control. Given the current limitations in experimental resolution, in silico study appears most promising for providing molecular insight into such causes of toxicity.

When investigating the interaction of hydrophobic C60 with bio-membranes, the role of collectivity is eminent in the growth of C60 aggregates (nano- or n-C60) in an aqueous environment prior to uptake and in the induced membrane remodelling during uptake itself, while these aggregates are thought to disperse into tiny clusters and even individual C60 in the hydrophobic interior of the membrane. In general, NP binding and uptake characteristics were found to vary considerable in the size range 3–50 nm, comparable or somewhat larger than the membrane thickness. The average size of stable n-C60 assemblies in water is at the higher end of this range, and the shapes that n-C60 adopt by spontaneous aggregation is generally far from spherical, introducing additional and less studied factors in NP-membrane binding. Yet, computational studies have focussed on the binding of tiny n-C60 clusters to flat membranes in a liquid state. Here, we systematically consider the role of lipid packing on uptake and dispersion of n-C60 of more realistic dimensions, in order to provide new insight in the role of size, membrane curvature and lipid domain formation in n-C60 binding to and insertion into actual membranes. Since PC lipids like DPPC or DOPC are primary components of the various types of biomembranes that are macroscopically in a liquid state at body temperature, we do adopt the standard choice for a single-component model membrane. In our case, this component is DPPC, the majority lipid (50%) in a lung membrane and already considered in many computational studies of C60 binding. The versatile role that cholesterol plays in multi-component membranes is beyond the scope of the current investigation, but we will shortly discuss its anticipated role in the Discussion section.

The choice for pure DPPC, providing a membrane in the gel state below 41.3 °C, requires some additional explanation. While coarsening seriously extends the computationally accessible time and length scales of an in silico investigation, any map from fine to coarser is formally state dependent, meaning that exact reproduction of thermodynamic properties for states other than the mapped one is not guaranteed. Earlier studies employing the popular Martini coarse-grained (CG) DPPC lipid representation showed that the gel-to-liquid disordered (L_β-to-L_α) phase transition is recovered at the CG level, but for a transition temperature T_g that shifts down as much as 20 °C, to 22 ± 5 °C, compared to the experimental value.⁹ Although such mismatches typically have little impact if handled with care, it should be understood that, as a rule of thumb, the membrane state in CG will never exactly match the AAMD or the experimental state. For DPPC, the temperature at which the CG simulation takes place does matter, as even the actual phase will differ in the temperature range 27–41 °C between experiment/atomistic simulation and CG simulation. The alternative reasoning, i.e. that a precise separation of enthalpic and entropic factors to the free energy usually cannot be established when projecting to a coarser representation, may suggest that the delicate balance between enthalpic gain and the entropic costs for bending lipids around monomeric C60 or around n-C60, which to some extent determines the stability of n-C60 inside a particular membrane, will be perturbed by switching to another, in this case coarsened, representation. One may, however, conclude that this balance is at least qualitatively retained at the CG level from a recent Martini study of heterogeneous lipid membranes, which shows the expected preference of n-C60 for domains of poly-unsaturated (and thus most flexible) lipids.¹⁰ In the same study, C60 was found to repel cholesterol, which confirms the enhanced partitioning of cholesterol into ordered domains within a mixed ordered/disordered liquid membrane upon addition of C60 observed previously, and shows that C60 stabilises raft formation.¹¹ We will exploit this information – cholesterol and C60 compete in terms of packing around lipids – to anticipate the role of cholestrerol in the Discussion section.

A factor previously identified important for the binding of small molecules or NPs is the characteristics of so-called lipid defects, expressed in the defect constant or average defect size π. Defects relate to lipid packing and are therefore inherently sensitive to temperature. This dependence of π was previously analysed for DPPC membranes using atomistic MD, and π was found to increase linearly with T for liquid membranes, whereas π experiences a sudden drop upon crossing T_g into the gel phase.¹² One may thus conclude that one way to alleviate issues with AA-to-CG mapping is to tune the simulation temperature, a thermodynamic parameter. It illustrates that setting up a computational study of small NP binding that focusses on the extraction of quantitative data or state-dependent mechanisms requires careful thinking ahead. Here, we set up systems to analyse local curvature and phase separation – relevant for many passive processes and transport phenomena, including the regulation and modulation of signalling events, trafficking, recognition, and more – separately via biased and unbiased coarse-grained (CG) molecular dynamic (MD) simulations. Tensionless flat membranes and membranes that buckle under applied steady strain are used to consider the role of (strong) curvature. Temperatures well above and below the (CG) transition temperature are employed to consider well-defined membrane phases, gel (L_β) and liquid (L_α). In addition, we analyse the role of C60 aggregate (n-C60) size by considering two C60 concentrations. Since we concentrate on passive binding processes, we will further neglect the role of membrane-bound proteins.

While the study of strong curvature on (n-)C60 uptake and stability is rare, the role of phase behaviour as well as the perturbation of mechanical properties due to uptake have received some attention both experimentally and by simulation. Analysing vesicles formed by different lipids at the same temperature, an increase of the Young's modulus in gel DPPC or DSPC membranes was found upon addition of (n-)C60 during vesicle preparation, but not in their liquid DMPC counterparts.⁵ In pulmonary membranes, the outer lipid mono-layer is in the liquid ordered state and coexists with curved fluid (liquid disordered) domains that provide elasticity to expand and compress during respiration.¹³ Lipid concentration gradients and membrane stresses along the membrane do not only underlie such phase coexistence, but they can also induce strong local curvatures, varying rigidity and permeabilities, all of which play a role in the diverse functions of membranes. The 20 μs simulation trajectories for the eight considered systems will be analysed via a number of C60 and membrane attributes. Membrane defect characteristics are extracted via a recently developed procedure that is applicable for arbitrary membrane shapes and able to distinguish between defects in flat and curved domains.¹⁴ Binding and (dis)assembly kinetics of n-C60 are monitored via the lipid-C60 contact fractions and the total number of monomeric C60s. One would like to analyse the changes of the local membrane structure that are induced by binding and insertion of (n-)C60, but extraction of such information is far from straightforward. In particular, the usual quantities like the area per lipid (APL), membrane thickness, orientational tail order parameter P₂, or lipid diffusion rates are ensemble values, and their extraction for curved membranes can still pose a challenge. We take a more direct route and apply a recently developed machine learning (MLLPA) approach that is capable of determining whether individual lipids belong to a gel or liquid phase for arbitrary membrane shapes with good accuracy. Projecting the instantaneous conformations of individual lipids that are used for this classification back to the membrane provides direct access to such local information.

2 Material and methods

All simulations were performed with the GROMACS 2021 package,¹⁵ using a Martini 2 representation for C60 and DPPC.¹⁶ We selected this slightly older version of Martini because a proper parameterisation of C60 in Martini 3 was not available at the time of this study. As mentioned, we considered lipid membranes consisting of Dipalmitoylphosphatidylcholine (DPPC) lipid only, generating initial membrane configurations using the INSANE building tool.¹⁷ In addition to the standard flat geometry, we considered buckled membranes to investigate the role of local curvature. The initial buckled membrane was generated using the BUMPy tool.¹⁸ To transform the flat coordinates to buckled membrane, we used a radius of 7.5 nm and length 20.0 nm to maintain the similar number of lipids to flat membrane. We considered system with two different C60-to-lipid ratios in the same simulation volume, which leads to different concentrations and sizes of the n-C60 aggregates that are formed in the aqueous phase prior to absorption. Most studies consider low concentrations of C60 or disperse C60 in close proximity of the membrane, rendering the role of n-C60 NP size rather unexplored. The cutoffs for electrostatic and van der Waals interactions were set at 1.4 nm with dielectric constant at ε_r = 15. Long-range electrostatics interactions were evaluated using a Particle mesh Ewald method.¹⁹

For all flat bilayer simulations, the pressure in the npT production run was controlled at 1 bar using a Parinello–Rahman barostat²⁰ via semi-isotropic pressure coupling and a compressibility of β = 3 × 10⁻⁴ bar⁻¹. Buckled membranes were simulated using semi-isotropic pressure coupling, with the long axis aligning with the z-dimension, using a reference pressure of 1 bar only in the z-dimensions to preserve the curvature during simulations. Periodic boundary conditions were employed throughout, and simulations were performed using a v-rescale thermostat²¹ at either 285 K (gel membrane) or 325 K (liquid disordered membrane) with the solvent and lipid bilayer coupled separately to avoid issues with equilibration.²² All reference system i.e. lipid membrane without C60, were energy minimized using the standard steepest descent algorithm,²³ followed by a short 10 ns equilibration run. Final production simulations were then performed for 10 μs in all cases. For each reference membrane with C60, the initial conformations was similarly minimised using the standard steepest descent algorithm.²³ This was further followed by two small equilibration intervals of 10 nanoseconds (ns), with a 1 femtosecond (fs) time step and another 10 ns using a 10 fs time step with C60 position constraint to their initial position. This was followed by a 5 ns run after releasing all the constraints on the system. Finally, production runs of 20 microseconds (μs) were performed for all systems, with a time step of 20 fs, and the last 1 μs used for analysis. Mean square displacements (MSD) as a function of time were calculated using the built-in GROMACS routines. If reported, values for strongly curved membranes are and should be evaluated using methods that account for non-flat geometries. The membrane mean curvature with and without C60 was determined by the program g_lomepro.²⁴ Curvature plots use data extracted from the last 1 μs of the trajectory, and we applied a bandpass filter with an upper limit of 2.0 nm⁻¹ to accurately capture the mean curvature for the highly curved membrane. We calculated the C60-DPPC contact fraction using,²⁵

Here, C represent number of contacts. It should be noted that only the phosphate (PO₄) bead was considered for the DPPC lipid in this analysis, and the threshold distance was set to 0.8 nm. To quantify the contacts in the flat and buckled regions, lipids were classified based on their z-coordinates at 10 ns intervals throughout the 20 μs production run. At each time point, lipids were assigned to either the flat or buckled region according to their z-position. Subsequently, the number of contacts between C60 and the lipids in each region were calculated as previously described.

2.1 Free-energy calculation

To determine the free energy associated with C60 insertion or dimerisation inside the lipid phase, we calculated potentials of mean force (PMF) using an Umbrella Sampling (US)²⁶ technique. The PMF for the insertion of a C60 into the lipid bilayer from the solvent phase was calculated using the distance between the center of mass (COM) positions of the C60 and the lipid bilayer as a reaction coordinate, in agreement with standard practice. The validity of this choice is confirmed by the observation that the perturbation of lipid conformation stemming from the insertion of a single C60 into the membrane is negligible. Note that the response to n-C60 insertion will differ, meaning that our findings for single C60 cannot simply be extrapolated. Initially, the C60 was manually inserted 2.0 nm above the equilibrated lipid membrane and subsequently solvated with CG water beads of which 5–10% were antifreeze water beads to avoid unphysical effects. We applied a force constant of K = 500 kJ mol⁻¹ nm⁻² and an equidistant spacing along the reaction coordinate of 0.125 nm, giving rise to 35 configurational windows. PMFs for C60 dimerisation inside the lipid membrane were obtained for T = 285 (gel phase) and 325 K (liquid disordered phase), keeping the two transverse C60 (z) positions fixed. A total of 35 configurations were generated using the COM distance between the two C60s as a reaction coordinate. Here, each con-figuration window is simulated for 2.5 μs for all the free energy calculations. Subsequently, in all cases, the standard weighted histogram analysis method (WHAM)^27,28 was applied to extract the free energy profile. Errors associated with the PMF determination were analyzed with the bootstrap routine that is available in GROMACS.

2.2 Analysing lipid phases using MLLPA

Characterizing small, localized perturbations in the lipid matrix using conventional membrane descriptors – such as average curvature, membrane thickness, area per lipid, or the tail order parameter P₂ – is challenging, since subtle effects are often masked by ensemble-averaged properties of the entire membrane, which itself must be sufficiently large to minimize boundary artifacts. Focussing on local properties, we benefit from a recently developed machine-learning-assisted lipid phase analysis (MLLPA) method that analyses membrane phases via individual molecular conformations.²⁹ Besides ensemble values, such as the liquid and gel percentages, MLLPA enables the automated identification of lipids that are sensitive to the (n-)C60 kinetics in simulation shapshots. We trained MLLPA using our C60-free data for the pure gel and liquid membranes. We particularly used only a minority percentage (20%) of the lipid conformations present in each flat membrane. It should be noted that lipid orientations play no role in MLLPA, since they are aligned during the procedure, and that the overall membrane curvature only plays a role in neighbor selection for the Voronoi tesselation protocol,²⁹ which is dealt with by proper instruction of the MDAnalysis routines that underlie MLLPA. The trained MLLPA was subsequently used to predict the percentages of fluid/gel in all four pure DPPC membranes. Compared to the averaged 86.8 ± 0.7% accuracy previously observed for a CG Martini representation of DPPC,²⁹ here we find a 89.1% fluid at 325 K (flat liquid membrane) and 93.8% gel at 285 K (flat gel membrane), see Fig. 2C (average 91.45%), corroborating that MLLPA provides proper predictions for the membrane states at the considered temperatures. Predictions for the pure buckled membranes are provided in the Results section. For completeness, we note that the issues of MLLPA in distinguishing between the ordered liquid phase (L_o) and the gel phase (L_β) are not relevant for this study.²⁹ The results of MLLPA analysis that is alternatively trained on the buckled pure DPPC membranes, see Fig. S4, corroborate that membrane curvature indeed does not play a role in this procedure.

2.3 Characterisation of membrane defects

Defect characteristics are extracted using a recently developed Coarse-Grained to Triangular Surface (CG2TS) algorithm implemented in a Python 3.7 code that builds upon the open source libraries MDAnalysis,³⁰ Numpy, SKLearn, Open3D,³¹ and Optimesh. Here, the normalized frequency of finding a defect p(A) of size A (in Ų) extracted from simulation is fitted to the mono-exponential decay function b·exp(A/π), rendering a defect constant or average defect size π in units of Ų. Defect constants π for a variety of membranes and examples of standard fitting diagnostics can be found in the original study, and show that the defect constant π is sensitive to several thermodynamic parameters and internal membrane properties, including membrane curvature and lipid phase.¹⁴ Since the procedure projects particular atom or CG bead volumes onto a 2D grid to extract their overlap, the calculated π values are resolution and force field dependent.¹² Using Martini 3.0, it was previously shown that a defect constant π = 17.5 obtained for a flat liquid disordered DPPC membrane at 325 K reduces to π = 8.8 for a flat liquid ordered DPPC : cholesterol (58 : 42) membrane at 303 K and to π = 3.7 for a flat gel DPPC membrane at 285 K, illustrating both the condensing effect of cholesterol and the difference between a liquid ordered and a gel membrane. Here, we analyse the defect characteristics using the last 1 μs of the simulation trajectories, and only compare π values determined within this study. As a result of considering long trajectories, the standard deviations of the extracted π in this study are negligible.

3 Results

3.1 Free energy of insertion and dimerization: the role of lipid order

Hydrophobic C60 will spontaneously aggregate in an aqueous environment, to reduce their exposure to water, and form stable, irregular n-C60 aggregates of mesoscopic dimensions. Since the process of C60 aggregation varies with experimental conditions, the n-C60 shape, stability and size cannot be predicted, and we consider individual C60s in the free energy analysis, in agreement with standard practice. With a size of less than 1 nm diameter, they are, unlike n-C60, able to spontaneously permeate into lipid membranes. After insertion, they exhibit a tendency to position themselves in regions of high density in order to maximize the dispersion interactions with the surroundings. As a result, PMFs display a minimum away from the low-density interstitial space,³² showing that the perturbation of the lipid matrix by a single C60 does not lead to a significant entropic loss. This position of the PMF minimum, d_min, relates to the maximum density on average and thus to a particular conformational sampling of the lipid over time, meaning that it is sensitive to the acyl tail length (N_a) and the number of unsaturated bonds (b_u) per tail. To understand how different lipid states that are present in biological structures such as the lung membrane affect insertion, we calculated the PMF for a single C60 in the gel and liquid (liquid disordered) phase, see Fig. 1A for a schematic overview. Although a liquid ordered membrane state differs from a gel, we consider the gel and liquid states of a pure DPPC membrane representative for the liquid ordered and disordered domain in the lung membrane, see Methods and discussion section for details. The average membrane thickness, as determined from the positions of PO₄ beads in Martini, is 4.0 nm at 325 K and 4.8 nm at 285 K.

Fig. 1 (A) Cartoon highlighting lung surfactant membrane structures as a complex of multiple parallel lipid mono- and bilayers connected by strongly curved domains. (Top and bottom) Zooms of (A) illustrating the different phases, i.e. liquid disordered (L_d) and ordered (L_o), and curvatures that coexist within the structure. The flat monolayer closest to the air is in the ordered phase. (B) PMFs for C60 insertion into flat DPPC lipid bilayer in the gel (285 K, blue) and liquid (325 K, red) phase, as a function of their COM separation distance. (C) PMFs for C60–C60 interaction inside a DPPC lipid bilayer in the gel (285 K, blue) and liquid (325 K, red) phase, as a function of the separation distance. The quasi two-dimensional nature of these PMFs stems from our choice to keep the transverse coordinate of the C60 NPs, related to their depth within the membrane, fixed. The grey shaded regions in both (B) and (C) correspond to statistical uncertainty determined by a bootstrapping technique of GROMACS. Significantly longer equilibration times for T = 285 K give rise to an increased uncertainty compared to T = 325 K, albeit that this increase is modest.

As expected, the PMF for the liquid phase at 325 K, red line in Fig. 1B and Fig. S1B, aligns well with published Martini CG results for the same temperature.³³ In particular, this PMF reproduces prominent features of the atomistic PMF of Qiao et al.¹ such as the shape and the position of the minimum, but also the small energy barrier around 3 nm. From the position of the lipid head beads, at 2.05 nm, one may conclude that this barrier can be attributed to the hydrophobic C60 starting to loose their hydration shell. The passing of the domain rich of hydrophilic head groups, however, is not visible in the PMF, a finding that agrees well with the idea that permeation is a barrier-free process that takes place via spontaneously formed short-lived membrane defects of an appropriate size.

Next, we consider the insertion of C60 in a more ordered gel phase at 285 K. Lowering T below the gel transition temperature renders the DPPC membrane slightly more hydrophilic on average, because lipid defects, which are known to regulate the exposure of the hydrophobic interior of the membrane, significantly reduce in size. Consequently, one may conclude that the barrier for insertion into the membrane interior increases due to the denser hydrophilic domain that C60 need to cross.¹⁴ Moreover, the entropic costs of accommodating C60 within the DPPC matrix may rise, owing to a tighter lipid packing in the gel phase. The calculated PMF, see the blue curve in Fig. 1B, indeed shows a new but tiny barrier around 2.5 nm, i.e. close to the head group domain, signalling the reduced defect size. Overall, the insertion energy is, however, found to increase, see Fig. 1B, by ΔG = −1.38 kcal mol⁻¹. Apparently the perturbation of the lipid matrix by a single C60s can be easily accommodated even in the gel phase, and the enhanced dispersion interactions due to a higher local density dominate the insertion energy. The location of the minimum, which is the same in the gel as in the liquid phase, can be explained by the fact that the change in conformational sampling with temperature is limited for a fully saturated lipid such as DPPC.

To evaluate whether C60s remain aggregated within a DPPC membrane or disassemble into single C60 after entry, we calculated the PMF for dimerization for the two membrane states. As with the C60 position, the stability of C60 aggregates upon incorporation into the membrane has not been entirely resolved. In particular, C60 generally enters the membrane as an aggregate, in line with the dimerisation energy of about 3.6 kcal mol⁻¹ in water that was identified by atomistic calculations and reproduced on the CG level.³⁴ Nano-sized aggregates of C60, i.e. n-C60, were found to be stabilized in water by electrostatic repulsion as a result of the negative surface potential that builds up on colloidal aggregates, an effect that recently has been attributed to charge polarisation.³⁵ Previously, it was shown that adding C60 to a POPC bilayer environment reduces their propensity to dimerise compared to octanol,³⁶ which was attributed to enthalpic contributions, in particular the weakening of dispersion interactions that is the result of aggregation within a lipid matrix.³⁶ As also discussed by the authors, decomposing a PMF into enthalpic and entropic contributions suffers from errors, as this free energy is generally a result of adding two large and opposing terms,³⁷ making it hard to quantify uncertainty.

Our PMF profile for a pair of C60s as a function of the distance between their COM in a liquid membrane of DPPC, see curve for T = 325 K in Fig. 1C and Fig. S1D, is very similar to the one obtained for POPC at 313 K in Barnoud et al.³⁶ and at 300 K by Lavagna et al.³⁸ A free energy gain in the first minimum, when C60 are brought into contact, is almost absent, and even smaller than the tiny −1.5 kJ mol⁻¹ determined for POPC.³⁶ There is an absolute minimum of ΔG = −0.56 kcal mol⁻¹ at 1.54 nm, corresponding to a C60 pair separated by a DPPC lipid, and this and the dimeric case are set apart by a barrier of 3.52 kcal mol⁻¹. These results suggest that an aggregated C60 state is less likely stable upon insertion into a DPPC membrane, in agreement with what was already concluded for POPC.³⁶ In the gel phase, at 285 K, the locations of PMF minima are slightly displaced compared to the liquid membrane, reflecting the change in lipid packing. Interestingly, the stability of C60 dimers in the gel phase is significantly enhanced, with the minimum of ΔG = −4.01 kcal mol⁻¹ when C60 are at close proximity (1.03 nm). The depth of this well is even more pronounced than for octane in the study of Barnoud et al.,³⁶ which features the highest amount of aggregated C60s of all considered media for unconstrained simulations with increasing C60 concentration. Yet, the increased stability of dimers at lower temperature appears inconsistent with the previously identified significance of enthalpic contributions due to the lipid matrix.³⁶ As we performed our PMF calculations with the greatest care, we conclude that subtleties in the considered approaches apparently do not allow for simple extrapolation.³⁷ Based on our free-energy calculations, the gel phase of DPPC is expected to favor stable C60 aggregates within the lipid bilayer, while the liquid DPPC phase promotes disassembly into monomeric C60 within the membrane.

3.2 Curvature and lipid packing in the absence of C60

First, to ascertain that our computational treatment indeed deals with the two different phases, we consider the fluidity of purely lipidic flat and curved membranes. Curvature does play a role in the diffusion of individual lipids and membrane-bound macromolecules like proteins, and separating the role of curvature and of the projection used for extracting lipid trajectories on a curved geometry is far from trivial, both in experiments and continuum descriptions, which may lead to anomalous values.^39,40 For flat membranes, we calculated the lateral diffusion coefficients D for DPPC from the slope of the mean square displacement (MSD) versus time using the standard procedure implemented in Gromacs, see Fig. 2B, which is applicable for this geometry. We found D^liquid_DPPC = 6 × 10⁻⁷ cm² s⁻¹ for the liquid phase and D^gel_DPPC = 2 × 10⁻⁹ cm² s⁻¹ for the gel phase, in good agreement with values obtained by experimental measurements and for an earlier version of CG Martini.⁴¹ The lateral diffusion rate determined by fluorescence experiments in the gel phase, in the order of 10⁻¹¹ cm² s⁻¹,⁴² is lower, but these measurements were performed for supported bilayers, and the chemical nature of the support was found to play a role on the lipid kinetics. Overall, the two orders of magnitude difference in diffusion rates going from the liquid to the gel phase is to be expected. On a (strongly) curved geometry like our buckled membrane, geodetic instead of Cartesian distances should be considered. Application of an appropriate Vertex-oriented Triangle Propagation algorithm implemented in CurD⁴³ showed that lipid mobility in curved structures is most sensitive to the headgroup packing.

Fig. 2 (A) Simulation snapshots of pure DPPC membranes for the four considered setups: flat membranes in the gel (285 K) and liquid (325 K) phase, and buckled membranes in the two phases. Lipids are shown in a stick representation: hydrophobic groups are coloured distinctly for the four cases, hydrophilic groups are grey. (B) Mean-squared displacement (MSD) versus time for DPPC lipids in the flat membranes, obtained using the standard projection implemented in Gromacs. Curves were obtained for the last 2 μs of a 20 μs production run, averaged over all lipids in the both the leaflets. (C) Percentages obtained by machine learning (ML) of lipids belonging to a gel or liquid phase for the four considered setups. The ML analysis considers individual conformations to decide whether the lipid is in the liquid or gel phase. (D) Accumulated defect distribution for all membranes setups. (E) The defect constant π and the orientational order parameter P₂ with respect to the director orthogonal to the membrane for all setups. The defect constant π was obtained from fitting the defect distribution versus defect area shown in panel D to a single exponent.¹⁴ Analysis of individual cases was performed using the last 1 μs of a 10 μs production run, and obtained as an average over all lipids in both leaflets.

We selected an alternative approach for analysing the role of curvature in the lipid packing characteristics. Machine-learning-assisted lipid phase analysis (MLLPA) was recently developed to analyse membrane phases on the basis of individual molecular conformations that are subsequently clustered into lipid domains,²⁹ see Method section for details. While lipid phase domains and dynamics are not necessarily one-to-one related, as exemplified by the existence of liquid ordered membrane domains or rafts with overall liquid-like dynamics, the situation for single component membranes is quite simple. Curvature-related changes in lipid diffusion were only observed for lipid mixtures, and attributed to the sorting of more flexible lipids, i.e. lipids that prefer disorder, to curved parts of the membrane.⁴⁴ For liquid (L_α) membranes of pure POPC, which is incapable of phase separation in the considered temperature range, curvature was found not to notably change diffusion rates. Analysing the buckled gel membranes using MLLPA, we find that the fluid phase increases more than twice (from 6.2 to 13.2%) in size, while the buckled liquid membrane shows fluid/gel percentages that are clearly insensitive to curvature. The latter correlates well with the earlier result that lipid diffusion does not change with curvature when considering single-lipid (POPC) membranes in the liquid phase.⁴⁴

In addition to the chain statistics, we also considered lipid membrane defects, as they provide entry points for C60 and affect the binding kinetics by tuning the apparent hydrophobicity of the entire membrane. Packing defect constants or average defect sizes π for our setups were determined by extracting and fitting the defect size distributions using a recently developed procedure that is applicable to membrane of arbitrary shape, see Methods section for detail.¹⁴ As before, the characteristic defect size π in the gel is significantly reduced compared to the liquid phase. This is consistent with lipid defect formation being a transient and local effect due to spatially uncorrelated fluctuations. The π values determined for our flat pure DPPC membranes are π = 12.8 ± 0.3 Ų in the liquid phase and π = 3.5 ± 0.1 Ų for the gel, signifying the reduced conformational sampling in the gel phase (Fig. 2D). Upon buckling, the defect constants increase to π = 17.5 ± 0.3 Ų for a DPPC liquid and π = 19.4 ± 0.5 Ų for a gel, consistent with earlier findings and showing that the order is reversed. We will use these values as reference for the C60 inserted case. While the defect constant mildly increases with curvature in the liquid phase, in line with an earlier identified linear relation with inverse radius,⁴⁵ this increase is much more dramatic in the gel phase. Previously, this finding was related to π for the entire membrane being dominated by domains having a large π, with lipid packing defects in the positively curved region of the gel membrane, where packing is frustrated, becoming large.¹⁴ The orientational order parameter P₂ shows an opposing effect, demonstrating higher order in flat membrane at 285 K and a significant reduction of lipid order for buckled membrane at the same temperature. The π value in flat parts of the buckled gel membrane is consistent with that of an entirely flat membrane.¹⁴ Overall, both MLLPA and lipid defect analysis indicate that curvature plays a critical role in the gel phase, perturbing chain statistics and increasing the defect number and their size, suggesting that the effect on defect-dependent or defect-driven mechanisms should be most critical to a DPPC membrane in the gel phase.

Spatially resolved curvature information, complementing Fig. 2A, is generated in the form of two-dimensional color maps depicting the mean curvature for all the systems. To assess the stability of the imposed curvature throughout the simulation, curvature profiles time-averaged over 1 μs intervals, along with the corresponding values averaged over the final 5 μs, are presented in Fig. S2 and S3. The curvature plots clearly shows the differences in the initial stages of the considered membrane systems with respect to their phases and morphology. For flat membranes, the liquid phase exhibits slightly lower curvature compared to the gel phase, while the buckled membrane in the gel phase displays the highest degree of deformation due to its buckled structure. Although one would expect the membrane to become more stiff in the gel phase, i.e. reducing curvature fluctuations, there is some evidence that the situation is more complex.⁴⁶

3.3 Nano-C60 insertion depends upon the lipid membrane phase

One could wonder how relevant the standard free energy analysis performed for introducing (individual) C60s in a membrane actually is, knowing that there is a considerable driving force for C60 aggregation in water. Indeed, the ∼4 kcal mol⁻¹ gain for dimerisation in water^47,48 confirms that the stable nano-aggregates of C60 identified in experiments are more realistic representatives of the actual situation prior to binding. Yet, the gain of ∼20 kcal mol⁻¹ for the membrane-inserted C60 state determined in the previous section tells us that – at least energetically – individual C60s prefer to insert themselves into the membrane rather than into an already formed aggregate in solution. The free energy gain associated with the membrane insertion of aggregate C60s, i.e. n-C60, may be expected to vary with cluster size and shape, but experimentally n-C60 has been identified to spontaneously translocate. While the computational study of insertion energies for aggregates is complicated by the extent of the size and shape design space, in addition to n-C60 (in)stability within the membrane, a free energy gain was indeed observed for the insertion of a single carbon-based spherical NP of 3.3 nm radius into a liquid POPC membrane.³⁸ Consequently, the initial C60 concentration in the solvent phase and their distance to the membrane are important simulation variables, since they set the size to which the C60 aggregate can grow in water and thus the degree to which computational studies mimic the actual pre-aggregation process in silico. Most simulation setups, however, balance the degree of solvation, i.e. the number of waters per lipid, needed to stabilize a bilayer and minimize the computational load, therefore positioning C60 close to the membrane, e.g. 1–2 nm away, which promotes membrane binding and suppresses the aggregation of C60 into larger n-C60 in solution. As a result, experimental aggregate sizes are effectively ruled out in simulations.

Our computational procedure is carefully tuned to provide larger n-C60 aggregates in solution prior to membrane binding and insertion, see details in the Methods section. As suggested by the dimerisation energy in solution, all C60 are seen to aggregate quickly into n-C60 of varying degrees of polymerisation (n) during the first 100 ns of simulation. Growth proceeds irreversibly, meaning that the size of these initial n-C60 aggregates increases depending on available C60 and is dictated by the proximity of each C60 to the closest nucleus, and thus by the C60 starting distribution, hence our choice to perform two simulations for different C60 concentrations. With time, n-C60s can be seen to grow further by fusing with other aggregates. Persistent shape irregularity shows that remodelling into an energetically more preferred spherical shape by C60 surface hopping, a mechanism that underlies the construction of equilibrated metal–oxide NPs, for instance by the Wulff construction,⁴⁹ does not take place, at least not on these time scales. With time, these irregular n-C60s are found to diffuse to the bilayer interface and bind, see Fig. 3. In the remainder, we make a distinction between aggregate NPs like n-C60, which form by assembly of smaller NPs but may also disassemble into smaller NPs under particular conditions, and standard NPs of constant size.

Fig. 3 Simulation snapshots after 20 μs of CG MD simulation of pure DPPC membranes with 640 C60 for the four considered setups: (A) flat membrane at 285 K, (B) flat membrane at 325 K (C) buckle membrane at 285 K, and (D) buckle membrane at 325 K. The NC3, PO4, GL1 and GL2 CG beads in DPPC lipid are shown as blue stickes (flat, gel), red sticks (flat, liquid), green sticks (buckle, gel), orange sticks (buckle, liquid) and all other atom types are omitted for clarity. The CG beads in the C60 are shown as grey spheres.

Analysing the trajectories of the unbiased simulations, see simulation snapshots in Fig. 3 and movies in the SI and evolution of aggregation number and contact fraction in Fig. 4, we observe distinct pathways of n-C60 insertion. The initial stage – nucleation and growth of aggregates from monomeric C60 that are dispersed in the solvent – should be rather indifferent of the membrane state, since aggregation takes place well away from it. Growth is regulated by the separation between mono- or multimeric C60, which act as seeds for further aggregation, and their diffusion rates, with seed diffusion slowing down with size or aggregation number n. Consequently, one can imagine that both the initial C60 concentration and temperature are important parameters for aggregate growth prior to membrane binding, with a higher concentration or an elevated temperature providing accelerated growth. The initial stages of the curves that evaluate the percentage of C60 monomers with time, see inset in Fig. 4, confirm that concentration is an important parameter, since a higher concentration for all cases leads to an increased loss of C60 monomers. Yet, these curves also show that temperature does not play a determining role in reducing the monomer survival rate during aggregate growth in solution. After this stage, which can be short-lived, the C60 solution gets depleted of smaller seeds. Next, there is a competition between further coarsening in solution and the removal of seeds for growth from solution by entering the membrane.

Fig. 4 Percentage of free C60 monomers versus simulation time in μs for four considered systems: (A) 320 and (B) 640 C60s. Inset represent the early stages of free monomer composition for all the cases, capturing the initial aggregation behavior. Contact fraction between C60 and DPPC for the four considered systems: (C) 320 and (D) 640 C60s. Here, a C60 “monomer” is defined as two C60 molecules separated by at least 1.08 nm. Colors relate to the following membrane state: flat gel (blue), flat liquid (orange), buckled gel (green) and buckled liquid (red).

The subsequent stage, namely the interaction of an aggregated nanoparticle (NP) with the membrane, can be compared to the corresponding case of individual NPs. The mechanism of binding and insertion of impartible NPs is known to vary considerably with size and hydrophobicity, ranging from spontaneous insertion for all smaller (radius ∼1 nm) NPs like C60, to membrane remodelling via wrapping for more hydrophobic NPs of sizes that are comparable to or slightly larger than the membrane thickness, while less hydrophobic NPs of this size range can still enter via insertion.⁵⁰ Material-based simulations for even larger NPs, i.e. between 10 and 50 nm, are currently computationally intractable even at the CG level, but continuum models suggest that the wrapping mechanism becomes more dominant.⁵¹ Wrapping, which may be only partial or full, depending on the energy of adhesion, can eventually lead to NP uptake into the membrane.⁵⁰

The role of shape irregularity has remained somewhat unexplored, but it may become a factor when the characteristic dimensions of a particular NP irregularity become compatible with the average size of lipid defects. The switching between binding and insertion mechanisms, and also the resulting equilibrium state, depends on a delicate balance of factors pertaining to both the membrane and the NP. For the NP, size and hydrophobicity, as reflected in the adhesion energy density, are key factors; for an aggregate NP, the ability to reorder or even disassemble into smaller parts are additional factors. These factors should be offset to membrane features such as the detailed lipid structure and overall membrane elasticity. Particularly wrapping implies the formation of high local curvature, which generally comes with an energetic penalty that relates to the lipid type. Commensurability between membrane defect structure and NP features may generate suitable hydrophobic pockets for enhanced binding. As an example, lipid defects were previously shown to induce a secondary structure transition in the hydrophobic domain of a membrane-active protein, by exposing a commensurate hydrophobic pocket for binding.⁵² Defects may also affect the overall hydrophilicity of the head group domain that poses an energetic barrier for insertion of hydrophobic NPs.⁵³

Binding and insertion for flat membranes. Starting with a flat gel membrane, i.e. the case that resembles the monolayer present next to the gas interface in the lungs, aggregate NPs (or n-C60) start to form in the solvent phase within the first 200 ns. A few monomeric C60 that are present in the solvent very close to the membrane are quickly absorbed to the interface and, with time, overcome the barrier for penetration, in line with the calculated PMF for this case. Translocation into the membrane interior is fast enough to enable these monomers to escape from the aggregation process that continues to take place in solution. In the solvent phase, n-C60 remain partially or entirely exposed to other n-C60 due to their slow permeation rate into the membrane, and continue to grow until only a few larger absorbed ones remain.

The absence of induced curvature via wrapping for all n-C60s suggests that the adhesion energy density for this material is rather low or that the non-spherical shape of the n-C60 plays a role. Adhesion is clearly not strong enough to overcome the penalty for locally curving the gel membrane, which is typically associated with a higher bending rigidity than the liquid membrane as well as smaller defect sizes on average. Plateaus in the contract fraction, see Fig. 4C, show that C60 aggregates remain intact upon partial incorporation into the membrane, in line with the PMF for dimer formation at this temperature discussed earlier. Close inspection shows that only C60 present in n-C60 at sharp edges can match the small defect size of hydrophobic pockets and are seen to disassemble from the n-C60 to which they are attached, finding a position within the membrane. This is reflected in a tiny slope of the monomer percentage. The low probability of such events, combined with slow transport within a gel membrane, hampers clustering of C60 monomers inside the membrane on the time scale of the simulation, despite the tendency to cluster observed in the earlier calculated PMF. Re-aggregation of these monomers to the absorbed n-C60 is even observed. As mentioned, prolonged (partial) exposure of these absorbed n-C60 to the solvent phase is seen to promote the growth of n-C60 in that phase by merging of smaller ones, see Fig. 3A.

For the flat liquid membrane, i.e. at a higher temperature of T = 325 K, solvated n-C60s are seen to quickly and spontaneously insert themselves fully into the membrane. Removal from solution affects further growth of n-C60 in solution, and the aggregates that are formed prior to membrane uptake are generally smaller than at 285 K, increasing the probability that at least one of the dimensions of the inserted n-C60s is commensurate with the membrane thickness. Apart from the entry pores or pockets, significant local membrane deformation are not observed. On a longer time scale and only after full insertion, n-C60s start to disassemble into C60 monomers, the rate of which is seen to vary with the total C60 concentration in the simulation volume, see Fig. 4. As this rate relates to the ability of the membrane to store individual C60s, the total loading should be of importance as well as the barrier for C60 flip flops between the two leaflets, which is in the order of 5–10 kcal mol⁻¹. In both considered cases, the loading is high.

In short, in the absence of curvature, n-C60 insertion is found to correlate well with the lipid defect characteristics, in full agreement with the anticipated role of transient hydrophobic pockets in uptake.⁵² For the gel membrane, associated with a tiny defect constant π, aggregates bind but uptake is absent apart from very few individual C60 or tiny n-C60. For a liquid DPPC membrane, π is much increased and n-C60 insert spontaneously. Interestingly, the migration of lipids out of the membrane plane is never observed, while such a mechanism was previously identified during poration of a liquid membrane by a carbon nanotube to shield these NPs from water contacts by a lipid monolayer.⁵⁴ In spite of a different penalty on bending, both gel and liquid membranes are seen to respond to n-C60 binding by enhanced membrane oscillations, see movies in the SI, at a larger scale. We did not analyse the characteristics, but it is reminiscent of capillary waves. The absence of strong local curvature in the liquid membrane may be explained by an earlier observation that bound proteins are able to lower the barrier for insertion via the recruitment of lipid defects from their direct environment.⁵² Apparently, this route is also available for n-C60 uptake by a liquid membrane, but not for the gel membrane.

Binding and insertion for curved membranes. For our the buckled gel membrane, i.e. a curved membrane at 285 K, n-C60s aggregate with time in the solvent phase to bigger clusters, with only a few monomers binding to and inserting into the membrane. Once a single n-C60 is formed, it binds to the domain of positive curvature in the membrane, inserting themselves partially, see SI movie and Fig. 3C. The finding that part of this aggregate can spontaneously insert itself into the curved membrane, in contrast to the situation for the flat gel membrane, without significantly perturbing the membrane, is attributed to the increased defect size in these curved domains.¹⁴ Also the overall change of the membrane shape induced by n-C60 uptake remains minimal, with almost no lipids being pulled out, while the n-C60 does not fall apart in individual C60 inside the membrane. The n-C60 shape inside the membrane is stable, suggesting little lipid re-organisation, since the external membrane-bound part of the aggregate limits the overall diffusion of the entire aggregate. As a result, the membrane-inserted C60 concentration is and remains high where aggregates are inserting, i.e. in curved domains. On the timescale of the simulation (20 μs), C60s are not further dispersed within the membrane.

Albeit less striking, also for the liquid buckled membrane at T = 325 K there is a preference for binding in the strongly positively curved membrane domains, where defects are larger on average, see SI movie. Yet, in contrast to the gel case, aggregates insert themselves entirely into the membrane with time. Since the disassembly of n-C60 into monomers is again a slower process than n-C60 uptake, the membrane shape and curvature may be perturbed by the uptake of n-C60, especially for larger ones. This perturbation is enhanced at the location of entry, thus in the curved domains, and recruitment of surrounding defects is again suggested as a mechanism that facilitates this. In particular, significant membrane remodeling around the aggregate, by lipids that diffuse out to form an encapsulating monolayer, is absent in most cases. Only larger n-C60 that are incommensurate with the membrane thickness upon insertion become encapsulated by a lipid monolayer, very similar to the insertion mechanism for 5 nm diameter NPs in our earlier study.⁵⁰ This finding supports the conclusion that direct insertion via defects is responsible for n-C60 uptake, rather than a mechanism based on wrapping.⁵⁰

As before, the total C60 concentration is seen to play a role. Whereas n-C60 uptake for the lowest concentration takes place without significant remodelling of the buckle, close inspection for the highest C60 concentration shows that curvature is significantly enhanced in the already curved parts where several n-C60 enter the membrane. In particular, the high C60 loading in the curved parts induces additional local curvature and pushes the bilayers at either side of the buckle into fusion and the formation of a membrane structure reminiscent of a vesicle with a very tiny diameter, see Fig. 5. Once the process of n-C60 disassembly sets in, monomers start to diffuse to less populated parts of the bilayer, reducing the local stress. As a result, curvature decreases.

Fig. 5 Simulation snapshots at different times T for the buckled membrane in the liquid state (325 K) with 20 mol% C60. To clarify the induced membrane curvature, snapshots are shown with (left, A) and without (right, B) C60. Lipids are shown in a stick representation: hydrophobic groups are coloured orange while hydrophilic groups are coloured white. The CG beads in the C60 are shown as grey spheres.

The role of curvature for the two states of the buckled membrane is further quantitatively analysed via the number of C60-PO4 (relating to DPPC) contacts in the curved and flat parts of the buckled membrane over time, see Fig. 6. Simulation snapshot showing the curved (blue) and flat (red) regions of buckled membrane systems at the two concentrations and temperatures are shown in Fig. S5. Overall, this analysis confirms our earlier conclusion that n-C60 prefer to bind to curved domains of both gel and liquid membranes, where defects are larger in size on average. This is most striking for the gel state. A comparison of the data for the two C60 concentrations and the buckled gel membrane, see Fig. 6B (10 mol% C60) and D (20 mol% C60), suggests that the binding and aggregation histories are important. In particular, for the lower (10 mol%) C60 concentration, the initial preference of n-C60 for the curved parts of the buckled membrane is replaced by a preference for the flat parts at later times. This finding confirms that membrane binding is preferred over aggregate growth in solution.

Fig. 6 Number of contacts between C60 and DPPC (PO4) groups in the curved and flat regions of buckled membrane systems at varying concentrations and temperatures: (A) 10 mol% C60 at 325 K, (B) 10 mol% C60 at 285 K, (C) 20 mol% C60 at 325 K, and (D) 20 mol% C60 at 285 K. Analysis of each cases was performed over the whole 20 μs production run, and obtained as an average over all lipids in curved and flat region in both leaflets.

In particular, full occupation of the curved part by inserted n-C60 drives the remaining solvated n-C60 to unoccupied parts of the membrane, which is necessarily in the flat part. For the higher (20 mol%) C60 concentration, aggregates are larger in size on average prior to binding, and their number is thus significantly reduced, ruling out this phenomenon. The balance is clearly quite subtle. It should be noted that in the gel case, n-C60 are stable both in solution and within the membrane, meaning that they can only grow. For the buckled membrane in the liquid state, however, and for the lowest C60 concentration, the differences in binding to curved/flat parts appear consistent with the relative areal difference, suggesting that there is no clear preference. In this case, n-C60 are unstable within the membrane, albeit that dissolution is a rather slow process. For the higher concentration of 20 mol% C60, i.e. larger n-C60 in solution on average, there is an initial (weak) preference for the curved parts, compare Fig. 6B (gel) and C (liquid), confirming our observation that there is a slight curvature preference also for the liquid case. Despite defects playing a much less significant role in the liquid membrane, we conclude that exclusion of binding sites resulting from the competition between uptake and dissolution leads to a weak curvature preference even when all defects are sufficiently large.

The uptake of n-C60 by membranes with varying and/or strongly curvature is rather unexplored, both owing to limitations in the experimental resolution and to restrictions in the volume that computational methods are able to address. A CGMD study of bicelles with constant negative/positive curvature apart from the edges identified a preferred C60 penetration in regions of positive mean curvature, both for symmetric and asymmetric compositions.⁵⁵ After insertion, C60 distribute almost evenly in the curved regions, but they strongly avoid regions with extreme positive or negative curvatures and do not induce curvature themselves up to very high loading. Yet, C60 were found to significantly increase the order of lipid tails with which they are in direct contact, leaving the order of other lipids unaffected. Experiments of Zupanc et al. showed that n-C60 aggregates can modify the average mean curvature of a liposome and even destabilize them via rupture.⁵⁶ This was confirmed by a study of Bouropoulos et al., who identified shape distortion and a general increase in liposome size.⁵⁷ Electrostatics was also found to play a role, with n-C60 disrupting positively charged small unilamellar vesicles (SUV) but not their negatively charged counterparts.⁵⁸ The diameter of the liposomes, and consequently also the range of curvatures, that can be evaluated in silico is restricted if one wants to avoid artifacts stemming from the (periodic) boundary conditions, particularly because a substantial water domain is needed to solvate the C60 and to allow them to aggregate prior to binding. Zhang et al. used Martini 2 to consider n-C60 binding to DPPC liposomes of 16.5 nm diameter.⁵⁹ Their PMF, again calculated for a single C60, was found substantially modulated by curvature, which reflects the asymmetric perturbation of the dispersion interactions stemming from the distribution of lipid tails in the two membrane leaflets, meaning that also the location of the C60 position is (asymmetrically) perturbed with respect to the midplane. For higher concentrations of C60, the DPPC vesicle was found to deform into an oblate shape by the binding of n-C60, introducing edges of higher curvature. Moreover, n-C60 were found to slowly disassemble inside the membrane, consistent with a liquid membrane simulated using Martini 2 for DPPC at T = 298 K.⁵⁹ Owing to the substantial computational requirements, the selected water content was very modest, and only the initial stages (200 ns) of the interactions between the liposome and n-C60 were simulated.⁵⁹ The observed increase of local curvature upon n-C60 insertion is in line with existing experimental and computational findings.^56,57,59 In some cases, even rupture was identified. As before, the membrane state and the final distribution of C60 within the membrane depends on the maximum C60 loading in the membrane. Nevertheless, at T = 325 K, the C60 distribution appears rather homogeneous in both liquid membrane states, see Fig. 3.

3.4 Quantifying n-C60 binding and dispersion

The earlier discussed competition between C60 and cholesterol in lipid packing suggest that their role is related. The three basic properties of cholesterol in a liquid bilayer – even distribution with curvature,⁶⁰ a large molar faction that can be loaded while maintaining membrane stability, and increasing the lipid tail order⁶¹ – have indeed previously also been observed for small NPs like C60.⁵⁵ In addition, cholesterol contributes to making the membrane thicker and more impenetrable, it regulates membrane fluidity by softening the main (gel/liquid crystalline) phase transition, and changes transport properties of proteins within membranes. A quantitative analysis of adding C60 to the lipid organisation thus makes sense, particularly after the C60 have dispersed within flat and curved membranes. Consequently, it makes sense to relate those findings to data for the gel case. A true comparison of the role of fillers like C60 and special lipids like cholesterol requires more detailed investigation and is only discussed in the Discussion section.

First, to understand how C60 disperse within the membrane, we have monitored the percentage of monomers during simulation within the entire volume. Fig. 4A and B clearly show that C60 fully dissolve inside the liquid membrane, in line with earlier findings,⁵⁹ at a rate that is sensitive to the total C60 loading, as discussed in some detail before. Moreover, we observe that the dissolution rate is generally reduced by curvature. Although dissolution hardly takes place in the gel membrane, also in this case the loading fraction plays a (minor) role, and the dissolution rate is slightly enhanced by curvature. This suggests that the perturbation to the lipid conformations by curvature very weakly contributes to solubilizing C60 on longer time scales. The contact fraction correlates well with the n-C60 dissolution kinetics for liquid membranes, and suggests that dissolution in the gel membrane does not correlate with membrane binding of n-C60, Fig. 4C and D.

While defects clearly play a role in n-C60 binding from solution, a finding that is supported by the significant increase in average defect size π with curvature both in the liquid and in the gel membrane, see Fig. 7A and B, our next interest lies in understanding how C60s affect the average defect size at longer time scales that are important for binding to C60-loaded membranes. Here, it should be noted that the average defect size is determined by statistical analysis of the defect size distribution, which requires statistically significant data – a large membrane area and long time sequences, see examples of distributions in Fig. 2D without C60 and Fig. 7 with C60 – to properly capture larger defects. Defect constants are therefore always the result of averaging and instantaneous values have no meaning.

Fig. 7 Analysis of the defect constant for all considered membranes after C60 incorporation: (A) accumulated defect distribution for all membranes with 320 C60 (10 mol%), (B) accumulated defect distribution for all membranes with 640 C60 (20 mol%), (C) defect constants π derived an exponential fit of the defect distributions, for flat and curved parts in the buckled gel membrane at T = 285 K, (D) fitted defect constants π for flat and curved parts in the buckled liquid membrane at T = 325 K.

Fig. 7A and B shows the defect size distributions in case when 320 (top left) or 640 C60 (top right) are added to the four states of the membranes. Fig. 7C and D enable a comparison of the average defect size π of the flat and curved domains in the buckled gel and liquid membrane upon different C60 loading. The values for buckled membranes without C60, utmost left in both panels, show that the increase in defect size is limited to curved domains, as π for flat domains are equal to π for entirely flat membranes, see Fig. 2E, in line with earlier findings.¹⁴ For the liquid membrane in which all C60s are sufficiently dispersed, the average defect size in the flat parts remains roughly the same with increasing C60 loading. In curved parts, however, π-values are seen to increase, reflecting the role of C60 in perturbing the conformational sampling of lipids under curvature. Since the loading is high, C60 spread through the membrane and both negative and positive curvatures appear affected; one should notice that the restricted surface area of negatively curved parts may lead to some undersampling, leading to mild imprecision in the determination of π. In defect analysis, fillers like C60 are considered as voids, whereas cholesterol is explicitly accounted for in the procedure.¹⁴ For cholesterol, the condensing effect in liquid POPC and mixed membranes was previously considered by atomistic simulation, and the average defect size π in flat membranes was seen to mildly reduce upon the addition of cholesterol in small quantities.^62,63 Surprisingly, a private and unpublished study showed a weak but opposite (increasing) effect on π using the coarser Martini description of the POPC/cholesterol system, a finding that reflects the subtleties of capturing packing effects at the coarser level.⁶⁴ Here, the average defect size in both the flat and curved parts of liquid membranes can be seen to increase with C60 loading, see Fig. 7D. This increase is most prominent in the curved domains, where lipid conformations are already affected by curvature.

It makes sense to correlate these findings with information about chain conformation statics via the earlier mentioned MLLPA method, which determines the ‘effective’ fluidity of an entire membrane via instantaneous conformations of individual lipids. In the absence of C60, we find that the effective fluidity of liquid membranes is rather insensitive to curvature, while curvature renders gel membranes clearly more fluid, see Fig. 8 and Table S1. Adding C60, a flat liquid membrane becomes effectively more fluid for a low (10 mol%) percentage and about equally fluid for a higher loading fraction (20 mol%). For the buckled membrane, the fluidity is conserved (10 mol%) and reduced (20 mol%). A similar trend is observed in the area per lipids, see Table S1.

Fig. 8 Lipid ratio analysed to be part of the fluid phase for three different C60 ratios X = 0, 0.1 or 0.2 in a DPPC bilayer. The considered systems were bilayers at 285 and 325 K both in the flat and buckled morphology. Fluidity was determined from instantaneous lipid conformations using the MLLPA method. MLLPA was used as is, meaning that it was trained using 20% of the data for gel and liquid phases in pure flat DPPC membranes.

Combining defect and fluidity analysis for liquid membranes, i.e. T = 325 K, we may draw two important conclusions. For a flat membrane, the nestling of C60 inbetween lipids increases the average defect size, but it has no significant effect on the lipid tail order or it even reduces that order. For curved membranes, the percentage of gel-like (stretched) conformations associated with increased tail order increases notably, a finding that should be attributed to the curved domains. This illustrates the role of C60 in lipid stretching upon curvature. As a consequence also the average defect size in curved domains increases significantly, in line with earlier findings.

For the gel membrane at T = 285 K, n-C60 disassembly of any significance is absent and n-C60 bind at the interface of flat membranes or insert partially in the curved regions of the buckled membrane. As a result, the defect (Fig. 7C) and fluidity analysis (Fig. 8) is dominated by membrane deformations that are a consequence of n-C60 binding or insertion, such as a local distortion of the lipid order (upon binding) or the formation of very large defects or holes for insertion, see Fig. 7C. As the binding or insertion kinetics appear static on the time scale of the simulation, this analysis is specific to the setup and does not provide additional information.

The spatially resolved curvature color maps for 20 mol% C60 at the end of the simulation trajectories visualise membrane deformations resulting from the long-term interactions between n-C60 and the lipid membrane, see Fig. 9. In the flat liquid membrane, solvated n-C60 molecules rapidly and spontaneously insert themselves fully into the membrane, and induce no significant long-lived membrane undulations (Fig. 9B). For the flat membrane in the gel phase, the stability of the n-C60s bound at the interface generate substantial undulations, see Fig. 9A, the amplitude of which is two-fold increased compared to the flat gel membrane without C60. In the case of the buckled gel membrane, the n-C60 inserts partially in the positively curved region and induces significant local membrane deformation, see Fig. 9C. Although we observe a significant increase in curvature upon the insertion of n-C60 in curved domains of the liquid buckled membrane, see the time series in Fig. 5, the mean curvature at the end of the simulation is equivalent to the situation without C60, compare Fig. S2D and Fig. 9D. Apparently n-C60 induces only transient curvature in the already curved liquid membrane, while they permanently deform the membrane in the curved gel phase.

Fig. 9 Mean curvature H = ½(κ₁ + κ₂), with κ_i = 1/r_i and r_i the radius of the matching sphere, for the four different membrane setups in the presence of 640 C60 (20 mol%): (A) flat gel, (B) flat liquid, (C) buckled gel, and (D) buckled liquid. In all cases, data was collected for an equilibrated membrane and analysed using the standard tool for extracting membrane curvature, “g_lomepro” with bandpass filter of q_high = 2.0 nm⁻¹. Shown curvatures were determined by averaging over the last 1 μs of a 20 μs production run.

4 Discussion

4.1 Lipid detail versus C60 insertion depth

The finding that the position of the PMF minimum is roughly equal for the flat gel and liquid DPPC membranes, d_min = 1.1 nm, is an interesting one, as it suggest that the dispersion interactions acting on C60 do not change significantly with phase. Only the depth and steepness of the attractive well around d_min can be see to change mildly. Meanwhile, d_min is sensitive to the force field and lipid specificity, denoted by (N_a : b_u), with N_a the acyl tail length and b_u the number of unsaturated bonds per tail. We shortly review literature for DPPC and other lipids. AAMD has provided d_min for DMPC (14 : 0) membranes of 0.6 nm at 310 K,³² 1.1 nm for DPPC (16 : 0) at 325 K,¹ and 0.8 nm for POPC (16–18 : 1) at 298 K,³⁴ suggesting that the equilibrium position is at least (weakly) dependent on lipid specificity. Using the CG Martini representation developed by Monticelli, C60 was found to position itself in a pure POPC membrane much further away from the membrane center, namely at d_min = 1.15 nm, than the previously identified 0.8 nm obtained by AAMD. However, since the membrane thickness is also slightly increased at the Martini level, and several key features of the atomistic PMF within the membrane match the CG PMF features quite well, such as the slope and free energy differences, the agreement was concluded to be satisfactory.³⁴ PMF calculations for DPPC at 325 K using an earlier coarse representation of C60 provided an ill-defined minimum around 1.0 nm (ref. 65) and an insertion energy or free energy difference between a membrane-inserted and solvated state of C60 of about 100 kJ mol⁻¹ (≈24 kcal mol⁻¹). The reference atomistic results of Qiao et al.¹ for DPPC display a much more pronounced minimum at 1.1 nm, and an insertion energy that is roughly 2.5 times smaller. Surprisingly, a later publication using the same CG Martini identified a much more pronounced PMF minimum around 1.1 nm and an insertion energy of slightly less than 20 kcal mol⁻¹,³³ suggesting that they used the updated CG representation for C60,³⁴ i.e. the one that we rely on, rather than the original one. The CGMD-derived value of about 20 kcal mol⁻¹ is consistent with insertion energies derived from experimentally C60 solubilities, between −12 and −45 kcal mol⁻¹, and the result of a polarizable continuum model (−24 kcal mol⁻¹) for the free energy difference between solvation of C60 fullerene in water and in n-hexane.⁶⁶ The finding of a considerably reduced insertion energy in atomistic calculations was previously attributed to issues with the considered force field to accurately capture solvation and partitioning.⁶⁵ Finally, for a POPC membrane at 300 K, gradually increasing the size of C60 whilst maintaining its chemical nature gives rise to a broading of the attractive well of the PMF and a minimum that is d_min ≈0.75 nm for F64 (1.87 nm diameter) to d_min ≈1.1 nm for F216 (3.31 nm diameter).³⁸ The most interesting finding, however, is that size tunes the stability of the NP aggregate within the POPC membrane, from little or no aggregation for nanoparticles with a diameter below 2 nm, more substantial aggregation in an aligned geometry for nanoparticles with intermediate diameter (3 nm), and strong aggregation into 3D clusters for the largest nanoparticles (5 nm).

4.2 Membranes preparation procedure for the different states

The transition temperature T_m between the liquid crystalline and gel membrane phase is determined by the lipid type, in particular the length of the acyl tails, the degree of tail saturation, and the charge and bulkiness of the head group. Attractive van der Waals interactions between adjacent lipids render lipids with longer hydrocarbon tails more likely to be in the gel phase. Unsaturated bonds change the chain statistics compared to fully saturated lipids, and generally gives rise to a reduced T_m, but also the position of double bond in the acyl chain can be important. While lipid packing is sensitive to the size of the headgroup, membrane fluidity also relates to the hydrogen bonding network close to and between headgroups, and internal or external factors, such a reduced local pH and higher membrane hydration, may change T_m for the same lipid composition. An extensive overview of all factor that play a role for PC lipids can be found here.⁶⁷ Membranes of saturated PC lipids like DPPC were experimentally identified to undergo three phase transitions with increasing temperature after equilibration at a low temperature. Starting from a lamellar crystalline subgel (L_c), the membrane transforms into a lamellar gel phase (L_β′) with tilted orientations via a so-called subtransition, then a rippled gel (P_β) phase that is a pre-transition to the lamellar liquid crystalline (L_α) phase, which develops from a P_β phase as a result of a highly cooperative process. The latter is the main or chain order/disorder transition. The ripple phase P_β, which is characterized by an alternation of interdigitated or disordered lipids and more regularly packed lipids and believed to be formed to relieve packing frustrations particular to some lipids, is either stable or metastable depending on the thermal history of the sample.⁶⁸ This sensitivity to experimental conditions is reflected in the difficulty of obtaining P_β phases by molecular simulation.⁶⁹ For the coarsened representation of Martini, the stability of lipid interdigitation is a particular issue that generally goes against finding such a phase.⁷⁰ Experimentally, the diffusion coefficient D was considered as an experimental descriptor for the membrane state, and a more than three orders of magnitude difference between the gel (L_β, T < 300 K) and liquid (L_α, T > 315 K) phase was measured.⁷¹ For 300 K < T < 315 K, the diffusion coefficient increases quasi-linearly with temperature, except for a plateau observed for a range between 306 and 312 K that was attributed to the ripple phase. The measured temperature for the main phase transition of DPPC is T_m = 314 K.⁷² Coarse-grained Martini 2 simulations for DPPC do reproduce the liquid L_α to gel L_β transition, but for a reduced temperature T^sim_m = 295 ± 5 K,^9,41 a discrepancy that may be attributed to the CG force field, amongst other to the state-dependence of the coarsening procedure. Our lower value of T = 285 K was selected with this reduced T_m in mind, meaning that we computationally consider two distinct membrane phases: a gel at T = 285 K and a liquid at T = 325 K. In addition to shifting the T_m, coarsening gives rise to enhanced diffusion rates because the thermal energy that drives the system is distributed over a reduced number of degrees of freedom. At 323 K, a four times speed-up was numerically found for the lateral diffusion of DPPC lipid compared to experimentally measured diffusion rates.¹⁶

4.3 Extrapolating model membrane findings to biomembranes

The feature shared by most biological membranes is the differentiation into finite and often dynamic domains that allow the membrane as a whole to function as an adaptive and sensitive barrier that allows transport. The often very diverse lipid content collectively enables such distinct emergent physico-chemical properties via partitioning, curvature formation, and other types of local restructuring, a process that originates from the instabilities that are imposed by thermodynamic forces on an entire (homogeneous) membrane. A well-known but yet debated example of such a phenomenon is the formation of liquid-ordered membrane domains or rafts that are enriched in particular lipids and that are anticipated to play a key role in protein binding.⁷³ Varying curvature was identified in various organelles. One of its suggested functions is to modulate the apparent hydrohobicity of the lipid membrane through the average size of lipid defects, a phenomenon that was proposed to play a role in molecular curvature recognition.⁷⁴ Computational studies of lipid membranes, however, face several inherent constraints, one of which is the choice of boundary conditions that typically enforce periodic structures. The study of curvature is particularly affected by this restriction of structural variability that can be explored. Analogous to cell membranes, lung membranes are both (strongly) curved and laterally heterogeneous, with more rigid domains (rich in saturated lipids/cholesterol) close to the interface with air and more fluid domains (rich of unsaturated lipids) in curved regions that are subject to compression, see Fig. 1A.

It is particularly interesting to investigate how distinct features of the lung membrane correlate with the simulated C60 insertion characteristics. An interesting final question is how well the different domains in the actual lung membrane, which is composed of DPPC, DOPC, POPC, and cholesterol in a ratio of 5 : 2 : 2 : 1,⁵³ are represented by our four model membranes of pure DPPC. The flat gel membrane was selected to represent the flat liquid ordered domain that is likely primarily composed of fully saturated DPPC and cholesterol. While adding small fractions of cholesterol is know to mildly perturb the defect characteristics of a phospholipid fraction, the mixture making up a liquid ordered has a defect constant π that is significantly lower than for a liquid disordered membrane but higher than for a gel phase. This may thus increase the insertion probability of n-C60 into the membrane, via defects of sufficient size. Yet, concentrating on the n-C60 stability within the membrane that is more important for the long term fate of the C60s, it should be understood that cholesterol decreases the tendency of C60 to disperse inside the membrane by competing with C60 for packing. As a result, we think that our gel membrane is a proper model for the actual binding. On the other hand, defects in the fluid phase are rather insensitive to lipid specificity and will contain only a minority of cholesterol, so the effect of considering a model membrane is minor. Finally, the previous study of Lavagna et al.,³⁸ identified a strong curvature sensitivity/selectivity for larger C60-like NPs in a POPC membrane, generalising the findings of the current study.

4.4 Additional insight

Combining the current results with information of binding mechanisms for individual NPs of larger size provides additional insight. The first one is the suggestion that n-C60 are not hydrophobic enough to be wrapped by the membrane, meaning that the gain in adhesion energy does not compensate the costs of locally strongly curving/deforming the membrane. Second, the spontaneously formed n-C60s in our setup are also not large enough to match the scale of low energy fluctuations of the membrane, which are always present and may induce (partial) wrapping for large NP. In reality, however, n-C60 are larger than the ones considered in our simulations, i.e. in the order of 100 nm or larger, well beyond the scales of CG simulation, and one may wonder if the possibility of wrapping should be considered after all. On the other hand, most of the studied mechanisms are for purely spherical NPs, whereas the natural irregularity of n-C60 shapes is likely to be important. In particular, with the presence of reduced dimensions due to shape irregularity, defects will remain important for binding and insertion phenomena. We thus expect our current setup to be sufficient for the investigation of general mechanisms.

5 Conclusions

We have considered the role of the membrane state and membrane curvature in the binding and subsequent uptake of hydrophobic nano-aggregates, i.e. nanoparticles (NPs) that are formed by the self-assembly of smaller NPs. As hydrophilic NPs are unlikely to spontaneously bind to a membrane, hydrophobicity is a shared property amongst membrane-interacting NPs, rendering the formation of NP aggregates in water a very common phenomena.

Here, we have focussed on the special case of nano-aggregates formed by C60 (called n-C60), which are individually so tiny that they are able to permeate the membrane spontaneously and do not to perturb it. This choice is inspired by the omni-presence of such chemicals in our surroundings and even in our own bodies. On the other hand, the membrane state (ordered and liquid disordered, or, in short, gel and liquid) and curvature are selected as discerning membrane features, because of the observation that such properties stand out in lipid barriers such as the lung membrane. A major complication for computational research of these mechanisms is the competition between n-C60 growth in solution – depending on C60 concentration, and generally rather slow – and the uptake of n-C60 by the membrane – depending on membrane phase and rather fast for a liquid membrane – favoring uptake and limiting n-C60 growth, whereas n-C60 in reality have already grown to mesoscopic dimensions prior to binding. Earlier work for NPs has shown that the binding and the insertion mechanism notably changes with size.

In this study, we have considered existing coarse-grained representations and simulation techniques to reach n-C60 and membrane dimensions that are needed for realistic simulation of insertion mechanisms. Focussing on the energetics, we find that single C60 insert spontaneously in both gel and liquid membranes. Similar analysis, however, shows that inserted C60s are dispersed in a fluid membrane and clustered in a gel membrane. Turning to C60 aggregates, which spontaneously form in water, we thus recognize that there is an important additional factor to be considered for the insertion mechanism.

Comparing the uptake mechanism for four different membrane states, gel or liquid, flat or curved, allows us to discriminate between the different factors. Focussing on flat membranes, we find that the fluid membrane is able to take up the entire n-C60s via a direct insertion mechanism, involving lipid defects and their recruiting from surroundings, after which the n-C60 disassemble at a slower time scale and disperse throughout the membrane. In flat gel membranes, on the other hand, n-C60 can only bind the membrane without insertion.

Introducing curvature, we find that n-C60 can suddenly insert themselves in the strongly curved part of the gel membrane, as a result of much inflated defect sizes, but they remain aggregated and do not further diffuse as a whole. Also for the liquid membrane we find increased binding to high curved domains, and insertion of large n-C60 may induce even higher local curvature, which relaxes once the C60 start to disperse throughout the membrane. For such a curved liquid membrane, we also observe the most pronounced effect of C60 on the lipid packing: particularly in the curved parts, C60 appear to stretch lipids, thereby generating more gel-like behaviour and much larger lipid defects than in the flat domains.

Conflicts of interest

There are no conflicts to declare.

Data availability

Trajectories movies and input file including final configuration files for all the cases with 10 mol% and 20 mol% C60 are made available via Zenodo (DOI: https://doi.org/10.5281/zenodo.15240543).

Supplementary information (SI) is available. See DOI: https://doi.org/10.1039/d5nr02822a.

References

1. R. Qiao, A. P. Roberts, A. S. Mount, S. J. Klaine and P. C. Ke, Nano Lett., 2007, 7, 614–619.
2. H. Johnston, G. Hutchison, F. Christensen, K. Aschberger and V. Stone, Toxicol. Sci., 2010, 114, 162–182.
3. N. Tagmatarchis and H. Shinohara, Mini-Rev. Med. Chem., 2001, 1, 339–348.
4. N. Tagmatarchis and H. Shinohara, Particuology, 2023, 73, 26–36.
5. J. Zhou, D. Liang and S. Contera, Nanoscale, 2015, 7, 17102–17108.
6. J. Wong-Ekkabut, S. Baoukina, W. Triampo, I.-M. Tang, D. Tieleman and L. Monticelli, Nat. Nanotechnol., 2008, 3, 363–368.
7. E. Castro, A. Garcia, G. Zavala and L. Echegoyen, J. Mater. Chem. B, 2017, 5, 6523–6535.
8. S. Sharma, L. Chiang and M. Hamblin, Nanomedicine, 2011, 6, 1813–1825.
9. S. Jaschonek, M. Cascella, J. Gauss, G. Diezemann and G. Milano, Biochem. Biophys. Res. Commun., 2018, 498, 327–333.
10. N. Nisoh, V. Jarerattanachat, M. Karttunen and J. Wong-Ekkabut, Biomolecules, 2022, 12, 639.
11. J. Barnoud, G. Rossi, S. J. Marrink and L. Monticelli, PLoS Comput. Biol., 2014, 10, e1003873.
12. R. Gautier, A. Bacle, M. Tiberti, P. Fuchs, S. Vanni and B. Antonny, Biophys. J., 2018, 3, 436–444.
13. S. Baoukina and D. Tieleman, Biochim. Biophys. Acta, 2016, 1858, 2431–2440.
14. R. W. van der Pol, B. W. Brinkmann and G. A. Sevink, J. Chem. Theory Comput., 2024, 20, 2888–2900.
15. D. Van Der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark and H. J. Berendsen, J. Comput. Chem., 2005, 26, 1701–1718.
16. S. J. Marrink, H. J. Risselada, S. Yefimov, D. P. Tieleman and A. H. De Vries, J. Phys. Chem. B, 2007, 111, 7812–7824.
17. T. A. Wassenaar, H. I. Ingólfsson, R. A. Bockmann, D. P. Tieleman and S. J. Marrink, J. Chem. Theory Comput., 2015, 11, 2144–2155.
18. K. J. Boyd and E. R. May, J. Chem. Theory Comput., 2018, 14, 6642–6652.
19. U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen, J. Chem. Phys., 1995, 103, 8577–8593.
20. M. Parrinello and A. Rahman, J. Appl. Phys., 1981, 52, 7182–7190.
21. G. Bussi and M. Parrinello, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 75, 056707.
22. T. Carpenter, C. Lopez, C. Neale, C. Montour, H. Ingolfsson, F. Di Natale, F. Lightstone and S. Gnanakaran, J. Chem. Theory Comput., 2018, 14, 6050–6062.
23. J. Fliege and B. F. Svaiter, Math. Methods Oper. Res., 2000, 51, 479–494.
24. V. Gapsys, B. L. de Groot and R. Briones, J. Comput.-Aided Mol. Des., 2013, 27, 845–858.
25. Y. Liu, A. de Vries, W. Pezeshkian and S. Marrink, J. Chem. Theory Comput., 2021, 17, 5876–5884.
26. G. M. Torrie and J. P. Valleau, J. Comput. Phys., 1977, 23, 187–199.
27. B. Roux, Comput. Phys. Commun., 1995, 91, 275–282.
28. J. S. Hub, B. L. De Groot and D. van der Spoel, J. Chem. Theory Comput., 2010, 6, 3713–3720.
29. V. Walter, C. Ruscher, O. Benzerara and F. Thalmann, J. Comput. Chem., 2021, 42, 930–943.
30. N. Michaud-Agrawal, E. J. Denning, T. B. Woolf and O. Beckstein, J. Comput. Chem., 2011, 32, 2319–2327.
31. Q.-Y. Zhou, J. Park and V. Koltun, arXiv, 2018, preprint, arXiv:1801.09847,  DOI:10.48550/arXiv.1801.09847.
32. L. Li, H. Davande, D. Bedrov and G. Smith, J. Phys. Chem. B, 2007, 111, 4067–4072.
33. D. Mohammadyani, H. Modarress, A. To and A. Amani, Braz. J. Phys., 2014, 44, 1–7.
34. L. Monticelli, J. Chem. Theory Comput., 2012, 8, 1370–1378.
35. J. Choi, S. Snow, J.-H. Kim and S. Jang, Environ. Sci. Technol., 2015, 49, 1529–1536.
36. J. Barnoud, G. Rossi and L. Monticelli, Phys. Rev. Lett., 2014, 112, 068102.
37. J. L. MacCallum and D. P. Tieleman, J. Am. Chem. Soc., 2006, 128, 125–130.
38. E. Lavagna, J. Barnoud, G. Rossi and L. Monticelli, Nanoscale, 2020, 12, 9452–9461.
39. J. Faraudo, J. Chem. Phys., 2002, 116, 5831–5841.
40. B. Fabian and M. Javanainen, J. Phys. Chem. Lett., 2024, 15, 3214–3220.
41. S. J. Marrink, J. Risselada and A. E. Mark, Chem. Phys. Lipids, 2005, 135, 223–244.
42. F. Harb, A. Simon and B. Tinland, Eur. Phys. J. E, 2013, 36, 1–7.
43. B. Fábián and M. Javanainen, J. Phys. Chem. Lett., 2024, 15, 3214–3220.
44. X. Woodward, M. Javanainen, B. Fabian and C. Kelly, Biophys. J., 2023, 122, 2203–2215.
45. S. Vanni, H. Hirose, H. Barelli, B. Antonny and R. Gautier, Nat. Commun., 2014, 5, 4916.
46. P. Diggins IV, Z. A. McDargh and M. Deserno, J. Am. Chem. Soc., 2015, 137, 12752–12755.
47. H. Kim, D. Bedrov and G. D. Smith, J. Chem. Theory Comput., 2008, 4, 335–340.
48. L. Li, D. Bedrov and G. D. Smith, J. Chem. Phys., 2005, 123, 204504.
49. A. Roosen, R. McCormack and W. Carter, Comput. Mater. Sci., 1998, 11, 16–26.
50. A. Singhal and G. A. Sevink, Nanomaterials, 2022, 12, 3859.
51. M. Deserno and T. Bickel, EPL, 2003, 62, 767.
52. H. Cui, E. Lyman and G. Voth, Biophys. J., 2011, 100, 1271–1279.
53. A. Singhal and G. A. Sevink, Nanoscale Adv., 2021, 3, 6635–6648.
54. E. Sarukhanyan, A. De Nicola, D. Roccatano, T. Kawakatsu and G. Milano, Chem. Phys. Lett., 2014, 595–596, 156–166.
55. Y. Cherniavskyi, C. Ramseyerb and S. Yesylevskyy, Phys. Chem. Chem. Phys., 2016, 18, 278.
56. J. Zupanc, D. Drobne, B. Drasler, J. Valant, A. Iglic, V. Kralj-Iglic, D. Makovec, M. Rappolt, B. Sartori and K. Kogej, Carbon, 2012, 50, 1170–1178.
57. N. Bouropoulos, O. Katsamenis, P. Cox, S. Norman, P. Kallinteri, M. Favretto, S. N. Yannopoulos, A. Bakandritsos and D. Fatouros, J. Phys. Chem. C, 2012, 116, 3867–3874.
58. Y. Zeng, Q. Wang, Q. Zhang and W. Jiang, RSC Adv., 2018, 8, 9841.
59. J. Zhang, X. Zhao and Q. H. Liu, RSC Adv., 2016, 6, 90388–90396.
60. E. Cino and D. Tieleman, Biophys. J., 2022, 121, 2060–2068.
61. C. Hofsäβ, E. Lindahl and O. Edholm, Biophys. J., 2003, 84, 2192–2206.
62. S. Sikdar, M. Banerjee and S. Vemparala, Soft Matter, 2021, 17, 7963.
63. Z. Qi, M. Wan, J. Zhang and Z. Li, J. Phys. Chem. B, 2023, 127, 1956–1964.
64. L. Borges-Araujo, A. C. Borges-Araujo, T. N. Ozturk, D. P. Ramirez-Echemendia, B. Fabian, T. S. Carpenter, S. Thallmair, D. Barnoud, H. I. Ingolfsson, G. Hummer, D. P. Tieleman, S. J. Marrink, P. C. T. Souza and M. N. Melo, J. Chem. Theory Comput., 2023, 18, 7387–7404.
65. J. Wong-Ekkabut, S. Baoukina, W. Triampo, I.-M. Tang, D. P. Tieleman and L. Monticelli, Nat. Nanotechnol., 2008, 3, 363–368.
66. D. Bedrov, G. D. Smith, H. Davande and L. Li, J. Phys. Chem. B, 2008, 112, 2078–2084.
67. R. Koynova and M. Caffrey, Biochim. Biophys. Acta, 1998, 1376, 91–145.
68. T. Kaasgaard, C. Leidy, J. Crowe, O. Mouritsen and K. Jørgensen, Biophys. J., 2003, 85, 350–360.
69. V. Walter, C. Ruscher, A. Gola, C. Marques, O. Benzerara and F. Thalmann, Biochim. Biophys. Acta, Biomembr., 2021, 1863, 183714.
70. P. Sharma, R. Desikan and K. Ayappa, J. Phys. Chem. B, 2021, 125, 6587–6599.
71. F. Harb, A. Simon and B. Tinland, Eur. Phys. J. E, 2013, 36, 140.
72. R. Koynova and M. Caffrey, Biochim. Biophys. Acta, Rev. Biomembr., 1998, 1376, 91–145.
73. P. Niemela, S. Ollila, M. Hyvönen, M. Karttunen and I. Vattulainen, PLoS Comput. Biol., 2007, 3, e34.
74. N. van Hilten, J. Methorst, N. Verwei and H. J. Risselada, Sci. Adv., 2023, 9, eade8839.

---