Ions are involved in multiple biological processes and may exist bound to biomolecules or may be associated with their surface. Although the presence of ions in nucleic acids has traditionally gained more interest, ion-protein interactions, often with a marked dependency on pH, are beginning to gather attention. Here we present a detailed analysis on the binding and distribution of ions around β-lactoglobulin using a constant-pH MD (CpHMD) method, at a pH range 3-8, and compare it with the more traditional Poisson-Boltzmann (PB) model and the existing experimental data. Most analyses used ion concentration maps built around the protein, obtained from either the CpHMD simulations or PB calculations. The requirements of approximate charge neutrality and ionic strength equal to bulk, imposed on the MD box, imply that the absolute value of the ion excess should be half the protein charge, which is in agreement with experimental observation on other proteins (Proc. Natl. Acad. Sci. U.S.A. 2021, 118, e2015879118) and lends support to this protocol. In addition, the protein total charge (including territorially bound ions) estimated with MD is in excellent agreement with electrophoretic measurements. Overall, the CpHMD simulations show good agreement with the nonlinear form of the PB (NLPB) model but not with its linear form, which involves a theoretical inconsistency in the calculation of the concentration maps. In several analyses, the observed pH-dependent trends for the counterions and co-ions are those generally expected, and the ion concentration maps correctly converge to the bulk ionic strength as one moves away from the protein. Despite the overall similarity, the CpHMD and NLPB approaches show some discrepancies when analyzed in more detail, which may be related to an apparent overestimation of counterion excess and underestimation of co-ion exclusion by the NLPB model, particularly at short distances from the protein.