A coarse-grained model of linear polyfunctional weak charged biopolymers was implemented, formed of different proportions of acid-base groups resembling the composition of humic substances. These substances are mainly present in dissolved organic matter in natural water. The influence of electrostatic interactions computing methods, factors concerning the structure of the chain, different functional groups, and the ionic strength on polyelectrolytes were studied. Langevin dynamics with constant pH simulations were performed using the ESPResSO package and the Python-based Molecule Builder for ESPResSo (pyMBE) library. The coverage was fitted to a polyfunctional Frumkin isotherm, with a mean-field interaction between charged beads. The composition of the chain affects the charge while ionic strength affects both the charge and the radius of gyration. Additionally, the parameters intrinsic to the polyelectrolyte model were well reproduced by fitting the polyfunctional Frumkin isotherm. In contrast, the non-intrinsic parameters depended on the ionic strength. The method developed and applied to a polyfunctional polypeptide model, that resembles a humic acid, will be very useful for characterizing biopolymers with several acid-base functional groups, where their structure, the composition of the different functional groups, and the determination of the main intrinsic proton binding constants and their proportion are not exactly known.