Lipid Nanoparticle Database for Understanding Structure and Designing Nucleic Acid Delivery

We have developed a web tool for LNPDB available at https://lnpdb.molcube.com↗. The interactive website allows users to view and search the database (Fig. 1c). Users can search for specific LNPs by properties such as library source, atomic characteristics of ionizable lipids, types of helper lipid, cholesterol, and PEG lipid, molar ratio ranges, and experimental properties. Alternatively, users can search for LNPs by ionizable lipid structure or sub-structure, either by selecting from a list of head, linker, and tail groups, or by using a chemical structure drawing sketch tool. The search functionality built into LNPDB allows researchers to systematically analyze the current LNP landscape and identify any underexplored regions of chemical space for potential lipid innovation. Additionally, researchers can deposit their own LNP structure-function data using the standardized template provided on the website. This helps ensure that LNPDB grows over time as the LNP field evolves.

Beyond global diversity patterns, LNPDB features diverse compositional and experimental features. As shown in Fig. 2b, the molecular weight of the 12,845 unique ionizable lipids ranges from 201.31 to 3984.45, with a mean of 864.24 ± 393.12. The number of nitrogens present in each ionizable lipid ranges from 1 to 17, with a mean of 2.54 ± 1.41. There are 12,000 ionizable lipids with only tertiary nitrogen(s); 6437 with both tertiary and secondary nitrogens; 299 with tertiary, secondary, and primary nitrogens; 677 with only secondary nitrogen(s); and 8 with both secondary and primary nitrogens.

As for the other components besides the ionizable lipid, among the 19,528 LNP formulations, the ionizable lipid-to-nucleic acid mass ratio is most often set to 10, with a range from 0.86 to 44.58 and mean of 10.78 ± 4.45 (Fig. 2c). The distribution of helper lipid type consists of DOPE (47.7%), DSPC (28.3%), DOTAP (12.6%), none (2.9%), MDOA (2.9%), DDAB (1.8%), 14:0 PA (1.8%), and 18:0 PG (1.8%). The distribution of PEG lipid type consists of DMG-PEG2000 (57.1%), DMPE-PEG2000 (32.4%), unreported (9.2%), none (0.9%), DSG-PEG2000 (0.1%), ALC-0159 (0.1%), DMG-C-PEG2000 (0.1%), C8-Ceramide-PEG2000 (0.1%), C16-Ceramide-PEG2000 (0.1%), DSPE-PEG2000 (0.1%), DMG-PEG5000 (0.01%), DPG-PEG2000 (0.01%), DPG-PEG5000 (0.01%), DSG-PEG5000 (0.01%), DOG-PEG2000 (0.01%), and DOG-PEG5000 (0.01%).

For experimental features, the distribution of nucleic acid cargo consists of mRNA (74.7%), siRNA (19.2%), and pDNA (6.1%) (Fig. 2d and Supplementary Fig. 1b). With respect to the type of cargo encoded by the nucleic acid, the distribution consists of firefly luciferase (90.9%), DNA barcode (3.4%), peptide barcode (2.1%), human erythropoietin (1.4%), Factor VII (0.6%), green fluorescent protein (0.6%), and renilla luciferase (0.3%). The primary delivery target involves in vitro (78.4%), lung epithelium (9.6%), liver (4.6%), muscle (2.5%), spleen (1.3%), multiorgan (1.0%), heart (0.5%), lung (0.5%), and kidney (0.5%). For the preparation method, 93.8% of LNPs were handmixed; the remaining 6.2% were prepared via microfluidics. The readout methods report luminescence (75.2%), discretized luminescence (15.7%), protein abundance (3.9%), cellular uptake (2.5%), editing efficiency (0.7%), LRP6 knockdown (0.6%), diameter (0.5%), zeta potential (0.5%), and percent hemolysis (0.5%). Luminescence measurements quantify LNP delivery performance by reporting the level of nucleic acid transfection in target in vitro or in vivo systems. Additional summary statistics are shown in Supplementary Fig. 1.

In a prior study¹⁷, we introduced lipid optimization using neural networks (LiON), a deep learning model for learning ionizable lipid design towards predicting LNP delivery efficacy. LiON uses deep message-passing neural networks (D-MPNNs)⁶⁵ to learn representations of ionizable lipid structure, while additional formulation details such as component ratios and experimental context are appended as auxiliary features to guide prediction.

The dataset used to train this original version of LiON involved 8727 LNP formulations¹⁷. Building on this foundation, LNPDB expands the dataset more than twofold by incorporating an additional 10,801 LNP formulations, bringing the total to 19,528. Beyond scale, LNPDB captures a more diverse and descriptive set of features for each formulation. The newly added formulations broaden the diversity of ionizable lipids and also place greater emphasis on varying the types and ratios of the other three LNP components. Moreover, as detailed in the next subsection, unlike the original dataset for LiON, LNPDB includes MD–ready CHARMM force field files for all constituent lipids, introducing a new, physics-based modeling modality altogether for assessing LNP structure-function relationships.

With our improved LiON model trained on LNPDB, we next sought to compare the predictive performance of our model with another published LNP deep learning model, AGILE¹⁶. We find that LiON trained on LNPDB achieves significantly better predictive performance compared to AGILE for the 4 held-out test sets evaluated across different delivery targets (Fig. 3b and Supplementary Fig. 3b). Model performance was assessed using five train/validation splits with a fixed held-out test set, with mean correlation coefficients and standard deviations reported across folds. LiON trained on LNPDB has 16-fold more training data compared to AGILE, potentially enabling it to learn a broader range of structure-function relationships and generalize better to unseen data. Altogether, LNPDB supports an improved deep learning model for predicting LNP delivery performance. Importantly, LNPDB establishes a framework of training data for the continued development of next-generation deep learning models for LNP design. Moreover, given that our results demonstrate significant variation in model accuracy across libraries, future research can leverage the training data of LNPDB to design alternative models that may be better suited for the specific LNP design strategy (e.g., helper lipid optimization) or delivery target (e.g., in vivo muscle) under investigation.

By providing CHARMM force field topology and parameter files for each lipid in 19,528 LNP formulations, LNPDB supports all-atom MD simulations of the full dataset, as well as any new formulations constructed from its constituent lipids. LNPDB represents a substantial advancement over existing CHARMM-GUI resources, which were limited to ionizable lipids comprising only 6 different head group types and 5 different tail group types⁶⁶. Moreover, new lipids uploaded to LNPDB will be automatically parametrized for MD simulation as well.

MD offers a complementary alternative to machine learning models for understanding the structure-function relationships of LNPs, differing not only in the type of data produced but also in the way the data is generated. MD simulations yield three-dimensional, temporal structural information for all constituent lipids—not just the ionizable lipid—capturing both type and ratio. For small organic molecules such as those comprising LNPs, the conformational dynamics and interactions with neighboring molecules may play an outsized role in function, making MD especially informative. The importance of MD is amplified by the limited accessibility of three-dimensional structural data for lipids, as the experimental techniques SAXS^21,22 and cryo-EM^22,23 are low-throughput, low-resolution, and cost-prohibitive. Furthermore, unlike machine learning approaches, MD does not require training data, which is particularly valuable given that, as shown in Fig. 2a, most new ionizable lipids lie outside the distribution of previously characterized structures, potentially limiting machine learning generalization.

For a given LNP formulation, two bilayer conditions were simulated: fully-neutral ionizable lipids and half-neutral, half-protonated ionizable lipids. These conditions represent, respectively, the neutral pH prior to cellular uptake and the early endosome environment (pH ~ 6.5) where roughly 50% of ionizable lipids would be protonated, assuming a pKa of 6.5⁷⁰. All-atom simulations (N = 134; 77 fully-neutral, 57 half-protonated) were run using OpenMM⁷¹ for 1.5 µs to allow for bilayer equilibration, which we observe occurring before 1 µs (Supplementary Fig. 6; see “Methods”). Any bilayer that did not remain intact—namely, through ionizable lipids escaping from the membrane—was terminated early. Additional details on simulation conditions are provided in Supplementary Table 2.

Snapshots of the final frame and density profiles of select LNP bilayer systems are shown in Fig. 4b. In all systems, protonated ionizable lipids are generally oriented with their head groups exposed at the membrane–water interface, whereas neutral ionizable lipids exhibited more variable behaviors—some remained at the surface (Fig. 4b center), while others were buried within the hydrophobic core (Fig. 4b left). We observed that certain simulated bilayers from the LM_2019 LNP library⁷ were unstable, with ionizable lipids dissociating from the membrane over the course of the simulation (Fig. 4b, right). Notably, we find that simulated bilayer stability is positively associated with experimental transfection, indicating in silico membrane behavior could be a useful screening criterion for the delivery potential of candidate LNPs (Fig. 4c).

Next, we aimed to analyze the simulated bilayers for additional structural features that may correlate with experimental delivery efficacy. One structural feature that we sought to quantify was CPP, which has been used to relate lipid shape to phase behavior and is hypothesized to influence endosomal escape efficiency^22,25,26,69. We computed CPP values averaged across ionizable lipids for each LNP bilayer according to two different approaches: one based on volume (CPP_V) and the other based on radii of gyration (CPP_Rg) (Fig. 4d; see “Methods”). CPP values < 1 indicate a cone shape (i.e., narrower at the hydrophobic tail region buried in the membrane than at the hydrophilic head group exposed to the aqueous interface). CPP values close to 1 indicate a cylindrical shape. CPP values > 1 indicate an inverted cone shape, which has been implicated in promoting an inverse hexagonal (H_II) phase, promoting membrane fusion and endosomal escape^22,69. Our CPP values derived from MD simulations demonstrate comparable relative differences consistent with experimental measurements⁷² for two ionizable lipids (Supplementary Fig. 7), with neutral forms having higher CPP values than their protonated counterparts (Supplementary Figs. 7 and 8), a finding which aligns with prior SAXS experiments²⁶.

We next analyzed the subset of N = 34 different LNP formulations from the LM_2019 library⁷ that formed stable bilayers during the simulations to assess if their ionizable lipid CPP values correlate with experimental delivery performance. We find that ionizable lipid CPP significantly predicts LNP performance. For protonated ionizable lipids, the CPP_V approach based on volume yields a Pearson correlation of 0.530 with delivery performance. This association strengthened when stratified by amine group: amine 12 (r = 0.761), amine 2 (r = 0.798), amine 3 (r = 0.501) (Fig. 4e top). The alternative CPP_Rg approach based on radii of gyration of the protonated ionizable lipids also shows a robust (r = 0.546) correlation with performance, with similarly improved associations when analyzed by amine group: amine 12 (r = 0.745), amine 2 (r = 0.760), amine 3 (r = 0.420) (Fig. 4f top). Compared to protonated ionizable lipids, neutral ionizable lipids in fully-neutral systems demonstrate comparably strong correlations between CPP and delivery performance (Fig. 4e–f bottom). For both CPP approaches, we also find that the significant correlation with performance holds when analyzing neutral ionizable lipids in the half-protonated systems (Supplementary Fig. 9). Moreover, when we focused our analyses on the subset of LNPs with mean CPP values greater than 1—corresponding theoretically to a transition to negative curvature—the correlative performance improved for both the CPP_V method (protonated: overall r = 0.723; neutral: overall r = 0.680) and CPP_Rg method (protonated: overall r = 0.621; neutral: overall r = 0.646) (Supplementary Fig. 10). Some correlations for amine 3 did not reach statistical significance, likely due to its limited representation of LNPs with CPP > 1. Importantly, overall, these MD-derived correlations with experimental delivery are greater than those of the LiON deep learning model for the LM_2019 fully held-out test set (r = 0.104) (Fig. 3b), underscoring the potential of MD as an alternative data-efficient modality for uncovering structure-function relationships. Moreover, we assessed whether the inclusion of PEG lipid or reducing the temperature to 298 K in simulations affected CPP and found no significant effect (Supplementary Fig. 11).

Next, we measured additional structural features of the simulated bilayers: membrane thickness, torque density, and compressibility (see “Methods”). Torque density values of the fully-neutral systems are positively associated with delivery performance (r = 0.461); however, high variance in this metric limits the strength of this conclusion. The remaining features for the simulated LNP formulations did not have any significant correlations with performance (Supplementary Fig. 12). Interestingly, CPP_V variance also significantly predicts delivery performance, suggesting that greater ionizable lipid polymorphism could allow for more effective delivery, potentially due to increased capacity to accommodate more inverse-conical lipid geometries (Supplementary Fig. 13).

To collect the data for the 19,528 LNPs featured in this initial version of LNPDB, we followed the same method as introduced in our prior study¹⁷. In summary, publications were selected from the literature based on the presence of large screening datasets, primarily focused on ionizable lipids, to allow for meaningful within-dataset comparisons. Additional publications were selected to broaden the representation of helper lipids, cholesterols, and PEG lipids. SMILES were created for each lipid for each publication. Functional data—most commonly delivery performance—were extracted from published heatmaps and bar plots by digitizing the figures and interpolating values based on either the color scale (heatmaps) or bar height (bar plots) as defined in the accompanying legends. Because delivery values are often reported on different scales across studies and modalities, for each publication and for each delivery context (e.g., in vivo or in vitro within the same publication), functional delivery data were standardized to have a mean 0 and a standard deviation 1. When raw luminescence values spanned several orders of magnitude, they were first log-transformed prior to standardization. Predictive performance results were separated by publication, as datapoints were treated as directly comparable within individual screens, but not necessarily across different screens or assay modalities. Moreover, standardization was applied to prevent overemphasis of any single dataset. The deep learning models used in this study do train across datasets to maximize the diversity of trainable data, though we recognize that the ability of data from one study to inform structure-function relationships in another is limited by inherent differences in experimental protocols, measurement modalities, and assay sensitivities across studies that may introduce systematic biases. LNPDB introduces experimental condition variables (e.g., solvents, dose) towards bridging studies, but these additions can only partially mitigate the systematic differences across laboratories and experimental setups (see “Discussion”).

A total of 269 ionizable lipids from the commercial supplier BroadPharm were also included in LNPDB, representing the full set available on the vendor’s website as of June 1, 2024.

LNPDB standardizes an encoding strategy for LNPs based on three general classes of features: composition, performance, and simulation. Additional organizational descriptors include LNP ID, experiment ID, formulation ID, publication link, and publication PubMed ID.

Composition features include ionizable lipid name, ionizable lipid SMILES, ionizable lipid amine name, ionizable lipid amine SMILES, ionizable lipid linker name, ionizable lipid linker SMILES, ionizable lipid tail 1 name, ionizable lipid tail 1 SMILES, ionizable lipid tail 2 name, ionizable lipid tail 2 SMILES, ionizable lipid tail 3 name, ionizable lipid tail 3 SMILES, ionizable lipid tail 4 name, ionizable lipid tail 4 SMILES, ionizable lipid molar ratio, ionizable lipid-to-nucleic acid mass ratio, helper lipid name, helper lipid SMILES, helper lipid molar ratio, cholesterol name, cholesterol SMILES, cholesterol molar ratio, PEG lipid name, PEG lipid SMILES, PEG lipid molar ratio, fifth component lipid name, fifth component lipid SMILES, fifth component molar ratio, aqueous buffer, and dialysis buffer.

Experimental features include mixing preparation method, model (i.e., in vitro or in vivo), model system (e.g., HeLa), model target (e.g., lung), route of administration, cargo, cargo type (i.e., encoded protein), nucleic acid dose, readout method, batching approach, and readout value (i.e., most commonly a measure of delivery performance).

For simulation features, CHARMM topology and parameter files for all lipid components were generated (see “All-atom molecular dynamics simulations”). Ionizable lipids were modeled in both neutral and +1e protonated states. A majority of LNPs (13,097) in LNPDB include ionizable lipids that have more than one nitrogen, often with several plausible protonation sites. For simplicity, LNPDB assigns one representative +1e protonated state per ionizable lipid. To select the nitrogen for protonation for each ionizable lipid, the following rule-based decision tree was applied. If the lipid contained only a single nitrogen, that nitrogen was protonated. If multiple nitrogens were present, the nitrogen with the highest priority was protonated based on the following hierarchy: tertiary amine, secondary amine, primary amine, imidazole, pyridine, tertiary aromatic amine, secondary aromatic amine, and primary aromatic amine. Groups comprising amide or sulfonamide structures and quaternary nitrogens were excluded. If multiple candidates of the same class were found, the nitrogen closest to the molecular periphery of the ionizable lipid head—defined as having the greatest graph eccentricity (i.e., the longest existing graph distance to a terminal atom)—was selected. For specific cases within the KZ_2016 dataset involving ionizable lipids with tail amines, the most centrally-located candidate (i.e., the lowest average squared distance to all other atoms) was chosen. Once selected, the nitrogen was protonated by assigning a + 1 formal charge and adjusting the SMILES accordingly. Future research is warranted to explore more accurate, dynamic protonation conditions⁷⁶.

The CHARMM topology and parameter files for all lipid components were generated as follows: ionizable lipids were modeled in both neutral and +1e protonated states, with parameters assigned manually via CHARMM force field analogy mapping and supplemented by the CGenFF workflow⁷⁷. We used the standard CHARMM force field definitions⁷⁸ for helper lipids, cholesterol, and PEG lipids.

Similar to the CHARMM-GUI Membrane Builder that supports diverse lipid types^66,79, LNPDB uses the latest version of the CHARMM C36 additive force field parameters. The CHARMM force field⁷⁷ is designed to allow a modular, building-block approach to create force fields for molecules composed of components (blocks) similar to the ones already parametrized. Many topologies and parameters of lipids and carbohydrates in the latest version of the CHARMM C36 force field were generated using this building-block approach; and the generated force fields were further validated by comparing simulations with experimental data. Accordingly, we used a similar building-block approach for the ionizable lipids in LNPDB. We have not yet seen abnormal all-atom simulation behavior (e.g., lipid flip-flop); however, further force field optimization is recommended for specific ionizable lipids if abnormal behavior is observed.

All bilayer systems were assembled using the standalone MolCube Membrane Builder—a commercial software application analogous to the CHARMM-GUI Membrane Builder tool⁸⁰—in membrane-only mode using approximately 100 lipids per leaflet and solvated with TIP3P water⁸¹ and 0.15 M NaCl.

Following the six-step equilibration procedure outlined in the CHARMM-GUI Membrane Builder protocol^80,82, NVT (constant particle number, volume, and temperature) simulations were conducted at 310 K (i.e., temperature of cells treated with LNPs) with strong harmonic positional restraints on lipid heavy atoms and dihedral restraints on ionizable head groups. The restraint force constants were gradually reduced to zero over the six equilibration steps for gradual membrane relaxation. Subsequently, unrestrained NPT (constant particle number, pressure, and temperature) production runs were conducted at 310 K and 1 bar for 1.5 µs using OpenMM with a 4 fs time-step enabled by hydrogen-mass repartitioning (HMR)^83,84. Temperature was maintained via a Langevin thermostat (collision frequency 1 ps⁻¹) and pressure via a semi-isotropic Monte Carlo barostat (coupling interval 0.4 ps)^85,86. Bonds involving hydrogen were constrained with SHAKE⁸⁷; van der Waals interactions were force-switched off between 10–12 Å⁸⁸, and long-range electrostatics were treated by the Particle-Mesh Ewald method with a 12 Å real-space cutoff⁸⁹. For ten bilayer systems, we additionally assessed whether the inclusion of PEG lipid or setting the temperature to 298 K (i.e., temperature of LNP synthesis) in simulations affected CPP and found no significant effect (Supplementary Fig. 11). Information for each simulation system is summarized in Supplementary Table 2. Note that the protonated systems have larger system sizes along the z-axis than their neutral counterparts because the protonated systems require more bulk region to fully solvate the system.

We ran all MD simulations on 48 NVIDIA RTX A5000 GPUs in parallel. Across the neutral systems, the average throughput was 363.86 ± 14.93 ns/day; across the protonated systems, it was 319.11 ± 11.65 ns/day. For the LM_2019 bilayers with a run duration of 1.5 μs, this corresponds to 4.1 days per neutral system and 4.7 days per protonated system. The full LM_2019 simulation batch, run on 48 GPUs (one per simulation), completed in about 10 days. Subsequent CPP calculations for all systems, executed on 384 CPU cores in parallel, finished in 5 h, with comparable compute times for CPP_V and CPP_Rg.

LNPDB has been developed as a RESTful application built on Django REST Framework and React, with a PostgreSQL backend, together with RDKit for native chemical structure storage and search. Lipid molecules are represented as molecular graphs and indexed via RDKit fingerprints, enabling fast SQL-level structure and substructure queries. User contributions are handled through a CSV upload portal requiring citation metadata. Each submission triggers curator review; data are validated, standardized, and then ingested into PostgreSQL. Approved entries become searchable and visible in all table views and interactive plots via the same API, ensuring seamless integration of new LNP formulations.

UMAP visualizations were created of the high-dimensional embedding landscapes of LNPs and ionizable lipids (Fig. 2a). The embedding landscape for ionizable lipids is represented by the top ten principal components (PCs) of Morgan fingerprints⁹⁰ (1024 bits, radius of 3) and the top ten PCs of Mordred descriptors⁹¹. The embedding landscape for LNPs is represented by the same axes as those for ionizable lipids, plus additional dimensions for molar ratios and the top five Morgan fingerprint PCs and top five Mordred descriptor PCs for helper lipids.

UMAP visualizations were also created of the 300-dimensional embedding landscape (i.e., fingerprints) from the LiON model of LNP formulations (Supplementary Figs.and). Fingerprints were extracted from the penultimate linear layer of the LiON model’s feedforward neural network trained on LNPDB. 2b 4a

As shown in Fig. 3a, we evaluated the predictive performance of LiON (lipid optimization using neural networks)¹⁷, which is based on the message-passing neural network architecture of chemprop⁶⁵. To compare models trained on LNPDB versus the original dataset from our prior study, we computed Spearman correlation values between predicted and experimental delivery outcomes. This evaluation was performed on held-out test sets using a 70–15–15% train-validation-test split by amine, consistent with the approach used in our prior study¹⁷. Datasets shared between LNPDB and the original dataset are compared.

As shown in Fig. 3b, we evaluated the predictive performance of LiON trained on LNPDB compared to another deep learning model, AGILE¹⁶. To perform this analysis, each of the four datasets was fully held out as an external test set. For each held-out set, models were trained using five cross-validation folds with 80%-20% train-test splits on the remaining data. Spearman correlation values between predicted and experimental delivery performance for each fully held-out dataset were computed. For LiON trained on LNPDB, holding out an entire dataset reduced the training sample size. In contrast, AGILE maintained its original training size of 1200 LNPs, as the held-out datasets did not overlap with its training data. To run the AGILE model, the GitHub repository provided in the study was referenced, and the HeLa transfection data was used for training¹⁶.

To assess whether MD simulations could provide meaningful correlations with experimental delivery performance (Fig. 4b–e), we selected N = 54 LNP formulations from a prior study (LM_2019) in LNPDB, which introduced an isocyanide-mediated three-component reaction approach for ionizable lipids⁷. For the sake of modeling, we randomly selected to model the subset of LM_2019 LNPs that contain ionizable lipids with amines A12, A2, or A3; isocyanides Iso5 or Iso9; and any alkyl ketone⁷. PEG lipids were excluded from simulations, as they are typically shed prior to endosomal escape⁶⁸, the key bottleneck for effective delivery⁶⁹, and the physiological context that we aimed to model. This subset, drawn from a single combinatorial ionizable lipid library, was chosen as a representative example of systematic lipid library design commonly employed in the field, while keeping the scope feasible within computational limits.

Additional simulations were conducted for illustrative purposes (Figs. 1b and 4a, b) that contain PEG lipid or the common control ionizable lipids of DLin-MC3-DMA, SM-102, and ALC-0315. Details of all bilayer simulations analyzed in this study are provided in Supplementary Table 2.

To quantify the spatial distribution of lipid components along the membrane normal (z-axis) as shown in Fig. 4b, density profiles were computed from the final 500 ns of each 1.5 μs MD trajectory. At each frame, atomic coordinates were re-centered to have the membrane center of mass be at z = 0. For each lipid molecule, a single representative atom was used to track z-position over time: the hydroxyl oxygen for cholesterol, the phosphorus atom for helper lipids, and the nitrogen atom on the ionizable head group for ionizable lipids. These atom positions were binned along the z-axis to generate one-dimensional density histograms for each lipid type. The profiles were averaged across frames and normalized by bin width to obtain continuous density distributions, reflecting the vertical organization of each component within the bilayer.

To quantify ionizable lipid shape, we computed CPP values. For a given stable LNP bilayer simulation, we analyzed each timestep among the final 500 ns of the 1.5 μs trajectory. For each ionizable lipid, atom subsets corresponding to the head group and tail were defined. For the LM_2019⁷ LNP bilayers, atoms in the amine and isocyanide groups—both found to be generally positioned at the membrane–water interface—were assigned to the head group; all other atoms (i.e., alkyl ketones) were assigned to the tail group. In line with its conventional formula^25,92, we calculated CPP based on volume as follows:1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathrm{CPP}}}_{V}=\frac{V}{{a}_{0}{l}_{c}}$$\end{document}CPPV=Va0lcwhere \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V$$\end{document}V is the volume formed by the ionizable lipid, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${a}_{0}$$\end{document}a0 is the surface area of the head group at the water–membrane interface, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${l}_{c}$$\end{document}lc is the average distance between the head group and tail ends (Fig. 4d). The head group area \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${a}_{0}$$\end{document}a0 was computed as the cross-sectional area of the circle formed in the membrane plane.2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${a}_{0}=\pi {r}_{{\mathrm{head}}}^{2}$$\end{document}a0=πrhead2where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${r}_{{\mathrm{head}}}$$\end{document}rhead is the head group radius, calculated as half the maximum pairwise distance between head atoms in the membrane plane. The tail radius \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${r}_{{\mathrm{tail}}}$$\end{document}rtail was computed similarly using terminal tail atoms. The lipid volume \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V$$\end{document}V was estimated by modeling the molecule as a truncated cone.3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V=\frac{1}{3}\pi {l}_{c}\left({r}_{{\mathrm{head}}}^{2}+{r}_{{\mathrm{head}}}{r}_{{\mathrm{tail}}}+{r}_{{\mathrm{tail}}}^{2}\right)$$\end{document}V=13πlcrhead2+rheadrtail+rtail2

These geometrical parameters were computed frame-by-frame for each ionizable lipid molecule across the trajectory to calculate mean CPP_V values. Standard error of the mean (SEM) was also calculated. We observe that CPP_V values exhibit greater variability across lipids and time steps compared to CPP_Rg. To focus on uncertainty in the central tendency of CPP_V across lipid molecules, SEM was used in place of standard deviation for CPP_V plots.

Lipids with CPP_V > 1 exhibit an inverse cone shape, favoring negative curvature, whereas those with CPP_V < 1 exhibit a cone shape, favoring positive curvature.

We also quantified CPP for ionizable lipids using an approach based on radii of gyration, CPP_Rg. We similarly analyzed the final 500 ns of each 1.5 μs trajectory of stable bilayers. For each ionizable lipid, atom subsets corresponding to the head group and tail were defined, and their centers of mass, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${R}_{{\mathrm{COM}}}^{{\mathrm{head}}}$$\end{document}RCOMhead and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${R}_{{\mathrm{COM}}}^{{\mathrm{tail}}}$$\end{document}RCOMtail were computed at every frame.4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${R}_{{\mathrm{COM}}}^{{\mathrm{mid}}}=\frac{1}{2}\left({R}_{{\mathrm{COM}}}^{{\mathrm{head}}}\,+\,{R}_{{\mathrm{COM}}}^{{\mathrm{tail}}}\right)$$\end{document}RCOMmid=12RCOMhead+RCOMtailfrom all atomic coordinates. We then translated each lipid to the origin by subtracting the midpoint. The orientation vector was computed as5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v=\,{R}_{{\mathrm{COM}}}^{{\mathrm{head}}}-\,{R}_{{\mathrm{COM}}}^{{\mathrm{tail}}}$$\end{document}v=RCOMhead−RCOMtailand aligned to the membrane normal \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{z}=[0,\,0,\,1]$$\end{document}z^=[0,0,1] by rotating the coordinate set through the angle between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v$$\end{document}v and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{z}$$\end{document}z^. With all lipids consistently oriented, the radius of gyration in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${xy}$$\end{document}xy plane was computed as6\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${R}_{g}\,=\,\sqrt{\frac{1}{M}{\sum }_{i}\,{m}_{i}\left({x}_{i}^{2}+{y}_{i}^{2}\right)}$$\end{document}Rg=1M∑imixi2+yi2where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${m}_{i}$$\end{document}mi and (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${x}_{i},{y}_{i}$$\end{document}xi,yi) are the mass and coordinates of atom \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i,$$\end{document}i, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M=\,{\sum }_{i}{m}_{i}$$\end{document}M=∑imi. We recorded the average \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${R}_{g}$$\end{document}Rg values for tail and head atoms across all frames—denoted \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${R}_{g}^{{\mathrm{tail}}}$$\end{document}Rgtail and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${R}_{g}^{{\mathrm{head}}}$$\end{document}Rghead—and computed CPP_Rg as7\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathrm{CPP}}}_{{Rg}}=\frac{{R}_{g}^{{\mathrm{tail}}}}{{R}_{g}^{{\mathrm{head}}}}$$\end{document}CPPRg=RgtailRghead

Lipids with CPP_Rg > 1 exhibit an inverse cone shape favoring negative curvature, whereas those with CPP_Rg < 1 exhibit a cone shape favoring positive curvature. Standard deviation was directly computed from the distribution of CPP_Rg values over the analysis window (single-block averaging).

The bilayer thicknessis defined as the instantaneous difference between the average-coordinates of phosphate atoms in the upper and lower leaflets: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$({d}_{B})$$\end{document} ( ) d B \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z$$\end{document} z 8 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${d}_{B}\left(t\right)=\, \langle {z}_{P,{{\mathrm{upper}}}}\left(t\right)\rangle \,-\, \langle {z}_{P,{{\mathrm{lower}}}}\left(t\right)\rangle$$\end{document} d B z P upper , z P lower , t t t = ⟨ ⟩ − ⟨ ⟩

We averagedover the final 500 ns of each 1.5 µs trajectory to yield a single representativeper system. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${d}_{B}(t)$$\end{document} d B ( ) t \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${d}_{B}$$\end{document} d B

Monolayer torque density () was calculated from the first moment of the lateral pressure profile, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tau }_{{\mathrm{mean}}}$$\end{document} τ mean \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p\left(z\right)={p}_{T}\left(z\right)-{p}_{N}(z)$$\end{document} p z z = − p T p N ( ) z 9 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${p}_{T}\left(z\right)=\,\frac{{p}_{{xx}}\left(z\right)+{p}_{{yy}}\left(z\right)}{2},\,{p}_{N}\left(z\right)=\,{p}_{{xx}}\left(z\right)$$\end{document} p T p N p x x z z z = , = p x x p y y z z + 2

Pressure profiles (0.2 Å bins) were integrated to give leaflet torques:and averaged as 10 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tau }_{{\mathrm{upper}}}={\int }_{0}^{\frac{{L}_{z}}{2}}{zp}\left(z\right){dz},\,{\tau }_{{l}{o}{w}{e}{r}}={\int }_{\frac{{-L}_{z}}{2}}^{0}{zp}\left(z\right){dz}$$\end{document} τ upper τ l o w e r = , = ∫ 0 L z 2 ∫ − L z 2 0 z p d z z p d z z z 11 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tau }_{{\mathrm{mean}}}\,=\frac{{\tau }_{{\mathrm{upper}}}+{\,\tau }_{{\mathrm{lower}}}\,}{2}$$\end{document} τ mean = τ upper τ lower + 2

Pressure and torque calculations were performed on velocity- and position-recoupled trajectories over the final 500 ns to ensure full equilibration. As context, for a stress-free symmetric bilayer, the monolayer torque density (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document}τ) is related to the bending modulus (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${k}_{c}$$\end{document}kc) and spontaneous curvature (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${c}_{0}$$\end{document}c0) as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau=\,{k}_{c}{c}_{0}$$\end{document}τ=kcc0. A monolayer with positive curvature is convex to the head group side, while negative curvature is concave⁶⁶.

The area compressibility moduluswas determined from fluctuations in the instantaneous projected bilayer area:whereis the Boltzmann constant,the simulation temperature (310 K). Standard deviation was directly calculated from the distribution ofvalues over the same 500 ns analysis block. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${K}_{A}$$\end{document} K A \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A(t)$$\end{document} A ( ) t \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${k}_{B}$$\end{document} k B \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T$$\end{document} T \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A(t)$$\end{document} A ( ) t 12 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${K}_{A}=\,\frac{{k}_{B}T\,\langle A\rangle }{\langle {\left(A\left(t\right)-\,\langle A\rangle \right)}^{2}\rangle }$$\end{document} K A = k B T A ⟨ ⟩ ⟨ ⟩ A A t − ⟨ ⟩ 2

Further information on research design is available in thelinked to this article. Nature Portfolio Reporting Summary