The identification of the enzymatic profile that achieves a maximal production rate of a given metabolite is an important problem in the biotechnological industry, especially if there is a limit on the number of enzymatic modulations allowed. The intrinsic nonlinear behavior of metabolic processes enforces the use of kinetic models, such as the generalized mass action (GMA) models, giving rise to nonconvex MINLP formulations with multiple local solutions. In this paper, we introduce a customized spatial branch-and-bound strategy devised to solve efficiently these particular problems to global optimality. A tight MILP-based relaxation of the original nonconvex MINLP is constructed by means of supporting hyperplanes and piecewise linear underestimators. The overall solution procedure is expedited through the use of bound tightening techniques and a special type of cutting plane. The capabilities of the proposed strategy are tested through its application to the maximization of the citric acid production in Aspergillus niger. We also provide a numerical comparison of our algorithm with the commercial package BARON and an outer approximation-based method earlier proposed by the authors.