In the context of heterogeneous catalysis, computational methods are relatively mature for characterizing the transition states of targeted reactions, performing microkinetic modeling on established reaction networks, and using descriptor-based methods for improving catalysts. However, a central challenge in characterizing new catalytic interfaces lies in establishing the kinetically relevant reaction network, which is often based on intuition and can be time-consuming and error prone to characterize ad hoc. Indeed, even seemingly simple heterogeneous reactions, like methane activation on metal oxide surfaces, can be decomposed into multiple elementary steps. Furthermore, catalytic cycles can involve many intermediates or even open ended reactant lists such that brute force enumeration and characterization are infeasible. Such examples include the oxidative coupling of methane and olefin oligomerization, each involving the forma-tion/dissociation of long carbon backbones as intermediates and an open set of olefins as potentially adsorbed reactants. For these reasons, the recent advent of automated reaction prediction approaches is potentially promising for elucidating reaction networks involving heterogeneous interfaces. These methods can be categorized on the basis of whether the potential energy surface (PES) is explored in detail to locate transition states or whether the reaction networks are enumerated using a closed set of reaction templates. The latter class includes packages like Network Generation (NetGen), and Reaction Mechanism Generator (RMG); however, due to the reliance on established reaction templates, this strategy is less relevant to characterizing exploratory catalysts where such data is typically absent. In contrast, methods that directly explore the PES circumvent this limitation, at least in principle. This class includes several approaches that are under active development, including the artificial force induced reaction (AFIR) method, stochastic surface walking reaction sampling (SSW), the ZStucture method from the Zimmerman group, and Yet Another Reaction Program (YARP), our recently developed methodology. All of these approaches are intrinsically more expensive than template-based methods because they sample the PES (e.g., using quantum chemistry calculations), which has been a major obstacle to applying them to inhomogeneous systems in an exploratory context. For example, both SSW and AFIR have been applied to successfully (re)discover the relatively simple heterogeneous water-gas shift reaction occurring at a copper surface. Nevertheless, this required millions of density functional theory (DFT) gradient calls, despite the small reactive system sizes. As the number of elementary steps grows and the complexity increases, a highly efficient reaction exploration scheme becomes even more indispensable to mitigate computational costs. Moreover, heterogeneous applications have several other technical obstacles to applying automated approaches that were designed for molecular systems. These include the larger system sizes that are typical of surface models; the occurrence of spectator atoms that do not participate in reactions but nevertheless play important non-covalent or structural roles in the reaction pathways; and the use of periodic versus molecular models of the reacting systems. The optimal manner of addressing these obstacles are all outstanding research questions.