In this perspective, we discuss how these computational approaches -- physics-based modeling and ML-powered data science -- can collectively drive breakthroughs in LNP research. By integrating mechanistic insights with predictive data-driven models, computational studies hold the potential to guide rational LNP design, improve therapeutic efficacy, and ultimately expand the possibilities of RNA-based medicines.
Physics-based modeling refers to the use of molecular-level simulation techniques grounded in physical laws (e.g., Newtonian or statistical mechanics) to investigate the behavior, structure, and dynamics of biomolecular systems such as lipid nanoparticles. Physics-based modeling of lipid nanoparticles is a rapidly developing field, especially driven by recent advances in multiscale modeling and high-performance computing techniques. Complementary to experimental efforts for LNP formulation and characterization, physics-based modeling is expected to offer molecular-level insight into the LNP structure and interactions, essential to connect LNP composition to their activities, which ultimately provides predictive power to guide LNP design. An increasing number of publications have begun to demonstrate the effectiveness of physics-based modeling in explaining experimental observations, the self-assembly process of LNPs, and interactions with various biomolecules under different conditions. The goal of LNP physics-based modeling will be to provide accurate, high-throughput, structure-based virtual screening for LNP development and, hopefully, reduce the experimental time and cost and the need for extensive tests of composition variations. Herein, we provide a brief review of current approaches and their limitations in the physics-based modeling of LNPs, including all-atom and CG-MD, and CFD simulations, along with forward-seeing perspectives on future directions for advancement.
MD is a family of computational techniques that model the time-dependent behavior of atoms and molecules by numerically solving Newton's equations of motion. It has been widely used in physics, chemistry, biochemistry, and related areas to connect the microscopic structures of molecules to their collective or macroscopic properties, which enables the computational investigation of systems ranging from simple argon liquid to complex biological systems like coronaviruses. A primer text is available for readers who are new to MD. More specifically, all-atom (AA) MD is a well-established technology for simulating lipid membranes and membrane-protein interactions, with numerous applications primarily aimed at enhancing our understanding of membrane dynamics, membrane remodeling processes, and membrane proteins. Recently, AA-MD models have also been used to examine the structure and dynamics of LNPs, although accurately modeling the protonation states of ionizable lipids in various membrane environments relevant to LNPs remains challenging. Importantly, the protonation states of ionizable lipids in LNPs -- factors that affect the overall charge of the LNPs as well as their interactions with biological systems -- are often environment-dependent when the pK values of ionizable sites are near the pH of the solution. This can significantly influence the overall charge and interactions of an LNP with cells and surrounding biological media (e.g., proteins binding to an LNP as part of the biocorona). Due to this environment-dependent nature of ionizable sites, the protonation states can also be affected by specific manufacturing conditions (such as the type of dialyzing buffer used during LNP production, which is known to influence the transfection efficiency of LNPs) and the types and concentrations of helper lipids surrounding a particular ionizable lipid. To address these challenges, it is essential to utilize more precise, constant pH molecular dynamics (CpHMD) models. Notably, a scalable CpHMD model has been reported, which performs at comparable speeds to standard MD models. This method implements l-dynamics based on the linear interpolation of partial charges between protonated and deprotonated states of appropriately parameterized ionizable sites. The additional computational cost associated with parameterization is offset by the substantial increase in performance, which allows for hundreds of ionizable sites to be modeled simultaneously. We anticipate that these models will effectively capture environment-dependent effects within LNPs, similar to how they can model the protonation states of peptides and permeation enhancers integrating into membranes during oral peptide absorption. Very recently, scalable CpHMD models have been implemented for LNP modeling and were shown to accurately reproduce the apparent pK values for different LNP formulations (mean average error (MAE) = 0.5 pK units) in which pH-dependent structures are observed.
Overall, a key strength of atomistic adaptive membrane models is their accuracy in capturing complex supramolecular interactions, such as the hydrophobic effect, which dictates membrane self-assembly. Entropy plays a significant role in these molecular interactions among various lipid components within the membrane, as well as in the interactions at the membrane-solvent interface. However, a major challenge associated with AA-MD models is their relatively high computational cost due to the need to treat all the atoms in the system explicitly, particularly the solvent molecules, which often represent more than 70% of the total atoms present. Some of these challenges can be addressed by establishing reduced model systems, such as bilayer or multilamellar membrane models combined with periodic boundary conditions to approximate larger lipid nanoparticle (LNP) structures. Furthermore, enhanced sampling techniques -- including umbrella sampling, metadynamics, replica exchange MD, steered MD, and biased MD -- can be employed to model events occurring on timescales that exceed the current capabilities of AA models. These advanced sampling techniques are specifically designed to improve the sampling of rare events during MD simulations, which would otherwise be extremely difficult to observe within the limited timeframes that can be simulated with classical MD. We anticipate that this enhancement will ultimately allow AA-MD simulations to model rare events crucial for LNP function. This includes but is not limited to membrane reorganization processes that occur during LNP manufacturing or the endosomal release of LNP-encapsulated RNA from endosomes.
Nevertheless, each collective variable (CV) sampled using enhanced sampling methods incurs significant additional computational costs. This limitation restricts the number of CVs that can be efficiently sampled. Furthermore, defining reasonable CVs for enhanced sampling often requires a hypothesis about a molecular mechanism, which makes the simulation outcomes dependent on these initial assumptions. This dependency can hinder the exploration of the potential energy surface for CVs that aren't well-represented in the selected set for enhanced sampling. To address this issue, it is essential to develop new multiscale computational techniques that can better bridge models at different resolutions hierarchically, enabling the exploration of systems over larger time and spatial scales without sacrificing the accuracy of all-atom models. Machine learning (ML) and artificial intelligence (AI) will be crucial in these efforts, facilitating effective feature representation and linking various models for coarse-graining and back-mapping tasks.
CG-MD is a simulation approach in which groups of atoms are represented by simplified interaction sites, allowing for the modeling of larger systems and longer timescales compared to all-atom MD simulations. MD simulations of coarse-grained (CG) models help understand the detailed molecular structures and mechanisms of LNPs, which are often difficult to characterize experimentally. Unlike AA models, there is a variety of CG models, ranging from the highly CG/low-resolution ones (e.g., 1 to 3 CG sites per lipid) to relatively fine-grained/high-resolution ones (e.g., over 6 sites per lipid). In the popular Martini-CG model, a typical lipid is represented by around 10-15 CG sites per lipid, with the key principle being a "four-to-one mapping" where ~4 heavy atoms are represented by a single CG site. The number of CG sites per lipid can vary slightly depending on the lipid structure, which can result in heterogeneity in the CG model and the resulting dynamics. The fine-grained CG models like Martini-CG retain essential chemical details of LNP and greatly facilitate parameterization and back mapping to AA models, which are useful to simulate LNPs with different lipid and nucleic acid compositions. Further reducing the model resolution, the highly CG models are useful to simulate LNPs on more relevant temporal and spatial scales, and thus suitable to study the LNP self-assembly, size dependence, mechanical properties, etc. The highly CG models are also limited in the chemical details and complexities, and their parameters are often not transferrable, which requires significant efforts to develop and validate such models. However, many tools have been developed to automate CG model construction and parameterization.
Given the pros and cons of AA and CG models, hierarchical simulations (Fig. 2) that combine multiple models seamlessly may be a way to get the best of both models, allowing for AA accuracy and CG efficiency. Current hierarchical simulations have been categorized by how information is transferred between different resolutions -- in serial or in parallel. (i) The serial multiscale method carries out modeling at different resolutions in sequence, which takes advantage of sampling efficiency at lower resolutions and detailed accuracy at higher resolutions. For instance, one can start modeling from the least detailed model and ultimately obtain a fully atomic model. This so-called top-down modeling is promising to simulate complex systems like LNPs. (ii) The parallel multiscale methods include two different classes. The "hybrid resolution" methods combine AA or united-atom (UA) models of a given subsystem of interest with a CG representation of the environment. New parameters, however, are often needed to account for the cross interactions between two resolutions. In short, these hierarchical methods can be useful to study LNPs, but many key issues, such as transformations between multiple resolutions, sampling effectiveness, and simulation protocol optimization, still need to be studied systematically to advance their applications to systematic LNP simulation and, eventually, LNP development.
In the synthesis of LNPs, achieving rapid and uniform mixing is crucial for producing particles with well-defined sizes and high encapsulation efficiency. To produce LNPs with low polydispersity via antisolvent precipitation, the process requires mixing times on the order of 100 ms. Research indicates that confined impinging jet mixers (CIJMs) and multi-inlet vortex mixers (MIVMs) are effective for facilitating rapid solvent exchange and nanoprecipitation. CFD simulations can be used to better understrand fluid flow and mixing dynamics in different mixers.
Microfluidic mixing has played a key role in the self-assembly of LNPs at the lab scale. A key challenge in these systems is achieving efficient mixing at low Reynolds numbers, where turbulence is largely absent, making diffusion the dominant transport mechanism. Diffusion-based self-assembly is impractical due to its slow timescales, making hydrodynamic mixing essential for rapid nucleation and controlled growth. Staggered herringbone mixers have been shown to produce monodisperse LNPs, but their low throughput presents a challenge. While parallelization can increase throughput, it also adds complexity and cost to the system. Higher throughput LNP production can occur using inertial micromixers at higher flow rates. Among these microfluidic mixers, Dean vortex-based micromixers are suggested for LNP manufacturing due to their ability to maintain efficient mixing at high throughput. Dean vortex-based micromixers use curved microchannels to generate transverse rotational flows, known as Dean vortices. These vortices arise due to flow instabilities in curved geometries and actively moving fluid between different regions of the channel, enhancing mixing even at low Reynolds numbers. This passive design offers effective mixing without complex structures. There is a critical transition regime in these devices, which influences the optimal flow conditions for LNP formation. For achieving LNPs with optimal encapsulation efficiency, charge, and monodispersity, it is crucial to operate above this transition regime, as performance is compromised when operating within or below it. These insights highlights the importance of computational fluid dynamics to define the physical parameters necessary for consistent LNP quality.
CFD has been instrumental in analyzing and optimizing mixing, providing insights into flow behavior and mixing efficiency . Various passive micromixer designs have been developed to enhance mixing performance, including split-and-recombine (SAR) micromixers, staggered herringbone mixers, and Dean vortex-based mixers. These designs enhance mixing by stretching and folding fluid layers, thereby increasing the interfacial surface area available for diffusion.
Large Eddy Simulations (LES) and Direct Numerical Simulations (DNS) have been extensively used to investigate turbulence-driven mixing in these systems, understanding the role of self-sustained oscillations and flow structures on mixing uniformity. Studies on confined impinging jet mixers (CIJMs) suggest that turbulent structures impact mixing and encapsulation efficiency.
Computational studies can be used to evaluate mixing dynamics for different micromixer designs. High-fidelity CFD simulations provide a detailed understanding of fluid dynamics, mixing efficiency, and nanoparticle size, complementing experimental measurements. Computational approaches enable researchers to investigate a broad range of design parameters, flow conditions, and geometric modifications saving time and reducing costs. By systematically examining the effects of flow regimes, e.g., Reynolds number (the ratio inertial to viscous forces), chaotic flow structures, and turbulence-driven mixing, these studies can help optimize mixing platforms for enhanced nanoparticle properties, encapsulation efficiency, and scalability.