The present article discusses some aspects concerning the inclusion of proton isomerism in simulations of the global protonation equilibrium of protein molecules. In the context of continuum electrostatic methods, the usual basis for these simulations, this isomerism can be treated as a coexistence of tautomeric forms in equilibrium in a rigid structure; furthermore, it can be formally extended to nontitrable sites with proton isomerism, such as alcohol groups and water molecules. We follow the previously adopted approach of transforming the real system of tautomeric sites into a thermodynamically equivalent one of nontautomeric pseudosites, establishing a proper relation between the two systems. The necessary energetic and entropic modifications of model compound pK(a) values are also discussed. Additionally, we discuss the new entropy term, named tautomeric entropy, that results from the explicit inclusion of tautomerism in the simulations and how it can be computed together with the occupational entropy. Simulations using tautomerism were done for hen egg white lysozyme (HEWL) using a simple set of tautomers at dihedral energy minima. A very good overall prediction of pK(a) values was obtained, presumably the best in the literature for HEWL, using a high value for the dielectric constant assigned to the protein region, is an element of (p). The explicit inclusion of water molecules treated under the extended tautomer formalism further improved the prediction, in contrast with previous works using rigid water molecules. In all calculations performed, the region with is an element of (p) approximate to 20 is shown the to be the optimal one. Some aspects of the somewhat controversial issue of the "proper" is an element of (p), value are also discussed.