  We present the latest development of the disk gravitational instability and
fragmentation model, originally introduced by us to explain episodic accretion
bursts in the early stages of star formation. Using our numerical hydrodynamics
model with improved disk thermal balance and star-disk interaction, we computed
the evolution of protostellar disks formed from the gravitational collapse of
prestellar cores. In agreement with our previous studies, we find that cores of
higher initial mass and angular momentum produce disks that are more favorable
to gravitational instability and fragmentation, while a higher background
irradiation and magnetic fields moderate the disk tendency to fragment. The
protostellar accretion in our models is time-variable, thanks to the nonlinear
interaction between different spiral modes in the gravitationally unstable
disk, and can undergo episodic bursts when fragments migrate onto the star
owing to the gravitational interaction with other fragments or spiral arms.
Most bursts occur in the partly embedded Class I phase, with a smaller fraction
taking place in the deeply embedded Class 0 phase and a few possible bursts in
the optically visible Class II phase. The average burst duration and mean
luminosity are found to be in good agreement with those inferred from
observations of FU-Orionis-type eruptions. The model predicts the existence of
two types of bursts: the isolated ones, showing well-defined luminosity peaks
separated with prolonged periods (~ 10^4 yr) of quiescent accretion, and
clustered ones, demonstrating several bursts occurring one after another during
just a few hundred years. Finally, we estimate that 40\%--70\% of the
star-forming cores can display bursts after forming a star-disk system.
