C–H activation in alkanes poses a major challenge in chemistry due to the inert character of this bond, manifesting the necessity of improved catalysts. Although various metal–oxo complexes are known to facilitate alkane hydroxylation, probing the mechanistic nature of the reaction is difficult due to the extremely fast rebound of the radical intermediate in the postulated oxygen-rebound pathway. Automated reaction mechanism discovery methods, such as the artificial force induced reaction (AFIR) method, enable the efficient exploration of both expected and unexpected reaction pathways, revealing the reaction mechanism. Here, we employed this approach combined with density-functional theory (DFT) to investigate the structure and reactivity of heterometallic cubane complexes similar to the oxygen-evolving complex of photosystem II. For a series of M(O)Co3O4 cubanes, where M(O) is a terminal oxo with M = Cr, Mn, Fe, Mo, Tc, Ru and Rh, we first computed the stability of the possible spin states and the radical (i.e., oxyl) character of the M(O) moiety as a measure of their potential activity in the catalytic hydroxylation of alkanes. DFT calculations on these reactions promoted by Ru(O)Co3O4 and Fe(O)Co3O4 suggest that the latter promotes the hydroxylation of methane with a rate-determining H-abstraction barrier of 24.6 kcal mol-1. The moderate height of this barrier, together with the low cost and low toxicity of iron and cobalt, suggest that the Fe(O)Co3O4 cubane is a promising candidate for the catalytic oxidation of methane to methanol. AFIR calculations showed that the oxygen-rebound pathway yields the lowest-energy profile, thus validating this mechanism for the hydroxylation of alkanes by heterometallic cubanes. Furthermore, unexpected intermediates in which the methyl radical couples with either the metal center or the bridging-oxo ligands were also observed.