The distribution of gene expression across cells is controlled by gene regulatory networks (GRN) that determine cell type identity, direct developmental trajectories, and respond to perturbations. While many methods exist for detecting changes in mean expression from single cell RNA-sequencing data, detecting changes in the variability and covariability of gene expression between groups of cells remain challenging. Here we introduce memento, a method that implements a generative model for estimating the mean, residual variance, and correlation of gene expression and a resampling strategy to detect differences in these moment parameters scalable to millions of cells and hundreds of samples. We apply memento to study the interferon GRN in human tracheal epithelial cells (HTEC) and peripheral blood mononuclear cells (PBMC). The analysis of expression variability from resting cells revealed canonical interferon stimulated genes (ISGs) to be amongst the most variable indicative of tonic interferon signaling in both HTECs and PBMCs. Upon recombinant interferon stimulation, cells exhibit decreased variability of ISGs that is associated with the topology of the GRN and the variability of canonical regulators. The analysis of gene correlations confirmed that the canonical interferon GRN is already active in resting cells. In response to recombinant interferon, upregulation of co-activators is associated with expression of non-canonical ISGs enriched for antigen processing and presentation genes. memento enables quantitative comparisons of gene expression distributions from single-cell RNA-sequencing data to reveal distinct modes of regulation for canonical and non-canonical ISGs. It can be applied to other large-scale datasets to detect gene expression changes imparted by disease states, perturbations, or genetic variants.