Multiplex networks are a common modeling framework for interconnected systems and multimodal data, yet we still lack fundamental insights for how multiplexity affects stochastic processes. We introduce a novel ``Markov chains of Markov chains'' model called multiplex Markov chains (MMCs) such that with probably $(1-ω)\in [0,1]$ random walkers remain in the same layer and follow (layer-specific) intralayer Markov chains, whereas with probability $ω$ they move to different layers following (node-specific) interlayer Markov chains. One main finding is the identification of multiplex convection, whereby a stationary distribution exhibits circulating flows that involve multiple layers. Convection cycles are well understood in fluids, but are insufficiently explored on networks. Our experiments reveal that one mechanism for convection is the existence of imbalances for the (intralayer) degrees of nodes in different layers. To gain further insight, we employ spectral perturbation theory to characterize the stationary distribution for the limits of small and large $ω$, and we show that MMCs inherently exhibit optimality for intermediate $ω$ in terms of their convergence rate and the extent of convection. As an application, we conduct an MMC-based analysis of brain-activity data, finding MMCs to differ between healthy persons and those with Alzheimer's disease. Overall, our work suggests MMCs and convection as two important new directions for network-related research.