EXPERIMENTAL General. Calculations were performed using a previously reported version of GROMACS 2021 (patched with PLUMED 61 that incorporates scalable -dynamics CpHMD methodology implemented with the phbuidler tool. Lipids and mRNA were simulated using the CHARMM36 forcefield water was modeled using the TIP3P water-model. In this scalable CpHMD model, dynamic protonation is modeled by interpolating between atomic charges of both  states as detailed by Groenhof and Hess. An initial steepest descent minimization was performed for 5,000 steps followed by five steps of equilibration. The first two equilibration steps -NVT ensemblewas performed using v-rescale coupling to maintain 310 K with a time constant of 0.5 ps. The first NVT ensemble utilized 0.5 fs timestep for 10 ps while the second implemented a 2fs timestep for 100 ps. The last two equilibration steps -NPT ensemblewas similarly performed with a using c-rescale to isotropically maintain a 1.0 bar pressure with a compression constant of 4.5 x 10 -5 bar -1 . The timestep and length of NPT ensemble equilibrations matched those used for NVT ensembles. Short-ranged nonbonded interactions were held to a 1.2 nm cutoff distance with electorstatic interactions calculated using the PME method and constraints imposed with the LINCS algorithm. Lastly, before full length production simulations were performed, short (100 ps) production runs were done using a short timestep (0.5 ps). Membrane simulations, without mRNA payload, were performed for 10 ns using a 2 fs timestep, while self-assembly simulations were performed for 100 ns with a 2 fs timestep. No biases, besides the pH-dependent potentials implemented by the CpHMD model to simulate different pH environments were implemented on production runs.