The evolution of Bose-Einstein condensates is amply described by the time-dependent GrossPitaevskii mean-field theory which assumes all bosons to reside in a single time-dependent oneparticle state throughout the propagation process. In this work, we go beyond mean-field and develop an essentially-exact many-body theory for the propagation of the time-dependent Schrödinger equation of N interacting identical bosons. In our theory, the time-dependent many-boson wavefunction is written as a sum of permanents assembled from orthogonal one-particle functions, or orbitals, where both the expansion coefficients and the permanents (orbitals) themselves are timedependent and fully determined according to a standard time-dependent variational principle. By employing either the usual Lagrangian formulation or the Dirac-Frenkel variational principle we arrive at two sets of coupled equations-of-motion, one for the orbitals and one for the expansion coefficients. The first set comprises of first-order differential equations in time and non-linear integro-differential equations in position space, whereas the second set consists of first-order differential equations with time-dependent coefficients. We call our theory multi-configurational timedependent Hartree for bosons, or MCTDHB(M ), where M specifies the number of time-dependent orbitals used to construct the permanents. Numerical implementation of the theory is reported and illustrative numerical examples of many-body dynamics of trapped Bose-Einstein condensates are provided and discussed.