Theoretical titration curves at different simulation time points, fit using the Henderson-Hasselbalch equation, of a prepacked DLin-MC3-DMA membrane -1ns, black; 3ns, red; and 5ns, blue (error represent the standard error for Henderson-Hasselbalch fits). (B) A representative structure of the simulated lipid bilayer system at pH = 6.0 after 5 ns of simulation time, visualized with an additional periodic image along the +x axis for clarity. The reported apparent pKa for this formulation (pKa TNS = 6.57) indicates that there should roughly be a 1:1 distribution of protonated to deprotonated lipids. While after a very short simulation time (1 ns), the computed apparent pKa value (5.12 ± 0.05) is not too far from the experimental apparent pKa value, further equilibration of the system renders the apparent pKa value of the system much more acidic for this ionizable lipid, as the structure rearranges into a kinetically trapped, unphysical state. In fact, after 5 ns, the representative structures at pH = 6 and pH = 7 only contain deprotonated lipids. Overall, this example shows that the local pKa values of the ionizable head groups are dictated by the nature of their environment, and poor correlation with experimental data, in this case indicates that pre-packed membrane conformations are not necessarily representative of the actual internal, pH-dependent structures of the LNP formulation used in the corresponding experiments. In turn, we found that selfassembly simulations (Figure ) conducted in different pH environments, which are able to more accurately model pH-dependent structures of LNP interiors without initial structural bias, provide (Table ) significantly better (MAE = 0.5 pKa units) agreement with experimental apparent pKa values for a dataset of five distinct ionizable lipids, which also includes the Dlin-MC3-DMA ionizable lipid. Previous computational efforts to predict and model effective pKa values of LNP formulations have been mostly focused on machine learning (Figure ) or classical, fixed-protonation state computational models (Figure ). While machine learning models can predict many influential properties of LNPs, such as apparent pKa, size, and zeta potential (ZP), such methods provide limited physical insight into fundamental structurefunction relationships that are confined by the quality and amount of the data provided during training. At the same time, fixedprotonation state models of LNP bilayers, performed by the Grossfield, Shelley, and Böckmann groups, were able to successfully highlight the differences in lipid aggregation as the distribution of ionized lipids changes. This fundamental work illustrates the propensity for neutral lipids to aggregate in the membrane in simulated acidic conditions and form bilayers in simulated basic conditions. Nevertheless, while some insight about the preferred position of charge vs. neutral ionizable lipids can be gained from fixed protonation-state models, fixed protonation-state methods provide a limited ability to calculate pKa values in different environments and/or broadly capture environment-dependent changes in protonation states that can take place upon conformational/environmental changes. This work now expands on the prior fixed-protonation state models by incorporating and validating dynamic protonation state MD methods and mRNA to model the self-assembly of five distinct LNP formulations. Our work advances upon a recently reported, scalable continuous constant pH MD (CpHMD) methodology for GROMACS 2021, which has very recently been implemented to model environment-dependent changes in ionizable lipid charge. Notably, preliminary work by the Böckmann group (Figure ) demonstrated pH-dependent charge distribution of a model lipid bilayer using standard CpHMD parameterization methods . In contrast, this work (Figure ) now expands CpHMD validation to a larger test set , which provides accurate experimental data for multiple ionizable lipids in an environment mimicking real LNP-delivery systems (including in the presence of mRNA), which enabled us to carefully validate theoretical apparent pKa values with experimental data. Furthermore, our work also incorporates, for the first time, mRNA and pHdependent self-assembly with CpHMD models to more accurately account for larger, pH-dependent structural changes that are commonly observed in LNP systems. In our hands, this methodology provides better agreement with experimental apparent pKa values than a simple lipid bilayer system, which resulted (Figure ) in predicted apparent pKa that were vastly different from the corresponding experimentally measured pKa values for some ionizable lipids. Last but not least, our work provides a new parametrization method for ionizable lipids, generated using a Hamiltonian replica exchange (HREX) methodology (which enhances the parametrization accuracy for challenging-to-parametrize ionizable lipids with slow torsional rotations). Overall, our new methodology provides good accuracy for modeling the dynamic protonation of ionizable lipids coexisting in different LNP microenvironments with CpHMD models, coupled with the influence of an mRNA payload with physics-based insight on the interior environment of an LNP. and Figures .