The kinetic modeling of the hydrolysis of cheese whey proteins using chymotrypsin immobilized on agarose gel particles is studied here. The classical approach to describe the progress of proteolytic reactions uses the “degree of hydrolysis” as the model input and the concentration of “hydrolysable bonds” as a pseudo-substrate. Since hydrolysis profiles usually follow Michaelis-Menten functional forms, their explicit algebraic equations are frequently applied to predict the overall rate of hydrolysis of peptide bonds. However, when the molecules of substrate are large, the enzyme-substrate molar ratio increases, and the pseudo-steady state formulation may be no longer valid. On the other hand, for immobilized enzyme diffusion delays may have to be considered. Our results indicate that a Michaelis-Menten equation, derived from the Briggs-Haldane pseudo-steady state approach, could not predict the kinetics of the reactions here considered. These models presented systematic deviations from experimental data even when free enzyme or low load chymotrypsin-agarose derivatives were used (regions where diffusive effects are negligible). Relaxing the pseudo-steady state assumptions, i.e. admitting that the dynamics of the intermediate pseudo-species will no longer be instantaneous, one obtains a seven-parameter model which was able to predict the experimental “hydrolysable bonds” profiles with remarkable accuracy. The model results in a set of stiff ordinary differential equations, which could be solved with appropriate numerical methods. In order to avoid sub-optimal sets of parameters, a hybrid search procedure was used for model fitting and parameter estimation, which combined global and local search methods: Simulated Annealing and Levenberg-Marquardt algorithms.