We have calculated the electrostatic potential and interaction energies of ionizable groups and analyzed the response of the protein environment to redox changes in Paracoccus denitrificans cytochrome c oxidase by using a continuum dielectric model and finite difference technique. Subsequent Monte Carlo sampling of protonation states enabled us to calculate the titration curves of all protonatable groups in the enzyme complex. Inclusion of a model membrane allowed us to restrict the calculations to the functionally essential subunits I and II. Some residues were calculated to have complex titration curves, as a result of strong electrostatic coupling, desolvation, and dipolar interactions. Around the heme a3-CuB binuclear center, we have identified a cluster of 18 strongly interacting residues that account for most of the proton uptake linked to electron transfer. This was calculated to be between 0.7 and 1.1 H+ per electron, depending on the redox transition considered. A hydroxide ion bound to CuB was determined to become protonated to form water upon transfer of the first electron to the binuclear site. The bulk of the protonation changes linked to further reduction of the heme a3-CuB center was calculated to be due to proton uptake by the interacting cluster and Glu(II-78). Upon formation of the three-electron reduced state (P1), His325, modeled in an alternative orientation away from CuB, was determined to become protonated. The agreement of these results with experiment and their relevance in the light of possible mechanisms of redox-coupled proton transfer are discussed.