We chose to work on the DEP-BChE system, from which the answer, once obtained, may be extended onto all other states. We performed enhanced-sampling QM/MM MD simulations, allowing protons to move freely within h-bonding networks we called "proton depots" (Figure ). All possible protonation states can be divided into groups that share the same amount of protons in both depots. Such a setup allowed us to study these states systematically. For each such group, we then modeled proton rearrangements in dynamics. We kept the positions of all protein heavy atoms close to the reference, thus biasing the whole protocol to output the probabilities of locations of protons given their number and a particular active site geometry. We note that due to these restraints, the absolute energy values in this setup are not directly interpretable, and we are interested only in relative energy differences between end states. Given the scale of the experiment, for the QM part we chose a fast semiempirical DFTB3 method that was previously shown to give reasonable results for biochemical systems. Each depot may realistically hold from 0 to 2 protons. More than 2 protons start to violate the saturation of h-bonding networks, while 0 protons could mean undersaturation. States within these 9 groups differ by exact locations of protons. To be able to model and distinguish them, we designed unique sets of collective variables (Figure , 3). For example, groups with no protons in depot 1 are all undersaturated, and the water molecule (HOH1117 in 1XLW, Fig2A) may freely choose its two interacting partners from three possible (Fig3, A-C). All relevant states were distinguishable with up to two tailored variables.