pH is a crucial physicochemical property that affects most biomolecules. Changes in protonation equilibrium of susceptible sites will modify the electrostatic environment and, consequently, have an effect on the molecular structure, stability, and catalysis. However, the protonation behavior of pH sensitive biomolecules is difficult to study using experimental techniques and can strongly benefit from using computational approaches. In this context, we have successfully studied several systems using the stochastic constant-pH molecular dynamics (CpHMD) method. In these studies, we were able to obtain titration curves for proteins, membranes, and peptides at the membrane water interface. In the later case, it was observed that, when the titrable groups are deeply inserted in the membrane, the conformational / protonation sampling becomes very limited. In this project, we extended the stochastic CpHMD method to introduce enhanced protonation sampling. We implemented a pH replica exchange scheme and applied it to ethylenediamine, a simple molecule with two strongly coupled macroscopic pKa values, and to hen egg white lysozyme (HEWL), a typical test system for pKa prediction methods. In the future, we will use this method to study challenging pH dependent phenomena in complex biological systems.