phbuilder: A Tool for Efficiently Setting up Constant pH Molecular Dynamics Simulations in GROMACS.

Constant pH molecular dynamics (MD) is a powerful technique that allows the protonation state of residues to change dynamically, thereby enabling the study of pH dependence in a manner that has not been possible before. Recently, a constant pH implementation was incorporated into the GROMACS MD package. Although this implementation provides good accuracy and performance, manual modification and the preparation of simulation input files are required, which can be complicated, tedious, and prone to errors. To simplify and automate the setup process, we present phbuilder, a tool that automatically prepares constant pH MD simulations for GROMACS by modifying the input structure and topology as well as generating the necessary parameter files. phbuilder can prepare constant pH simulations from both initial structures and existing simulation systems, and it also provides functionality for performing titrations and single-site parametrizations of new titratable group types. The tool is freely available at www.gitlab.com/gromacs-constantph. We anticipate that phbuilder will make constant pH simulations easier to set up, thereby making them more accessible to the GROMACS user community.

Best Practices in Constant pH MD Simulations: Accuracy and Sampling.

Various approaches have been proposed to include the effect of pH in molecular dynamics (MD) simulations. Among these, the λ-dynamics approach proposed by Brooks and co-workers [Kong, X.; Brooks III, C. L. J. Chem. Phys. 1996, 105, 2414-2423] can be performed with little computational overhead and hfor each typeence be used to routinely perform MD simulations at microsecond time scales, as shown in the accompanying paper [Aho, N. et al. J. Chem. Theory Comput. 2022, DOI: 10.1021/acs.jctc.2c00516]. At such time scales, however, the accuracy of the molecular mechanics force field and the parametrization becomes critical. Here, we address these issues and provide the community with guidelines on how to set up and perform long time scale constant pH MD simulations. We found that barriers associated with the torsions of side chains in the CHARMM36m force field are too high for reaching convergence in constant pH MD simulations on microsecond time scales. To avoid the high computational cost of extending the sampling, we propose small modifications to the force field to selectively reduce the torsional barriers. We demonstrate that with such modifications we obtain converged distributions of both protonation and torsional degrees of freedom and hence consistent pKa estimates, while the sampling of the overall configurational space accessible to proteins is unaffected as compared to normal MD simulations. We also show that the results of constant pH MD depend on the accuracy of the correction potentials. While these potentials are typically obtained by fitting a low-order polynomial to calculated free energy profiles, we find that higher order fits are essential to provide accurate and consistent results. By resolving problems in accuracy and sampling, the work described in this and the accompanying paper paves the way to the widespread application of constant pH MD beyond pKa prediction.

Scalable Constant pH Molecular Dynamics in GROMACS.

Molecular dynamics (MD) computer simulations are used routinely to compute atomistic trajectories of complex systems. Systems are simulated in various ensembles, depending on the experimental conditions one aims to mimic. While constant energy, temperature, volume, and pressure are rather straightforward to model, pH, which is an equally important parameter in experiments, is more difficult to account for in simulations. Although a constant pH algorithm based on the λ-dynamics approach by Brooks and co-workers [Kong, X.; Brooks III, C. L. J. Chem. Phys.1996, 105, 2414-2423] was implemented in a fork of the GROMACS molecular dynamics program, uptake has been rather limited, presumably due to the poor scaling of that code with respect to the number of titratable sites. To overcome this limitation, we implemented an alternative scheme for interpolating the Hamiltonians of the protonation states that makes the constant pH molecular dynamics simulations almost as fast as a normal MD simulation with GROMACS. In addition, we implemented a simpler scheme, called multisite representation, for modeling side chains with multiple titratable sites, such as imidazole rings. This scheme, which is based on constraining the sum of the λ-coordinates, not only reduces the complexity associated with parametrizing the intramolecular interactions between the sites but also is easily extendable to other molecules with multiple titratable sites. With the combination of a more efficient interpolation scheme and multisite representation of titratable groups, we anticipate a rapid uptake of constant pH molecular dynamics simulations within the GROMACS user community.

A spontaneous mitonuclear epistasis converging on Rieske Fe-S protein exacerbates complex III deficiency in mice.

We previously observed an unexpected fivefold (35 vs. 200 days) difference in the survival of respiratory chain complex III (CIII) deficient Bcs1lp.S78G mice between two congenic backgrounds. Here, we identify a spontaneous homoplasmic mtDNA variant (m.G14904A, mt-Cybp.D254N), affecting the CIII subunit cytochrome b (MT-CYB), in the background with short survival. We utilize maternal inheritance of mtDNA to confirm this as the causative variant and show that it further decreases the low CIII activity in Bcs1lp.S78G tissues to below survival threshold by 35 days of age. Molecular dynamics simulations predict D254N to restrict the flexibility of MT-CYB ef loop, potentially affecting RISP dynamics. In Rhodobacter cytochrome bc1 complex the equivalent substitution causes a kinetics defect with longer occupancy of RISP head domain towards the quinol oxidation site. These findings represent a unique case of spontaneous mitonuclear epistasis and highlight the role of mtDNA variation as modifier of mitochondrial disease phenotypes.

Exploring fast proton transfer events associated with lateral proton diffusion on the surface of membranes.

Proton diffusion (PD) across biological membranes is a fundamental process in many biological systems, and much experimental and theoretical effort has been employed for deciphering it. Here, we report on a spectroscopic probe, which can be tightly tethered to the membrane, for following fast (nanosecond) proton transfer events on the surface of membranes. Our probe is composed of a photoacid that serves as our light-induced proton source for the initiation of the PD process. We use our probe to follow PD, and its pH dependence, on the surface of lipid vesicles composed of a zwitterionic headgroup, a negative headgroup, a headgroup that is composed only from the negative phosphate group, or a positive headgroup without the phosphate group. We reveal that the PD kinetic parameters are highly sensitive to the nature of the lipid headgroup, ranging from a fast lateral diffusion at some membranes to the escape of protons from surface to bulk (and vice versa) at others. By referring to existing theoretical models for membrane PD, we found that while some of our results confirm the quasi-equilibrium model, other results are in line with the nonequilibrium model.

How to minimize dye-induced perturbations while studying biomembrane structure and dynamics: PEG linkers as a rational alternative.

Organic dye-tagged lipid analogs are essential for many fluorescence-based investigations of complex membrane structures, especially when using advanced microscopy approaches. However, lipid analogs may interfere with membrane structure and dynamics, and it is not obvious that the properties of lipid analogs would match those of non-labeled host lipids. In this work, we bridged atomistic simulations with super-resolution imaging experiments and biomimetic membranes to assess the performance of commonly used sphingomyelin-based lipid analogs. The objective was to compare, on equal footing, the relative strengths and weaknesses of acyl chain labeling, headgroup labeling, and labeling based on poly-ethyl-glycol (PEG) linkers in determining biomembrane properties. We observed that the most appropriate strategy to minimize dye-induced membrane perturbations and to allow consideration of Brownian-like diffusion in liquid-ordered membrane environments is to decouple the dye from a membrane by a PEG linker attached to a lipid headgroup. Yet, while the use of PEG linkers may sound a rational and even an obvious approach to explore membrane dynamics, the results also suggest that the dyes exploiting PEG linkers interfere with molecular interactions and their dynamics. Overall, the results highlight the great care needed when using fluorescent lipid analogs, in particular accurate controls.