For some simple colloidal systems, one can use the Poisson Boltzmann theory with the charge regulation boundary condition to study the degree of protonation. However, this approach breaks down for suspensions containing salts with multivalent counterions or when studying flexible molecules -such as proteins, polyelectrolytes, polyampholytes etc., whose three-dimensional conformation is directly determined by the protonation state of the macromolecules. To overcome this difficulty, in 1992 Reed & Reed (RR) proposed a constant pH (cpH) Monte-Carlo simulation method, to calculate the titration curves for such systems. Originally the RR-cpH algorithm was applied to systems with implicit ions -the monomers interact with each other through the screened Coulomb potentialsoon after, however, extensions to systems with explicit ions were proposed. To preserve the overall charge neutrality in such approaches, a protonation move is combined with an arbitrary deletion of one hydronium ion inside the cell and a deprotonation move with an arbitrary insertion of a hydronium ion at a random position inside the cell. Such arbitrary insertion/deletion of ions inside the simulation cell, however, violates the detailed balance condition, leading to incorrect results. The problem with the RR-cpH algorithm, when applied to explicit ions, was pointed out some 20 years ago by Labbez & Jönsson, but unfortunately went unnoticed inside the physical chemistry and soft matter communities, where RR-cpH became very popular to treat various systems in which charge regulation plays important role. An alternative canonical approach to calculate the titration curves was proposed by Ullner and Woodward. These authors used a closed spherical simulation cell with a fixed number of protons. They then used Widom insertion method to obtain the chemical potential inside the system for a given protonation state. Contrary to a common misconception, RR-cpH algorithm is not a canonical algorithm. For isolated systems the total number of protons is conserved and the titration algorithm must use the correct canonical acceptance probabilities for the protonation/deprotonation moves, which do not explicitly depend on the pH inside the simulation cell. In a canonical (isolated) system, pH is not a control parameter -it can not be specified a priori, but only be calculated a posteriori, after the simulation has fully equilibrated and the number of protonated groups has been determined. Therefore, if one wants to modify the RR-cpH algorithm to properly include the explicit ions, one must take into account that the simulation cell is implicitly in contact with an external reservoir of acid and salt. A protonation/deprotonation move should then be combined with a simultaneous insertion/deletion of an anion from the reservoir, in accordance with the correct grand canonical acceptance probability. Furthermore, since in such simulation polyelectrolyte is confined within the cell, while ions can freely exchange between the system and the reservoir, the two will be at different electrostatic potentials. This is known as the "Donnan potential", ϕ D .