The present study deals with the mathematical modeling of crosslinking kinetics of polymer–phenol conjugates mediated by the Horseradish Peroxidase (HRP)-hydrogen peroxide (H2O2) initiation system. More specifically, a dynamic Monte Carlo (MC) kinetic model is developed to quantify the effects of crosslinking conditions (i.e., polymer concentration, degree of phenol substitution and HRP and H2O2 concentrations) on the gelation onset time; evolution of molecular weight distribution and number and weight average molecular weights of the crosslinkable polymer chains and gel fraction. It is shown that the MC kinetic model can faithfully describe the crosslinking kinetics of a finite sample of crosslinkable polymer chains with time, providing detailed molecular information for the crosslinkable system before and after the gelation point. The MC model is validated using experimental measurements on the crosslinking of a tyramine modified Hyaluronic Acid (HA-Tyr) polymer solution reported in the literature. Based on the rubber elasticity theory and the MC results, the dynamic evolution of hydrogel viscoelastic and molecular properties (i.e., number average molecular weight between crosslinks, Mc, and hydrogel mesh size, ξ) are calculated.