Stochastic models that describe interacting processes, such as stochastic automata networks, feature a dimensionality that grows exponentially with the number of processes. This state space explosion severely impairs the use of standard methods for the numerical analysis of such Markov chains. In this work, we discuss the approximation of solutions by matrix product states or, equivalently, by tensor train decompositions. Two classes of algorithms based on this low-rank decomposition are proposed, using either iterative truncation or alternating optimization. Our approach significantly extends existing approaches based on product form solutions and can, in principle, attain arbitrarily high accuracy. Numerical experiments demonstrate that the newly proposed algorithms are particularly well suited to deal with pairwise neighbor interactions.