Systems biology describes cellular phenotypes as properties that emerge from the complex interactions of individual system components. Little is known about how these interactions have affected the evolution of metabolic enzymes. Here, we combine genome-scale metabolic modeling with population genetics models to simulate the evolution of enzyme turnover numbers (kcats) from a theoretical ancestor with inefficient enzymes. This systems view of biochemical evolution reveals strong epistatic interactions between metabolic genes that shape evolutionary trajectories and influence the magnitude of evolved kcats. Diminishing returns epistasis prevents enzymes from developing higher kcats in all reactions and keeps the organism far from the potential fitness optimum. Multifunctional enzymes cause synergistic epistasis that slows down adaptation. The resulting fitness landscape allows kcat evolution to be convergent. Predicted kcat parameters show a significant correlation with experimental data, validating our modeling approach. Our analysis reveals how evolutionary forces shape modern kcatS and the whole of metabolism.