Unraveling of the density matrix of a non-stationary quantum state as a function of a few observables provides a complementary view of quantum dynamics. Here we pay special attention to the entangled electron-nuclear motion. We have recently developed a practical way to identify the minimal set of the dominant observables that govern the quantal dynamics even in the case of strong non-adiabatic effects and large anharmonicity [Komarova et al., J. Chem. Phys. 155, 204110 (2021)]. Fast convergence in the number of the dominant contributions is achieved when instead of the density matrix we describe the time-evolution of the surprisal, the logarithm of the density operator. In the present work, we illustrate the efficiency of the proposed approach using an example of the early time dynamics in pyrazine in a Hilbert space accounting for up to four vibrational normal modes, {Q10a, Q6a, Q1 , Q9a}, and two coupled electronic states, the optically dark and the bright states. Dynamics in four-dimensional configurational space involve 19600 vibronic eigenstates. Our results reveal that the rate of the ultrafast population decay as well as the shape of the nuclear wave packets in 2D, accounting only for {Q10a, Q6a} normal modes, are accurately captured with only six dominant time-independent observables in the surprisal. Extension of the dynamics to three- and four-dimensional vibrational subspace requires only five additional constraints. The time-evolution of a quantum state in 4D vibrational space on two electronic states is thus compacted to only 11 time-dependent coefficients of these observables.