The Role of Blood Oxygen Level Dependent Signal Variability in Pediatrics: A Systematic Review

1. Introduction

Variability in the blood-oxygen-level-dependent (BOLD) signal has emerged as a metric with potential clinical relevance. It is no longer viewed as simply “noise” from confounding events during functional magnetic resonance imaging (fMRI) [1,2,3]. At its most inclusive, BOLD signal variability, hereafter referred to as BOLD SV, is a measure of the intraindividual change of the measured BOLD signal, a proxy for neural activity. BOLD SV has been associated with age and cognitive function over the lifespan [4], as well as clinical symptoms in eating disorders [5], attention deficit hyperactivity disorder (ADHD) [6], or 22q11.2 deletion syndrome [2]. Despite the growth of BOLD SV being utilized in pediatric research, there remains inconsistency around its definition and the metrics used to characterize variability. Which metrics exist and which should be used is unclear, rendering BOLD SV challenging to apply or interpret in a standardized manner, especially in a clinical setting. For example, cortical morphology metrics like cortical thickness confound BOLD SV measurements, but this is not consistently demonstrated or accounted for across all metrics [7].

Overall, younger individuals are reported to be more variable in neural processing than older populations. Specifically, an inverted U-shaped trend of cognitive performance and whole brain variability level has related to age. At the same, BOLD SV of individual brain regions presents differing trends across the lifespan [8]. Understanding pediatric BOLD SV may provide insight into critical neuro-developmental processes including maturation of neurotransmitter systems, pruning and neuroplasticity, myelination and white matter integrity, and functional network changes [9,10,11]. For example, higher BOLD SV in medial prefrontal areas comprising the default mode network (DMN) has been shown to positively correlate with ADHD symptom severity [6]. Importantly, infancy and adolescence are unique periods of brain development in which early screening and surveillance can mitigate neurodevelopmental issues. Diagnosing neurological and developmental disorders early on in a child’s life can result in earlier identification and therefore improved chances for intervention [5]. BOLD SV has the potential to be established as a neurodevelopmental biomarker, and contribute to the diagnosis, prognosis, and treatment of neurological disorders in pediatric development [4].

As pediatric BOLD SV is relatively novel, there is a need to provide a foundational framework that gives researchers an entry point into engaging with the topic. This includes clarifying the definition of BOLD SV through identifying and categorizing the various metrics used to measure BOLD SV, and how it has been utilized in the literature.

2. Materials and Methods

2.1. Operational Definitions

BOLD SV in this review was defined as any measure of intra-individual variance in the BOLD signal. The BOLD signal was acquired signal correlated with changes to blood flow and blood oxygenation to localized regions of the brain [12]. The signal also had to characterize the flow of oxygenated hemoglobin being used to support neuronal activity [12]. This definition was chosen to ensure a comprehensive set of all definitions, which are presently not well defined in the BOLD SV literature.

2.2. Article Search Strategy

This systematic review was conducted using the preferred reporting items for a systematic review and meta-analysis (PRISMA) [13]. The electronic literature search was conducted in November 2021 by using the MEDLINE (2003 to 2021), Ovid Healthstar (2003-2021), Psychinfo (2003-2021) and EMBASE (2003 to 2021) through the Ovid platform, as well as the Scopus (2003 to 2021) database. Google scholar was also searched, although no additional citations were captured. The same strategy was used for each database was used to search controlled vocabulary and keywords. These Key terms were: BOLD Signal AND (Variability OR standard deviation OR mean successive difference) AND (Paediatric OR Adolescent OR youth OR Infant. Forward and backward citation tracing were conducted to search for articles that may not captured in electronic databases or Google Scholar.

2.3. Study Inclusion and Exclusion

Articles included used BOLD SV to measure intra-individual variability, or variability changes across the lifespan. Measures of BOLD SV used to adjust for signal noise and other confounders were excluded. All studies possessed participants in their final study sample that were under the age of 18, even if individuals over the age of 18 were included in the study. The inclusion criteria permitted both studies that used an observational or experimental design. Only results relevant to pediatric samples were included in the final extraction. Review papers, conference abstracts and letters to editor case reports and case studies, non-peer reviewed studies, populations over 18, non-human studies, non-English papers, editorials, and any study without enough data (i.e., did not identify age) were all excluded.

2.4. Study Selection and Quality Assessment

A two-stage screening process was conducted to identify relevant articles. All articles were first identified through the electronic database searches and imported into the systematic review management system known as Covidence. Duplicates of captured articles were then removed. Abstract and Title screening was completed by two authors (MD, DRP) who screened in an inclusively, to avoid the removal of potentially relevant articles. Then, a full-text review was independently conducted by two authors (MD, DRP). The review team collaboratively reviewed full-text articles and resolved conflicts. Risk of bias (ROB) and quality assessment of the identified studies were performed using the Downs and Black [14].

2.5. Data Synthesis

Mean and standard deviations (SD) were used to describe the study population’s age and the number of control and diseased patients. In addition, the number and percentage of analysis conducted in studies, and the total number of male and female participants were also calculated from all studies included. Given the heterogeneity of the identified BOLD SV definition subgroups, and of the pediatric populations studied, a meta-analysis was not performed. A data extraction spreadsheet was developed to capture information pertaining to the definitions, study characteristics, sample characteristics, patient characteristics, variability metrics, scales, independent variables, and results of each article. Tables were constructed to summarize study characteristics, BOLD SV metrics and associated findings, and study objectives. Figures were constructed to summarize results of the search via the PRISMA flowchart, summarize overall findings of Metric type and prevalence, and to summarize significant findings in the literature.

3. Results

3.1. Study Sample

Of the five databases searched, Psychinfo contained 20 results, Healthstar had 41, Medline had 46, Embase had 96 and Scopus with 191. 185 unique studies were identified, 17 studies were included. The PRISMA flowchart, depicts the number of studies included and excluded at each screening step (Figure 1). This study included ten different metrics of BOLD SV, a total of 3258 participants and 2869 total pediatrics (1796M/1324F). Metrics classification and associated findings can be found in Table 2. Of the 17 studies, 8 used deviation from average BOLD signal [8,15,16,17,18,19,20,21], 4 used correlational measures of BOLD SV [9,22,23,24], 7 used signal complexity [15,19,21,25,26,27,28] and 1 used the structure of the HRF [29] (Table 1).

3.2. Study Characteristics

As presented in Table 1, all studies were observational. Of these 11/17, (65%) are cross-sectional studies and 6/17 (35%) are case control. 11/17 of the studies exclusively included children (those under the age of 18) and 6/17 included a mixed population that included more than 1 child in their study population.

The statistical analyses performed in the included studies were heterogenous, with the most common being partial least squares analysis (n = 4, 24%) [11,15,16,17], followed by temporal fractal analysis (n = 2, 12%) [25,26] and principle component analysis (PCA) (n = 1, 5.9%) [27]. Age-based variability trends (n= 5, 29.4%) [8,17,18,27,29], behavioural and psychological task performance (n = 2 , 12%) [16,22], and mental disorders (n = 5, 29%) [15,17,19,23,25], neurological disorders (e.g. Autism, ADHD, schizophrenia/psychotic symptoms, and anorexia) (n = 3, 18%) [26,28,29], and genetic conditions (n = 3, 18%) [17,18,24], were examined in relation to BOLD SV. Appendix Table A1 indicates the various objectives of the included articles.

3.3. BOLD SV Metrics (stopped editing here)

Variability metrics were grouped into four categories: BOLD signal deviation from the mean, measures of BOLD SV derived through a correlational analysis of signal, measures of BOLD signal complexity, and measures that utilize characteristics of the hemodynamic response function (HRF). Table 2 presents the descriptions of each variability metric identified in the review.

3.4. Findings Associated with Deviation from the Average BOLD Signal

3.4.1. Standard Deviation of the BOLD Signal (BOLD_SD)

Five studies utilized BOLD_SD (n=5, 72.5%) [16,17,18,20,21] quantifying an average amount of deviation from the typical BOLD signal (Table 2). Most studies had an objective of characterizing variability differences between healthy or typically developing controls and those with a mental disorder, neurological issue or genetic condition (n = 2, 40%) [17,18]. The next most common objective was identifying lifespan and age based patterns in variability (n= 2, 40%) [17,18]. Only one of these studies looked at structural changes across the lifespan (n=1, 20%)8, while the rest looked at network-based changes (n=4, 80%) [17,18,20,21].

A longitudinal study of white matter structure in healthy children, and a cross-sectional study on 22q11.2 deletion syndrome, assessed aging and its effect on BOLD_SD of the brain and age-related variability patterns, respectively. A global association of increasing BOLD_SD with age, particularly in the frontal gyrus, supramarginal gyrus, middle temporal gyrus, and superior parietal lobule was reported [20]. BOLD SV and white matter micro and macro structure metrics like white matter volume, mean fractional anisotropy (FA) and mean diffusivity (MD) measured at younger ages were predictive of BOLD SV at older ages [20]. BOLD_SD associations with these macro and micro structural alterations changed over the lifespan and across various regions throughout the brain.

Two studies focussed on 22q11.2 deletion syndrome, a genetic disorder commonly associated with schizophrenia, and its relationship to BOLD_SD [2,18]. Strong positive psychotic symptoms (PS+) were associated with aberrant age-relationships and concurrently saw BOLD SV increase in visual regions and decrease in the cortices of the prefrontal and orbitofrontal regions of the brain [18]. Both of the studies identified elevated and reduced BOLD_SD across different brain regions, often being lower in regions of the DMN (medial prefrontal cortex, posterior cingulate cortex (PCC), and lateral parietal cortex) [2,18]. Notably, the lack of association between age and BOLD_SD was identified in the dACC or DMN of patients with high psychotic symptom scores (PS+) [2] or schizophrenia, which was typical of healthy controls [2,18] and populations with less severe psychotic symptoms (PS-) and the 22q11.2 deletion [2]. This resulted in globally reduced variability in the dACC region of children with PS+ when compared to children with PS [2].

A task-based study measured internally directed creative cognition using a future simulation task, and an alternate uses task (AUT) also has correlations with BOLD_SD. Performance, which acts as an index of creativity, was negatively correlated with the BOLD_SD [16].

3.4.2. Mean Successive Squared Difference (MSSD)

MSSD is another way deviation from the average signal is measured (Table 1). Three studies utilized MSSD as a metric of BOLD SV [8,15,19] (n = 3, 27.5%). Although heterogeneous, studies in this section focused on diseased based findings. One study assessed BOLD SV’s relationship to autism spectrum disorder (ASD) [15], and another assessed lifespan related trends in various networks across the brain [8]. The third study looked at recovered anorexia patients, but found no significant findings using the metric [19].

A cross-sectional study of ASD and typically developing individuals used MSSD as a metric to quantify BOLD SV. Variability increased linearly in the SN nodes (anterior insula) and the ventral temporal cortex and decrease across subcortical, visual, sensorimotor, DMN, and central executive network (CEN) regions [8,15]. When MSSD was used in a population of children with ASD, positive correlations between MSSD of the BOLD time series and GE in structural networks are present. Brain regions that had positive correlations with GE also had a negative correlation with behavioural severity scores such as the social responsiveness scale (SRS) [15].

3.5. Findings Associated with Correlational Measures of BOLD SV

The following section includes all metrics that attained variability measures through correlation-based methods (n = 4, 24%). This includes temporal variability, GLM derived measures of variability and the BOLD% signal change. Supplementary Table S2 includes the definitions of these metrics.

3.5.1. Temporal Variability

To characterize temporal variability in the functional regions, temporal variability was used in 1 study which included subjects with mental disorders and healthy controls (n=1, 25%) [23]. In typically developing children, age-related trends demonstrated significant increases of temporal variability across the inhibition network, from childhood to adulthood [23]. When characterizing relationships in mental disorders, children with schizophrenia had decreased variability in DMN regions associated with increased activity and connectivity compared to typically developing patients. This decreased variability was also associated with neurocognitive symptoms characteristic of schizophrenia [23]. Increased variability was seen in subcortical regions (thalamus, putamen, and pallidum) in these patients, while children with ADHD, saw increased variability in regions of the DSN and decreased variability in subcortical regions. Variability levels were not the same in the DMN regions of typically developing children with autism and ADHD, with the medial frontal areas mainly affected in ASD, and the posterior cingulate in ADHD [23]. Importantly, regions with the highest variability in controls, (i.e trans modal areas) have lower levels of variability in disorders [23]. Those areas in controls with the lowest variability, like primary sensory regions, are more prevalent in mental disorders [23].

Half of the regions showing significant variability changes in the three disorders are in the top 10% of regions with the highest or lowest variability in controls [23]. Negative correlations between the variability of a brain region and its BOLD activity were found. Low DMN variability was consistently identified alongside strong functional connectivity (FC) within the DMN during resting state fMRI [23].

3.5.2. Multilinear and GLM Derived Variance Measurement

Two studies used multilinear models or GLM derived variance measures (n = 2, 50%) [22,24]. A reading skill task used to investigate deviations from the mean BOLD signal by measuring the standard deviation of a beta series representing mean activation. In the left inferior frontal gyrus pars triangularis, the SD of the series appears to account for additional variance in reading skill, measured as task performance [22].

Only one study measured BOLD SV as the % change in signal. This study was task-based, studying a “Go/No-Go” behavioural task in children and assessed neural correlates of inhibitory control and genetic variation in the DA receptor gene. Having 7-repeats in the variable number of tandem repeats (VNTR) of the DA receptor gene DRD4 (7R+). Having 7 repeats in the VNTR region of DRD4 (7R+) has been associated with psychiatric disorders that present self-regulation issues such as ADHD [24]. Lower variability was found in those with 7R+ when compared to 7R- groups during successfully inhibited prepotent motor response. This was observed in two regions located in the prefrontal cortex, one in the cerebellum and one in the occipital lobe [24]. There were no differences in behavioral performance of the “Go/No-Go” task and correlations between task-related BOLD response were not observed during the task [24].

3.5.3. Difference of Residuals

One cross-sectional study used the difference of residuals metric (n=1, 25%) [9]. This variability metric compares the difference in variability between the two residual models of observed and expected BOLD response. BOLD SV in the inhibition network, was reported as lower in children than adults during a successful stopping task [9].

3.6. Findings Associated with Signal Complexity

Signal complexity (as defined in Table 2) was also used to measure the variability in the BOLD signal. Signal complexity is also described as unpredictability of a signal over its time series [30]. Seven studies were identified that utilized this metric (n=7, 41%) [15,19,21,25,26,27,28].

3.6.1. Entropy/Sample Entropy

One way it is estimated is by using an entropy metric such as sample entropy (SE). SE is used to identify repetitive patterns in a time series, and the degree of regularity of patterns of activation observed [30].

Only 1 study used an entropy-based metric of BOLD SV in populations of children with and without autism spectrum disorder (n=1, 14%) [15]. Distributed brain regions showed increases in MSSD and entropy from childhood through adolescence and positive correlations between entropy, general efficiency (GE) and age in both ASD and typically developing groups [15]. Negative correlations with SRS scores and entropy [15]. Lower levels of sample entropy are seen in ASD individuals during social and non-social tasks [15].

3.6.2. Fractal Dimensionality

Two studies used fractal dimension, obtained by fractal analysis, as a measure of the complexity derived from hurst exponents (n = 2, 29%) [25,26]. FD is a statistical measure of how completely a fractal appears to fill the space in the geometric sense. When used for signals, it can become a metric of structural complexity across a given time domain [31].

Reduced signal complexity in ASD participants with respect to controls in the amygdala, the vermis, the basal ganglia, and the hippocampus. Decreases were correlated with autism diagnostic interview-revised (ADI-R) and autism diagnostic observation schedule (ADOS) scores [25]. The nucleus accumbens and the caudate head showed significantly reduced fractal dimension. Regions of the cerebellum in the ASD cohort showed significant reduced FD, specifically in the vermis with mild correlations with the ADIRRB and ADOSRRB metrics [25].

A cross-sectional study of children with mTBI utilized FD as a metric of BOLD SV. There were 11 brain regions where FD significantly decreased for mTBI patients, including the caudate nucleus and nucleus accumbens [26]. The FD also decreased for mTBI patients when compared with the uninjured control group in both these areas.

3.6.3. Power Based Metrics

Power or spectral density-based variability metrics are an index of the signal amplitude of sinusoidal oscillations within and across frequencies over a time series [10]. Signal which demonstrates scale-free behavior, require fractal-like self-similarity in a spatial or temporal scale to use a power law measure of complexity [21].

Only one study used this metric to quantify complexity in healthy participants (n = 1, 14%) [21]. Increases in complexity were found throughout the whole brain during adolescence and early adulthood, excluding the DMN and attention control networks. Complexity did not change with age in a subset of gray matter regions and dorsal attention networks [21]. Decreases in complexity were observed in other regions of the brain, but the largest reductions occurred in the subcortical gray nuclei [21]. A strong positive correlation between local connectivity (ReHo) and complexity in endogenous brain activity fluctuations was also identified [21]. White matter and areas of gray matter with lower local connectivity exhibited more randomness in their BOLD fluctuations. Basal ganglia, thalami, and spinocerebellum showed relatively lower complexity than would be predicted from ReHo [21].

3.6.4. Fractional amplitude of low-frequency fluctuation (fALFF)

Three studies utilized fALFF as a metric of complexity (n = 3, 43%) [19,27,28]. One cross-sectional study found that fALFF metrics were associated with age, with 5.2% of variability in age attributed to the complexity. This distinguished areas of the DMN and salience network from occipital, temporal, superior parietal, and pre/post central gyral regions. The age-associated fALFF component was also distinguishable from the posterior from anterior cortical regions. DMN anterior regions showed a more pronounced decline, relative to non-DMN posterior regions [27].

Two studies identified used fALFF in populations with mental disorders or neurological conditions [19,28]. It was found that in children with nocturnal enuresis (NE), fALFF was higher in right insula and in the typical spectral band. Regional Homogeneity (ReHo) increased in the left insula and the right thalamus in NE children, and fALFF was increased in the right insula in NE patients [28]. In the slow-5 frequency band fALFF increased in the superior cerebellum, superior temporal gyrus in those with NE. The fALFF in slow-2 were mainly seen in white matter and appeared negative in other bands [28]. In anorexia nervosa patients and healthy controls, values indicated alterations in the temporal gyrus and cerebellum of recovered anorexic patients. Group differences in fALFF were also observed in the cerebellum, specifically in the vermis [19].

3.7. Findings Associated with Characteristics of the Hemodynamic Response Function (HRF)

The hemodynamic response function (HRF) describes behaviour of the BOLD response over time by measuring the change in the HRF with respect to time. A single cross-sectional study uses the “time to peak” in the function, as well as the overall shape, as a metric of variability.

A study of pediatric epilepsy patients discovered that the shape of the HRF changes in children, specifically the amplitude, decreased significantly with greater EEG spike frequency in epileptic patients. When looking at age related trends in the these patients, the intrasubject variability of the amplitude of the HRF does not vary significantly across the age groups of epileptic children [29]. In epilepsy cases, authors reported differences in age could not be distinguished using time to peak or amplitude based metrics of the HRF [29].

4. Discussion

4.1. Summary of Evidence

4.1.1. Metrics

Seventeen studies and a total of 10 unique metrics of BOLD SV were included in this review. These metrics were categorized into 4 types of variability measures: 8 used deviation from average BOLD signal, 4 used correlational measures of BOLD SV, 7 used signal complexity and 1 used the structure of the HRF. In addition, only 7 studies included healthy controls (HCs) exclusively, while 10 included patients with neurological, psychiatric, or genetic disorders.

4.1.2. The Inverted U Trend and BOLD SV

Global increases in whole brain BOLD SV were associated with aging in 3 of the 4 metric types in pediatric populations. This literature has competing perspectives on age-based variability trends. Typically, an inverted U like pattern of variability is reported over the lifespan, suggesting that variability is functionally related to cognitive performance, and increases through childhood to adulthood and decreases in older age [32]. Multiple studies included in this review demonstrate that in both resting state and task-based protocols, whole brain variability increased from 0-18 years of age [15,17,18,25]. This inverted U may relate to development of cognitive capacity, dynamic range and therefore efficiency that increases in childhood, peaks at young adulthood and declines in older age [32]. One cross-sectional study that was identified proposed that lifespan-based trends in BOLD SV are better characterized by networks that simultaneously increase and decrease in variance over the lifespan. They report resting state fMRI data which indicates regions of the SN nodes (anterior insula) and the ventral temporal cortex increase while there are decreases across subcortical, visual, sensorimotor, DMN, and CEN regions [8].

Though not established yet in pediatrics, the inverted U trends of cognitive performance and variability over the lifespan, also follows dopamine signalling strength [34]. Cognition appears to change as DA signaling strength becomes too low or too high, and may play a role in the change in variability throughout the lifespan [34,35]. It should be noted that higher levels of BOLD SV has been specifically associated with elevated levels of cognitive flexibility [32]. In younger and older adults, suboptimal dopamine synthesis capacity is also associated with reduced cognitive flexibility, in addition to reduced BOLD SV [34,37].

4.1.3. BOLD SV Trends in Mental and Neurological Conditions

BOLD SV trends were identified in individuals with ASD, ADHD, schizophrenia, epilepsy, NE, recovered anorexia mTBI and VNTR 7R deletion syndrome. In all disorders, atypical variability trends were uniquely observed across a variety of brain regions and compared to health controls. Regions associated with mental disorders, brain injury and higher symptom severity were also associated with changes in variability.

Interestingly lower variability of DMN specifically was seen in those with schizophrenia, while increased DMN variability in Autism/ADHD were identified, reporting opposite relationship [23]. Schizophrenia and ASD’s/ADHD’s opposite variability trends in similar regions, may be a result of their differing function in the social development patterns associated with the two diseases [33]. In the identified mental disorders, variability trends were inconsistent across regions, with some showing increases, aligning with the global variability trend, while others showed decreases or no change, breaking the trend. Elevated and reduced variability levels compared to healthy controls were consistent across all studies with a focus on mental disorders and neurological conditions, often in affected regions (Figure 3).

4.1.3. Recommendations for Clinical Applications:

Although the literature around BOLD SV does not currently support its use as a clinical biomarker at the moment, it has the potential to become an important one in the future. This includes the use of BOLD SV as a biomarker for normal neurodevelopment and neurodevelopmental issues such as in mental disorders like schizophrenia. Reduced or elevated levels of BOLD SV may also act as a target for treatment in pediatric populations with dopaminergic drugs, especially for those with poor cognitive function and low BOLD SV.

Given global variability’s associations with age in childhood development, there are benefits to producing standardized thresholds of BOLD SV in a typical healthy patient across the different neurodevelopmental milestones in different populations (healthy, diseased, injured etc.). This can be for associated brain regions, networks or across the whole brain. This would allow for the establishment of BOLD SV as a biomarker of neurodevelopment and neurological conditions. To achieve this, consistent use of BOLD SV metrics in the literature and the recruitment of larger pediatric cohorts with low risk of bias is vital. For BOLD SV to be a useful clinical biomarker, researchers must produce normative and non-normative distributions of BOLD SV for various regions of the brain in healthy populations and those effected by different neurodevelopmental disorders respectively.

4.2. Future Directions and Limitation

Given the plethora of metrics identified, future studies should seek to implement multiple variability metrics into their analysis to validate their findings are consistent across metrics and add to the existing knowledge on variability metrics and their limitations. More work must be done in this field to establish or standardize these metrics in a logical way that highlights the strengths and weaknesses of each from a signal processing perspective. This will help create a common framework for their use in a variety of situations as more high-quality evidence is procured.

A complimentary risk of bias analysis was conducted (Table A2), identifying that included articles displayed a risk of bias for external and internal validity overall. Using the Modified Downs and Black checklist identified total scores of 64% for reporting bias, 37% for external validity, 42% internal validity bias, 37% for internal validity confounding and a 11% score for power (Appendix Table A2). An overall score of 47% was obtained from all articles (Appendix Table A2). Criteria that were not included in studies were given a score of 0 unless otherwise indicated. In addition, 7 studies included healthy controls while the rest were cross-sectional and included none. Of the 17 studies, 10 (n=10, 59%) studies included patients with neurological, psychiatric, or genetic disorders. Many of the findings identified come from cross-sectional and case control studies with a risk of bias. These studies are both observational and non-randomized, making it difficult to make etiologic or casual statements regarding risk factors effect on variability. Higher quality evidence with a lower risk of bias will be important to the future of this promising and developing field in order to validate present trends in the literature [34,38]. Once studies on BOLD SV present higher quality results in terms of typical and atypical pediatric populations, BOLD SV can be used as an important biomarker for neurodevelopment.