The efficiency of molecular dynamics is limited by the time step that can be used to integrate equations of motion, which is dictated by the highest frequency motion in the system. Multiple time step (MTS) integrators alleviate this issue by decomposing the forces acting on the particles into “fast” and “slow” components, which can then be integrated using different time steps. In ab initio MTS, an inexpensive, low level electronic structure method can be used to integrate the fast components, while its difference with an expensive but accurate high level method is used for the slow components. In this work, we present a machine learning-enhanced multiple time step (ML-MTS) method for performing accurate Born–Oppenheimer molecular dynamics at significantly reduced computational cost. We propose two alternative ML-MTS schemes, which invoke different timescale separations and result in stable and accurate trajectories. In the first scheme, ML force estimates bypass the need for a high level calculation, resulting in speedups of two orders of magnitude over standard velocity Verlet (VV) integration using a hybrid exchange–correlation functional. In the second scheme, we keep the high level calculation for the slow component and an ML correction is applied to the fast component, allowing a fourfold increase in time step compared to modern ab initio MTS algorithms without any loss of stability, thus yielding speedups up to almost an order of magnitude over straightforward VV.