While from the perspective of statistical mechanics the behaviour of coupled titratable sites is well understood, the subsequent challenge arises: What is the best way to extract information about changes in residue protonation using existing computational approaches? Molecular dynamics (MD) simulations offer an attractive solution to this question by providing efficient sampling of well defined thermodynamic ensembles. Over the years, several approaches have emerged to quantify pK a shifts based on MD simulations, among the most popular are constant-pH simulations (CpHMD) and alchemical free energy calculations. While CpHMD explicitly models protonation changes dynamically over the course of a simulation, an alchemical approach considers discrete protonation states and computes the free energies between them (see Appendix A for additional details). Here, we focus on the latter approach and demonstrate how to connect double free energy differences with a statistical mechanics description of population changes in coupled residue protonations.