The boundary element method (BEM) in combination with continuum electrostatics is employed to compute the acid-dissociation constants K-a (pK(a)= -log K-a) of titrating sites in proteins. The boundary element method determines the electrostatic potential from the solution of two coupled integral equations, valid on a triangulated surface enclosing the macromolecule. First the intrinsic pK(2)(intr) is computed from the shift of the pK(a) of a model compound containing the titrating site upon its transfer from solution to the protein; during this step interactions between titrating sites are neglected. Subsequently, interactions between titrating sites are included by means of a Monte Carlo scheme to sample protonation states of the protein. A convenient vector-matrix formulation in terms of the BBM is given which allows the use of a single atom or a detailed charge model (or a combination of both) to describe the titrating sites. The method has been applied to four proteins : bovine pancreatic trypsin inhibitor, calbindin, lysozyme, and ovomucoid third-domain. Different choices for the dielectric constant of proteins ranging from 4 to 78.5 were investigated in a systematic fashion. Comparisons are made with pK(a) values calculated by the finite difference method and those determined experimentally for these four proteins. Our results indicate that accurate pK(a) values are obtained with the BEM when a dielectric constant for the protein of 20 or higher is used, For calbindin, different choices for the ionic strength were considered and comparison was made with pK(a) values obtained experimentally and from a simulation model using explicit ions for the solvent.