Here we present a new implementation and extension of the stochastic titration method which makes it possible to perform MD simulations at constant pH and reduction potential. The method was applied to the redox titration of cytochrome c(3) from Desulfovibrio vulgaris Hildenborough, and a major finding of this study was that the method showed a better performance when the protein region is assigned a high dielectric constant. This dependence on the value of the protein dielectric constant was not found in previous constant-pH MD simulations and is attributed to excessively high heme-heme interactions at low dielectric constants. The simulations revealed strong coupling between hemes in close proximity, and we also showed how these couplings can be used to estimate the sensibility of the heme reductions to small pH changes.