We discuss the multi-configurationtime-dependent Hartree (MCTDH) method for the approximation of the time-dependent Schrödinger equation in quantum molecular dynamics. This method approximates the high-dimensional nuclearwave function by a linear combination of products of functions dependingonly on a single degree of freedom. The equations of motion, obtained via the Dirac-Frenkel time-dependent variational principle, consist of a coupled system of low-dimensionalnonlinear partial differential equations and ordinary differential equations. We show that, with a smooth and bounded potential, the MCTDH equations are well-posed and retain high-order Sobolev regularity globally in time, that is,as long as the density matrices appearing in the method formulation remain invertible. In particular, the solutions are regular enough to ensure local quasi-optimality of the approximation and to admit an efficient numerical treatment.