An approach to modeling the protonation equilibria in molecular systems with conformational flexibility, applied to pK(a) prediction in proteins, is formulated. This approach is based on the widely used concept of model system, usually identified with a free amino acid (model compound). However, analysis of the influence of conformational variability on the observed pK(a) for the free amino acids Points to the necessity of reducing the model system to only the set of atoms defining an the essentially nonflexible group (model group). A method for parametrizing the corresponding pK(a)(model), based on the hybrid Molecular Dynamics and Poisson-Boltzmann model, is presented. The results are tested on two proteins: hen egg white lysozyme and turkey ovomucoid third domain. Comparison between the predicted and experimental pK(a) values for these systems allows one to examine the problem of the magnitude of the dielectric constant to be used for the protein in this kind of modeling. The proposed methodology may be extended to any system involving proton exchange.