Abstract Motivation Molecular profiling techniques have evolved to single-cell assays, where dense molecular profiles are screened simultaneously for each cell in a population. High-throughput single-cell experiments from a heterogeneous population of cells can be experimentally and computationally sorted as a sequence of samples pseudo-temporally ordered samples. The analysis of these datasets, comprising a large number of samples, has the potential to uncover the dynamics of the underlying regulatory programmes. Results We present a novel approach for modelling and inferring gene regulatory networks from high-throughput time series and pseudo-temporally sorted single-cell data. Our method is based on a first-order autoregressive moving-average model and it infers the gene regulatory network within a variational Bayesian framework. We validate our method with synthetic data and we apply it to single cell qPCR and RNA-Seq data for mouse embryonic cells and hematopoietic cells in zebra fish. Availability and implementation The method presented in this article is available at https://github.com/mscastillo/GRNVBEM.