Ionizable groups of proteins play major roles in almost all biological processes, including catalysis, proton transport and folding/misfolding. Understanding these roles requires evaluating both interactions of the ionized groups and the energetics of the ionization process itself. Particularly for enzyme-catalyzed reactions, the pKa values of specific residues and their catalytically competent protonation states are of primary importance as they are the key for their adapted role in catalysis. I will present our recent work inspired by our long-term goal of developing efficient combined quantum mechanics/molecular mechanics (QM/MM) models to understand how carbohydrate-active enzymes work.1,2 As the stringent tests for the developed QM/MM methods, pKa shift calculations are carried out based on the thermodynamic integration approach for titritable groups along the catalytic process. I will focus on the characterization of the catalytically competent protonation states of catalytic residues in selected glycoside hydrolases including sialidase/neuraminidase and β-N-acetylglucosaminidases3 . The combined QM/MM simulations rationalize the molecular origins of the substantial pKa shifts for the key residues during catalysis and highlight the importance of the electrostatic interactions and desolvation effects induced by substrate binding and glycosylation. The simulations also provide the opportunity to characterize the important intermediate states and offer structural and electronic information of the species along the reaction path, which is of great value for rational design of novel inhibitors as potential therapeutics.