We formulate a time-dependent density functional theory for the coupled dynamics of electrons and nuclei that goes beyond the Born-Oppenheimer (BO) approximation. We prove that the time-dependent marginal nuclear probability density |(R,t)|^2, the conditional electronic density n_R(r,t), and the current density J_R(r,t) are sufficient to uniquely determine the full time-evolving electron-nuclear wave function, and thus the dynamics of all observables. Moreover, we propose two time-dependent Kohn-Sham schemes which reproduce the exact conditional electronic density and current density and the exact N-body nuclear density. The remaining task is to look for functional approximations for the Kohn-Sham exchange-correlation scalar and vector potentials. Using a model driven proton transfer system, we numerically demonstrate that the adiabatic extension of a beyond-BO ground state functional captures the dominant nonadiabatic effects in the regime of slow driving.