The genetic and molecular approaches to heterosis usually do not rely on any model of the genotype-phenotype relationship. From the generalization of Kacser and Burns' biochemical model for dominance and epistasis to networks with several variable enzymes, we hypothesized that metabolic heterosis could be observed because the response of the flux towards enzyme activities and/or concentrations follows a multi-dimensional hyperbolic-like relationship. To corroborate this, we used the values of systemic parameters accounting for the kinetic behaviour of four enzymes of the upstream part of glycolysis, and simulated genetic variability by varying in silico enzyme concentrations. Then we "crossed" virtual parents to get 1,000 hybrids, and showed that best-parent heterosis was frequently observed. The decomposition of the flux value into genetic effects, with the help of a novel multilocus epistasis index, revealed that antagonistic additive-by-additive epistasis effects play the major role in this framework of the genotype-phenotype relationship. This result is consistent with various observations in quantitative and evolutionary genetics, and provides a model unifying the genetic effects underlying heterosis.