Abstract The molecular equilibria involved in the second and third macroscopic deprotonation processes of Cys have been theoretically characterized at B3LYP/6-31++G** computing level. The role of solvent was analyzed by using the supermolecule (up to six discrete water molecules), the continuum, and the hybrid supermolecule-continuum models. Also, a novel approach to reveal the solvation effect of the bulk water was employed. Calculations performed with the supermolecule or the Onsager models overestimate absolute p K as, whereas the PCM continuum model yields data much closer to the experimental values. The supermolecule-PCM approach estimates the p K a values for the amino group much better than for the thiol group.