A new method is presented for simulating the simultaneous binding equilibrium of electrons and protons on protein molecules, which makes it possible to study the full equilibrium thermodynamics of redox and protonation processes, including electron-proton coupling. The simulations using this method reflect directly the pH and electrostatic potential of the environment, thus providing a much closer and realistic connection with experimental parameters than do usual methods. By ignoring the full binding equilibrium, calculations usually overlook the twofold effect that binding fluctuations have on the behavior of redox proteins: first, they affect the energy of the system by creating partially occupied sites; second, they affect its entropy by introducing an additional empty/occupied site disorder (here named occupational entropy). The proposed method is applied to cytochrome c3 of Desulfovibrio vulgaris Hildenborough to study its redox properties and electron-proton coupling (redox-Bohr effect), using a continuum electrostatic method based on the linear Poisson-Boltzmann equation. Unlike previous studies using other methods, the full reduction order of the four hemes at physiological pH is successfully predicted. The sites more strongly involved in the redox-Bohr effect are identified by analysis of their titration curves/surfaces and the shifts of their midpoint redox potentials and pKa values. Site-site couplings are analyzed using statistical correlations, a method much more realistic than the usual analysis based on direct interactions. The site found to be more strongly involved in the redox-Bohr effect is propionate D of heme I, in agreement with previous studies; other likely candidates are His67, the N-terminus, and propionate D of heme IV. Even though the present study is limited to equilibrium conditions, the possible role of binding fluctuations in the concerted transfer of protons and electrons under nonequilibrium conditions is also discussed. The occupational entropy contributions to midpoint redox potentials and pKa values are computed and shown to be significant.