Simulating the non-perturbative and non-Markovian dynamics of open quantum systems is a very challenging many body problem, due to the need to evolve both the system and its environments on an equal footing. Tensor network and matrix product states (MPS) have emerged as powerful tools for open system models, but the numerical resources required to treat finite-temperature environments grow extremely rapidly and limit their applications. In this study we use time-dependent variational evolution of MPS to expore the striking theory of Tamescelli et al. [1] that shows how finite-temperature open dyanmics can be obtained from zero temperature, i.e. pure wave function, simulations. Using this approach, we produce a benchmark data set for the dynamics of the Ohmic spin-boson model across a wide range of coupling and temperatures, and also present detailed analysis of the numerical costs of simulating non-equilibrium steady states, such as those emerging from the non-perturbative coupling of a qubit to baths at different temperatures. Despite ever growing resource requirements, we find that converged non-perturbative results can be obtained, and we discuss a number of recent ideas and numerical techniques that should allow wide application of MPS to complex open quantum systems.