Redox reactions are crucial to biological processes that protect organisms against oxidative stress. Metalloenzymes, such as cytochrome c peroxidases which reduce excess hydrogen peroxide into water in the periplasm of multiple bacterial organisms, play a key role in detoxification mechanisms. While accurate computational tools can be used to simulate ground state redox potentials in biomolecules, adapting such approaches to properly describe redox reactions in transition metal complexes, particularly in hemes in heterogeneous protein environments, remains a significant challenge. Here we present the results of polarizable hybrid QM/MM studies of the reduction potentials of two heme sites in the cytochrome c peroxidase of Nitrosomonas europaea. The simulated redox potential of the catalytic site Low Potential (LP) is in good agreement with the experiment, while for the High Potential (HP) heme the computational estimate significantly overestimate the experimental value. We have found that environment polarization shifts the computed value of the redox potential of the catalytic LP heme by 1.3 V, while it does not affect that of the non-catalytic High Potential (HP) heme. We demonstrate that it is necessary to account for mutual polarization of heme site and the protein environment when describing redox processes, particularly those that involve more charged heme sites. We have explored the role of various factors such as heme geometries, axial ligands, propionate side chains, and electrostatic field of the protein in tuning the redox potentials of hemes in NeCcP. The fluctuations in computed vertical ionization and electron attachment energies are predominantly affected by fluctuations in the electrostatic field of the environment but not by fluctuations in heme geometries. We attribute the difference in computed LP and HP heme reduction potentials of 0.05 V and 1.15 V, respectively, to different axial ligands and electrostatic interactions of the hemes with the protein environment. / 2023-06-30T00:00:00Z