  We extend the twin-space formulation of the hierarchical equations of motion
approach in combination with the matrix product state representation
(introduced in J. Chem. Phys. 150, 234102, [2019]) to nonequilibrium scenarios
where the open quantum system is coupled to a hybrid fermionic and bosonic
environment. The key ideas used in the extension are a reformulation of the
hierarchical equations of motion for the auxiliary density matrices into a
time-dependent Schr\"odinger-like equation for an augmented multi-dimensional
wave function as well as a tensor decomposition into a product of low-rank
matrices. The new approach facilitates accurate simulations of non-equilibrium
quantum dynamics in larger and more complex open quantum systems. The
performance of the method is demonstrated for a model of a molecular junction
exhibiting current-induced mode-selective vibrational excitation.
