We investigate the possibility of extracting the probability distribution of the effective environmental tuning and coupling modes during the nonadiabatic relaxation through a conical intersection. Dynamics are dealt with an open quantum system master equation by partitioning a multistate electronic subsystem out of all the nuclear vibrators. This is an alternative to the more usual partition retaining the tuning and coupling modes of a conical intersection in the active subsystem coupled to a residual bath. The minimal partition of the electronic system generally leads to highly structured spectral densities for both vibrational baths and requires a strongly nonperturbative non-Markovian master equation, treated here by the hierarchical equations of motion (HEOMs). We extend-for a two-bath situation-the procedure proposed by Shi et al. [J. Chem. Phys. 140, 134106 (2014)], whereby the information contained in the auxiliary HEOM matrices is exploited in order to derive the nuclear dissipative wave packet, i.e., the statistical distribution of the displacement of the two tuning and coupling collective coordinates in each electronic state and the coherence. This allows us to visualize the distribution, all along the nonadiabatic decay. We explore a large parameter space for a symmetrical conical intersection model and a symmetrical initial Franck-Condon preparation. Some parameters could be controlled by external fields, while others are molecule dependent and could be designed by molecular engineering. We illustrate the relation between the strongly coupled electronic and bath dynamics together with a geometric measure of non-Markovianity.