At the molecular level, a major factor involved in controlling not only the loading efficiency but ultimately also the biodistribution and transfection efficiency is an LNP's pHdependent charge distribution. This crucial property is set by the local pKa values of the ionizable lipids, which can be drastically influenced by the local environment/conformations of the ionizable lipids. Furthermore, the distribution of charged and uncharged ionizable lipids of an LNP also governs its internal molecular structures in different pH environments, such as the formation of mesophases in acidic environments as well as its interactions with biomacromolecules (either internally as the cargo for delivery, or on the surface, as the biocorona). For all of these reasons, the apparent pKa of an ionizable lipid in an LNP is not constant throughout an LNP, but rather dependent on the detailed local environment of the ionizable lipid, such as (i) the conformation/position of the ionizable lipid within the LNP, (ii) the clustering/distribution of helper lipids and/or RNA around the ionizable lipid, (iii) the overall charge of the LNP, and (iv) the presence of proteins on the surface of an LNP (i.e. the biocorona). Yet, these crucial environment-dependent effects are very difficult to capture with current experimental methods (which are mostly designed to simply resolve the overall, apparent pKa of ionizable lipids, averaged over a whole LNP). Therefore, to advance the fundamental understanding of the complex structure-activity relationships governing LNP formulations, new computational methods are needed that can accurately capture the local protonation states of ionizable lipids in different environments.