Quantitative Neuroimaging Meets Normative Modelling: The Last Mile for Precision Medicine Applications

1. Introduction

Since the introduction of X-ray techniques in the late nineteenth century, it quickly became evident that medical imaging would drastically transform our approach to studying and understanding the brain (Narang and Jha, 2017). Today, after more than a century of use, neuroimaging methods have become part of the clinical routine for the diagnosis and monitoring of many brain disorders, as well as tools to define interventions, assist neurosurgery, and inform patient treatment plans. With the advancement of technology and scientific knowledge, neuroimaging methods have significantly evolved (Mishra and Singh, 2010), adding to their capacity to pinpoint the spatial localization of both structural and functional alterations, quantitative information on the magnitude of these aberrant changes. Quantitative neuroimaging refers to the use of imaging techniques, such as magnetic resonance imaging (MRI), positron emission tomography (PET), single-photon emission computed tomography (SPECT), computed tomography (CT), and others, in combination with computational methods to quantify various aspects of brain structure and function (Bremner, 2005). Examples of these include brain morphometry (i.e., quantification of the size, shape, and volume of different brain structures) (Mechelli et al., 2005), brain connectivity (i.e., analysis of the structural and functional connections between different regions of the brain) (Jirsa and McIntosh, 2007), brain metabolism (i.e., measurement of brain glucose or oxygen consumption using techniques such as PET) (Guedj et al., 2022; Sokoloff et al., 1977), and molecular profiles (i.e., quantification of neurochemicals, neurotransmitters, or their receptors using techniques like magnetic resonance spectroscopy and molecular imaging) (Gunn et al., 2015).

Quantitative neuroimaging has played a crucial role in understanding various neurological and psychiatric disorders, as well as in studying normal brain development and aging (Cole and Franke, 2017). It has enabled researchers to objectively assess changes in the brain associated with different conditions and to investigate the effects of novel and existing interventions. Most of these findings follow traditional population-level neuroimaging studies, exemplified by cross-sectional analyses, aimed at identifying differences between average or representative subjects from distinct groups, such as healthy controls versus individuals with a specific brain disorder. While these approaches have proven valuable for gleaning insights into disease mechanisms and therapeutic strategies, they often overlook a fundamental aspect of the disease: its inherent variability across the population. By exclusively examining population averages, cross-sectional analyses make the implicit assumptions that the measured parameters manifest homogeneously across all individuals and that the disease exerts a uniform effect on all patients. In practice, these assumptions frequently do not align with reality (Kapur et al., 2012; Meyer-Lindenberg, 2023).

Recent years have witnessed a notable pivot toward embracing more personalised approaches in the realms of disease diagnosis and treatment (Verdi et al., 2021; Williams and Whitfield Gabrieli, 2025), often defined as “the tailoring of a medical treatment to the individual characteristics of a subject, reducing the need for unnecessary testing and therapies” (Ginsburg and Phillips, 2018). As a result, the landscape of neuroimaging research is undergoing a significant paradigm shift from the conventional “population-level approach” in favour of a more personalised perspective, with implications not only for research but also for clinical applications (Brammer, 2009). One pioneering approach at the forefront of the neuroimaging field is normative modelling (NM). NM is an emerging and data-driven statistical framework that tries to overcome some limitations of traditional cross-sectional analysis. By describing the neuroanatomical variation of brain imaging between subjects and quantifying the individual’s deviation from an expected pattern learned from a reference population, NM aims to identify subject-specific neuroimaging biomarkers of disease state.

In this work, we showcase the current body of work in the field of NM by presenting the practical steps of NM implementation, the main methodologies adopted for model estimation, the evaluation of model fit, and scoring metrics commonly used as measures of deviation from normality. We additionally summarise the state of the art in the most common neuroimaging modalities, including structural MRI (sMRI), functional MRI (fMRI), diffusion MRI (dMRI), and molecular imaging (PET and SPECT). Lastly, we discuss limitations and potential future development for the field.

2. The NM Framework in Neuroimaging Studies

NM always begins with the estimation of reference standards against which individual patients are compared. The rationale behind NM methods is that, to fully characterise the value of any biological measure in disease, one must first understand how it is distributed within a reference healthy population. Once the “normal” range is established, it becomes possible to quantify whether the same measure is altered in the presence of a condition, and if so, the magnitude of this alteration.

In the field of neuroimaging, NM offers the significant advantage of modelling the between-subject variability of any imaging-based biomarker, while preserving the spatial organization of brain measures, without assuming that the manifestation of brain abnormalities – whether structural, functional, or molecular – is homogeneous across individuals for any given disease (Marquand et al., 2016a, 2019). NM implementation requires neuroimaging measures from a reference cohort (generally healthy individuals) and a target cohort (e.g., a patient population). The former is used to model a statistical relationship between the biomarker of interest and a predefined set of covariates of interest. The latter is used for model application to estimate individual maps of deviation from normality.

More specifically, the NM framework involves several steps, as suggested by (Marquand et al., 2019) and summarised in Figure 1: 1) the study design, 2) the statistical model definition and identification, 3) the model validation, and 4) the application to the target cohort. When successfully applied, NM returns both the brain average pattern of the biomarker of interest and its centiles of variation, as a function of a set of covariates of interest. These centiles are used to describe the variability of the reference population and make statistical inferences about the significance of the deviation exhibited by a target subject. The concept is analogous to that of growth charts for pediatric populations, where a newborn’s normal body development (e.g., height, weight, and head circumference) is described as a function of age (Cole, 2012; Rutherford et al., 2022a), and the centiles of variation provide a measure of the expected variability and allow the detection of potential abnormalities.

2.1. Study Design: Data Collection, Harmonization, and Covariates Selection

The study design step is where operative choices are made and necessary data are identified, gathered, and prepared for the analysis.

First, the relevant biological features for the biomarker of interest are computed from the neuroimaging data. This could correspond to some morphological brain measures, as well as some functional ones. The raw data are rarely used for NM in their original forms, but rather processed, quality-checked, and harmonised with each other, to identify from the neuroimaging measures the signal of interest and to remove possible confounding effects. Neuroimaging data can be used with different levels of spatial granularity. This can be at the voxel or vertex level, providing the highest spatial resolution, or at the region-of-interest (ROI) level, which sacrifices spatial resolution for a higher signal-to-noise ratio. If ROI-based NM is adopted, the investigator must also select an appropriate atlas for brain parcellation, which should be consistent across all data. Once the NM is estimated, brain parcellation cannot be changed and must be consistent for both training and application.

One important choice regards the identification of the reference and target cohorts: the reference cohort should represent a normative population, often consisting of healthy individuals or matched controls for a particular class of patients; the target cohort usually represents a group of individuals with a specific disease or condition of interest for the investigation. Particularly critical is the identification of the reference cohort, which should be a large sample that captures a wide range of variability within the whole normative population. As a rule of thumb, the larger the sample, the more accurate the description of population variability, but so far, an optimal strategy for sample size calculations in the context of NM has not been defined.

Independent of the modality, the multicentric nature of large-scale neuroimaging initiatives inevitably induces pooling data from different imaging sites and/or scanners (Fortin et al., 2018), creating data inconsistency and batch effects (similarly to what happens in genomics). This adds to the NM analysis an extra level of complexity, as effects of interest might be hidden by differences in acquisitions (Marzi et al., 2024). Several ways have been proposed to either harmonise the data or account for the site effect when estimating the model (Bayer et al., 2022, 2021; Kia et al., 2022). The simplest approach would consist of the estimation of a different model for each site, thus without the need to handle site differences. However, this approach is unable to take advantage of the richness of available data and can be prone to overfitting when the sample size of each batch is small. Amongst the alternative solutions are the Combat-family methods (Fortin et al., 2018; Pinto et al., 2025). These methods allow for data harmonization by removing the variance associated with batch effects and preserving a priori known sources of variation in data (accounted for in a design matrix) and unknown sources of variation that are not correlated with batch effects. However, they suffer from the risk of removing unknown biologically relevant factors or partially preserving batch effects because of correlation with other covariates of interest. Additionally, it requires a re-harmonization of the data and re-training of the model any time a new scanner/batch is included. On the other side, hierarchical models overcome this issue by modelling instead of removing any side effects. More specifically, these models account for site differences by including an additional random effect. Random effects differ from traditional covariates (fixed effects) as they assume that not all the observations have been observed and allow for the approximation of the effect of a new observation. The effect of a new scanner can be estimated directly in the model without the need for harmonization with other methods and retraining the model. Hierarchical models allow for federated learning approaches in normative modelling, allowing for learning, updating, and adaptation of model parameters on decentralised neuroimaging data (Kia et al., 2022).

Another critical choice in the design of a NM study is related to the selection of biological, clinical, and demographic covariates to use as independent regressors in the statistical model. These variables should be able to explain a “significant” amount of the response variance, which ultimately is defined by the neuroimaging biomarker variability across individuals. On the other hand, it is important to consider the practical availability of these covariates both for the reference and target cohorts. While variables like age or sex are easily accessible, other potential regressors for neuroimaging data, like polygenic risk scores (D’Ambrosio et al., 2022; Janahi et al., 2022) or neuropsychological assessments (Bigler, 1994), may be more rarely available or inconsistent across individuals. In general, as for neuroimaging measures, data collection for predictors should be standardised and homogenous for both reference and target samples.

While most NM applications aim to characterise normative variability in a neuroimaging biomarker as a function of biological, clinical, and demographic covariates, a peculiar class of NM is represented by brain age models. Brain age models focus on predicting individual brain age based on brain imaging data and computing the brain age gap, i.e., the difference between predicted brain age and actual chronological age (Chen et al., 2022a, 2022b; Cropley et al., 2021; Erus et al., 2015). In this approach, we have an inverse problem compared to traditional NM applications in neuroimaging, with the brain imaging biomarkers adopted as model predictors and subjects’ age as model response.

2.2. Model Definition and Identification

This step involves selecting the best model(s) fitting the data at hand and estimating its parameters on the training set. When deciding on the type of modelling, three main families of methodologies can be considered: frequentist, Bayesian, and Deep Learning (DL) methods.

Frequentist Methods

Frequentist methods are so-called for their use of frequentist or inferential statistics. This class of methods can be divided into two sub-groups.

The first group focuses on the mean effect only. Several models have been used in literature belonging to this group, including simple linear regression (LR) or polynomial regression (Coupé et al., 2022; Romascano et al., 2024; Zhang et al., 2023), support vector regression (SVR) (Cropley et al., 2021; Erus et al., 2015), and generalised linear model (GLM) (Ji et al., 2023; Zhang et al., 2022). LR is the simplest method and assumes a linear relationship between the neuroimaging variable and predictors. The point-estimate parameters of the predictors (i.e., mean parameters) are generally estimated via maximum likelihood estimation (MLE). SVR is a type of supervised machine learning algorithm used for regression tasks, where the goal is to predict a continuous numerical output. Analogously to Support Vector Machines (SVM), by employing a kernelized approach, SVR maps input data into a higher-dimensional space, enabling the algorithm to identify an optimal hyperplane that minimises prediction error. The algorithm considers a user-defined tolerance for deviations from the actual values and ignores the errors of any training data close to the model prediction within this tolerance (Erus et al., 2015; Smola and Schölkopf, 2004). Lastly, GLM is a statistical framework that extends the ordinary linear regression model to accommodate various types of response variables, including binary, count, and categorical data. GLMs model the relationship between the mean of the response variable and the predictors through a link function (e.g., logit, probit, and log), which transforms the mean of the response variable to the scale of the linear predictor and specifies a probability distribution for the response variable that belongs to the exponential family of distributions (Wüthrich and Merz, 2023). By design, this class of methods estimates only the mean of the response variable and assesses individual deviations from normality by computing model residuals, i.e., the difference between observed and predicted responses. These approaches are quite appealing for their simplicity and permit obtaining group-level inferences. However, they do not quantify the uncertainty and variation of the biomarker across the whole population, and this must be used by post-hoc regressing the residuals with the true response variable (Erus et al., 2015; Marquand et al., 2019).

This limitation is overcome by other frequentist methods that directly measure the centiles of variation within the population. An example is the quantile regression (Bedford et al., 2024; Lin et al., 2023; Lv et al., 2021; Wong et al., 2023), which estimates the relationship between specific quantiles of the data distribution and a set of covariates (Huizinga et al., 2018; Koenker and Hallock, 2001). Unlike traditional LR, which models the mean, quantile regression focuses on one or more quantiles of interest (e.g., median, quartiles). Another example is given by the Nandaraya-Watson (NW) estimator, a non-parametric kernel regression method that estimates the conditional probability of the response from the set of predictors via computation of a weighted average of the observed values of the response with weights depending on the distance between the points (Elad et al., 2021). More recently, the Generalized Additive Model of Location, Scale, and Shape (GAMLSS) model (Bethlehem et al., 2022; Dinga et al., 2021), has been proposed in neuroimaging studies. GAMLSS extends traditional GLM by allowing not only the mean (i.e., location), but also the scale, skewness, and kurtosis (i.e., shape parameters) of the response variable’s distribution to be modelled as functions of selected covariates. As a generalized additive model, it also overcomes the assumption of linearity of GLM models.

Notably, alongside standard NM approaches, few studies uniquely relied on correlation analysis (e.g., by exploiting Pearson’s based similarity index) to identify deviations from the typical pattern of brain inter-regions anatomical similarity or functional connections (Antoniades et al., 2021; Baldwin et al., 2022; Doucet et al., 2020; Yang et al., 2020).

Bayesian Methods

Bayesian methods represent the bulk of NM approaches. Among them, the most used approaches are Bayesian Linear Regression (BLR) and Gaussian Process Regression (GPR).

BLR is a Bayesian parametric estimation approach that extends traditional LR (Baldwin and Larson, 2017; Corrigan et al., 2024; Schulz et al., 2018) by assuming that both the response variable and the model parameters come from a Gaussian probability distribution. Model estimation will thus reduce to the definition of probability distributions (rather than the point estimates) of these parameters. By adopting Bayesian inference, BLR enables the inclusion of prior knowledge and expectations about model parameters and responses. According to Bayes’ rule, once a set of data is collected, the probability distribution for the model parameters (i.e., the posterior distribution) is computed as the product of the likelihood function, representing the probability of the data given the parameters, and the prior probability. Conceptually, starting from an initial estimate provided by the prior, the model is updated as new evidence is gathered from the data. The resulting probability distribution over model parameters, informed by both the data and the prior, can then be used to estimate the predictive distribution of the response for a new set of predictors (Harnett et al., 2024; Kasper et al., 2024; Savage et al., 2024). It is worth noting that BLR assumes a linear relationship between the response variable and covariates, an assumption that is rarely satisfied in real-world scenarios. One possible solution to this problem involves using a non-linear function to transform the set of inputs into a new, linearised feature space. This approach offers high flexibility to model functions of any shape, but it requires identifying the optimal mapping from numerous available options, which can be challenging in practice. Similarly, extensions of the BLR model, such as B-spline BLR (Shao et al., 2024) and Warped BLR (Lotter et al., 2024; Meijer et al., 2024; Rutherford et al., 2023), have been developed to handle non-Gaussian data distributions. Specifically, these models employ basis functions to map the data onto a Gaussian space (Fraza et al., 2021).

One flexible data-driven Bayesian alternative method is represented by GPR (Blok et al., 2022; Liu et al., 2024; Sun et al., 2023), which allows addressing some limitations of BLR approaches, such as the linearity of response (Marquand et al., 2016a, 2019; Schulz et al., 2018; Wang, 2023). The estimation concept in GPR is similar to BLR, but the mapping is chosen implicitly, allowing the data to determine the shape of the regressing function. Unlike BLR, which determines a distribution over the coefficients, GPR adopts Bayesian inference to directly define a probability distribution over the regressing function, which is modelled as a Gaussian Process (GP). A GP is a statistical process that defines a distribution over functions so that any set of points follows a joint multivariate Gaussian distribution.

The first stage in GPR approaches is to establish a prior GP that defines the prior probability distribution over the set of interpolating functions. This prior is determined by a mean function, usually assumed to be zero, and a covariance function (kernel), that specifies the shape, smoothness, and uncertainty of the function. After establishing the prior, Bayes’ rule is applied to incorporate the knowledge gained from the training data and restrict the GP distribution according to the observed data points. By conditioning the joint Gaussian prior distribution on the collected observations, the posterior distribution over the regressing function can be determined and used to obtain a predictive distribution on new input data. For a more comprehensive and detailed description of the GPR models refer to (Rasmussen and Williams, 2006).

Bayesian methods directly provide a full predictive distribution of the response variable (Ge et al., 2024b), with the mean function describing the expected trajectory and the variances indicating predictive uncertainties. The more data that is available, the more accurate the predictive distribution becomes; conversely, a smaller amount of data results in a wider spread of the posterior distribution.

Bayesian approaches can additionally be estimated as hierarchical models. These models adopt random effects as additional covariates, which differ from the other covariates (referred to as fixed effects), as they assume that not all observations have been observed. The inclusion of random effects allows the approximation of the effect of a new, unseen observation on the fly. This is particularly useful when modelling site effects, since the effect of a new scanner can be estimated directly in the model without the need for harmonization with other methods (Bayer et al., 2022, 2021; Kia et al., 2020).

DL Methods

DL methods represent an emerging family of applications within the context of NM, overcoming some limitations of traditional statistical approaches. In neuroimaging, the study of neurological or psychiatric disorders involves analysing complex, multivariate data with highly variable spatial patterns that standard statistical approaches often fail to model (Lopez Pinaya et al., 2020). To address these challenges, recent studies have proposed DL methods to detect brain abnormalities starting from MRI images (Kumar et al., 2023a, 2023b; Lawry Aguila et al., 2022; Pinaya et al., 2021, 2019). These approaches often use autoencoders to identify patterns relating the response and the covariates in the data. Autoencoders are neural network (NN) architectures composed of two blocks, an encoder and a decoder, trained with unsupervised strategies (Bank et al., 2021). Specifically, the encoder compresses the input into a low-dimensional representation, and the decoder reconstructs the original image from this representation. The objective is to produce a reconstructed image that closely resembles the original.

Autoencoders are particularly suited for applications in NM, and two main types used are variational autoencoders (Feng et al., 2024a, 2024b; Kumar et al., 2023a, 2023b; Lawry Aguila et al., 2022; Wang et al., 2022) and adversarial autoencoders (Pinaya et al., 2021). In DL applications, the contribution of the covariates is integrated into the NM in the form of a loss function for the NNs (Pinaya et al., 2019). Once the training on the reference data is completed, deviation maps are calculated from the neuroimaging data of a new individual as reconstruction error, i.e., the difference between the reconstructed and the original brain image (Kumar et al., 2023b). Individual deviation maps are usually normalised with respect to the mean and variance of the reference deviation map, assuming a minimal reconstruction error for the reference cohort (Kumar et al., 2023b).

2.3. Model Evaluation

After the parameter estimation, the performance of NM is evaluated in terms of data fit and model accuracy to predict the response variable. This step is particularly crucial: if the reference data are not adequately modelled, the estimated subject’s abnormality scores in the target population may be biased. To ensure the generalizability of the model, the evaluation must be performed on unseen data, either in-sample, through k-fold cross-validation (mostly with k=5 or k=10) (Floris et al., 2021; Ge et al., 2024b; Wolfers et al., 2018), or, ideally, out-of-sample on an independent cohort.

Practically, the assessment of model performance reduces to the evaluation of a set of performance metrics measuring the ability of the NM to fit the data. Several alternatives have been proposed in the literature. A detailed description of the metrics can be found in (Dinga et al., 2021) where they are subdivided into two subgroups: measures of goodness of fit of the central tendency and measures of total goodness. Measures of goodness of fit of the central tendency assess how accurately the centre of the distribution (i.e., the mean) is predicted. This is crucial, as poor fit can lead to systematic under- or over-estimation of the predicted response. The most common metrics include the explained variance (EV, or alternatively the coefficient of determination, R²), standardised mean squared error (SMSE), and Rho coefficient between true and predicted responses (Rutherford et al., 2022b). In addition to these, the mean absolute error (MAE) has been used in some works as a goodness-of-fit metric (Dinga et al., 2021; Macleod et al., 2020). On the other hand, measures of total goodness of fit are sensitive to all possible misspecifications of the model, including location, scale, and shape of the predicted distribution. These include measures based on the maximization of model likelihood, such as the mean standardised log-loss (MSLL), Validation Global Deviance (VGD), and the Logarithmic Score (log-score).

In addition to these validation metrics, a few articles also analysed the beta coefficients defining the weight of each regressor and resulting from model fit (Cropley et al., 2021; Floris et al., 2021; Wolfers et al., 2021, 2018). Plotting the estimated beta, especially at the voxel level, may pose challenges due to potential distributions over multiple parameters. Nevertheless, it can offer a rapid and basic method to assess the importance of each regressing covariate on the response variance, allowing a better interpretation and evaluation of the estimated model.

Independent of the specific approach and metric adopted, the evaluation of model performance metrics is a fundamental step in the application of NM since these metrics provide a benchmark of NM performance, prior to its application in the patient cohort. For this reason, they should always be reported before the analysis of any testing cohort. Moreover, when several competing statistical models are estimated, this step allows for an objective and quantitative model selection strategy.

2.4. Model Application to the Target Cohort

The final step of NM is its application: once the model has been defined and trained, given the set of covariates characterizing the target subject, the predictive response for a given target subject – based on their covariates – is quantified, either as a point estimate or as a predictive distribution. Once the predicted and true responses are available, one can measure the discrepancy between them through the computation of a deviation score, which will quantify the distance between the subject’s measures of interest and their expected normal range (Marquand et al., 2016a, 2019).

The most used deviation metric is the Z-score, defined as the ratio between the prediction error (i.e., the bias between true and predicted values of the response) and the model variability (i.e., the standard deviation of the predicted response) (Rutherford et al., 2022b). The Z-score enables evaluation of prediction error while accounting for differences in predictive uncertainty and supports probabilistic inferences about which individuals significantly deviate from the normative pattern (Marquand et al., 2016a). Depending on the granularity at which the neuroimaging biomarker and model have been estimated, Z-scores can be used to generate a deviation map, also known as a Normative Probability Map (NPM), for each subject, either at voxel or ROI level (Marquand et al., 2016a).

Practically, applying a symmetric threshold to Z-scores, defining the range of normality, allows the identification of extreme deviations from the norm (Floris et al., 2021; Janssen et al., 2021; Meijer et al., 2024; Wolfers et al., 2021, 2018). These extreme deviations are interpreted as significant alterations of the biomarker under investigation. The threshold choice reflects the desired significance level (e.g., a Z-score of 2.6 corresponds to a p-value of about 0.01, while a Z-score of 2 corresponds to a p-value of 0.05) and the choice to account or not for multiple comparison.

Participant-level summary statistics of deviation from the normative pattern can also be obtained by aggregating Z-scores across brain regions. This can be performed by averaging all Z-scores or calculating the proportion of extreme deviations relative to the total. Depending on the study goal, either signed (positive or negative) or unsigned (absolute) deviations may be analysed. Importantly, the interpretation of these deviations depends on various study choices, such as imaging modality and Z-score threshold.

An alternative to Z-score thresholding is the extreme value theory, which focuses on the statistical modelling and analysis of the tail ends of a distribution, where rare events occur. These events may correspond to exceptionally large or small values, depending on the context (Gomes and Guillou, 2015). Extreme value statistics is widely used in several fields like environmental sciences, economics, and healthcare to model low-probability, high-impact outcomes (Cotter, 2006; Hamdi et al., 2021; Thomas et al., 2016). This approach can be applied to NM when large enough sample sizes are available, fitting a generalised extreme value (GEV) distribution to the tails of the Z-score distribution of the reference population. This approach offers more accurate modelling of rare deviations and can enhance the detection of atypicality, particularly in psychiatric research (Fraza et al., 2022).

While the Z-score remains the most common deviation metric, alternative measures are also employed. For example, the MAE is often used in DL approaches to quantify reconstruction error, which tends to increase as individual data deviates more from the normative pattern that the trained network has learned (Feng et al., 2024b, 2024a). Additionally, effect size metrics such as Cliff’s delta (Pinaya et al., 2021, 2019) and other specific deviation scores are also used.

3. NM Applications to Neuroimaging

The current section provides an overview of NM applications across major neuroimaging modalities, including sMRI, fMRI, dMRI, and molecular imaging. Literature searches were conducted in November 2024, resulting in the identification of a total of 131 studies from an initial pool of 173. Full details on the search query and inclusion/exclusion criteria are provided in the Supplementary Materials. Selected papers included NM methodological contributions as well as studies presenting NM applications to human brain imaging biomarkers for the study of different clinical or physiological conditions. Specifically, 110 applicative studies were included (Supplementary Table S1). Summary information for these studies is visually illustrated in Figure 2.

At first glance, it appears clear that NM is gaining momentum across various neuroimaging modalities, with an increase of 67% in the number of NM applications to neuroimaging data in 2024 (Figure 2a) with 44 published studies. This approach has been widely adopted in the study of neurodevelopmental conditions, including autism and attention deficit hyperactivity disorder (ADHD) (Bayer et al., 2021; Bedford et al., 2024; Bethlehem et al., 2019; Duan et al., 2024; Floris et al., 2021; Nunes et al., 2020; Pinaya et al., 2019; Rutherford et al., 2022a; Shan et al., 2022; Wolfers et al., 2020; Zabihi et al., 2020, 2019), as well as psychiatric disorders, specifically within the schizophrenia and psychosis spectrum (Berthet et al., 2024; García-San-Martín et al., 2024; Janssen et al., 2021; Kia et al., 2022, 2020; Pinaya et al., 2019; Rutherford et al., 2022a; Segal et al., 2023; Wolfers et al., 2021, 2020, 2018), inherently characterised by high heterogeneity in brain structure and function. A consistent number of studies have also applied NM to neurocognitive disorders (Guo et al., 2024; Kia et al., 2022, 2020; Pinaya et al., 2021) including Alzheimer’s disease, mild cognitive impairment, and Lewy body dementia, mood disorders (Fang et al., 2024; Kia et al., 2022, 2020; Rutherford et al., 2022a; Wolfers et al., 2021, 2018) such as bipolar disorder and depression, as well as to physiological processes such as aging (Bethlehem et al., 2022; Di Biase et al., 2023; Rutherford et al., 2022a). Fewer applications have been reported in the context of pediatric and perinatal conditions (Cromb et al., 2024; Dimitrova et al., 2021; Gaiser et al., 2024), obsessive-compulsive disorders (Han et al., 2023b; Segal et al., 2023, 2024), neurodegenerative diseases (Lawry Aguila et al., 2022) (Chen et al., 2022b) (Villalón-Reina et al., 2023) (Bhome et al., 2024), genetic conditions (Fukami-Gartner et al., 2023; Villalón-Reina et al., 2023, 2022), and other clinical populations (Figure 2b).

Concerning the design of the NM models (Figure 2c), around 60% of the studies employed a Bayesian model, with BLR, either hierarchical or not hierarchical, and GPR models being the most adopted. On the other hand, only a few studies adopted a hierarchical version of the GPR model. Several different frequentist approaches have also been applied (around 28% studies), ranging from log-log regressor (Reardon et al., 2018), and local polynomial regression (Bethlehem et al., 2019), to NW estimator (Elad et al., 2021) and quantile regression (Bedford et al., 2024; Lin et al., 2023; Lv et al., 2021; Wong et al., 2023). More complex approaches were also proposed by more recent works, such as the Rotation-Based Iterative Gaussianization (Snyder et al., 2021), which helps in handling non-Gaussianity of data, or the GAMLSS, for a complete characterization of the probability distribution of the response variable (Dinga et al., 2021; Ge et al., 2024b). Finally, only ~11% of the reviewed studies employed a DL approach for NM (Kumar et al., 2024, 2023a, 2023b; Lawry Aguila et al., 2022; Pinaya et al., 2021, 2019), with the majority of these studies preferring a variational autoencoder as architecture.

Independent of the specific application field, the model set up, and the neuroimaging modality involved, age and sex are generally the most common covariates, followed by data acquisition site and scanner type (Figure 2d).

Approximately 77% of the 78 studies employing multi-site data addressed possible batch effects arising from inter-site differences in scanner and acquisition protocols, either by harmonizing data with Combat, hierarchical modelling, inclusion of site as a predictor, or, more simply, by standardizing preprocessing pipelines to reduce scanner-related variability. A few recent studies have adopted more sophisticated approaches, such as transfer learning and dataset adaptation, to fine-tune pre-trained models on new target datasets (Figure 2e).

Regarding model validation, performance evaluation was explicitly reported only in 78 of the 110 reviewed studies.

The numerosity of reference cohorts ranges from less than 100 scans (Snyder et al., 2021) to over 50,000 (Rutherford et al., 2022a), with one study including more than 100,000 scans (Bethlehem et al., 2022). Data availability is generally highly dependent on the imaging modality chosen (Figure 3). More common modalities, such as sMRI, generally support larger datasets. In contrast, molecular neuroimaging techniques (e.g., PET) often involve smaller datasets, primarily due to higher acquisition costs and modality-specific limitations, such as radiation exposure. Overall, the establishment of international consortia (e.g., ENIGMA (Thompson et al., 2014), ADNI (Petersen et al., 2010)), open-source datasets containing acquisitions from both healthy and clinical populations (e.g., ADNI (Mueller et al., 2005; Petersen et al., 2010), OASIS (Marcus et al., 2007), ARWiBo (Frisoni et al., 2009)), and standardised data-sharing frameworks (e.g., BIDS and its derivatives (Gorgolewski et al., 2016; Holdgraf et al., 2019; Niso et al., 2018; Norgaard et al., 2021; Pernet et al., 2019)) have facilitated the aggregation of large-scale neuroimaging cohorts. Notably, efforts have recently led to the creation of the Normative Neuroimaging Library, including structural, diffusion, and functional MRI scans as well as comprehensive demographic information on around 1,900 healthy individuals, designed as a reference dataset to support the diagnosis and therapeutic development of neurological and psychiatric conditions (Gage et al., 2024).

Interestingly, a different number of quantitative metrics has been adopted for each neuroimaging modality (Figure 4), with sMRI being the best one to be characterised and molecular imaging being the one involving the lowest number of features (only 2 metrics). This trend reflects a greater number of NM applications to MRI data — particularly sMRI data — compared to molecular imaging.

3.1. sMRI

The majority of NM applications (88 out of 110 studies) in neuroimaging focus on sMRI data. The higher number of sMRI publications compared to other modalities reflects the greater availability of this data type, due to its relatively low cost and high adaptability across experimental settings. Applying NM to sMRI data aims to map physiological interindividual variability in brain morphology metrics as a function of age and sex. Among the most commonly used measures are regional cortical thickness (Alvarez et al., 2019; Bayer et al., 2021; Gaiser et al., 2024; Harnett et al., 2024), surface area, and cortical and subcortical volumes (Bedford et al., 2024; Bethlehem et al., 2022; ENIGMA Clinical High Risk for Psychosis Working Group et al., 2024), as well as gray matter volumes (García-San-Martín et al., 2024; Han et al., 2024, 2023b, 2023a), white matter volumes (Sampaio et al., 2024; Wolfers et al., 2020, 2018), and whole brain and intrcranial volumes (Ge et al., 2024a; Laidi et al., 2022; Lawry Aguila et al., 2022; Pinaya et al., 2021). NM approaches have been widely adopted in the study of psychiatric, neurodevelopmental, and mood disorders and have allowed the detection of morphological abnormalities in schizophrenia (Janssen et al., 2021; Kia et al., 2022, 2020; Pinaya et al., 2019; Rutherford et al., 2022a; Wolfers et al., 2021, 2020, 2018), psychosis spectum (Berthet et al., 2024; ENIGMA Clinical High Risk for Psychosis Working Group et al., 2024; García-San-Martín et al., 2024), depression (Fang et al., 2024; Kia et al., 2020; Rutherford et al., 2022a), bipolar disorder (Fang et al., 2024; Kia et al., 2022, 2020; Rutherford et al., 2022a; Wolfers et al., 2021, 2018), but also autism (Bayer et al., 2021; Bethlehem et al., 2019; Floris et al., 2021; Nunes et al., 2020; Pinaya et al., 2019; Rutherford et al., 2022a; Shan et al., 2022; Zabihi et al., 2020, 2019), and ADHD (Bedford et al., 2024; Kia et al., 2022, 2020; Rutherford et al., 2022a; Wolfers et al., 2020). NM provides indeed a powerful tool to account for the high interindividual variability in disease manifestation in these cohorts and detect morphological alterations linked to the pathology. Additionally, a few works also applied NM to study neurodegenerative conditions such as Alzheimer’s disease (Guo et al., 2024; Kia et al., 2022, 2020; Pinaya et al., 2021) or mild cognitive impairment (Kia et al., 2022, 2020).

Another key application of NM to sMRI data is to map brain development and aging through morphological markers (Bethlehem et al., 2022; Di Biase et al., 2023; Rutherford et al., 2022a). Similar to traditional growth charts, brain charts offer a powerful tool to track structural changes across the lifespan, providing normative reference points against which individual brain data can be compared. A notable example is the CentileBrain project, which applied NM, artificial intelligence, and machine learning techniques to brain morphometric data from over 37,000 scans of healthy controls, aged 3 to 90, collected across multiple sites within the ENIGMA consortium (Ge et al., 2024a; Muili et al., 2024). The resulting model establishes reference values for the size of various brain regions, enabling earlier and more accurate detection of atypical structural patterns. Commercial models are also currently available and approved for diagnostic clinical use. An example is icobrain, a CE-marked, FDA-cleared medical imaging software developed by icometrix for the automated analysis of MRI and CT scans, which provides tools for lesion and abnormality detection to adopt in the diagnosis of neurodegenerative disorders (Sima et al., 2024).

3.2. fMRI

Although less represented in the literature compared to sMRI applications (20 out of 110 studies), the use of the NM framework in the analysis of resting state and task fMRI has already demonstrated its ability to capture individual functional alterations across different clinical contexts, including schizophrenia (Lawn et al., 2024; Lin et al., 2024; Oliveira-Saraiva and Ferreira, 2023; Rutherford et al., 2023; Yamashita et al., 2018), bipolar disorder (Lawn et al., 2024; Oliveira-Saraiva and Ferreira, 2023), depression, autism (Jiang et al., 2024; Looden et al., 2022; Savage et al., 2024; Sun et al., 2024) e ADHD (Lawn et al., 2024; Oliveira-Saraiva and Ferreira, 2023; Savage et al., 2024), but also Alzheimer’s and cognitively impaired patients (Huo et al., 2025; Wang et al., 2022). Despite the high heterogeneity of methodologies adopted, most of these studies modelled a response variable derived from resting state functional connectivity analysis of fMRI data, either static or dynamic, including the functional connectivity matrix itself or related derivative metrics (Jiang et al., 2024; Jing et al., 2023; Liu et al., 2024; Looden et al., 2022). These studies aimed to investigate physiological and pathological deviations in normative inter-region connectivity networks (Huo et al., 2025; Jiang et al., 2017; Liu et al., 2024). Only a few papers applied NM to task fMRI data (Floris et al., 2024; Savage et al., 2024; Zabihi et al., 2024). Alongside more standard NM approaches, a few studies uniquely relied on correlation analysis to identify deviations from the typical pattern of brain inter-region connections. As an example, in the study by Yang and colleagues (Yang et al., 2020), voxel-wise fMRI time series were averaged across healthy controls to construct a normative spatio-temporal template of synchronised brain activity. A measure of individual patient deviation from this template was then derived via the computation of similarity maps describing the voxel-wise correlation between the BOLD time series of patients and the norm template, allowing to highlight areas of divergence from the normative pattern. In addition to the studies employing NM directly on fMRI data, one recent study investigated the application of fMRI data enriched by molecular targets (Lawn et al., 2024). In this study, Lawn and colleagues enriched fMRI data with six molecular targets using the REACT framework (Dipasquale et al., 2019). The authors developed region-wise models of functional connectivity related to the six molecular targets using a hierarchical BLR approach, with age and sex as fixed effects and site as the random effect. The constructed models were then applied to cohorts of patients with schizophrenia, ADHD, and bipolar disorder. The estimated NM, by characterizing neurodisorders as divergence from normative ageing, enhanced robustness of between-group similarities and deviation-symptom relationships compared to traditional approaches. Along with clinical research studies, fMRI data was finally adopted in a methodological NM study by Kia and Marquand (Kia and Marquand, 2018) aimed at exploring and validating a novel scalable multi-task GPR model on high-dimensional fMRI data.

3.3. dMRI

A consistent number of applications of the NM to dMRI data is reported in the literature (16 out of 110 studies). By allowing the mapping of water molecules diffusion within the brain tissue, diffusion imaging enables the visualization and delineation of white matter tracts, providing insights into the organization and integrity of white matter pathways, and thus the identification of possible alterations in white matter microstructure. The majority of NM applications to dMRI data aim to model brain regional fractional anisotropy (FA) and mean diffusivity (i.e., measures of the degree of anisotropy and average rate of water molecule diffusion, respectively), often in combination with structural measures derived from anatomical MR images, as a function of age, sex, and other covariates of interest. In the study by Liv and colleagues (Lv et al., 2021), for example, the quantile regression approach was adopted to map regional FA and CT measures as a function of subject age and sex and then applied to patients affected by schizophrenia, highlighting the strong heterogeneity across patients in anatomical loci presenting deviations. A few studies also applied the NM framework to the analysis of white matter fiber geometry derived from tractography. Feng and colleagues proposed a novel DL-based NM framework for analysing white matter microstructure and macrostructure in the context of Alzheimer’s disease, demonstrating its high potential to detect the effects of degenerative disease on white matter structure (Feng et al., 2024a, 2024b). Additionally, Shen and colleagues (Shen et al., 2024) proposed the Brain Representation via Individualized Deep Generative Embedding (BRIDGE) framework, which integrates Bayesian NM and deep generative AI, to model the neurotypical development of the brain structural connectome derived from dMRI tractography. This framework enables the characterization of age-related changes in connectivity and the identification of region-specific patterns of neurodivergence.

3.4. Molecular Imaging

So far, NM applications to the study of molecular imaging data have been notably limited. As far as we know, there are no examples of NM applications to SPECT data, and only recently, a few studies have investigated the validity of the NM framework in PET imaging. The first example is represented by the study of Giacomel and colleagues (Giacomel et al., 2025), who used [¹⁸F]FDOPA scans on healthy controls to estimate voxel-wise BLR models with age and sex used as covariates. These models were applied to a cohort of patients affected by first-episode or chronic psychosis, and the analysis of Z-score allowed the investigation of whole-brain abnormalities of the dopaminergic system in these patients. Comparisons with current state-of-the-art analyses (i.e., assessment of the striatal binding) for the specific tracer and disorder under investigation showed the ability of the NM framework to highlight patterns of alteration in psychosis, normally overlooked in studies with [¹⁸F]FDOPA (i.e., non-striatal regions). Another recent study by Kumar and colleagues (Kumar et al., 2023a) combined standardised uptake values ratio (SUVr) maps from [¹⁸F]-AV45 Amyloid PET images with sMRI data using a multimodal variational autoencoder to identify functional and anatomical deviations in patients with Alzheimer’s disease with respect to a healthy condition. After training the model on cognitively impaired subjects and out-of-sample validation, the model was applied to individuals at various stages of Alzheimer’s disease. The deviation of each Alzheimer’s disease subject from the latent distribution of healthy controls was found to be sensitive to different Alzheimer’s disease stages and significantly associated with individual cognition status. In another recent preprint, Kumar also applied a DL approach to model tau SUVr data (Kumar et al., 2024).

4. Open-Source Tools to Implement Normative Modelling in Neuroimaging Studies

Implementing NM in neuroimaging is not easy, due to the multivariate nature of neuroimaging data and the datasets size. Most studies presented in the literature implement NM with in-house solutions. Nevertheless, several open-source tools have been developed to support and facilitate the application of the NM framework. These tools strongly contribute to promoting NM application in clinical research studies by facilitating its methodological implementation.

• Predictive Clinical Neuroscience Toolkit (PCNtoolkit) (Rutherford et al., 2022b): Python package of pipelines for NM estimation and application to neuroimaging data; based on Bayesian approaches, it is highly flexible in handling multi-site data and complex relationships between response and covariates (https://github.com/amarquand/PCNtoolkit).
• Predictive Clinical Neuroscience Portal (PCN Portal) (Barkema et al., 2023): online platform integrated with the PCNtoolkit providing a user-friendly interface where researchers can access pre-trained normative models estimated from large datasets, upload their data, and obtain normative models without requiring local computational resources or programming expertise (https://github.com/predictive-clinical-neuroscience/PCNportal?tab=readme-ov-file).
• PyNM (Harvey and Dumas, 2022): lightweight Python NM implementation designed to facilitate the implementation of NM operative steps (model estimation, comparison, validation, etc) for users with low technical expertise, providing both a Python API and a command-line interface and a comprehensive tutorial (https://github.com/ppsp-team/PyNM?tab=readme-ov-file#readme).
• CentileBrain (Ge et al., 2024a): web portal providing pre-trained models of typical brain trajectories estimated on large multi-site ENIGMA data for their application to users’ datasets (https://centilebrain.org).

5. Discussion

In this work, we made the case for the use of NM in neuroimaging. The goal of this methodological framework is to model the population variability of neuroimaging biomarkers of interest, such as cortical thickness, responses to functional stimuli, or orientation of brain fibers, and to use it to interpret meaningful alterations when applied to a single patient.

This methodology has been developed and applied to neuroimaging with the intent of overcoming some of the limitations imposed by traditional case-control studies and clustering methods. These methods rely on average representations of subject groups and are the most predominant methods used for the investigation of group differences and within-cohort heterogeneity. However, they both still present some shortcomings. Case-control studies have the limitation of necessitating a priori-defined classes (e.g., healthy controls and schizophrenia patients) and disregarding, or considering as noise or measurement error, the naturally occurring variability of the response variable in the population. This is problematic for clinical domains, such as psychiatry, where different disorders present with a set of heterogeneous and overlapping symptoms and are diagnosed solely from the manifestation of these (Segal et al., 2025). Clustering methods, on the other hand, used for the stratification of patients into subgroups, require the user to select, through appropriate methods, the optimal number of clusters and, once this number is determined, the clustering algorithms will necessarily converge to it (Marquand et al., 2016b). NM relies on a variety of methods to explicitly model the inter-subject variability of the response of interest within the population, without the need for a predefined assignment of subjects of the target cohort to clinical subtypes. Moreover, this methodology allows the detection of distinct patterns of abnormality without requiring a consistent pattern of deviation across individuals (Han et al., 2023b; Huo et al., 2025; Janssen et al., 2025; Lv et al., 2021; Segal et al., 2024, 2023). By doing so, NM addresses key limitations of traditional neuroimaging analysis methods, providing insight into whether a single individual deviates from the normative range, as well as the extent and spatial characteristics of that deviation. This marks a significant shift from group-level to subject-specific statistical analysis, enhancing the clinical utility of neuroimaging biomarkers. Indeed, quantitative neuroimaging biomarkers, such as DAT scans for Parkinson’s disease diagnosis and amyloid imaging as companion diagnostics for newly approved Alzheimer’s disease therapies, rely fundamentally on normative reference values (Klunk et al., 2015; Varrone et al., 2013). Without these benchmarks, their clinical application would not be feasible.

On the other hand, the estimation of NM faces many methodological challenges that hinder its applicability and portability across sites. Firstly, it is worth mentioning the general high computational demand required for the estimation of the NM for the reference cohort, especially in Bayesian methods, which necessitate running multiple chains of Markov Chain Monte Carlo (MCMC), and DL methods, where additional computing power from GPUs is mandatory. Additionally, the NM framework suffers from the need for big datasets for the estimation of the model. Despite the efforts to investigate the sensitivity of model performance to varying training sample sizes (Bozek et al., 2023), explicit, standardised, and comprehensive guidelines for defining the minimum sample size required to ensure reliable estimation of normative models are still lacking.

Generally, to increase the reference cohort size for NM estimation, it is common practice to combine data from multiple studies. However, this often introduces additional sources of variability due to differences in scanners, acquisition protocols, preprocessing pipelines, and inclusion criteria or different definitions of the “healthy” cohorts. The latter, in particular, represents a critical aspect as NM studies typically rely on matched controls from previous cross-sectional studies, where the controls are generally defined merely by the absence of the specific clinical condition of interest and are not necessarily healthy. It is therefore crucial to clearly define the inclusion criteria and, when possible, acquire ad hoc data to establish well-characterised reference cohorts. In general, the variability included by pooling data from different centres does not necessarily represent a limitation, as it increases the heterogeneity of the sample. However, it requires that the complexity of the study is increased by adding some sort of harmonization strategy, either in the model (i.e., hierarchical models) or as an additional preliminary step before modelling (e.g., with the Combat family of methods). Regardless of the specific model adopted, most of these studies on multi-site data included some sort of data harmonization or accounted for the site as a model predictor or random effect. In some cases, specific preprocessing steps for data standardization were adopted.

So far, most NM studies in neuroimaging have used Bayesian approaches, particularly GPR (Fraza et al., 2021), for model estimation. Bayesian approaches provide a probability distribution over the response variable; consequently, they directly allow both the prediction of the response variable, given a new input, and the estimation of the uncertainty of the prediction. They also allow the user to include any a priori information during model estimation. Among Bayesian approaches, GPR has proven particularly powerful as a flexible data-driven method allowing the mapping of different relationships (both linear and non-linear) between the response and covariates, without any a priori assumption on the shape of the regressing function. However, although suitable for smaller datasets (up to a few thousand participants), these approaches do not scale well to larger cohorts (Fraza et al., 2021; Rutherford et al., 2022b). A novel BLR approach with likelihood warping was specifically proposed to address major limitations of widely used GPR methods, namely the poor scalability to large datasets and model non-Gaussian predictive distributions (Fraza et al., 2021). Moreover, Xu et al. (Xu et al., 2021) argued that this type of model alone may not capture the full variability required for accurate normative model estimation.

Recently, DL has emerged as a powerful alternative for the application of NM in neuroimaging. In recent years, DL methods have been increasingly used across various neuroimaging tasks, supported by the growing availability of open datasets and advances in computational power. However, the number of studies applying DL specifically to NM remains relatively limited. DL offers several advantages in the NM context. As “black box” models, DL methods do not require prior assumptions about the underlying system or populations under study. Moreover, recent work has demonstrated the potential of DL to construct multimodal models that incorporate multiple neuroimaging modalities. For example, Kumar and colleagues proposed a multimodal DL-based NM framework combining sMRI and PET data to quantify deviations in Alzheimer’s disease patients (Kumar et al., 2023a). Despite their strong learning capacity, DL approaches face challenges, particularly regarding interpretability and generalizability. Unlike GPR, which may not scale efficiently with large training datasets, DL models can fully exploit large sample sizes. However, in some cases, these methods have underperformed compared to more traditional frequentist approaches in detecting abnormal brain structural patterns (Pinaya et al., 2021, 2019).

Despite methodological complexity and operative issues, NM has been successfully applied to the study of different anatomical, functional, and molecular neuroimaging markers across a wide range of psychiatric diseases (especially in the schizophrenia and psychosis spectrum), neurodevelopmental disorders, but also mood disorders, neurocognitive diseases, and others. This highlights the flexibility and disease-agnostic nature of this approach, reinforcing its potential as a powerful tool for advancing personalised medicine in brain disorders.

6. Conclusion

The NM framework has recently emerged as a pioneering approach in neuroimaging studies of brain disorders. This methodology overcomes the limitations of traditional population-level studies by enabling the detection of subject-level deviations from the normal condition, without requiring a shared pattern of abnormalities across multiple individuals. While future research should prioritise the development of well-defined normative cohorts tailored to different neuroimaging techniques, the rapid expansion of NM demonstrates that this approach is no longer confined to experimental medicine but is increasingly paving the way for the integration of quantitative neuroimaging into clinical practice.