MOLECULAR SIMULATION METHODS FOR pH RESPONSIVE MOLECULES

Abstract

pH is an important parameter in many chemical systems, as it determines the protonation state of an ionisable site in a molecule and thus affects the structure, dynamics and function of the molecule in a solution. Although the Henderson-Hasselbalch equation is widely used to determine the ionisation state of a molecule at a given pH, it is based on several approximations and do not account for structural changes at the molecular scale and the local environment of the molecule. Molecular simulations have the potential to bridge this gap and thus use as a predictive tool to determine the pH-effect in novel molecules. However, in conventional molecular simulations, the protonation state of a system is fixed (and not the pH) and cannot adapt to the local environment. In this study, we have developed a more realistic and thermodynamically rigorous scheme, where the protonation state of an ionisable site in a molecule is allowed to change by continuous protonation-deprotonation processes, as will be expected at a given pH condition. This method is based on a ! −dynamics approach, where we introduce an additional dimensionless degree of freedom (“particle”), !, for every ionisable site indicating its protonation state (! ≈ 0 for fully protonated and ! ≈ 1 for fully deprotonated). This ! particle is propagated in time by the Newton’s equation of motion using the interpolated forces between protonated and deprotonated states. We treat each ionisable site as a mixed state – linear combination between deprotonated and protonated states; free protons are not handled explicitly. This method relies on precalculated empirical functions for each ionisable group (reference free energy simulations), which are obtained by running multiple conventional molecular dynamics (MD) simulations on our system. These constant pH simulations are computationally intensive as the energetics of charge changes upon protonation and deprotonation must be rigorously modelled and such simulations must sample large number of protonation states to give reproducible results. For performing the simulations, we first developed and tested our own library of MD code in C++ and then incorporated it as a patch of an open-source MD software, GROMACS in order to make it more versatile and computationally efficient. Simulations were first performed on a simple and weak acid, hydrogen fluoride, to test this method and do a parametric study. This was followed by the simulations for a relatively complex molecule, acrylic acid, to test the accuracy of this method. The pKa value predicted by this method is are found to have good agreement with experimental results. The method is also able to reproduce the titration curves for these acids. Finally, we performed some preliminary 2 simulations of polyacrylic acid, which could not be completed due to time constraints. Nonetheless, this study provides a strong foundation and a detailed parametric study to simulate more complex molecules.