Abstract Multidimensional fitness landscapes provide insights into the molecular basis of laboratory and natural evolution. Yet such efforts are rare and focus only on limited protein families and a single enzyme trait, with little concern about the relationship between protein epistasis and conformational dynamics. Here, we report the first multiparametric fitness landscape for a cytochrome P450 monooxygenase that was engineered for the regio- and stereoselective hydroxylation of a steroid. We developed a computational program to automatically quantify non-additive effects among all possible mutational pathways, finding pervasive cooperative sign and magnitude epistasis on multiple catalytic traits. By using quantum mechanics and molecular dynamics simulations, we show that these effects are modulated by long-range interactions in loops, helices and beta-strands that gate the substrate access channel allowing for optimal catalysis. Our work highlights the importance of conformational dynamics on epistasis in an enzyme involved in secondary metabolism and offers lessons for engineering P450s.