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

Background: As pediatric BOLD Signal Variability (SV) analysis is relatively novel, there is a need to provide a foundational framework that gives researchers an entry point into engaging with the topic. This begins with clarifying the definition of BOLD signal variability by identifying and categorizing the various metrics utilized to measure BOLD SV. Methods: A systematic review of the literature was conducted. Inclusion criteria were restricted to studies utilizing any metric of BOLD SV and with individuals younger than 18 in the study population. The definition of BOLD SV was any measure of intra-individual variability in the BOLD signal. Five databases were searched: Psychinfo, Healthstar, MEDLINE, Embase, and Scopus. Results: A total of 17 observational studies, including male (n = 1796) and female (n = 1324) pediatric participants were included. Eight studies quantified variability as the amount of deviation from the average BOLD signal, seven used complexity-based metrics, three used correlation measures of variability, and one used the structure of the hemodynamic response function. In this study, 10 methods of quantifying signal variability were identified. Associations and trends in BOLD SV were commonly found with age, factors specific to mental and/or neurological disorders such as attention deficit disorder, epilepsy, psychotic symptoms, and performance on psychological and behavioral tasks. Conclusions: BOLD SV is a potential biomarker of neurodevelopmental and neurological conditions and symptom severity in mental disorders for defined pediatric populations. Studies that establish clinical trends and identify the mechanisms underlying BOLD SV with a low risk of bias are needed before clinical applications can be utilized by physicians.


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 [1], as well as clinical symptoms in eating disorders [2], attention deficit hyperactivity disorder (ADHD) [3], or 22q11.2 deletion syndrome [4]. Although it is true that other physiological pulsations can confound the true neural activation found in the BOLD signal, the predictability of these associations should decrease after regression of these confounding post-analysis techniques such as independent component analysis (ICA) denoising and RETROspective Image CORrection (RETROICOR) [1,5]. This allows for the isolation of neural activation effects from other external and physiological contributions that confound the association and improve the signal-to-noise ratio [1].
Overall, younger individuals are reported to be more variable in neural processing than older populations. Older populations report reduced BOLD SV during aging, primarily in subcortical regions. In addition, correlations have been made between poorer performance and reduced BOLD SD in these regions [6]. Although BOLD SV of individual brain regions presents differing trends across the lifespan, an inverted U-shaped trend of cognitive performance and whole brain variability level is observed over the lifespan. Specific trends in the pediatric period have yet to be investigated [6][7][8][9][10].
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 such as cortical thickness confound BOLD SV measurements, but this is not consistently demonstrated or accounted for across all metrics [11].
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 [7,8,10,12]. 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 [3]. 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 [2]. 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 [1].
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. BOLD SV is not yet clearly defined, nor is it measured using a single metric of variability deemed more effective in clinical and research environments. This systematic review attempts to clarify the definition of BOLD SV by identifying and categorizing the various metrics used to measure BOLD SV, and how each metric has been utilized in the literature.

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 [13]. The signal also had to characterize the flow of oxygenated hemoglobin being used to support neuronal activity [13]. This definition was chosen to ensure a comprehensive set of all definitions, which are presently not well defined in the BOLD SV literature.

Article Search Strategy
The preferred reporting items for a systematic review and meta-analysis (PRISMA) guidelines were used to conduct the systematic review [14]. The electronic literature search was conducted in November 2021 by using MEDLINE (2003 to 2021), Ovid Healthstar (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020)(2021), Psychinfo (2003Psychinfo ( -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 to search for 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 was conducted to search for articles that may not have been captured in electronic databases or Google Scholar.

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.

Study Selection and Quality Assessment
A two-stage screening process was conducted to identify relevant articles. All articles were first identified through 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 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 was then conducted using the Downs and Black checklist [15].

Data Synthesis
Mean and standard deviations (SD) described 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 the results of the search via the PRISMA flowchart, summarize overall findings of Metric type and prevalence, and summarize significant findings in the literature.

Study Sample
Of the five databases searched, Psychinfo contained 20 results, Healthstar had 41, Medline had 46, Embase had 96, and Scopus had A total of 185 unique studies were identified, and 17 studies were included. The PRISMA flowchart depicts how many studies were included or excluded at each screening step ( Figure 1). This study included 10 different metrics of BOLD SV, a total of 3258 participants, and 2869 total pediatrics (1796 M/1324 F). Of the 17 studies, 8 used deviation from the average BOLD signal [4,10,[16][17][18][19][20], 4 used correlational measures of BOLD SV [7,[21][22][23], 7 used signal complexity [17,19,20,[24][25][26][27], and 1 used the structure of the HRF [28] (Table 1). Metrics classification and associated findings can be found in Table 2.  Epilepsy, case definition of epilepsy not explicit but EEG-fMRI data were ascertained in children who met the following criteria: (1) An indication for an anatomical scan on the basis of the necessity to investigate a lesion seen on a prior anatomical MRI scan or to diagnose them with epilepsy syndrome and exclude pathological changes. (2) Participants had frequent spikes (N 10 in 20 min) identified on EEG outside the scanner, without occurrence in bursts.
Evaluation of spontaneous regional brain activity in weight-recovered anorexia nervosa Abbreviations: BOLD, Blood oxygen level dependent; SV, signal variability; BOLD SD, blood oxygen level dependent standard deviation, GMV, gray matter volume; fALLF fractional amplitude of low-frequency fluctuations; ASD, autism spectrum disorder; DRD4, Dopamine Receptor D4; VNTR, Variable Number of Tandem Repeats; VNTR 7-repeats present, 7R+; No VNTR 7-repeats present, 7R−; ADHD, attention deficit hyperactive disorder; recAN weight-recovered anorexia; acAN, acute anorexia; fMRI, functional magnetic resonance imaging; resting state, rs; recovered anorexia, recAN; MSSD, mean square successive difference; mTBI, mild traumatic brain injury, NE, nocturnal enuresis. The BOLD time series were segmented into non-overlapping windows, a whole brain signal measure is obtained using Pearson correlation, and a region's variability is compared to others.
Lower variability of DMN in schizophrenia, and increased variability in Autism/ADHD. Changes in variability were closely related to symptom scores and in the 10% most variable regions. Variability increases with age in the inhibition network. More variability in the network was associated with less variability in behavioral performance. Low variability in the DMN was correlated with high FC. Lower variability in 7R+ when compared to 7R− when participants successfully inhibited a prepotent motor response. Primarily seen in the prefrontal cortex, occipital lobe, and cerebellum. The degree of regularity of these patterns of activation was also observed, with fewer complex signals being more random. Positive correlations were identified between entropy, GE, and age. Negative correlations with SRS severity scores and FD in social and non-social tasks, ADIR and ADOS. Grey matter rs-BOLD FD in mTBI patients had reduced FD. Power law exponents remained unchanged or decreased with age and are linearly related to ReHo, which covaried across subjects and gray matter regions. Grey matter rs-BOLD FD in mTBI patients had reduced FD. The fALFF increased with age, distinguishing posterior, and anterior regions.
Higher fALFF values in recAN patient's cerebellum and the inferior temporal gyrus compared to controls. The fALFF decreased in the right insula in children with NE. The ratio of the low-frequency power spectrum, specifically in the range of 0.01-0.08 Hz, to the entire signal frequency range.
Structure of Hemodynamic Response Function (Jacobs et al., 2008 [28]) HRF Structure Using the structure of the HRF, like peak time, amplitude or other signal characteristics not mentioned above.
Could not identify an age-specific HRF. Longer peak times of the HRF 0 to 2 yrs.

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. Here, 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 one child in their study population.

Study Characteristics
As presented in Table 1, all studies were observational. Of these 11/17, (65%) are crosssectional studies and 6/17 (35%) are case control. Here, 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 one child in their study population.

BOLD SV Metrics
Variability metrics were grouped into four categories: BOLD signal deviation from the mean, measures of BOLD SV derived through a correlational analysis of the 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.
A longitudinal study of white matter structure in healthy children, and a crosssectional 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 [10]. BOLD SV and white matter micro and macro structure metrics such as white matter volume, mean fractional anisotropy (FA), and mean diffusivity (MD) measured at younger ages were predictive of BOLD SV at older ages [10]. BOLD SD associations with these macro and micro structural alterations changed over the lifespan and across various regions throughout the brain.
Two studies focused on 22q11.2 deletion syndrome, a genetic disorder commonly associated with schizophrenia, and its relationship to BOLD SD [4,30]. 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 [4]. 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) [4,30]. 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+) [30] or schizophrenia, which was typical of healthy controls [4,30] and populations with less severe psychotic symptoms (PS-) and the 22q11.2 deletion [30]. This resulted in globally reduced variability in the dACC region of children with PS+ when compared to children with PS- [30].
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 [18].

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 [16,17,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) [17] and another assessed lifespan-related trends in various networks across the brain [16]. 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 [16,17]. When MSSD was used in a population of children with ASD, positive associations 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 behavioral severity scores such as the social responsiveness scale (SRS) [17].

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. Appendix A Table A2 includes the definitions of these metrics.

Temporal Variability
Temporal variability was assessed in one study, which included subjects with mental disorders and healthy controls (n = 1, 25%) [22]. In typically developing children, age-related trends demonstrated significant increases in temporal variability across the inhibition network, from childhood to adulthood [22]. When characterizing relationships in mental disorders, children with schizophrenia had decreased BOLD SV in DMN regions associated with higher activity and connectivity compared to typically developing patients. This decreased variability was also associated with neurocognitive symptoms characteristic of schizophrenia [22]. Increased variability was seen in subcortical regions (thalamus, putamen, and pallidum) in these patients, while children with ADHD, saw increased BOLD SV in regions of the DSN and decreased BOLD SV 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 [22]. Importantly, regions with the highest variability in controls, (i.e., trans modal areas) have lower levels of variability in disorders [22]. Those areas in controls with the lowest variability, such as primary sensory regions, are more prevalent in mental disorders [22].
According to a study, 50% of regions with significant changes in BOLD SV in the three disorders are in the top 10% of regions with the highest or lowest variability in controls [22]. Negative correlations between the variability of the signal in a brain region and its level of activity were found. Low DMN variability was consistently identified alongside strong functional connectivity (FC) within the DMN during resting state fMRI [22].

Multilinear and General Linear Model (GLM)-Derived Variance Measurement
Two studies used multilinear models or GLM-derived variance measures (n = 2, 50%) [21,23]. A reading skill task was 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 [21].
Only one study measured BOLD SV as the % change in signal. This study was taskbased, studying a "Go/No-Go" behavioral task in children and assessed neural factors associated with inhibitory control and genetic variation in the Dopamine (DA) receptor gene, having seven repeats in the variable number of tandem repeats (VNTR) of the DA receptor gene DRD4 (7R+). The presence of seven repeats in the VNTR region of DRD4 (7R+) is associated with psychiatric disorders that present self-regulation issues such as ADHD [23]. 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 [23]. There were no differences in behavioral performance of the "Go/No-Go" task and correlations between task-related BOLD responses were not observed during the task [23].

Difference of Residuals
One cross-sectional study used the difference of residuals metric (n = 1, 25%) [7]. 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 [7].

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 the unpredictability of a signal over its time series [31]. Seven studies were identified that utilized this metric (n = 7, 41%) [17,19,20,[24][25][26][27].

Entropy/Sample Entropy
Another way BOLD SV 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 [31].
Only one study used an entropy-based metric of BOLD SV in populations of children with and without autism spectrum disorder (n = 1, 14%) [17]. 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 [17]. Negative correlations with SRS scores and entropy [17]. Lower levels of sample entropy are seen in ASD individuals during social and non-social tasks [17].

Fractal Dimensionality
Two studies used fractal dimension, obtained by fractal analysis, as a measure of the complexity derived from hurst exponents (n = 2, 29%) [24,25]. 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 [32].
Reduced signal complexity was seen 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 [24]. The nucleus accumbens and the caudate head showed significantly reduced fractal dimension. Regions of the cerebellum in the ASD cohort showed significantly reduced FD, particularly in the vermis with mild correlations with the Autism Diagnostic Interview restricted and repetitive behaviors (ADIRRB) and Autism Diagnostic Observation Schedule Restricted and Repetitive Behaviors (ADOSRRB) metrics [24].
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 [25]. The FD also decreased for mTBI patients when compared with the uninjured control group in both these areas [11].

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 [8]. This signal, which demonstrates scale-free behavior, requires fractal-like self-similarity in a spatial or temporal scale to use a power law measure of complexity [20].
Only one study used this metric to quantify complexity in healthy participants (n = 1, 14%) [20]. 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 [20]. Decreases in complexity were observed in other regions of the brain, but the largest reductions occurred in the subcortical gray nuclei [20]. A strong positive correlation between local connectivity (ReHo) and complexity in endogenous brain activity fluctuations was also identified [20]. White matter and areas of gray matter with lower local connectivity exhibited more randomness in their BOLD fluctuations. In the basal ganglia, thalami, and spinocerebellum. relatively lower complexity than would be predicted from ReHo [20] was observed.

Fractional Amplitude of Low-Frequency Fluctuation (fALFF)
Three studies utilized fALFF as a metric of complexity (n = 3, 43%) [19,26,27]. One cross-sectional study found that fALFF metrics were associated with age, with 5.2% of the variability in age attributed to complexity. They distinguished areas of the DMN and salience network in occipital, temporal, superior parietal, and pre-and post-central gyral regions. The age-associated fALFF component was also distinguishable from the posterior from anterior cortical regions. Anterior regions of the DMN had a more evident decline compared to posterior regions not a part of the DMN [26].
Two studies identified used fALFF in populations with mental disorders or neurological conditions [19,27]. It was found that in children with nocturnal enuresis (NE), fALFF was higher in the right insula and in the typical spectral band. Regional Homogeneity (ReHo) rose in the left insula and the right thalamus in children with NE, and the right insula saw increased fALFF in NE patients [27]. In the slow-5 frequency band, fALFF increased in the superior cerebellum and superior temporal gyrus in those with NE. The fALFF in slow-2 was primarily seen in white matter and was observed to be negative in other bands [27]. In anorexia nervosa patients and healthy controls, values indicated alterations in the temporal gyrus and cerebellum of recovered anorexic patients. Between-group differences in fALFF were also observed in the cerebellum, specifically in the vermis [19].

Findings Associated with Characteristics of the Hemodynamic Response Function (HRF)
The hemodynamic response function (HRF) describes the behavior 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 decreases significantly with greater EEG spike frequency in epileptic patients. When looking at age-related trends in these patients, the intrasubject variability of the amplitude of the HRF does not vary significantly across the age groups of epileptic children [28]. In epilepsy cases, authors reported differences in age could not be distinguished using time-to-peak or amplitude-based metrics of the HRF [28].

Metric Utilization
A total of 17 studies and 10 unique metrics of BOLD SV were included in this review. These metrics were categorized into four types of variability measures: Eight used deviation from average BOLD signal, four used correlational measures of BOLD SV, seven used signal complexity, and one used the structure of the HRF. In addition, only seven studies included healthy controls (HCs) exclusively, while 10 included patients with neurological, psychiatric, or genetic disorders.
Presently, deviation-based and complexity-based metrics appear to be the most viable for clinical application and utilization in pediatric BOLD SV research, as the majority of studies included that used these metrics, produced significant results. In addition, complexity and deviation-based methods of quantifying signal variability have been previously established in other neuroscientific contexts that use biomedical signal analysis such as brain variability, heart rate variability, and variability in respiratory rates [33,34]. This attests to these metrics' promise as future biomarkers of BOLD SV. These metrics have been associated with neurophysiological, psychological, genetic, and developmental findings and have been presented in these studies. The pediatric BOLD SV literature specifically includes global and regional changes or associations of variability and neurostructural alterations, factors relating to mental and neurological disorders, genetic markers of disease, and psychological performance and aging ( Table 2).
Deviation-based metrics can be used in the BOLD SV literature to describe how much the magnitude of the BOLD signal changes in the y-axis compared to the baseline level of activation over the course of the time series. Deviation-based metrics appear to be more versatile, as they are used to assess pathological or non-pathologically related trends in neuropsychological development [4,18]. This includes developmental, pathological, and psychologically relevant findings in the BOLD SV literature ( Table 2). Complexity metrics best describe the degree of structure and information of signals and how associated they are with other functional networks [35]. Signals can have equal amounts of deviation from the average but still differ in complexity. Studies included in this review report on only developmental or age-related trends, neurological pathologies, and mental disorders ( Table 2) since aging and pathologies of the nervous system result in signal degradation over time [35][36][37]. Typically, the more complex a signal is, the less it is impacted by unhealthy pathologies or age-related complexity degradation, and the metric should be utilized with this in mind.
Non-deviation or non-complexity-based subgroups include HRF metrics and correlational BOLD SV metrics. HRF-based variation currently is not well established or justified as an effective BOLD SV metric, given the single study using the HRF variability metric did not show significant results [28]. Correlational metrics are highly specific to statistical models utilized in the study, and it is unclear how useful they will be unless they can be applied to clinical settings in a standardized manner. Unlike the HRF, however, they have produced significant results that have established pediatric BOLD SV trends [21,22].

The Inverted U Trend and BOLD SV
Global increases in whole brain BOLD SV were associated with aging in three of the four metric types in pediatric populations. The literature has competing perspectives on agebased 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 [8]. 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 [4,17,24,30]. This inverted U may relate to the development of cognitive capacity, dynamic range, and therefore, efficiency that increases in childhood, peaks in young adulthood, and declines in older age [8]. 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 [16].
Though not established yet in pediatrics, the inverted U trends of cognitive performance and variability over the lifespan also follow dopamine signaling strength [38].
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 [38,39]. It should be noted that higher levels of BOLD SV have been specifically associated with elevated levels of cognitive flexibility [8]. In younger and older adults, suboptimal dopamine synthesis capacity is also associated with reduced cognitive flexibility, in addition to reduced BOLD SV [38,40].

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 was identified, reporting opposite relationships [22]. Schizophrenia and ASD/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 [3,29,41]. 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). Interestingly, lower variability of DMN specifically was seen in those with schizophrenia, while increased DMN variability in Autism/ADHD was identified, reporting opposite relationships [22]. Schizophrenia and ASD/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 [3,29,41]. 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).

Recommendations for Clinical Applications:
Though it is likely that deviation-based or complexity-based metrics have the potential to be used as a clinical biomarker in the future, not enough information is known at this time to make specific assessments on which metrics if any should be used and applied to clinical environments at this stage. More studies conducted using these metrics will allow for a more comprehensive appraisal of BOLD SV metrics. For instance, each metric's ability to be used in combination with signal processing and statistical techniques to mitigate confounding from non-neuronal sources of BOLD SV should be considered before clinical utilization.
Although the variability in neuronal activity is a key component of BOLD SV, there

Recommendations for Clinical Applications
Though it is likely that deviation-based or complexity-based metrics have the potential to be used as a clinical biomarker in the future, not enough information is known at this time to make specific assessments on which metrics if any should be used and applied to clinical environments at this stage. More studies conducted using these metrics will allow for a more comprehensive appraisal of BOLD SV metrics. For instance, each metric's ability to be used in combination with signal processing and statistical techniques to mitigate confounding from non-neuronal sources of BOLD SV should be considered before clinical utilization.
Although the variability in neuronal activity is a key component of BOLD SV, there are other significant contributing sources to variability before processing, filtering, and statistical modeling, and ICA techniques are applied to the signal. Non-neuronal physiological sources of BOLD SV include brain hemodynamics influenced by heartbeat, respiration, and low-frequency oscillations (LFOs). LFOs in and of themselves are derived from a multitude of sources, potentially including Mayer waves, vasomotion from oscillations in vascular tone, CO 2 vasodilation, variations in heart rate or respiratory volumes, gastric oscillations seen using electrogastrograms, and aliased signals of cardiac and respiration due to long signal repetition times (TRs) [42].
These three categories of physiological factors contribute to BOLD SV due to their ability to impact fluctuations in oxygenated hemoglobin concentration in various regions of the brain. Heart rate, vasodilation vasoconstriction, respiratory rate, and respiratory depth over the course of the signal can all cause oxygenated hemoglobin concentration to change independent of brain region activation. This can result in the true associations or effects of BOLD SV being confounded by non-neuronal sources [43]. For instance, age or diseaserelated correlations with decreased BOLD SV could truly be a result of a decrease in resting cerebral blood flow, cerebral metabolic rate of O 2 consumption, and vascular reactivity from aging or disease-related causes. If this variability is not removed or controlled for in the analysis inferences from study conclusions cannot be validly applied to clinical practice [43,44].
These non-neuronal sources of signal can contribute between 20% and 70% of total BOLD SV prior to filtering, component analysis, or other correction steps [42]. Nonneuronal LFO variability is a particularly important target for controlling confounding. The LFO frequency range is where the neuronally related contributions to the BOLD signal can be found and is generally within the range of the low-frequency band (0-0.15 Hz) [45,46]. This is a consequence of the speed at which neuronal activation occurs relative to confounding physiological sources of BOLD SV [45,46]. Caution should be exercised when preparing raw signal for analysis. Approximately 30% of BOLD SV in grey matter is non-neuronal BOLD SV in the same 0-15 Hz frequency range that neuronal BOLD SV is found [42,45,46]. Each contributor to non-neuronal BOLD SV must be mitigated or removed from the signal through statistical modeling to control for these sources of variability. This should be carried out post spectral filtering of the BOLD signal, which already removes 10-15% of the variability in the BOLD signal derived from respiration and cardiac factors [43,47]. If certain metrics are better equipped to isolate neuronal BOLD SV in this frequency range, they should be favored in clinical and non-clinical practice.
Future areas of investigation into BOLD SV should continue to focus on developing BOLD SV into a biomarker of neurodevelopment and a risk factor for neurological issues and mental disorders such as schizophrenia would be a critical advancement of the pediatric BOLD SV literature. Reduced or elevated levels of BOLD SV in particular may identify the need for early intervention or treatment in pediatric populations [48]. In addition, 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 a low risk of bias is vital. For BOLD SV to be a useful clinical biomarker, researchers must produce normative and nonnormative distributions of BOLD SV for various regions of the brain in healthy populations and those affected by different neurodevelopmental disorders, respectively.

Future Directions and Limitations
Given the plethora of metrics identified, future studies should seek to implement multiple variability metrics into their analysis to validate that their findings are consistent, while simultaneously ensuring high-quality evidence is procured. Although it appears deviation-based and complexity-based metrics are most utilized, complexity-based metrics were primarily utilized in assessing associations between BOLD SV and either aging or developmental trends, or neurological pathologies or mental disorders. To verify if these metrics are more effective at characterizing variability in these populations, more work must be done to logically standardize these metrics in a way that highlights the strengths and weaknesses of each from a signal-processing perspective. Each metric should be utilized across diverse populations of patients so that findings may contribute to a future common framework for BOLD SV metric utilization.
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% for internal validity bias, 37% for internal validity confounding, and an 11% score for power (Table A2). An overall score of 47% was obtained from all articles (Table A2). Criteria that were not included in studies were given a score of 0 unless otherwise indicated. In addition, seven 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 [38,39]. 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.  Moment-to-Moment BOLD Signal Variability Reflects Regional Changes in Neural Flexibility across the Lifespan Determine regional variability trends across the lifespan using rs-fMRI data across ages 6-85. Identify brain variability trends in healthy brain development Mulligan et al., 2014 [23] Neural correlates of inhibitory control and functional genetic variation in the dopamine D4 receptor gene Demonstrate carriers of DRD4 (7R+) have reduced prefrontal inhibitory control than non-carriers (7R−) confirmed by locally reduced BOLD % signal change and Go/No-Go task performance.
Clinical: identify pathophysiology of genetic alterations (addition of 7 tandem repeats) of the DRD4 receptor gene Table A1. Cont.

Author and Year Title Objectives of Included Articles Application
Zhang et al., 2016 [22] Neural, electrophysiological, and anatomical basis of brain-network variability and its characteristic changes in mental disorders Characterize temporal variability of the functional architecture of pediatrics with mental disorders in associated regions and networks and establish functional-structural relationships.
Clinical: identify pathophysiology and symptom severity of ADHD, ASD, and schizophrenia and identify SV trends in typically developing patients Zöller et al., 2017 [29] Psychotic symptoms influence the development of anterior cingulate BOLD variability in 22q11.2 deletion syndrome Use the BOLD signal and BOLD SD to assess brain dynamics in patients with psychotic symptoms seen in schizophrenia.
Clinical: identify pathophysiology and symptom severity of schizophrenia Dona et al., 2017 [24] Temporal fractal analysis of the rs-BOLD signal identifies brain abnormalities in autism spectrum disorder Utilize a model-free complexity analysis based on the fractal dimension of the rs-BOLD signal to identify BOLD SV trends in children with ASD.
Identify pathophysiology and symptom severity ASD Jacobs et al., 2008 [28] Variability of the hemodynamic response as a function of age and frequency of epileptic discharge in children with epilepsy Identify age-related changes in the HRF after interictal epileptic discharge and identify the relationship between the number of spikes and a change in the HRF.
Clinical: Identify Pathophysiology of Epilepsy/interictal epileptic discharge in development Seidel et al., 2020 [19] Evaluation of spontaneous regional brain activity in weight-recovered anorexia nervosa Investigate intrinsic regional brain activity differences between recAN, acAN, and healthy controls using resting-state measures of fALFF, ReHo (as reported in acAN, MSSD, DC, and VHMC (as new measures of spontaneous regional brain activity).
Clinical: Identify long-term effects and pathophysiology of recAN and patients Anderson et al., 2013 [20] Complexity of low-frequency blood oxygen level-dependent fluctuations covaries with local connectivity Determine if BOLD fluctuations exhibit temporal complexity associated with local connectivity (ReHo), across gray matter regions utilizing power law behavior.
Identify associations between grey matter ReHo and temporal complexity Dona et al., 2017 [25] Fractal Analysis of Brain Blood Oxygenation Level Dependent (BOLD) Signals from Children with Mild Traumatic Brain Injury (mTBI) Utilize FD to quantify local neural activity in the brain as a metric of functional changes in mTBI patients.
Clinical: Identification of pathophysiology of mTBIs.

Author and Year Title Objectives of Included Articles Application
Wang et al., 2021 [10] The longitudinal relationship between BOLD signal variability changes and white matter maturation during early childhood Investigate longitudinal relationships between BOLD SV, age, and white matter structure in early childhood.
Identify variability trends between brain region SV, structure, and age.
Zheng et al., 2021 [27] Frequency-specific alterations of the resting-state BOLD signals in nocturnal enuresis: an fMRI Study Use associations of spontaneous brain activity and scores on urinary intention-related wakefulness in 5 frequency bands to detect abnormalities of activity in local brain regions of children with nocturnal enuresis.