Abstract Galacto-oligosaccharides (GOS) are highly valued prebiotics in the food industry with many significant health benefits. The enzymatic transgalactosylation of lactose to produce GOS, catalyzed by β-galactosidase, is a consolidated route in industry. Kinetic models of this reaction network presented in the literature usually adopt broad assumptions by considering only the formation of GOS with a glucose residue in the reaction, neglecting the presence of galactobiose (Glb), which is also a Galacto-oligosaccharide that has bifidogenic effects. Proposed models that take into account the presence of Glb and use lactose as a substrate do not give information about this molecule concentration during the reaction. In this paper, a kinetic model is proposed with nine adjustable parameters to characterize the formation of Glb, trisaccharides, tetrasaccharides and enzyme inactivation. A bootstrap method, based on the residuals obtained during the fitting procedure, is employed to calculate parameters’ confidence intervals. This model fits very well to experimental data and addresses many of the existing issues in the literature. As a result, this model may be useful for bioreactor optimization and process design, as well as process control and state estimation studies.