Protein internal ionizable groups can exhibit large shifts in pK(a) values. Although the environment and interaction changes have been extensively studied both experimentally and computationally, direct calculation of pK(a) values of these internal ionizable groups in explicit water is challenging due to energy barriers in solvent interaction and in conformational transition. The virtual mixture of multiple states (VMMS) method is a new approach designed to study chemical state equilibrium. This method constructs a virtual mixture of multiple chemical states in order to sample the conformational space of all states simultaneously and to avoid crossing energy barriers related to state transition. By applying VMMS to 25 variants of staphylococcal nuclease with lysine residues at internal positions, we obtained the pK(a) values of these lysine residues and investigated the physics underlining the pK(a) shifts. Our calculation results agree reasonably well with experimental measurements, validating the VMMS method for pK(a) calculation and providing molecular details of the protonation equilibrium for protein internal ionizable groups. Based on our analyses of protein conformation relaxation, lysine side chain flexibility, water penetration, and the micro environment, we conclude that the hydrophobicity of the microenvironment around the lysine side chain (which affects water penetration differently for different protonation states) plays an important role in the pK(a) shifts.