Recent work has highlighted the scale and ubiquity of subject variability in observa ons from func onal MRI data (fMRI). Furthermore, it is highly likely that errors in the es ma on of either the spa al presenta on of, or the coupling between, func onal regions can confound cross-subject analyses, making accurate and unbiased representa ons of func onal data essen al for interpre ng any downstream analyses.Here, we extend the framework of probabilis c func onal modes (PFMs) [Harrison et al. ] to capture cross-subject variability not only in the mode spa al maps, but also in the func onal coupling between modes and in mode amplitudes. A new implementa on 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. Using res ng-state data from , subjects collected as part of the Human Connectome Project [Van Essen et al. ], and an analysis of subjects in a variety of con nuous task-states [Kieliba et al.], we demonstrate how PFMs are able to capture, within a single model, a rich descrip on of how the spa o-temporal structure of res ng-state fMRI ac vity 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 assump ons, much more of the (behaviorally relevant) cross-subject variability in fMRI ac vity should be a ributed to the variability in spa al maps. This has fundamental implica ons for the interpreta on of cross-sec onal studies of func onal connec vity that do not capture cross-subject variability to the same extent as PFMs.