The stability of colloidal particles in aqueous suspensions is intrinsically connected with their surface charge density, which is controlled by the pH of solution. Similarly the activity of many biologically relevant proteins and polyelectrolytes is controlled by solution's pH and ionic strength . Quantitative understanding of charge regulation in such complex systems is, therefore, of a paramount importance in a wide range of industrial and medical applications. For some simple colloidal systems with a regular distribution of surface active groups, one can use the Poisson-Boltzmann theory with the charge regulation boundary condition to study the particle protonation state . However, this approach breaks down for suspensions containing multivalent counterions or when dealing with flexible molecules, such as proteins or polyelectrolytes, whose three dimensional conformation is intimately coupled with the protonation state of the molecule. For such systems one is forced to rely on computer simulations .