where ze is the charge and ψ is the mean electrostatic potential. Thus the partitioning of charge-regulating solutes is affected by the difference in ψ between coacervate and the supernatant in two different ways: (1) at fixed valency of the solute, it determines the electrostatic contribution to the excess chemical potential; (2) at pH which is not too far from pK A of the solute, valency of the solute can be switched as it is exchanged between the phases. While the first effect has been considered in the existing theories, the second one has not yet been analyzed theoretically, although there is clear experimental evidence that a change in pH dramatically affects solute partitioning. To bridge this knowledge gap between theory and experiment, various molecular models have been designed for use in molecular simulations, which have become a strong tool for the polymer physics. The prevalent approaches in particle-based simulations of polyelectrolyte complexes can be classified in three categories. First, simulations exploring only a single pair of oppositely charged chains in dilute conditions. These studies correctly identified the main thermodynamic driving forces of polyelectrolyte complexation but they provided very little information about the behaviour of the complexes in the bulk. Second, simulations using many chains in a large simulation cell, directly emulating the phase separation. In these simulations, a nanodroplet of polyelectrolyte complex is formed inside an otherwise empty box. Such simulations suffer from massive finite-size effects, caused by the explicit interface of the coacervate nanodroplet. Extrapolation of these results to the macroscopic systems is problematic not only because the interface contributes much less to the bulk properties but also because the bulk of the coacervate phase is poorly defined in such simulation. Third group of the simulations is based on the Gibbs ensemble, where one simulates coexistence of two bulk phases (PEC and supernatant) using two simulation cells, avoiding the need for simulating their interface. 102 Such a simulation typically assumes that there is no polymer in the supernatant phase, while simultaneously small ions and other solutes can be exchanged between both phases. This third group of simulations is well suited for predicting the phase stability of coacervates and partitioining of solutes between the phases. It avoids the drawbacks of the previous two, therefore, the only limitation of this group consists in the quality of the molecular model. Indeed, the simulations of the third group qualitatively reproduced the experimentally observed phase diagrams as a function of salt concentration and also the partitioning of monovalent ions. In addition, by tuning the model parameters, a quantitative agreement between such simulations and specific experiments could be achieved. Our method falls within this third category, however in contrast to the previous studies, it also includes the effects of pH and charge regulation. These effects have not yet been addressed in molecular simulations although they have been addressed by phenomenological thermodynamic models. Unlike preceding simulation studies, we have not designed our model to reproduce the properties of a specific experimental system. Instead, our model of polymers and solutes are constructed as a generic beadspring model, possessing only a small number of parameters which can be linked with solvent permittivity, size and valency of the monomers and of the small ions. Such a setting allows us to systematically explore the effect of these parameters on the phase stability of the polyelectrolyte complex and on the partitioning of small ions. Only in simulations these parameters can be varied independently, while keeping all other interactions constant, thereby allowing us to distinguish generic features of the phase separation and partitioning from those which are determined by the specific chemistry of the polymers and small ions.