We extend the twin-space formulation of the hierarchical equations of motion approach in combination with the matrix product state representation [R. Borrelli, 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ödinger-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.