To ensure a reliable operation of the rf systems in charged particle accelerators, understanding their interaction with the beam is crucial. With increasing beam intensities, beam loading becomes more important and can cause beam quality degradation. In addition, high intensities can bring rf amplifiers based on vacuum tubes to operate in a nonlinear regime with high anode current, which affects the beam loading compensation. This contribution provides a methodology for the mathematical modeling of the response of tetrode vacuum tube amplifiers in rf systems, which can be coupled with a particle tracking code to take the dynamic beam behavior into account. An example of a numerical implementation is shown for the case of the 80 MHz rf systems in the CERN proton synchrotron. The implementation was validated with low- and high-power measurements. Detailed studies have been performed in order to verify that the variation of the tube operating regime with beam intensity and cavity voltage was correctly reproduced. Finally, a beam instability driven by the impedance of the high-level rf systems has been observed experimentally and reproduced in simulations, demonstrating that the model captures the dynamic behavior of both, the beam and the rf systems.