In this thesis, we have presented and investigated the possible strategies for modelling complex heterogeneous catalytic systems in operando regimes. By introducing modern computational approaches to sample the potential energy surface of the catalytic active site, we have attempted to account for the reactive conditions, solvent presence, additives inclusion, and structural dynamics of the active site. There is growing spectroscopic and theoretical evidence of the critical role of the active site dynamics for the catalytic performance, advocating for the active site representation as an ensemble of possible isomers. Challenged by the complexity of the reactive environment and common heterogeneous catalysts, we strongly believe that addressing these factors and incorporating them explicitly into the model description will contribute to a more realistic representation of the catalytic system.