Recent work has highlighted the scale and ubiquity of subject variability in observations from functional MRI data (fMRI). Furthermore, it is highly likely that errors in the estimation of either the spatial presentation of, or the coupling between, functional regions can confound cross-subject analyses, making accurate and unbiased representations of functional data essential for interpreting any downstream analyses. Here, we extend the framework of probabilistic functional modes (PFMs) (Harrison et al., 2015) to capture cross-subject variability not only in the mode spatial maps, but also in the functional coupling between modes and in mode amplitudes. A new implementation of the inference now also allows for the analysis of modern, large-scale data sets, and the combined inference and analysis package, PROFUMO, is available from git.fmrib.ox.ac.uk/samh/profumo. A new implementation of the inference now also allows for the analysis of modern, large-scale data sets. Using simulated data, resting-state data from 1000 subjects collected as part of the Human Connectome Project (Van Essen et al., 2013), and an analysis of 14 subjects in a variety of continuous task-states (Kieliba et al., 2019), we demonstrate how PFMs are able to capture, within a single model, a rich description of how the spatio-temporal structure of resting-state fMRI activity varies across subjects. We also compare the new PFM model to the well established independent component analysis with dual regression (ICA-DR) pipeline. This reveals that, under PFM assumptions, much more of the (behaviorally relevant) cross-subject variability in fMRI activity should be attributed to the variability in spatial maps, and that after accounting for this functional coupling between modes primarily reflects current cognitive state. This has fundamental implications for the interpretation of cross-sectional studies of functional connectivity that do not capture cross-subject variability to the same extent as PFMs.