Mechanistic insights into interactions between ionizable lipid nanodroplets and biomembranes

Abstract

Delivery of RNA into cells using lipid nanoparticles (LNPs) has been a significant breakthrough in RNA-based medicine, with clinical applicability expanded through the use of ionizable lipids (ILs). These unique lipids can alter their charge state in response to pH changes, which is crucial for pH-triggered endosomal escape and effective lipid-mediated RNA delivery. In this study, we conducted a comprehensive set of molecular dynamics (MD) simulations to investigate interactions between IL-containing lipid nanodroplets (LNDs) and cell membrane models. Using an atomistic resolution model, we investigated the merging process of LNDs with cell membrane models under neutral conditions relevant to an intercellular environment and acidic pH conditions found in late endosomes. Our observations revealed that at neutral pH, LNDs merged with lipid membranes while preserving the bilayer structure. Under acidic conditions, the LNDs remained attached to the bilayer without fusing into the membranes. Importantly, the presence of ILs did not disrupt the original biomembrane structure during the simulation period. The MD simulations provided valuable atomistic insights into the mechanism of interaction between IL-containing nanodroplets and biomembranes, which could aid the rational design of ILs to develop more efficient LNPs for RNA therapies.

Communicated by Ramaswamy H. Sarma

Keywords:

• Ionizable lipids
• molecular dynamics
• membrane simulations
• lipid nanoparticles
• lipid mixing
• ionizable lipid nanoparticles

Introduction

RNA-based therapy is an emerging medical technology that has already proved its potential in vaccinations, as well as treatments for cancer and rare diseases (Hou et al., 2021; Jayaraman et al., 2012; Wang et al., 2020). Safe delivery of therapeutic RNA to its cellular target initially posed a significant barrier which was recently overcome by the encapsulation of RNA into lipid nanoparticles (LNPs) (Patel et al., 2019; Shin et al., 2018). This capability facilitated the distribution of over a billion doses of mRNA vaccines during the COVID-19 pandemic (COVID-19 Vaccine Doses Administered by Manufacturer, n.d.; Mathieu et al., 2021). However, despite significant progress in this field, many important questions remain unanswered. For example, fine details of the mechanism of cargo delivery and the nature of interactions of LNP lipids with cellular compartments, including biomembranes, have not yet been explored.

Use of ionizable lipids (ILs) represents a crucial breakthrough (Jayaraman et al., 2012; Semple et al., 2010) that has increased pharmaceutical applications of RNA-carrying LNPs. Initially, cationic lipids were considered due to their RNA encapsulation capability, but their practical use was hindered by their cytotoxicity (Kumari et al., 2020). Subsequently, they were replaced by ILs because the pH-dependent charge of the IL molecules decreases the LNP cytotoxicity (Cheng et al., 2022) whilst preserving the RNA encapsulation capacity. IL-based LNPs are assembled under acidic conditions, which favour electrostatic interactions between the positively charged (protonated) ILs and negatively charged RNA (Eygeris et al., 2022). When present in an intercellular medium, the system is neutralized and LNP ILs become uncharged (deprotonated). When uncharged LNPs encounter a cell membrane, they are encapsulated into endosomes (Patel et al., 2019). As endosome maturation progresses, the pH decreases and the ILs become reprotonated (Eygeris et al., 2022). The apparent pK_a of a LNP mixture is a tuneable property, with the peak RNA delivery efficiency usually obtained at pK_a 6-7 (Carrasco et al., 2021), resulting in an efficient response of LNPs to the pH change during endosome maturation and release of RNA into the cytosol (Eygeris et al., 2022).

Gaining atomistic insights into LNP structure and behaviour inside the human body is an important albeit very challenging task. Whereas experimental techniques have provided valuable information about the RNA delivery efficiency, assessment of the inner LNP structure remains an experimentally unexplored field and structural details of the interaction of LNPs with cellular membranes have yet to be investigated. So far, experimental techniques have provided a reliable description of the LNP size, which ranges from 50 to 300 nm in diameter, depending on the methods used for LNP preparation and size measurement (Malburet et al., 2022; Shepherd et al., 2023). Further, structural insights into LNP organization have been obtained from in vitro cryo-TEM experiments of either onion- (Kulkarni et al., 2018) or blebbed-phase (Brader et al., 2021) LNPs. Complementary to experiments, molecular dynamics (MD) simulations can provide unique insights into structural and dynamic features of very complex systems on a femtosecond timescale and with atomic resolution. MD simulations have been used to evaluate the lipid phases in COVID-19 vaccine mixtures (Paloncýová et al., 2021), structural organization of LNPs (Bruininks et al., 2020; Leung et al., 2012) and pH-dependent behaviour of LNPs (Paloncýová et al., 2023; Trollmann & Böckmann, 2022). Recently, simulations of lipid bilayers composed of increasing fractions of deprotonated ILs revealed formation of a separated IL phase at high IL concentration (Dehghani-Ghahnaviyeh et al., 2023). Whereas large models of LNPs with radii of dozens of nm focus on the overall arrangement, small systems of single nanometer scales can reveal the behaviour and function of individual molecules.

Small lipid bilayers can serve as good models for MD simulations of biomembranes and lipid-based systems, offering insights into the structure, mechanism and behaviour of these essential natural components. Precise representation of vital properties has been reported, such as interactions with small molecules (Paloncýová et al., 2014), peptides and proteins (Sandoval-Perez et al., 2017) and nanomaterials (Chen et al., 2016; Lazaratos et al., 2020). Moreover, previous studies of non-lamellar lipid phases (Marrink & Mark, 2004; Marrink & Tieleman, 2002; Ramezanpour et al., 2020) have provided representative lipid models with adequate and justifiable information even on limited size systems. Hence, lipid bilayer models can be used to unravel molecular mechanisms underlying critical biological processes that cannot be determined by traditional experimental techniques. In the present study, we used a lipid bilayer model to evaluate interactions of cell membrane and LNPs mimicking lipid mixtures in COVID-19 vaccines.

We carried out an extensive set of MD simulations to investigate interactions between LNP-mimicking lipid mixtures and mammalian cell membranes. As cell membranes are complex environments containing a variety of lipid species (Doktorova et al., 2020; Harayama & Riezman, 2018; Lorent et al., 2020), we explored a series of models of various phospholipid bilayers differing in lipid headgroup. To evaluate their interactions with ILs currently used in LNPs, we simulated the lipid bilayer merging with an uncharged or charged lipid nanodroplet (LND) mimicking the LNP lipid composition. We focused on describing the merging process and the effect of ILs and helper lipids (DSPC, CHL) in the LNP composition on the phospholipid bilayer structure and stability. We observed that neutral nanodroplets merged completely with membrane lipids and did not destabilize the membrane structure, whereas with the acidic form of ILs, the LND stayed connected to the lipid membrane but was not fully incorporated into it. These findings provide unique insights into the fate of ILs in an organism, which could aid further structural and functional optimalization of LNPs.

Methods

Composition of the systems

POPC (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine), POPE (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoethanolamine), POPG (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoglycerol) and POPS (1-palmitoyl-2-oleoyl-sn-glycero-3-phospho-L-serine) were used to build the membranes. The IL ALC-0315, [(4-hydroxybutyl)azanediyl]di(hexane-6,1-diyl) bis(2-hexyldecanoate), denoted A3D in the deprotonated form and A3P in the protonated form (Figure S1), and the helper lipids cholesterol (CHL) and 1,2-distearoyl-sn-glycero-3-phosphocholine (DSPC) were used to build the LNDs. Membranes of regular size were modelled as bilayers with 100 lipids in each leaflet, and LNDs comprised 100 lipids with a 50:40:10 molar ratio of IL:CHL:DSPC. To evaluate the effect of the size of the system, POPC lipid bilayers were also prepared with 200 or 300 POPC molecules per leaflet. The systems were neutralised and adjusted to the physiological concentration of 0.15 M by adding K⁺ and Cl⁻ ions.

Simulation protocol

The lipid membranes were assembled, minimised and subjected to a 200 ns equilibration run. The LNDs were assembled, minimised and subjected to a 100 ns equilibration run. Afterwards, both systems were combined and hydrated (see run details below, Figure 1), then an unbiased “free” 500 ns simulation of the LND-membrane system was performed. As the unbiased simulation did not lead to a spontaneous merging, the same inputs were then subjected to a 100 ns pulling simulation, during which the centre of mass (CoM) of the LND was pulled towards the CoM of the membrane along the z-axis. The systems were then subjected to a 1000 ns (or 1500 ns in the case of a POPE membrane with an uncharged LND) unbiased production run (see Figure 1 for a schematic of the process). Each of the pulling + production simulations were performed in three replicas.

Figure 1. Scheme of the simulation protocol. Membrane headgroups are depicted in orange (PC), light orange (PE), purple (PG) and green (PS), all lipid tails in grey, CHL in yellow, ILs in pink (deprotonated, neutral) and cyan (protonated, positively charged).

Software and force fields

All simulations were performed in the GROMACS 2018.1 software (Abraham et al., 2015; Berendsen et al., 1995; Hess et al., 2008; Lindahl et al., 2001; Páll et al., 2015; Pronk et al., 2013; Van Der Spoel et al., 2005). The forcefields used were Amber Lipid17 for the standard lipids (Gould et al., 2018), OPC (Izadi et al., 2014) for the water model and Joung − Cheatham (Joung & Cheatham, 2008) for the ions. The ILs used were parametrised in our previous work (Paloncýová et al., 2021) by combining parameters of Lipid17 (Gould et al., 2018) for the tails and GAFF (J. Wang et al., 2004) for the ionizable headgroup.

The pure membrane systems were assembled using the CHARMM-GUI (Jo et al., 2008; Lee et al., 2016) webpage (Jo et al., 2007, 2009; Lee et al., 2019, 2020; Wu et al., 2014), whereas the IL containing systems were assembled using an in-house bilayer builder, originally created by Gabin Fabre for (Paloncýová et al., 2014). The software used for analyses were as follows: GROMACS 2018.1 (for order parameter calculations), MDAnalysis 2.3.0 (Gowers et al., 2016; Michaud-Agrawal et al., 2011) (and its dependencies SciPy (Virtanen et al., 2020), NumPy (Harris et al., 2020) and Biopython (Chapman & Chang, 2000; Cock et al., 2009; Hamelryck & Manderick, 2003)) for density based and system-property based analyses, LiPyphilic 0.10.0 (Smith & Lorenz, 2021) (with freud (Ramasubramani et al., 2020)) for leaflet-based analyses and gromacs_LS v2016.3 (Torres-Sánchez et al., 2015, 2016; Vanegas et al., 2014; Vanegas & Arroyo, 2014) (with LAPACK (Anderson et al., 1999)) for lateral pressure calculations. Graphs were prepared using Matplotlib (Hunter, 2007). The Python-based analysis tools were run in Python 3.9.12 as standalone scripts or within the Jupyter Notebook (Beg et al., 2021). The visualisation software used were PyMol 2.5.5 (Schrödinger & LLC, 2015) and VMD 1.9.3 (Humphrey et al., 1996).

Simulation details

The pure membrane systems were subjected to steepest descent minimization and multi-step of ∼3 ns subsequent equilibration runs set at 310 K, in accordance with the procedure recommended by CHARMM-GUI. The systems were then subjected to a 200 ns production run.

The lipids used to build the LNDs were selected to mimic the composition of the Pfizer&BioNTech mRNA Comirnaty vaccine. Therefore, a bilayer of 50 ILs (ALC-0315), 40 CHLs and 10 DSPC molecules was created using our in-house script. We avoided using the PEGylated lipid ALC-0159 because its low content required including only two molecules per system, which would likely have generated simulation artefacts. The resulting bilayer was then subjected to the CHARMM-GUI equilibration process. Next, the LND was inserted into a larger box with at least 2 nm of solvent on all sides (corresponding to approximately 12 × 12 × 7 nm box size for the uncharged, and 10 × 10 × 8 nm for the charged system), re-equilibrated and simulated for 100 ns under isotropic conditions.

Each membrane and LND pair (after their respective 200 ns and 100 ns of simulation, Figure 1) were manually combined so that the LND lay above the membrane. The system was rehydrated and subjected to the same minimization and then a 50 ps simulation to settle the box. The systems were simulated for 500 ns without any bias. As the LNDs and membranes remained separate in these free simulations, the starting structures were subjected to 100 ns of pulling and 1000 ns of unrestrained simulation (350 ns for large POPC models, see the supplementary spreadsheet). The pulling was performed by applying a harmonic biasing force to decrease the distance along the z-axis by 0.05 nm per ns with a force constant of 10,000 kJ∙mol⁻¹∙nm⁻². As the lipid dynamics of the POPE systems with deprotonated IL was still ongoing 1000 ns post-pulling, this one system was prolonged to 1500 ns of free run. The simulations were performed with three replicas for each system for the regular size (100 lipids per leaflet) and two replicas for the larger POPC membrane models.

The simulation step was 2 fs and periodic boundary conditions were applied in all directions. The thermostat was set to 310 K according to Berendsen (Berendsen et al., 1984) for the equilibration and Nosé-Hoover (Hoover, 1985; Nosé & Klein, 1983) for the production runs with a temperature coupling constant of 1.0 ps. The pressure coupling was semi-isotropic (Nosé & Klein, 1983; Parrinello & Rahman, 1981) with a pressure coupling constant of 5.0 ps, pressure of 1 bar and compressibility of 4.5 × 10⁻⁵bar⁻¹, unless stated otherwise (see the supplementary Systems and Simulation Details spreadsheet). LINCS (Hess, 2008) restraints were applied to all H-bonds, and long-range electrostatics were calculated via the particle-mesh Ewald method (Essmann et al., 1995). All the simulation .mdp files are available in the supplementary dataset, and the full list of simulated systems, their abbreviations, amount of water/ion molecules and other simulation details can be found in the supplementary spreadsheet. Furthermore, all simulation protocols, starting/final structures, topologies, trajectories and scripts are available in the online Zenodo repository (doi: 10.5281/zenodo.10679359).

In this paper, the membrane side from which the LND was pulled (and to which the spike remained connected) will be called upper and the opposite one lower. The main parts of the systems are depicted in Figure 2D and E. Division of the system into four sections (above and below the membrane, lower and upper leaflet) was defined based on most populated position of the atoms of the membrane headgroups, averaged in 50 ns blocks (Figure 2D). The headgroups and tail groups used were defined by the Amber Lipid17 forcefield. Peaks were located using the find_peak function of the SciPy (version 1.9.3) Python package with height parameter set to 0.19 (g∙cm⁻³) and distance parameter set to 55 (histogram bins). Because the centre of membrane moved in the z-axis during the simulation, the centre of the membrane was defined as the average of the position of the two peaks, dynamically in 50 ns blocks. The lipids were assigned a position in the respective section based on the z-axis position of the CoM of each residue. As the ILs remained in the upper curve of the spike until its dissolution, the fraction of ILs above the membrane was used to describe the mixing process.

Figure 2. Outer and inner views of (A) deprotonated and (B) protonated LNDs (ILs in pink (deprotonated) and cyan (protonated), CHL in yellow, helper DSPC lipid headgroup in orange, tails in grey; (C) Structure of IL ALC-0315 (A3D) represented in uncharged (deprotonated) form; (D) Density profile of a starting structure of a deprotonated IL LND and POPC membrane, with mass densities of the headgroup shown in orange, tail in grey, CHL in yellow, IL in pink and water in light blue. The system was divided into four sections (above and below membrane, upper and lower leaflet) based on the position of the headgroup density peak; (E) Representation of parts of the system: LND (containing charged ILs) interacting with upper leaflet of a membrane forming a lipid spike connected to the membrane via the interface.

Lateral pressure profile calculations were performed on selected snapshots of the simulation by running a separate 1 ns simulation with 1 ps sampling and 2 nm cut-off treatment of electrostatics. 10,000 frames of this run were selected for analysis in 100 frame blocks, as recommended by the gromacs_LS manual. The pressure profiles were processed with a Gaussian filter with standard deviation of 0.2 nm and visualized in the z-direction.

For calculation of free energy profiles of the flip-flop process, POPC membranes were prepared with 36 lipids per leaflet. The membrane was equilibrated for 100 ns, then a single molecule of A3P, A3D or CHL was added at the top of the simulation box and another 100 ns free simulation was performed. From the position of the added molecule closest to the membrane centre, it was pulled at a rate of 0.01 nm∙ns⁻¹ until it reached the centre of the membrane. Both unbiased and biased trajectories were used to extract starting frames for umbrella simulations separated by 0.1 nm, with an umbrella potential of 2000 kJ∙mol⁻¹∙nm⁻¹ between the CoM of the membrane and a lipid head group (nitrogen atom in A3P and A3D and oxygen in CHL). Free energy profiles were extracted by the weighted histogram analysis method (Hub et al., 2010; Kumar et al., 1992) from 100 ns umbrella simulations, taking the first 50 ns as an equilibration period.

Results

pH determined the shape of LNDs

To describe the behaviour of the IL-containing mixture in both the neutral intercellular medium and acidic endosomal environment, LND models were built for uncharged and charged forms of the IL (denoted A3D and A3P, respectively). The pK_a of the IL was calculated as 10.0 (Eygeris et al., 2022), but the deprotonated form of the LND was considered as a model of the lipid mixture behaviour under neutral conditions (pH ∼7) as the apparent pK_a of the whole LNP for an experimentally relevant scale was lower, that is, 6.1 (P. Patel et al., 2021). MD simulations (100 ns) of the lipid droplets revealed that they remained stable with diameters of ∼6 nm and ∼7.5 nm for deprotonated and protonated IL forms, respectively. The shape of the LNDs depended on the IL protonation state (Figure 2). In the neutral form of the IL, the LND formed a sphere with IL distributed across the whole volume. In the charged form of the IL, relevant to the acidic endosomal environment, the LND was shaped like a nanodisc, with a CHL-DSPC-rich bilayer-like centre and the IL clustered close to the nanodisc edges, in accord with our previous observations (Paloncýová et al., 2021). In both cases, the DSPC molecules preferentially associated with CHL (Figure 2) and CHL molecules diffused between locally stable clusters via the disordered IL-rich phase.

Merging of LNDs with lipid bilayers required an external force

In unbiased 500-ns-long MD simulations, the LNDs did not spontaneously merge with the lipid bilayers (Figure S5) due to a kinetic barrier. Such a barrier preventing spontaneous lipid membrane merging has been described theoretically and observed in simulation (Bruininks et al., 2020; Leikin et al., 1987; Schneck et al., 2012), even though in large-scale experiments, merging would likely occur. However, the kinetic barrier can be overcome, for example, by dehydration of a contact region (Bruininks et al., 2020). In the present study, the CoM of the LND was pulled towards the CoM of the membrane for 100 ns. The pulling induced membrane deformation, which spontaneously relaxed within the first 50 ns of the following unbiased simulation (Figure 3). After release of the biasing force, the LND slid partially out of the membrane due to a high lateral pressure (Figure 3), forming a spike-like structure. In all cases, extensive lipid exchange between the spike and membrane occurred, as described later.

Figure 3. Time evolution of lateral pressure profiles of an uncharged (top) and charged (bottom) system with POPC lipid and respective simulation snapshots. Negative values in the pressure profiles correspond to attractive forces, whereas positive peaks correspond to repulsive regions. The time values indicate the time from the beginning of the pulling simulation. Therefore, 100 ns corresponds to the end of the pulling simulation, after which an unbiased simulation was run.

LNDs with neutral ILs dissolved in the lipid membrane

In the case of the uncharged form of the IL, the LND fully merged with the membrane’s bilayer lipids and the spike dissolved into the membrane. However, the dynamics of the process differed between the POPE bilayer and other systems. In the case of bilayers composed of POPC, POPG and POPS lipids, the spike merged with the membrane rapidly (Figures 3, S4, and S6a,b). In the three simulations with deprotonated IL, the LND had dissolved in the bilayers by 350, 400 and 300 ns (median values for the three replicas) after pulling for the POPC, POPG and POPS membranes, respectively. The LND dissolution was accompanied by rapid flip-flop of the IL molecules to the lower membrane leaflet (Figures 4 and 5). Similar results were observed for the larger POPC membranes containing 200 or 300 POPC molecules per leaflet (Figure S7). In POPE, the IL separated into a disordered phase residing in the middle of the membrane, resulting in inverse curvature of the upper leaflet, as documented before for PE lipids (Kaasgaard & Drummond, 2006) (Figure S6d). This was due to the smaller volume of the POPE lipid headgroup making a conical molecular shape that favours an inversely curved surface (Ding et al., 2015; Kaasgaard & Drummond, 2006). After 1300 ns of post-pulling (median value, Table S2), the IL content above the POPE membrane had decreased to under 10%, together with rapid accumulation of the IL in the lower leaflet, analogously to other simulations.

Figure 4. Fraction of IL content above the membrane and in the lower leaflet during simulations averaged over three replicas. The solid lines correspond to systems with deprotonated IL, whereas the dotted lines correspond to the protonated IL. The grey vertical line indicates the end of pulling.

Figure 5. Ratio of different lipids above the membrane (upper panel) and within the two membrane leaflets (middle and lower panel) averaged over three replicas. ML (grey) denotes the membrane lipid content, CHL (yellow), A3D (pink) the deprotonated and A3P (cyan) the protonated form of the IL.

While the uncharged IL possess a very small, low polar head group allowing its facile flip-flop between leaflets, no helper lipids from the LND flipped into the opposite membrane leaflet (Figure 5). A rare CHL flip-flop was observed with deprotonated ILs in at least one replica per system (Table S1). The flip-flop was facilitated by a highly disordered IL-rich phase created by IL accumulation within the membrane. To evaluate this phenomenon thermodynamically, we created a series of simulations with POPC membrane and a single lipid molecule (CHL, A3D, A3P) pulled inside the membrane. The extracted snapshots were used to calculate a free energy barrier for the flip-flop of ILs and CHL in a pure POPC membrane. While CHL had a 5.8 kcal/mol free energy barrier, A3D had a barrier of only 4.4 kcal/mol (Figure S8) and a further barrier decrease was expected in a disordered A3D-rich phase. Our simulations showed that the flip-flop of uncharged IL was more favourable than that of CHL and the system used it to alleviate the crowded upper membrane (Figure 3).

The flip-flop mechanism of the uncharged IL makes advantage of the high flexibility of the molecule. In the first step (Figure 6A), the IL head stretched towards the bulk water, while the deprotonated nitrogen was located under the main lipid phosphate (1.1 nm from the middle of the membrane, Figure S3). Next, the whole molecule moved to the area between the leaflets (Figure 6A1), where it remained for up to tens of nanoseconds. Afterwards, the tails oriented to the other leaflet (Figure 6A2), with the head group still in the membrane core. The membrane core was enriched by polar oxygens located in the middle of IL tails, as also seen in the small symmetric membrane prepared separately (Figure S2). Finally, the headgroup flopped so that it reached the bulk solution on the opposing side of the membrane, with the deprotonated nitrogen remaining under the main lipid phosphates (Figure 6A3).

Figure 6. Mechanism of diffusion of (A) uncharged IL (A3D) from the upper to lower leaflet, (B) charged IL (A3P) from an LND-induced spike to the membrane, (C) POPC from the membrane to the spike. The lipid of interest is shown by sticks, other molecules as balls. POPC headgroups are shown in orange, POPC tails in grey, charged IL in cyan/green, uncharged IL in light/bright pink, nitrogens of the lipid of interest in blue and oxygens of the lipid of interest in red. Numbers in the figures of the left column correspond to the position of ionizable nitrogen/POPC phosphorus atoms of the steps on the right.

Due to the different abilities of the IL and CHL to flip-flop, the effect of the LND lipids on the membrane itself was leaflet-dependent. In the deprotonated systems, where the entire LND dissolved in the membrane, the upper leaflet took up (almost) all the CHL molecules (Figure 5). This increased the deuterium order parameters of the membrane palmitoyl tails (Table S2) to the level of a corresponding POPX:CHL (65:35%) membrane (Table S3). Correspondingly, the area per lipid of the upper leaflet decreased by 15–20% (Tables S2 and S3). In contrast, the order parameters of the lower leaflet corresponded to a pure POPX membrane. The area per lipid of the lower leaflet was not significantly different from that of the pure membrane in the POPC system and was slightly higher in the other systems. Other leaflet properties, such as membrane thickness were not significantly affected by the presence of ILs (Table S2).

Mixing of deprotonated ILs with cell membrane phospholipids resulted in a bilayer that was consistent with current knowledge of IL behaviour. Previously, IL phase separation between membrane leaflets was observed for two-chained ILs DLin-KC2-DMA (Ramezanpour et al., 2019) and DLin-MC3-DMA (Ermilova & Swenson, 2020)). However, our systems contained branched ILs with only ∼16 mol% of IL in the final mixture, more akin to the very recent findings of Dehghani-Ghahnaviyeh et al. (Dehghani-Ghahnaviyeh et al., 2023). Under low concentrations (<20% of the IL), Lipid-5 (with one tail less than ALC-0315) was found to be localised both on the membrane surface and deeper in the tail and formed a separated phase only at high IL concentrations.

LND with charged ILs formed a spike connected to the lipid membrane

In the systems with the charged IL, relevant to the acidic environment in late endosomes, the lipid spike remained dynamically stable in all systems, with some differences in the inner spike conformations between systems with different headgroups. In the PC and PS systems, the spike retained a globular shape on the membrane, whereas in the PE and PG systems with smaller area per lipid, the spike connected across the periodic boundary conditions, forming a tube connected to the membrane via a highly disordered IL-rich phase (Figure S6). The conformation of the system allowed the dynamic exchange of lipids between the membrane and the spike. At least 10% of the IL and 20% of CHL diffused to the upper layer of the membrane, and oppositely at least 10% of the upper-layer membrane lipid diffused into the spike (Figure 5). However, no flip-flop of the charged IL was observed in any of the systems, in contrast to the behaviour under neutral conditions. The flip-flop free energy barrier was estimated to be as high as 15.2 kcal/mol. Due to the absence of the IL flip-flop, the asymmetric lateral pressure could not be equilibrated between the lipid leaflets and equilibrium was reached via the presence of a lipid spike (Figure 3). Even though the spike induced local deformation of the membrane surface, the overall structure remained intact.

The differences in the lipid mixing process in the neutral system stemmed from the increased hydrophilicity of charged nitrogen ILs. Concerning the mechanism of lipid mixing, two simultaneous processes were observed in the interfacial region: spontaneous translocation of ILs from the lipid spike to the bilayer and diffusion of membrane lipids in the opposite direction. Lipid translocation occurred after diffusion of the ionizable and membrane lipids towards the interfacial region (Figure 2E). The lipids adopted an extended conformation where their tails spread between the lipid bilayer and lipid spike at their interface (Figure 6B1 and C2), connecting both parts of the system. Subsequently, the remaining lipid tails from both lipids followed the first one (Figure 6B3 and C3), completing the translocation process. As a result, the protonated nitrogen of the IL resided on the level of the membrane lipid phosphates within the bilayer, whereas membrane lipid accumulated alongside CHL, forming a partially ordered conformation in the lipid spike.

The MD simulations provided an insight into the dynamics of interactions between IL-containing mixtures and lipid biomembrane models, along with an explanation of the mixing process. Furthermore, they revealed that flip-flop of IL molecules may underlie the mixing during their mutual interaction. Our models of LNDs with neutral ILs demonstrated that the LND lipids rapidly integrated with cell membranes, ultimately leading to complete dissolution of the LND into the lipid membrane. In contrast, the presence of charged ILs in the LND led to weaker interaction with the biomembrane facilitated only by the lipid spike on the LND-membrane interface.

Conclusions

We conducted MD simulations to investigate the mechanism of interaction between IL-containing lipid nanodroplets (LNDs) and mammalian membrane lipids. To explore individual stages of lipid-mediated delivery systems, we studied ILs in both their protonation states. In the deprotonated (uncharged) state, which is relevant under neutral pH conditions in the bloodstream, LNDs merged with and dissolved seamlessly into the membrane. The uncharged ILs readily flip-flopped across the membrane, indicating smooth integration with mammalian cell lipids. These findings indicate the high biocompatibility of neutral ILs, causing minimal disruption to the membrane structure. Conversely, in the protonated state, charged IL-containing nanodroplets did not merge into the membrane but instead induced membrane deformation in the form of a lipid spike. This phenomenon may be attributed to the inability of protonated ILs to flip-flop, affecting predominantly only one membrane leaflet. These atomistic insights into IL-membrane interactions offer great potential for advancing the research and rational design of ILs in the rapidly progressing field of RNA medicine. Furthermore, the in-depth understanding presented here opens new avenues for future exploration in IL-biomembrane simulations at larger scales, which could ultimately yield a comprehensive understanding of endosomal escape mechanisms.

Acknowledgements

We acknowledge IT4Innovation, Czech Republic for granting this project access to the LUMI supercomputer owned by the EuroHPC Joint Undertaking, hosted by CSC (Finland) and the LUMI consortium with the support of the Ministry of Education, Youth and Sports of the Czech Republic through e-INFRA CZ (grant 90254). COST Action CA21101 is also acknowledged.

Disclosure statement

MO has a share in a biosimulations company, InSiliBio, which focuses mostly on membrane permeation.

Data availability statement

The starting structures and topologies as well as trajectories, simulation details and scripts are available in the online Zenodo repository (https://doi.org/10.5281/zenodo.10679359).

Additional information

Funding

The work was supported by the ERDF/ESF project TECHSCALE (No. CZ.02.01.01/00/22_008/0004587) and by the European Union’s Horizon Europe research and innovation programme under grant agreement No 101070865 (MINIGRAPH).