The importance of accounting for residue coupling is multifaceted and particularly relevant in the context of enzymatic active sites that are often enriched in protonatable residues. The theoretical formalism to describe such couplings in polyprotic acids has been detailed by Ullman. Based on earlier work by Edsall and co-workers, Ullman defines equilibrium protonation constants for all microstates and derives partition functions that fully describe the thermodynamics of these systems. This approach was then applied in a protein context using continuum electrostatics calculations to resolve the free energies between microstates. Here, we demonstrate that such a framework can be readily extended to double free energy difference calculations based on molecular dynamics simulations and introduce a thermodynamic cycle approach that allows one to resolve the apparent pK a s of coupled residues. This demonstration and extension is particularly relevant in the context of alchemical free energy calculations, which give access only to such ∆∆G values.