Driving forces in the assembly of lipid nanoparticles containing mRNA revealed by molecular dynamics simulations at acidic and physiological pH

To simulate the experimental preparation conditions of LNPs, MD simulations were performed under two different external environmental pH, mimicking the pH variation observed in wet lab setups⁶. Initially, the LNPs were generated in a citrate buffer solution at pH 4.5, a condition that maintains the positively charged state of SM-102P, facilitating the assembly of the lipid components into nanoparticles. This configuration is denoted as the Positive system. Subsequently, PBS was introduced to adjust the pH to 7.4, inducing a charge transition in SM-102: SM-102 molecules on the outer surface of LNPs become neutral (SM-102N), while those in the inner region of the LNP, with a pH of 4.5, remains positively charged (SM-102P). This second LNP configuration is denoted as the Neutral system. To investigate the effect of the SM-102 protonation state on mRNA encapsulation, we also created a third system in which all SM-102 molecules are in the neutral state (SM-102N), referred to as the All-Neutral system. Note that in the following sections of this paper, SM-102 will refer to the molecule regardless of its charge state.

The preparation of the simulation system was carried out using Packmol¹⁴ with a box dimension of 160 × 160 × 160 ų. The number of molecules in the system was determined based on the ratio of protonated amine groups (N) of SM-102 at pH 4.5 to the number of negative charges from the phosphate groups (P) of the mRNA⁷. We used an N/P of 30:2 where one molecule of mRNA consisting of 21 nucleobases, namely 5'-UUCGUUGUCAAUGACGCUGCA-3', which has a total negative charge of − 20 on the phosphate groups. The number of SM-102 molecules required to achieve an N/P ratio of 30:2 is 300. The ratio of SM-102:Cholesterol: DSPC: DMG-PEG2000 used was 50:38.5:10:1.5, which was found to be effective based on literature^6,15,16. In this study, the LNP consists of 300 SM-102, 231 cholesterol, 60 DSPC, 9 DMG-PEG2000, and 1 mRNA. The number of SM-102P and SM-102N in the Neutral system is 62 and 238, respectively, obtained by adjusting the protonation state of SM-102P at a distance greater than 15 Å from mRNA. Table 1 presents the detailed composition of the simulation systems.

MD simulations were performed using Amber 22¹⁷. The force field used for mRNA was OL3¹⁸. LIPID21¹⁹ was used for cholesterol and DSPC. For H₂O and Na⁺, TIP3P²⁰ was applied. SM-102, DMG-PEG2000, and citrate were parameterized using GAFF2 ²¹. We performed RESP charge parameterization for SM-102 and DMG-PEG 2000 at the ESP-HF/6-31G(d)//MP2/aug-cc-pVDZ level of theory, according to the structure fragments in Figure S1. Structure optimization was conducted using Gamess US VERSION 30 SEP 2020 (R2)²² while ESP charge calculations were done with Gaussian 09²³. The parameters for citrate were obtained from our previous work²⁴.

Minimization was performed in four stages with a restraint of 10 kcal mol^{– 1} Å^{– 2}. The first and second stages of minimization used the steepest descent method for 10,000 cycles each. The first stage minimized mRNA, lipids, DMG-PEG2000, and ions, while the second stage was for H₂O. The third and fourth stages of minimization used MD simulations in the NVE ensemble for 1 ns each. The third stage was for mRNA, lipids, DMG-PEG2000, and ions, whereas the fourth stage was for H₂O.

System volume equilibration at 300 K was performed for 3 ns, with the last 2 ns of equilibration using the SHAKE algorithm for bonds involving hydrogen. System density equilibration was conducted for 10 ns with pressure regulation using the Berendsen barostat. During both equilibration stages, all system components except H₂O and ions were restrained with a force constant of 0.1 kcal mol^–1 Ų.

Production simulation was conducted for 1000 ns. The simulation was performed in triplicate from the third minimization stage to the production stage. Each replicate was run with a different random seed, potentially leading to variations across trajectories of the same system. Non-bonded interaction calculations were performed using the PME method with a 10 Å cutoff. Temperature regulation was carried out using the Berendsen thermostat, while pressure regulation used the Monte Carlo barostat. The analysis of the total 9000 ns simulation trajectory was performed using CPPTRAJ²⁵ particularly on the last 500-ns of MD trajectories unless otherwise specified. Interaction energy analyses were carried out through molecular mechanics with generalized Born and surface area (MMGBSA) calculations using MMPBSA.py²⁶.

We utilized Biovia Discovery Studio Visualizer v21.10.020298 (Dassault Systèmes, San Diego, USA) and Visual Molecular Dynamics (VMD) 1.9.4²⁷ for figure generation. The R packages of ggridge²⁸ ggplot2²⁹, gridExtra³⁰ tidyr³¹ ggforce³² and ggpubr³³ were used to make graphs on Jupyterlab 3.6.7 (Project Jupyter, Berkeley, CA, USA)³⁴ under an R programming language environment version 4.3.1 (R Foundation for Statistical Computing, Vienna, Austria)³⁵. Artworks were created using Inkscape 1.3 (0e150ed, 2023-07-21) (The Inkscape Project, Boston, MA, USA)³⁵.

In the first step, we analyzed the radius of gyration, computed as a scalar quantity representing the overall atomic distribution. To ensure reproducibility, three independent replicates were conducted for each system, yielding consistent results (Fig. 1A). Figure 1A illustrates the evolution of the radius of gyration for both individual lipid components and 'All lipids' (see the caption of Fig. 1) in two LNP systems: Positive and Neutral. In the Positive system, both individual lipid components and 'All lipids' display an initial rapid decrease in the radius of gyration (Fig. 1A), indicating a shift toward a more compact and ordered arrangement. This behavior suggests that, after the initial structural instability, the LNP components rearrange into a denser configuration. By contrast, the Neutral system exhibits a consistently tighter configuration for 'All lipids' throughout the simulation (44.09 ± 0.27 Å) compared to the Positive system (45.84 ± 0.93 Å), indicating reduced fluctuations and greater structural stability. Nonetheless, unlike other lipid components, cholesterol in the Neutral system exhibits an increased radius of gyration across all replicates, suggesting a potential structural perturbation that warrants further investigation to elucidate the underlying mechanisms.

An interesting aspect is that the gyration radius profiles for 'All lipids' closely mirror those of SM-102 in both LNP systems. This similarity suggests that the structural behavior of SM-102 plays a significant role in determining the overall compactness and stability of the LNP systems.

Since the radius of gyration analysis indicates the more compact LNPs in both Positive and Neutral systems throughout the MD simulations, we proceeded with a mass density analysis of the 'All Lipids' and SM-102 lipid components over various time ranges (Figure S2, Figure S3, Figure S4, and Figure S5). Mass density profiles were calculated by combining trajectories from all three replicates for each system to enhance statistical robustness in resolving collective structural organization and improve clarity in visualization. As shown in Fig. 1B (upper panel), the 'All Lipids' in the Positive system during the late simulation period (grey line) generally exhibit shorter distances from the mRNA mass center compared to the early period (0–10 ns, orange line) across all axes, suggesting a decrease in LNP size. Interestingly, the SM-102 lipid component shows similar mass density profiles (Fig. 1B middle panel, Figure S4) to those of 'All Lipids', supporting our previous notion that SM-102 significantly influences the overall compactness and stability of the LNP. Additionally, the mass density on the x-axis of the PEG2000 moiety likely shifts toward the center of the LNP (Fig. 1B lower panel, grey line).

Mass density profiles of 'All lipids' in the Neutral system may experiences gradual decreases in distance from the LNP center of mass over time, particularly along the x and z axes (Figure S3). Consequently, the mass density profiles in the late time range (last 500 ns, grey line) show smaller distances compared to those in the early time range (0–10 ns, red line) (Fig. 1C upper panel, the x and z axes). The SM-102 lipid component also follows a similar trend of mass density decreases like those of 'All Lipids' (Fig. 1C middle panel, Figure S5), which further supports our previous notion that SM-102 significantly influences the overall compactness and stability of the LNP. Moreover, the mass density of PEG2000 is more distributed toward the LNP (Fig. 1C lower panel).

The solvent accessible surface area (SASA) analysis may provide complementary insight into the compactness of LNPs. A lower SASA corresponds to reduced solvent accessibility and a more compact structure, whereas higher values reflect greater surface exposure³⁷. In the Positive system, the SASA profiles for 'All Lipids' (Fig. 1D, left panel) decrease sharply within the first 10 ns across all replicates, demonstrating rapid reorganization into a tightly packed configuration. This compaction aligns with the gyration radius trends (Fig. 1A). Like those in gyration radius analysis, SM-102 (Fig. 1D, left panel) have a similar SASA profiles to 'All Lipids'.

The 'All Lipids' in the Neutral system (Fig. 1D right panel) shows a gradual decrease of SASA values, indicating the LNP becoming more ordered, which agrees with our gyration radius analyses (Fig. 1A). SM-102 (Fig. 1D right panel), which combines positively (SM-102P) and neutrally charged SM-102 (SM-102N) molecules, similar SASA profiles to 'All Lipids'. Further analysis on R03, the moiety of SM-102P bearing the positive charge, reveals increased SASA values in the first 10-ns, followed by fluctuated SASA values with a steady trend (Figure S6). Meanwhile, 03N, which is the R03 with a neutral charge in SM-102N, exhibits similar SASA profiles to SM-102 and 'All Lipids'. Like the gyration radius analysis, the SASA of 'All Lipids' in the Neutral system (106,847.26 ± 3,953.77 Ų) indicates a more compact LNP structure compared to the Positive system (110,485.45 ± 6,396.65 Ų).

We hypothesized that SM-102 molecules drive LNP compactness and stability via electrostatic interactions between the negatively charged mRNA and the protonated R03 headgroup of SM-102 at pH 4.5. Therefore, we calculated radial distribution functions (RDFs) between the center of mass of mRNA (Figure S7) and the R03 headgroup in the Positive system using combined trajectories from all replicates to enhance statistical resolution. As illustrated in Fig. 2A (upper panel), the RDFs of mRNA and R03 are distributed quite evenly between distances of 4.35 and 95.65 Å during the early time interval (0–10 ns). In the late time interval, the RDF is split into two: one at 2.85–33.75 Å and the other at 33.75–88.25 Å, hereafter denoted as the first and the second coordination shells, respectively. The RDF is more concentrated at the first coordination shell closer to the mRNA molecule. A more detailed analysis over various time ranges (Figure S8) shows an increasing RDF near the mRNA, indicating that SM-102 molecules approach the mRNA. This movement supports our notion that charge-charge interactions between mRNA and R03 affect LNP compactness. Furthermore, the lower distribution between 55.00 and 88.25 Å in the last 500 ns interval compared to the 0–10 ns interval suggests a decrease in the size of the LNP in the Positive system.

Additionally, we computed the RDF between the citrate ions, with a charge of − 1 (citrate (− 1)), and the center mass of mRNA (citrate (− 1)-mRNA) over various time intervals (Figure S9). Interestingly, the RDF profiles at the 0–10 ns and last 500 ns time intervals (Fig. 2A, middle panel) exhibit similar patterns to those of R03 and mRNA (Fig. 2A, upper panel). The RDF of citrate (− 1)-mRNA near the center mass of mRNA increases over time. RDF profiles over stepwise time ranges (Figure S9) also show patterns like those of R03 and mRNA. This observation emerges a notion that citrate (− 1) follows the movement of SM-102 due to electrostatic interaction. Interestingly, Fig. 2A also shows an overlap of the second coordination shell of R03-mRNA and citrate (− 1)-mRNA, indicating the presence of citrate (− 1) on the surface of the LNP. Citrate (− 1) ions may pair with SM-102 and potentially compensate for the repulsive forces between the positively charged SM-102 molecules, thereby increasing ionic strength. Moreover, we analyzed the number of LNP components in each coordination shell (Table 2) and calculated RDFs between the R03 headgroup (specifically its oxygen and nitrogen atoms, denoted as O^R03 and N^R03) and the oxygen atom of citrate (− 1) (O^{citrate (−1)}) (Figure S10). The results confirm the presence of citrate (− 1) ions around R03 in the first and second coordination shells, with the highest RDF intensity at 2.65 and 2.85 Å for O^R03-O^{citrate (−1)} and N^R03-O^{citrate (−1)}, respectively. Citrate (− 1) neutralizes the surface of the LNP, potentially playing a role in LNP aggregation and resulting in larger particle sizes. This notion is consistent with the findings of Nakamura and coworkers that the use of a citrate buffer during LNP formation resulted in relatively larger particles compared to the use of acetate and lactate buffers³⁸.

Since the analyses of gyration radius, mass density, SASA, and RDF suggest that the LNP in the Positive system decreases in size over time, it is likely that water molecules inside the LNP are released. To further investigate this hypothesis, we calculated the RDF between water molecules and mRNA's center of mass (water-mRNA) over different time ranges. The results in Fig. 2A (lower panel) and Figure S11 describe the distance between the center mass of mRNA and water molecules not atoms to atoms (see Figure S7). Figure 2A (lower panel) indicates that the RDF from 0.00 to 50.00 Å in the last 500 ns time range is lower compared to the 0–10 ns range. However, the distribution from 50.00 to 80.00 Å in the last 500 ns is higher than that in the 0–10 ns range, suggesting that most water molecules inside the LNP are released. This is confirmed by the fluctuation in the number of water molecules in the first coordination shell throughout the simulation (Figure S12). The high peak close to 0 Å in the last 500-ns range (Fig. 2A, lower panel) suggests that some water molecules form intensive interactions with the mRNA molecule. Furthermore, the RDF from 30 to 50 Å declines over time, reflecting that the channels or clefts in the LNP of the Positive system become smaller, which suggests the LNP becomes more compact.

For the LNP in the Neutral system, we initially organized the R03 headgroup at the inner region, while the neutral charged of R03 headgroup (03N) at the outer region (see Methods). Figure 2B (upper panel) shows that the RDF of R03 and the center of mass of mRNA (R03-mRNA) remains constant. Meanwhile, the RDF of 03N and the center of mass of mRNA (03N-mRNA) (Fig. 2B) slightly moves toward the mRNA in the last 500 ns compared to 0–10 ns. It does not have a sharp peak, which may be due to the migration of some 03N molecules toward the inner region of the LNP, where the distribution is classified into the coordination shells at 13.45–75.85 Å. Moreover, the RDF between 55.00 and 77.45 Å decreases at a time range of last 500 ns compared to that of 0–10 ns (Fig. 2B). Additionally, overlapping RDF profiles of R03-mRNA in both systems in the last 500 ns more clearly indicate that the exterior part of LNP in the Neutral system is contracted (Figs. 2C and 3N-Positive) compared to the Positive system (Fig. 2C, R03-Positive).

Like the profile in the LNP Positive system, the RDF profile of citrate (− 1)-mRNA in the Neutral system follows that of R03-mRNA (Fig. 2B, middle panel). Citrate (− 1) ions are around R03 as shown in Figure S13 and Table 2, with the highest RDF intensity at 2.65 and 2.75 Å for O^R03-O^{citrate (−1)} and N^R03-O^{citrate (−1)}. Both RDFs of R03-mRNA and citrate (− 1)-mRNA slightly move away from the mRNA at the late time range. Meanwhile, the RDF of water-mRNA (Fig. 2B, lower panel; Figure S14) shows a decrease between 25 and 45 Å in the last 500 ns. This decline implies a possible reduction in the size of water permeation channels within the LNP, indicating increased structural compaction over time.

RDFs of citrate ions bearing charges of − 3 (citrate (− 3)) and the center mass of mRNA (citrate (− 3)-mRNA) show distributions of around 10.00 to 100.00 Å at the time range of 0–10 ns (Figure S14), suggesting that the diffusion of citrate (− 3) ions quickly occur at the early time range. Interestingly, at the time range of last 500 ns, the RDF is split to between around 5.85–36.35 Å and beyond 40.95 Å, leaving the distances between 36.45 and 40.85 Å having zero density. This observation suggests that some citrate (− 3) ions may move approaching SM-102P, with a positive charge, while others may move away from the outer surface of LNP. The analysis of the number of LNP components in each coordination shell (Table 2) shows the presence of a small number of citrate (− 3) ions in the first coordination shell. The majority of citrate (− 3) is located outside the LNP, as shown in Figure S15. Citrate (− 3) ions surrounding the outer surface coordinate with Na⁺ ions, which is confirmed by the RDF between Na⁺ and the oxygen atoms in citrate (− 3) (Na⁺-O^{citrate (−3)}) (Figure S16), with the highest RDF intensity at 2.35 Å. Citrate ions are believed to act as a chelators for Na⁺ ions, and thus, using citrate as a buffering agent may potentially reduce mRNA hydrolysis³⁹.

The analysis of the number of LNP components in each coordination shell (Table 2) in the Positive and Neutral systems show the presence of R03 around mRNA, supported by the RDF analysis between the R03 headgroup (N^R03 and O^R03) and the oxygen atoms of mRNA (O^mRNA) (Figure S17). The highest RDF density values for N^R03-O^mRNA and O^R03-O^mRNA in the Positive and Neutral systems are 2.75 and 2.65 Å, respectively. In contrast, in the system where all R03 molecules are deprotonated to 03N (All-Neutral system), no sharp peak is observed in the RDF curve between the 03N headgroup (N^03N and O^03N) and O^mRNA. Additionally, its density profiles are lower than those in the Positive and Neutral systems, indicating a minimal presence of 03N around mRNA (Figure S17). This suggests the crucial role of electrostatic force in the interaction between mRNA and SM-102.

The MMGBSA energy decomposition analyses (EDAs) between mRNA and all components in the Positive system are visualized in Fig. 3A. The figure indicates that R03 exhibits the strongest binding energy with mRNA among others. Further analysis (Fig. 3B) shows that the strong binding energy arises from the electrostatic energy term, supporting our hypothesis that SM-102 influences LNP compactness and stability through interactions between negatively charged mRNA and positively charged R03. Positive values of polar solvation energy in Fig. 3B suggest that both R03 and mRNA are favorably solvated by water molecules, which may explain water molecules remaining in the interior of LNP (Fig. 2A, lower panel).

EDA (Fig. 3A, citrate (− 1)) also exhibits intriguing result that citrate (− 1) ions also display attractive binding energy with some mRNA residues, which is unexpected since their net charges are negative. Visual inspection (Figure S18A) unveils that citrate (− 1) ions interact with phosphate and base residues of the mRNA via hydrogen bonds (H-bonds). These results may explain the RDF peak of citrate (− 1)-mRNA between distances of 2.95 and 10.55 Å at last 500 ns (Fig. 2A, middle panel). Further visual analysis (Figure S18B) reveals citrate (− 1) ions bridging between the R03 headgroups via H-bonds and electrostatic interactions. Moreover, the analysis of MMGBSA interaction energy terms shows that some citrate (− 1) ions form highly electrostatic interactions with R03, as portrayed in Figure S19. These results support our previous hypothesis (in the Radial Distribution Analyses section above) that citrate (− 1) ions may pair with R03 and potentially compensate for the repulsive forces between the positively charged SM-102 molecules, thereby increasing ionic strength. These results also address our other hypothesis, from RDF analyses, that citrate (− 1) ions follow the movement of SM-102 due to electrostatic interaction.

To gain a deeper understanding of the LNP self-assembly process, we calculated the diffusion constants (Fig. 4A, B, C, and Figure S20) and the linear interaction energy (LIE) (Fig. 5A, B, C, D, and E). R03 exhibits the highest diffusion constant during the initial 10 ns simulation period within the first coordination shell (Fig. 4A, green), corresponding to the movement of R03 towards mRNA in the Positive system. R03 and mRNA in the first coordination shell experience an electrostatic attractive force, as shown in the LIE plot (Fig. 5A, blue histogram), while those in the second coordination shell do not (data not shown). R03 also attracts some DSPC molecules, through phosphatidylcholine moieties via electrostatic interactions (Figure S21). From LIE (Fig. 5B and C, blue histogram) and EDA (Figure S22), we also found that SM-102P and other lipids in the Positive system attract each other through van der Waals (vdW) forces. In the LNP formation, vdW forces between lipids may play a crucial role in LNP self-assembly, which is consistent with the work reported by Maeki and coworkers⁴⁰. The self-assembly of lipid hydrophobic chains is a critical stage in LNP formation, driven by increased solvent polarity, due to ethanol dilution by an aqueous buffer, to reduce surface energy^40,41.

In the Neutral system, the profiles of EDA and MMGBSA interaction energy terms between mRNA residues and all LNP components (Figure S23A and B) follows a similar pattern to those in the Positive system (Fig. 3A and B). Like in the Positive system, electrostatic interactions significantly contribute to the total interaction energy between R03 and mRNA. Meanwhile, in diffusion constant, R03 exhibits low fluctuations (Fig. 4B). The high fluctuations in the diffusion constant of 03N during the first 100 ns of the simulation (Fig. 4C) are related to LNP contraction. This contraction is the result of van der Waals (vdW) strengthening among lipids due to the protonation change of R03 to its neutral form (03N) as the pH shifts to 7.4. In the Neutral system, vdW interactions of SM-102 and other lipid components in the second coordination shell are stronger than those in the Positive system (Fig. 5C, D, and E). Charge alteration to neutral causes 03N more hydrophobic than R03, as reflected by their dipole moments of 1.1763 and 6.8580 Debye, respectively, leading to a greater tendency for SM-102N to interact with other lipids to reduce surface energy. This enhanced lipid–lipid affinity may also explain the observed increase in the radius of gyration of cholesterol in the Neutral system (Fig. 1A). Notably, ethanol, which is commonly used in LNP production, can disrupt hydrophobic interactions⁴² thus potentially affecting LNP formation. Experiments confirmed that rapid ethanol dilution is crucial for producing small and uniform LNP sizes^40,43.

Our MD simulations primarily highlight the critical role of the protonation state of ionizable lipids in mRNA encapsulation and LNP formation. LIE calculations for the All-Neutral system, in which all R03 moiety are deprotonated to 03N, indicate the absence of electrostatic interactions between 03N and mRNA (Fig. 5A, red). This finding is consistent with RDF calculations between N^03N, O^03N, and O^mRNA (Figure S17). Consequently, in the absence of positively charged lipids, the mRNA remains unencapsulated (Fig. 6).

Ionic strength, pH, and solvent polarity are interrelated parameters that may modulate mRNA payload and LNP size. Elevated ionic strength is known to shift the protonation equilibrium of ionizable lipids toward their neutral form⁴⁵ and to decrease encapsulation efficiency in citrate‑based formulations^46,47. We hypothesize that higher pH and ionic strength in experimental formulations would increase the fraction of neutral SM‑102. Additionally, higher ionic strength may promote stronger interactions between mRNA and counterions such as Na⁺, reducing the efficiency of encapsulation and promoting aggregation^38,44,48. Meanwhile, adjusting the solvent polarity is crucial to determine the ratio of water to ethanol in LNP formation, which is suspected to be related to the proportion of SM-102P and SM-102N. Future simulations exploring variations in ionic strength, pH, and solvent polarity will be necessary to dissect their relative contributions to mRNA encapsulation.

In the Positive system, the LNP undergoes dramatic changes throughout the simulation (Fig. 7, upper panel). The initial configuration reveals water molecules within the lipid nanoparticle (LNP). During the first 10 ns of the simulation, these water molecules migrate out of the LNP through a large channel that formed, providing a pathway for their escape. This migration of water molecules may be due to strong electrostatic interactions between R03 and mRNA (Fig. 3B), squeezing the LNP assembly. As the simulation progressed, the channel is gradually closed, trapping some water molecules inside. This closing coincides with increases electrostatics interactions between R03 and mRNA (Fig. 3B) as well as the hydrophobic parts of SM-102 (R01 and R02) and other lipids (Figure S24). These interactions suggest that the reorganization within the LNP, driven by these forces, stabilized the structure and limited further water egress. The positive values of polar solvation energy suggest favorable solvation by water molecules, which may explain the retention of some water molecules in the interior of the LNP (Fig. 2A, lower panel, and Fig. 7). Additionally, some citrate (− 1) ions are adsorbed onto the surface of the LNP, likely due to pairing with R03 molecules residing on the LNP exterior (Figure S18B). Furthermore, as the simulation progresses, DMG-PEG2000 becomes increasingly buried within the LNP, while the PEG2000 segment folds (Fig. 7, upper panel, cyan), aligning with other studies that report PEG2000 adopting a prolate spheroidal⁴⁹.

In the Neutral system, MD production trajectory visualizations differ from those in the Positive system (Fig. 7). This difference may be attributed to the neutralization of R03 to 03N on the LNP exterior, which leads to stronger vdW interactions among lipids and enhances the integration of SM-102N with other lipids. Moreover, conversion of R03 to 03N causes citrate (− 3) ions desorbing from the LNP surface (Figure S15) and forming clusters with sodium ions as discussed in the Radial Distribution Analysis section.