To perform the simulations we use a cubic simulation box of side length L that contains either colloidal particles, polyelectrolyte, or protein molecules in a completely deprotonated state. In the present discussion we will use a primitive model, which treats water as a uniform dielectric continuum of Bjerrum length λ B = q 2 /k B T ϵ w = 7.2 Å, where q is the proton charge and ϵ w is the dielectric constant of water. There is, however, no difficulty to modify the algorithm to account for the explicit water or to combine it with a molecular dynamics simulation. The simulation cell also contains fully dissociated salt and acid ions -H 3 O + , Cl -, and Na + . We start the simulation with the number of H 3 O + equal to the number of negatively charged polyelectrolyte monomers. We then run the reactive canonical MC algorithm described above to obtain the equilibrium number of protonated groups and the number of free hydronium ions, see Fig. . Note that in a canonical simulation we do not have a direct access to the pH, which will be determined by the activity of hydronium ions in equilibrium. This can only be obtained using a separate Widom like particle insertion simulation that will allow us to probe the electrochemical potential of hydronium ions after the equilibrium has been established. To change pH inside the system, we can add a base such as NaOH. Since the spontaneous hydrolysis of water is so weak, addition of 1 base molecule will result in the formation of one water molecule, and an appearance of a Na + ion inside the simulation box. The net effect is, therefore, a replacement 1H + → 1Na + . We can then rerun the simulation to obtain the new protonation state. Repeating this process until all hydroniums are replaced by Na + , we can cover the full pH range.