Solution pH is a physicochemical property that has a key role in cellular regulation, and its impact at the molecular level is often difficult to study by experimental methods. In this context, several theoretical methods were developed to study pH effects in macromolecules. The stochastic titration constant-pH molecular dynamics method (CpHMD) was developed by coupling molecular sampling methods, which are appropriate to study the conformational ensemble of biomolecules, with continuum electrostatics approaches, which properly describe pH-dependent protonation states. However, in difficult cases, the protonation sampling can be too slow for the commonly accessible computational times. In this work, we combined a pH replica exchange scheme with this CpHMD method and explored several optimization strategies and possible limitations.