In this paper, a model reduction procedure is proposed for the simplification of biochemical reaction network models. The approach is capable of reducing ODE models where the right hand side of the equations contains polynomial and/or rational function terms. The method is based on a finite number of mixed integer quadratic programming (MIQP) steps where the objective function effectively measures the fit between the time functions of the selected concentrations of the original and the reduced models, and the integer variables keep track of the presence of individual reactions. The procedure also contains the re-estimation of rate coefficients in the reduced model to minimize the defined model error. Two examples taken from the literature illustrate the operation of the method.