Since their introduction more than four decades ago, molecular dynamics (MD) computer simulations have come of age. Thanks to improvements in computer hardware, algorithmic developments, as well as increased accuracy of force fields, MD simulation has evolved into a predictive technique that can complement experiment by providing atomistic insights into the dynamics of complex systems. While many experimental conditions can be modelled with good accuracy, the aqueous proton concentration, or pH, is typically accounted for indirectly by constraining the protonation states of titratable residues to their, presumed, most probable form at the start of the simulation. Because the electrostatic interactions depend critically on the protonation state of the residues, the pH affects the conformational ensemble. Conversely, because the conformation can influence the proton affinity of the residues, or pK a , a direct correlation exists between pH and conformational dynamics, which cannot be captured if protonation state is kept fixed in the simulation. To overcome this limitation in classical MD simulations and include the effect of pH on the conformational sampling directly, several solutions have been proposed in the last decades. These solutions can be roughly divided into a category that relies on discrete changes in protonation states, and a second category, in which a protonation state can change continuously. More recently, a third category that relies on the transfer of protonlike particles between titratable sites, including protein residues and solvent molecules, was proposed for the Martini force field. In the discrete constant pH approaches, the protonation state of a residue can change at regular intervals of the simulation according to a Metropolis Monte Carlo criterion. To avoid a low acceptance rate due to unfavourable solvent configurations, the Monte Carlo step is performed based on free energies calculated using the approximation of either an implicit solvent representation, or a short all-atom thermodynamic integration. Most continuous approaches for MD at constant pH are based on the λ-dynamics technique developed by Brooks and co-workers. A one-dimensional λ-coordinate with fictitious mass m λ is introduced for each titratable site, and the equations of motion for these additional degrees of freedom are integrated along with the Cartesian positions of the atoms. The λ-coordinate defines the protonation state of the residue: at λ = 0 the residue is protonated and interacts with the rest of the system as such, while at λ = 1 it is deprotonated.