26 September 2017 Mean and variability in functional brain activations differentially predict executive function in older adults: an investigation employing functional near-infrared spectroscopy
Author Affiliations +
Neurophotonics, 5(1), 011013 (2017). doi:10.1117/1.NPh.5.1.011013
Objective: although the preponderance of research on functional brain activity investigates mean group differences, mounting evidence suggests that variability in neural activity is beneficial for optimal central nervous system (CNS) function. Independent of mean signal estimates, recent findings have shown that neural variability diminishes with age and is positively associated with cognitive performance, underscoring its adaptive nature. The present investigation sought to employ functional near infrared spectroscopy (fNIRS) to derive two operationalizations of cerebral oxygenation, representing mean and variability [using standard deviation (SD)] in neural activity, and to specifically contrast these mean- and SD-oxyhemoglobin (HbO) estimates as predictors of cognitive function. Method: a total of 25 older adults (71 to 81 years of age) completed a test of cognitive interference (Multisource Interference Task) while undergoing fNIRS recording using a multichannel continuous-wave optical imaging system (TechEn CW6) over bilateral prefrontal cortex (PFC). Time-varying covariation models were employed to simultaneously estimate the within- and between-person effects of cerebral oxygenation on behavioral performance fluctuations. Results: mean effects were predominantly observed at the between-person level and suggest that greater concentrations of HbO are associated with slower and less accurate performance. Greater HbO variability at the between-person level was associated with slower performance, but was associated with faster performance at the within-person level. Conclusions: these findings are in keeping with assertions that mean and variability confer complementary (as opposed to redundant) sources of information regarding the effective functioning of a neural system and suggest that fNIRS is a viable methodology for capturing meaningful variance in the hemodynamic response that is characteristic of adaptive CNS function.
Halliday, Mulligan, Garrett, Schmidt, Hundza, Garcia-Barrera, Stawski, and MacDonald: Mean and variability in functional brain activations differentially predict executive function in older adults: an investigation employing functional near-infrared spectroscopy



Moment-to-moment variability in neural activity is an emerging area of research that shows promise for elucidating nuances of the human nervous system. In stark contrast to the number of studies examining cognitive behavioral and physiological variability, there is a paucity of research on variability in neural activity, despite the longstanding knowledge that neural variability is not merely noise, but rather a central feature of a stable and well-functioning neural system.12.3 In the present context, neural variability refers to within-person fluctuations in functional brain activity, with evidence to date primarily derived from functional magnetic resonance imaging (fMRI) and electroencephalogram (EEG) research. Neural variability exhibits an inverted U-shaped pattern, increasing through early life and declining through late life, with higher levels generally considered adaptive.2,3 Theoretically, increased neural variability is indicative of a more sophisticated neural system that can explore multiple functional states.3 Variability in neural activity may be functionally significant by facilitating a greater dynamic range of potential responses, according to Bayesian optimization and by enabling itinerant dynamics to avoid determinacy. Indeed, variability in brain activations may arguably represent the substrate for adaptive and stable neural function.3

Empirical support for the functional significance of neural variability includes examples from lifespan developmental phenomena as well as linkages to behavioral performance. McDonnell and Ward4 argue that neural networks are more robust when they are generated in the presence of greater noise (through “stochastic facilitation”), which is further supported by studies showing increased neural variability from infancy to early adulthood.56.7.8 In contrast, decreased variability in the functional magnetic resonance imaging blood-oxygen-level dependent (fMRI-BOLD) signal has been associated with increasing age and diminished behavioral performance in older adulthood.3,9,10 Recent investigations of variability in cerebral oxygenation, using either functional near-infrared spectroscopy (fNIRS) or fMRI, have shown positive associations with superior behavioral performance during measures of scene recognition11 and cognitive flexibility, but not cognitive stability,12 suggesting that neural variability during higher-order cognitive tasks is not only beneficial, but also construct specific. Further, Garrett and colleagues13 examined the impact of increasing cognitive demands on the modulation of brain signal variability. Increasing cognitive load was associated with broad (i.e., multiregion) increases in brain variability for younger and faster-performing adults, but comparatively fewer changes in brain variability for older and slower-performing adults. These age group differences, particularly under increasing cognitive load, were interpreted to reflect a neural system for the younger adults characterized by greater dynamic range between brain states (fixation versus cognitive demand) and an enhanced ability to efficiently process stimuli.

Neural variability also appears to be sensitive to developmental phenomena at both ends of the lifespan and to recovery following injury. EEG data from infants56.7.8 and fMRI data from older adult populations3,9,10 suggest a developmental trajectory of increasing-then-declining neural variability. Garrett and colleagues9 observed that older adults exhibited less neural variability than their younger counterparts, and direct comparison of standard deviation-versus mean-based BOLD patterns indicated that the former shared a five-fold stronger association with age. Further, increased neural variability has been associated with superior behavioral performance3,6,89.10,12,13 and better recovery from traumatic brain injury (TBI).14 The majority of studies examining neural variability in older adults have used fMRI to index the BOLD signal; however, given that the average sampling rate of fMRI methodology is comparatively slow (e.g., Garrett and colleagues3,9,10 used a repetition time (TR) of 2000 ms, resulting in 1 image every 2 s), these studies provide a relatively coarse estimate of BOLD variability reflecting sample-to-sample fluctuations. Further investigation of variability in cerebral oxygenation is therefore an important avenue of exploration, and may be better suited to neuroimaging methodologies indexing cerebral oxygenation that use faster sampling rates (e.g., fNIRS with signal sampling occurring once every 20 ms), to derive estimates that are statistically more precise. For example, for a single 30-s performance block, the SD estimate for a common fMRI TR of 2 s would be based upon only 15 samples (one every 2 s); in contrast, for the same length block within an fNIRS paradigm (at 50 Hz sampling), the within-person SD would be computed across 1500 individual samples (50  per second×30  s). As the hemodynamic signal is inherently noisy,15 and more importantly as moment-to-moment brain signal variability is hypothesized to also contain a durable characteristic signal independent of stochastic noise,1,3,16 the increased sampling density of fNIRS may be particularly useful for deriving precise and reliable estimates of neural variability.


Within- and Between-Person Analyses

Given that the majority of the extant neural variability literature is based on average between-person effects (e.g., does the older age group exhibit less variability than the younger age group?), it is of particular interest to investigate whether associations between behavior and neural function are also observed within persons (i.e., for any given individual, on occasions when neural variability is higher, is behavioral performance also better?). It has long been known that analyses executed at separate levels of nested-data hierarchies (i.e., the between- and within-person levels) do not necessarily yield equivalent results. Referred to as the ecological fallacy,17 results at the individual level may be of particular magnitude and direction, but when aggregated at the group level, can not only differ in pattern but may also be influenced by between-person confounds (e.g., age and cognitive status) that may obscure the effect of interest. For example, it is both conceptually and statistically possible for predictors to account for large proportions of variance within-persons, but to exhibit relatively little or no effect (or even opposite effects) when pooled across individuals.18 Similarly, research on ergodicity19 underscores the importance of examining within-person associations and suggests that (a) process-oriented phenomena (such as neural variability) are best examined within-persons over time and (b) a singular reliance on between-person averages to study such associations is incongruent with theoretical accounts. Therefore, despite being defined by the same outcome, the variation at between- and within-person levels of analysis may represent and be driven by completely different theoretical constructs (e.g., the reason person A performs better on a cognitive task relative to person B may be completely different from why both individuals each exhibit gains across the duration of a task).


Present Study

In the present study, we examined mean and variability operationalizations of the cortical hemodynamic response, at both the within- and between-person levels, for older adults performing an executive functioning task. Several key research aims were explored. First, functional estimates for mean and variability during an executive function task were derived. Relative to fMRI, the faster fNIRS sampling frequency arguably yields a more precise and reliable estimate of neural variability. Although several operationalizations of neural variability have been employed including a block normalized SD approach (fMRI literature) as well as multivariate multiscale entropy (MSE: EEG literature), these various approaches and the signals they index may not be equivalent.3 The findings reported in this investigation are based upon the normalized SD computation.9

Time-varying covariation models were then employed to estimate the within-person (at the individual level) and between-person (at the sample average level) associations for these two oxyhemoglobin (HbO) indicators on behavioral cognitive performance [Multisource Interference Task (MSIT) response latency and accuracy] across all performance blocks within the experiment. For variables such as functional activation that are studied over any longitudinal interval, an estimate for any given cross-sectional sample (e.g., block 1) comprises both between- and within-person sources of variance.20 Accordingly, failure to decompose variance into between- and within-person sources results in estimates that conflate the two—of particular concern to the extent that between-person variance for the phenomenon under study (e.g., the association between neural variability on cognitive performance) is distinct from patterns observed for within-person variance. Our investigation is particularly innovative as most investigations of neural variability are based upon between-person differences; here, we directly examine effects at both the between-person (i.e., effects that pertain to differences in level of neural variability between persons in relation to cognitive function) and within-person (i.e., effects that pertain to neural variability and how this couples with cognitive function within an individual) levels. As so few fMRI and fNIRS investigations have explored neural variability, with virtually all focusing exclusively on between-person differences, key foci of this objective were (a) to examine associations of HbO mean and variability with cognitive function at both the between- and within-person levels, (b) to document whether the patterns were similar in magnitude and direction of effect (an empirical question given the fact that between- and within-person findings need not be identical), and (c) to establish benchmark patterns for both the between- and within-person effects, representing patterns for unconfounded sources of variance, to be replicated in subsequent studies.

Finally, a related research objective pertains to understanding whether the complementary operationalizations (mean and standard deviation) of cortical hemoglobin concentrations were differentially modulated as a function of cognitive load (control versus interference conditions) on the MSIT task. Neural variability remains poorly understood for higher-order cognitive measures, such as executive function, and the modulation of neural variability has not been explored within persons. Consistent with some interpretations offered in the functional literature for the impact of load modulation on mean neural activations,2,13,21 we hypothesized that the more demanding interference condition would result in greater recruitment of neural tissue,2 but that the patterns would differ across operationalizations of the hemodynamic response as well by behavioral metric (accuracy versus latency) and cortical region.





This study was approved by the University of Victoria Human Research Ethics Board and was conducted in accordance with institutional guidelines. Data were collected from a narrow-age cohort 71- to 81-years of age (m=75.88 and SD=3.28) of 25 older adults (13 females and 12 males), 92% of whom were right-handed. Exclusionary criteria included self-report of (a) a physician-diagnosed major medical illness with residual motor or sensory deficits (e.g., Parkinson’s disease, stroke, heart disease, dementia, cancer, and brain tumor), (b) severe sensory impairment (e.g., difficulty reading newspaper-sized print, difficulty hearing a normal spoken conversation, and difficulty writing or pressing buttons), (c) drug or alcohol abuse, (d) history of inpatient psychiatric treatment, (e) significant cognitive impairment (i.e., Mini Mental State Examination score below 24), or (f) English as a second language.



Participants completed the MSIT, which is a computerized task of cognitive interference, and was designed to activate the cingulo-frontal-parietal cognitive attention network in subjects while undergoing functional neuroimaging.22,23 The task shares similarities with the Stroop, Eriksen Flanker, and Simon tasks and is suitable for the purposes of this investigation for notable reasons. Previous investigations using the MSIT have elucidated a reliable network of neuroanatomical correlates that are engaged during task performance,22,23 including several that are accessible by fNIRS methodology using a relatively basic array positioned over the forehead (see Fig. 1). The MSIT was designed with the assessment of neuropsychiatric populations in mind and is regularly employed in the studies of cognitive aging; however, it has yet to be used to examine variability in functional brain activity. The nature of the MSIT is such that behavioral performance can be readily yoked to neural activity, as the task demands remain relatively constant from trial to trial within a block of the experiment.

Fig. 1

The MSIT. Participants are presented with three numbers and indicate the value of the number that is different. The value and location are congruent during control and incongruent during the interference condition.


In the MSIT, participants are presented with an array of three numbers (ranging in value from 0 to 3), one of which is a different numerical value than the others. Using a serial response input device (Psychology Software Tools, Inc.) to ensure +/1  ms timing latency, participants respond to the value of the odd target as quickly as possible while remaining accurate, across a total of 15 trials within a 30-s block. Trial durations are fixed at 2000 ms, allowing behavioral responses to be time-locked to the corresponding samples within a neuroimaging recording. Participants begin with either the control (location and value of the target are congruent) or interference condition (location and value of the target are incongruent), and complete a total of four blocks for each of the conditions, which alternate (see Fig. 1). A rest block of 20 s separates each experimental block and a baseline of at least 90 s preceded the onset of the task. A measure of interference can be derived by comparing the easier (control) to the more demanding (interference) condition, across several outcome measures of interest, including accuracy (percent correct), response time (RT), and hemodynamic response.


Functional Near-Infrared Spectroscopy Recording

fNIRS data were recorded using a continuous-wave TechEn CW6 system (TechEn Inc., Milford, Massachusetts), with a sampling frequency of 50 Hz (corresponding to TR=20  ms) (For selected individuals (n=6), the frequency of fNIRS data acquisition was downsampled to 25 Hz. Despite this, sufficient individual samples were available to derive reliable estimates of HbO mean and variability, with analyses indicating no significant influences of these sampling differences on our estimates or pattern of results.), resulting in 100 images every 2 s. During computerized testing, the fNIRS system recorded cortical hemodynamic responses that were time-locked to events within the tasks. Participants wore custom-built fNIRS headgear consisting of an array positioned over prefrontal cortex (PFC) [Fig. 2(a)] containing 10 channels (8 at 3 cm separation, 2 at 1.5 cm separation) and were tested using wavelengths of 690 and 830 nm to index deoxyhemoglobin (HbR) and oxyhemoglobin (HbO), respectively. The array was designed to maximize the coverage of PFC, given the relevance of PFC areas during performance on the MSIT.

Fig. 2

(a) The array design showing the location of the sources (star), detectors (circle), channels (dashed line), and Brodmann’s areas (squares). (b) Sensitivity profile based on the Monte Carlo forward model (106 photons) for the fNIRS array.


The optical array was positioned relative to several 10 to 20 landmarks (Fpz, Fz, F7, and F8). 3-D coordinates of scalp reference and optode locations were obtained using a Polhemus Fastrak digitizer system (Polhemus, Colchester, Vermont), to perform probabilistic spatial registration.24,25 Following this procedure, Montreal Neurological Institute (MNI) coordinates were generated for the midpoint of each source–detector pair (i.e., channel) for each participant, as well as the average and composite standard deviation values across the group (Table 1). Lastly, we converted the MNI coordinates to Brodmann’s areas (BA) to ascertain macroanatomical labels using Talairach Client software.26 The lateral-most channels in both hemispheres recorded over BA 46, with all the remaining channels recording over BA 10 [Fig. 2(b)]. The 1.5 cm short-separation channels (B4 and C8) were not of interest for the present investigation, given the inability of these channels to capture information at the cortical surface or to facilitate regression of the cardiac signal. Therefore, they were dropped from subsequent analyses. Left hemisphere channels covered inferior frontal gyrus (A1) and middle frontal gyrus (A2, B2, and B3) and right hemisphere channels covered superior frontal gyrus (C5 and C6) and middle frontal gyrus (D6 and D7). To further ascertain whether the array facilitated adequate coverage from PFC regions of interest, the probabilistic path of the light photons using the Monte Carlo forward model (106 photons) was simulated, to derive a sensitivity matrix27 [Fig. 2(b)]. This model was based on the Colin27 atlas, which specifies the absorption properties of scalp, skull, cerebral spinal fluid, gray matter, and white matter. As is evident, the array captured information that is uniformly distributed across PFC.

Table 1

Group average MNI coordinates for the eight-channel array.


The methods undertaken as part of this study were well situated to investigate the within- and between-person associations between HbO variability and behavioral performance, as well as the modulation of neural variability by cognitive load. Although fNIRS is limited in its spatial resolution to cortical regions, the greatest age-related group differences in BOLD variability (between young and old) have emerged in relation to cortical regions.9 Further, the PFC in particular is linked with greater inconsistency in behavioral performance28 and is heavily implicated in executive functions (e.g., cognitive interference).29 Thus, the array used in this study is limited to the coverage of PFC, which is the area relevant to the phenomenon of interest.



Preprocessing of the fNIRS data was performed using Homer 2 software.30 After converting the raw wavelengths to optical density values, we corrected for motion using a wavelet transformation algorithm31 using an interquartile range of 0.1.32,33 Next, we applied bandpass filtering to correct for physiological noise using a high-pass filter value of 0.01 Hz and a low-pass filter value of 0.1 Hz. We then converted from optical density to hemoglobin concentrations by applying the modified Beer–Lambert law, and then exported for subsequent operationalizations. Figure 3 displays the group-averaged time-course data, sampled at 50 Hz.

Fig. 3

Group-averaged time-course data for the MSIT task, sampled at 50 Hz. The presence of each hemodynamic response function corresponds to a 30-s experimental block, separated by 20-s of rest.



Mean and Standard Deviation Operationalizations

Mean- and SD-based operationalizations of HbO were derived for the purpose of this study. Outliers were identified as values >3 SD from the sample mean across block estimates and were deleted pairwise, with subsequent modeling based on restricted maximum likelihood.

Mean estimates of HbO were estimated by aggregating across all samples contained within a given experimental condition. We employed a block design to index HbO within each condition. Specifically, we derived single estimates of HbO within a given block. This equated to eight segments for the MSIT (four control and four interference). Signal variance was estimated with block independent estimates, based on in-house software following the approach used by Garrett and colleagues.3,9,10 Percentage change from the onset value of a given block was computed for each sample, with the corresponding values subsequently averaged within a given block. Given our interest in examining the within-person time-varying covariation between HbO estimates (mean and SD) and behavioral performance (described further below), we did not normalize concatenated blocks to examine single SD estimates, condition wise. In similar fashion, the corresponding estimates for SD HbO were computed per block to facilitate estimation of the time-varying covariation models examining within-person associations between neural activity and cognitive performance.




Behavioral Data

MSIT scores were first screened for blocks with accuracy performance <50%. Given the nature of the task and potential executive functioning difficulties experienced by the participants, there were several blocks of interference results in which participants appeared to have reversed the criteria, responding to the location of the target instead of its value. These blocks were removed from both behavioral and fNIRS analyses (36 interference blocks removed from a total of 100 across all participants). Table 2 displays the demographic and behavioral data.

Table 2

Descriptive statistics for the demographic and behavioral data. Data for MSIT control and interference conditions are based on blocks with at least 50% accuracy and RT data are based on correct trials only.

Years of education17.402.771122
Control accuracy97.933.8985100
Control RT631.22128.80484.921073.85
Interference accuracy82.1313.4953100
Interference RT1142.09126.43988.201449.75


Examining Within- and Between-Person Associations of Cognitive Function with HbO Mean and Variability

Hierarchical linear and nonlinear modeling (HLM) 6.08 software was used to fit linear mixed models to examine the time-varying covariation (or “coupling”) between each operationalization of hemoglobin (i.e., mean and SD) and behavioral performance (i.e., MSIT accuracy and response time), with separate models for each fNIRS channel (i.e., A1, A2, B2, B3, C5, C6, D6, and D7). In order to avoid a type II error and incorrectly rejecting a true finding, we did not correct for multiple comparisons.34 Rather, given that the MSIT and variability operationalizations have not been used with fNIRS in older adults, our preference was to provide a full account of the results. Given our a priori hypotheses regarding expected directional effects, we employed one-tailed (p<0.05) tests for specific planned comparisons. To derive distinct between- and within-person estimates of cortical hemoglobin on behavioral performance, each person-mean operationalization of hemoglobin was centered, before entering it into the model.20 In person-mean centering, the person mean of the time-varying predictor is subtracted from the original time-varying predictor, such that the new time-varying predictor represents variation about one’s own mean level, and thus facilitates the partitioning of HbO-cognition associations into discrete and orthogonal between- and within-person estimates. The following equations outline the analyses conducted to examine the block-to-block covariation between hemoglobin and behavioral performance (i.e., accuracy and response latency):

where behavior represents the outcome measures of accuracy or response latency for person j and block i. Within-person variance is reflected in the level-1 residuals Var(eij), associated with within-person variability block to block. Between-person variance is reflected in the level-2 residuals, Var(U0j) and indicates the amount of variability in behavior that exists between-persons. The parameter estimate for γ10 represents the between-person effect of block on the cognitive outcomes (MSIT accuracy and RT), with parameter γ00 (the fixed intercept) reflecting the between-person average of the cognitive outcome for values of 0 on all predictors (e.g., the very first block for the person-mean-centered values of SD HbO). The parameter estimate for γ01 represents the test of between-person differences in the predictors (mean and SD HbO) on the cognitive outcomes. Specifically, for every unit increase in the person mean of mean or SD HbO, the mean of the cognitive outcome goes up (or down) by the value of the γ01 parameter estimate. Essentially, this parameter reflects the extent to which, for every unit increase in the person mean of HbO (based on either mean or SD computation), the mean of the behavioral performance outcome variable increases or decreases. In contrast, the γ20 parameter estimate represents the test of coupled within-person associations between the predictor and the outcome. Within any given individual, within-person coupling tests whether on occasions when mean or SD HbO is higher, the corresponding estimates of cognitive performance are lower (or higher). Essentially, this parameter reflects the extent to which deviations from an individual’s average amount of HbO across the eight experimental blocks (based on either mean or SD computation) are associated with differences in the same individual’s behavioral performance. The results for the between- versus within-person analyses do not have to be identical and can conceivably differ in magnitude of effect, direction of effect, or both.17,18,20

HbO mean: The within- and between-person associations of HbO mean on variation in behavioral performance were tested separately for each behavioral metric and fNIRS channel. Table 3 summarizes the effects. At the between-person level, greater HbO concentration in the control condition was associated with less accurate performance in two adjacent left hemisphere channels (A1: γ01=0.022, p=0.029; A2: γ01=0.018, p=0.012). At the within-person level, greater HbO concentration in the control condition was associated with faster performance in one right hemisphere channel (D6: γ20=23.14, p=0.019, pseudoR2=0.137), such that on blocks when participants recruited more HbO relative to their own mean, they tended to perform faster (23.14 ms faster for every μMol increase in mean HbO, relative to a given person’s average). Of the total within-person variation in response time, 13.7% was accounted for by changes in HbO concentration, with the pseudo-R2 estimate based upon the Snijders and Bosker computation.35,36 Greater HbO concentration in the interference condition was associated with less accurate performance in both left and right hemisphere channels (A1: γ01=0.104, p=0.015; D6: γ01=0.060, p=0.074) and with slower performance in two left hemisphere channels (A1: γ01=68.46, p=0.063; B3: γ01=59.55, p=0.085) at the between-person level. At the within-person level, greater HbO concentration was associated with less accurate performance in one right hemisphere channel (D7: γ20=0.060, p=0.063, pseudo-R2=0.049), such that on occasions when participants recruited more HbO relative to their own mean, they tended to perform less accurately (0.06% less accurate for every μMol increase in mean HbO, relative to a given person’s average). Of the total within-person variation in accuracy, 4.9% was accounted for by changes in HbO concentration.

Table 3

Two-level multilevel models examining between- and within-person associations for mean HbO with cognitive performance. Coefficients are reported for the between- (γ01) and within-subject (γ20) slope estimates. MSIT, multisource interference task; BS, between-subject; WS, within-subject; BA, Brodmann’s areas.

MSIT accuracyMSIT response time
ControlHbO-BS (γ01)HbO-WS (γ20)HbO-BS (γ01)HbO-WS (γ20)
BA 46 (A1)0.022**0.00514.6811.97
BA 10 (A2)0.018**0.0022.090.07
BA 10 (B2)0.0060.0242.990.96
BA 10 (B3)0.0070.0087.430.32
BA 10 (C5)0.0050.0028.999.95
BA 10 (C6)0.0020.00313.7811.20
BA 10 (D6)0.0080.0027.8823.14**
BA 46 (D7)0.0070.00421.009.95
InterferenceHbO-BS (γ01)HbO-WS (γ20)HbO-BS (γ01)HbO-WS (γ20)
BA 46 (A1)0.104**0.01668.46*9.23
BA 10 (A2)0.0240.00821.1522.09
BA 10 (B2)0.0030.01634.2721.77
BA 10 (B3)0.0250.01559.55*12.22
BA 10 (C5)0.0400.02034.2041.95
BA 10 (C6)0.0210.02244.9014.52
BA 10 (D6)0.060**0.04422.4190.48
BA 46 (D7)0.0140.060*17.4197.83


p<0.05, one-tailed.


p<0.05, two-tailed.

HbO SD: The within- and between-person associations of HbO SD on variation in behavioral performance were tested separately for each behavioral metric and fNIRS channel.

Table 4 summarizes the effects. At the between-person level, greater HbO variability in the control condition was associated with slower performance in two adjacent left hemisphere channels (A1: γ01=0.66, p=0.041; A2: γ01=0.39, p=0.095). At the within-person level, greater HbO variability in the control condition was associated with faster performance in one left hemisphere channel (A2: γ20=0.09, p=0.089, pseudo-R2=0.063), such that on blocks when participants were more variable in HbO relative to their own average variability, they tended to perform faster (0.09 ms faster for every μMol increase in SD HbO, relative to a given person’s average). Of the total within-person variation in response time, 6.3% was accounted for by changes in HbO variability. We also found that greater HbO variability in the interference condition was associated with more accurate performance in one right hemisphere channel (D6: γ20=29.04, p=0.087, pseudo-R2=0.243), such that on blocks when participants were more variable in HbO relative to their own mean, they tended to perform more accurately (0.003% more accurate for every μMol increase in SD HbO, relative to a given person’s average). Of the total within-person variation in accuracy, 24.3% was accounted for by changes in HbO variability.

Table 4

Two-level multilevel models examining between- and within-person associations for SD HbO with cognitive performance. Coefficients are reported for the between- (γ01) and within-subject (γ20) slope estimates. MSIT, multisource interference task; BS, between-subject; WS, within-subject; BA, Brodmann’s areas.

MSIT accuracyMSIT response time
ControlHbO-BS (γ01)HbO-WS (γ20)HbO-BS (γ01)HbO-WS (γ20)
BA 46 (A1)1.4240.1340.66**0.02
BA 10 (A2)0.7040.1840.39*0.09*
BA 10 (B2)0.6441.8540.140.00
BA 10 (B3)0.3540.0840.230.05
BA 10 (C5)0.0640.3840.050.08
BA 10 (C6)0.4440.0440.020.06
BA 10 (D6)0.0540.4240.120.02
BA 46 (D7)0.5540.6740.240.15
InterferenceHbO-BS (γ01)HbO-WS (γ20)HbO-BS (γ01)HbO-WS (γ20)
BA 46 (A1)2.8140.4540.020.03
BA 10 (A2)0.2641.0440.280.15
BA 10 (B2)2.1141.9240.120.11
BA 10 (B3)0.9541.1340.160.16
BA 10 (C5)5.5340.2640.210.15
BA 10 (C6)2.5640.1140.010.05
BA 10 (D6)0.76429.04*0.050.36
BA 46 (D7)9.4740.5940.080.23


p<0.05, one-tailed.


p<0.05, two-tailed.



This investigation sought to examine complementary operationalizations of the cortical hemodynamic response for older adults completing a measure of executive functioning. Neuroimaging data have historically been examined through central tendency computations, in order to derive what have been perceived as more stable estimates based upon block averaging of neural activity. Recent work has revisited previous assertions that neural variability is not merely noise9,10 and has shown that the variability inherent in neural activity conveys functional significance, including associations with brain maturation,56.7 brain senescence,3,9,10 behavioral performance,3,6,8,1011.12 and better recovery from TBI.14 The variability literature on cerebral oxygenation has been largely restricted to fMRI methodology, which is limited by temporal sampling resolution. Thus, more densely sampled profiles of variability for cerebral oxygenation remain virtually unexplored. Another impetus for the present investigation was to examine neural variability using fNIRS and to exploit its comparatively faster sampling frequency to derive variability estimates based upon a greater number of samples that may be more statistically reliable. Further, the ecological advantages of fNIRS (e.g., low cost, portability, and noninvasiveness) may be particularly advantageous for studying variability outside of highly controlled laboratory settings. Finally, although evidence continues to mount in favor of the functional significance of neural variability, it is less clear whether these effects are driven by within- or between-person factors. As emphasized earlier, effects at the between- versus within-person level of analysis may systematically differ due to unaccounted for confounds (e.g., age group differences in performance and individual differences in practice) or due to fundamental differences in the theoretical constructs indexed at each level of analysis (e.g., between-person differences in neural variability may reflect systemic differences in level of central nervous system function, whereas within-person fluctuations in neural activations may reflect more transient influences including stress or momentary lapses of attention). To the extent that variation between- versus within-persons reflects distinct underlying sources, the observed patterns for these distinct levels of analysis may differ in magnitude or direction of effect, or both.20 Thus, it is both conceptually and statistically feasible that as a predictor, neural variability may account for variance at the between- but not the within-person level (i.e., that it may be associated with superior behavioral performance, but be driven by between-person factors). Higher-order cognitive tasks may be more likely to exhibit practice effects relative to other cognitive constructs of interest, further underscoring the importance of separating within- and between-person sources of variance. Among the few fNIRS studies to examine variability, large between-group differences have been reported;37 however, no investigations of variability within-persons have been reported to our knowledge.

Using fNIRS, two operationalizations of HbO were derived, based on central tendency (mean) and variability (SD). Behavioral associations with each operationalization of HbO differed across task difficulty (control or interference), behavioral metric (accuracy or response latency), cortical area, and level of analysis (between- versus within-person). In general, greater amounts of mean HbO (i.e., higher mean activation) were associated with less accurate and slower response in both MSIT conditions, in lateral regions of PFC (including DLPFC). Greater amounts of mean HbO were associated with faster performance in the easier control condition, but with slower performance in the more challenging interference condition. One interpretation of this finding is that recruiting additional neural tissue may have been beneficial as the task remained relatively easy, but that as it became more challenging, additional recruitment of neural tissue was no longer able to compensate for the increased task demands.2 Greater variability in HbO at the between-person level showed some associations to slower performance, especially in the easier condition. At the within-person level, however, greater variability in HbO showed associations to greater accuracy and faster performance in both conditions, with these effects occurring predominantly in lateral areas of PFC. These within-person patterns replicate previous fMRI findings linking BOLD variability to faster and more accurate cognitive performance.20

On balance, these results are in keeping with fMRI findings that mean and variability in the HbO signal confer complementary sources of information.9 Although greater mean was associated with accuracy more so than response time, variability showed a trend for the opposite pattern and showed more reliable associations to response time. With regard to measurement reliability of rapidly changing internal dynamics, response time metrics may afford more reliable estimates relative to accuracy, given the sensitivity of the scales associated with each metric.3839.40 Accordingly, the greater proportion of HbO variability associations with response time suggests that neural variability may be effectively capturing moment-to-moment changes in internal dynamics. Of additional note is that these effects were observed at both between- and within-person levels. Although virtually all studies have tested neural variability–cognition associations at the population (group) rather than at individual level, the most theoretically-empirically matched test of the hypothesis (i.e., that on occasions when neural variability is higher, corresponding estimates of cognitive function are also higher) should be demonstrated within-persons.41 For the more conservative within-person test, greater HbO variability tended to couple with better behavioral performance, in keeping with previous claims that greater variability is adaptive.1, In contrast, the findings linking greater between-person variability to poorer cognitive function may be artifactual based upon population mean confounds, such as age differences.18 Although the coverage of cortical regions is limited by both methodology (fNIRS) and design (array covering regions of PFC only), these preliminary findings are consistent with previous results demonstrating that mean and variability are spatially distinct and are orthogonal in nature.9


Limitations and Future Directions

The observed effects that have been reported in this study are limited due to sample size constraints and duration of task (i.e., number of blocks), resulting in decreased power. To ascertain fully the extent to which each operationalization of HbO is driven by within- or between-person effects, a greater sample size is required; ideally, one in which the incidence of a potential cognitive impairment is apparent (e.g., mild cognitive impairment). Similarly, additional experimental blocks would allow for greater exploration of the within-person associations between HbO and behavior. Event-related designs will also allow for more precise yoking of neural variability with behavioral performance, with an increased number of blocks increasing statistical power for the coupling analyses. The use of short-separation channels would have allowed for a more precise estimate of high-frequency physiological artifacts (e.g., Mayer waves),42 and future studies would benefit from this approach. In using bandpass filtering with relatively conservative thresholds, the reported variability operationalization may result in underestimating the true relationship of neural variability with behavioral performance; future research might replicate these findings employing less conservative screening criteria. As previous results have shown patterns of mean- and SD-based computations of neural activity that are spatially and inferentially distinct,9 the discrepancies reported here between associations of each operationalization of HbO with behavioral performance are likely to become clearer with greater statistical power. The block normalized SD computation was conducted for each channel yielding a total of eight estimates; a distinct approach from that employed in fMRI studies to date. In these studies, a whole-brain estimate across all voxels and regions is derived, which is used subsequently during inferential statistical comparisons (e.g., between young and old adults). Future investigations may seek to examine regional differences in neural variability (e.g., between frontal and parietal cortices) using expanded head coverage, as well as whole-brain SD, to facilitate comparison with fMRI reports. Given the originality of neural variability operationalizations in the fNIRS literature, replication and extension of the reported findings will be essential. Multimodal fNIRS–fMRI recording in particular would allow for an investigation of the comparability of the neural variability signal derived from the same collection.



Variability operationalizations of neuroimaging data are emerging as complementary metrics to conventionally employed mean-based computations. Our results demonstrate that variability in HbO recorded using fNIRS is sensitive to behavioral performance and further substantiate claims that increased neural variability is adaptive. Notably, these associations were apparent at the within-person level, suggesting that they were not driven by between-person confounds such as age differences or individual differences in practice or learning. These patterns are consistent with suppositions that increased neural variability connotes moment-to-moment fluctuations indicative of a nimble, responsive, and dynamic system.3,43 Mean and variability operationalizations of HbO may be sensitive to different metrics of behavioral performance (accuracy or response time), with variability showing greater sensitivity for moment-to-moment changes in rapidly changing internal dynamics. Additional work is needed to further examine the associations between alternative operationalizations of HbO from tasks of varying complexity (e.g., examining the modulation of the HbO signal in tasks that vary in complexity and difficulty). The relatively high temporal sampling frequency of the current fNIRS systems places the methodology in good standing for such future work, as it facilitates deriving variability estimates that are relatively precise and reliable. Further, this line of research may represent a next frontier in facilitating our understanding of the complex interrelations among brain function, cognition, and age.


All authors declare no conflicts of interest.


This research was supported by a Canada Graduate Scholarship from the Canadian Institutes of Health Research (D. H.), and by grants from the Natural Sciences and Engineering Research Council of Canada (S. M. and M. G-B), and the National Institutes of Health/National Institute on Aging (R. S. and S. M., grant number R21AG045575). Further information about the study may be obtained by contacting Dr. MacDonald at smacd@uvic.ca.


1. L. R. Pinneo, “On noise in the nervous system,” Psychol. Rev. 73, 242–247 (1966). http://dx.doi.org/10.1037/h0023240 Google Scholar

2. C. Grady, “The cognitive neuroscience of ageing,” Nat. Rev. Neurosci. 13, 491–505 (2012).NRNAAN1471-003X http://dx.doi.org/10.1038/nrn3256 Google Scholar

3. D. D. Garrett, “Moment-to-moment brain signal variability: a next frontier in human brain mapping?” Neurosci. Biobehav. Rev. 37, 610–624 (2013).NBREDE0149-7634 http://dx.doi.org/10.1016/j.neubiorev.2013.02.015 Google Scholar

4. M. D. McDonnell and L. M. Ward, “The benefits of noise in neural systems: bridging theory and experiment,” Nat. Rev. Neurosci. 12, 415–426 (2011). http://dx.doi.org/10.1038/nrn3061 Google Scholar

5. S. Lippé, N. Kovacevic and A. R. McIntosh, “Differential maturation of brain signal complexity in the human auditory and visual system,” Front. Hum. Neurosci. 3, 48 (2009). http://dx.doi.org/10.3389/neuro.09.048.2009 Google Scholar

6. A. R. McIntosh, N. Kovacevic and R. J. Itier, “Increased brain signal variability accompanies lower behavioral variability in development,” PLoS Comput. Biol. 4, e1000106 (2008). http://dx.doi.org/10.1371/journal.pcbi.1000106 Google Scholar

7. V. A. Vakorin, S. Lippe and A. R. McIntosh, “Variability of brain signals processed locally transforms into higher connectivity with brain development,” J. Neurosci. 31, 6405–6413 (2011).JNRSDS0270-6474 http://dx.doi.org/10.1523/JNEUROSCI.3153-10.2011 Google Scholar

8. B. Mišić et al., “Brian noise is task dependent and region specific,” J. Neurophysiol. 104, 2667–2676 (2010). http://dx.doi.org/10.1152/jn.00648.2010 Google Scholar

9. D. D. Garrett et al., “Blood oxygen level- dependent signal variability is more than just noise,” J. Neurosci. 30, 4914–4921 (2010). http://dx.doi.org/10.1523/JNEUROSCI.5166-09.2010 Google Scholar

10. D. D. Garrett et al., “The importance of being variable,” J. Neurosci. 31, 4496–4503 (2011). http://dx.doi.org/10.1523/JNEUROSCI.5641-10.2011 Google Scholar

11. T. Angsuwatanakul, K. Iramina and B. Kaewkamnerdpong, “Brain complexity analysis of functional near infrared spectroscopy for working memory study,” in Biomedical Engineering Int. Conf. (BMEiCON), pp. 1–5 (2015). http://dx.doi.org/10.1109/BMEiCON.2015.7399531 Google Scholar

12. D. J. N. Armbruster-Genç, K. Ueltzhöffer and C. J. Fiebach, “Brain signal variability differentially affects cognitive flexibility and cognitive stability,” J. Neurosci. 36, 3978–3987 (2016). http://dx.doi.org/10.1523/JNEUROSCI.2517-14.2016 Google Scholar

13. D. D. Garrett et al., “The modulation of BOLD variability between cognitive states varies by age and processing speed,” Cereb. Cortex 23, 684–693 (2013). http://dx.doi.org/10.1093/cercor/bhs055 Google Scholar

14. A. Raja Beharelle et al., “Brain signal variability relates to stability of behavior after recovery from diffuse brain injury,” NeuroImage 60, 1528–1537 (2012).NEIMEF1053-8119 http://dx.doi.org/10.1016/j.neuroimage.2012.01.037 Google Scholar

15. P. Mohr and I. E. Nagel, “Variability in brain activity as an individual difference measure in neuroscience?” J. Neurosci. 30, 7755–7757 (2010). http://dx.doi.org/10.1523/JNEUROSCI.1560-10.2010 Google Scholar

16. C. L. Grady and D. D. Garrett, “Understanding variability in the BOLD signal and why it matters for aging,” Brain Imaging Behav. 8, 274–283 (2014). http://dx.doi.org/10.1007/s11682-013-9253-0 Google Scholar

17. W. S. Robinson, “Ecological correlations and the behaviour of individuals,” Sociol. Rev. 15, 351–357 (1950). Google Scholar

18. S. M. Hofer and M. J. Sliwinski, “Understanding ageing: an evaluation of research designs for assessing the interdependent of ageing-related changes,” Gerontology 47, 341–352 (2001).GERNDJ0304-324X http://dx.doi.org/10.1159/000052825 Google Scholar

19. P. C. M. Molenaar, “A manifesto on psychology as idiographic science: bringing the person back into scientific psychology, this time forever,” Meas. Interdiscip. Res. Perspect. 2, 201–218 (2004).MSRMDA1536-6367 http://dx.doi.org/10.1207/s15366359mea0204_1 Google Scholar

20. L. Hoffman and R. S. Stawski, “Persons as contexts: evaluating between-person and within-person effects in longitudinal analysis,” Res. Hum. Dev. 6(2–3), 97–120 (2009). http://dx.doi.org/10.1080/15427600902911189 Google Scholar

21. L. J. Otten and M. D. Rugg, “Task-dependency of the neural correlates of episodic encoding as measured by fMRI,” Cereb. Cortex 11, 1150–1160 (2001). http://dx.doi.org/10.1093/cercor/11.12.1150 Google Scholar

22. G. Bush et al., “The multi-source interference task: validation study with fMRI in individual subjects,” Mol. Psychiatry 8(1), 60–70 (2003). http://dx.doi.org/10.1038/sj.mp.4001217 Google Scholar

23. G. Bush and L. M. Shin, “The multi-source interference task: an fMRI task that reliably activates the cingulo-frontal-parietal cognitive/attention network,” Nat. Protoc. 1(1), 308–313 (2006).1754-2189 http://dx.doi.org/10.1038/nprot.2006.48 Google Scholar

24. A. K. Singh et al., “Spatial registration of multi-channel multi-subject fNIRS data to MNI space without MRI,” NeuroImage 27, 842–851 (2005).NEIMEF1053-8119 http://dx.doi.org/10.1016/j.neuroimage.2005.05.019 Google Scholar

25. D. Tsuzuki and I. Dan, “Spatial registration for functional near infrared spectroscopy: from channel position on the scalp to cortical location in individual and group analyses,” NeuroImage 85, 92–103 (2014).NEIMEF1053-8119 http://dx.doi.org/10.1016/j.neuroimage.2013.07.025 Google Scholar

26. J. L. Lancaster et al., “Automated talairach atlas labels for functional brain mapping,” Hum. Brain Mapp. 10, 120–131 (2000). http://dx.doi.org/10.1002/(ISSN)1097-0193 Google Scholar

27. C. M. Aasted et al., “Anatomical guidance for functional near-infrared spectroscopy: AtlasViewer tutorial,” Neurophotonics 2, 020801 (2015). http://dx.doi.org/10.1117/1.NPh.2.2.020801 Google Scholar

28. D. T. Stuss et al., “Staying on the job: the frontal lobes control individual performance variability,” Brain 126, 2363–2380 (2003).BRAIAK0006-8950 http://dx.doi.org/10.1093/brain/awg237 Google Scholar

29. E. K. Miller and J. D. Cohen, “An integrative theory of prefrontal cortex function,” Annu. Rev. Neurosci. 24, 167–202 (2001). http://dx.doi.org/10.1146/annurev.neuro.24.1.167 Google Scholar

30. T. J. Huppert et al., “HomerER: a review of time-series analysis methods for near-infrared spectroscopy of the brain,” Appl. Opt. 48, D280–D298 (2009).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.48.00D280 Google Scholar

31. B. Molavi and G. A. Dumont, “Wavelet-based motion artifact removal for functional near-infrared spectroscopy,” Physiol. Meas. 33(2), 259–270 (2012).PMEAE30967-3334 http://dx.doi.org/10.1088/0967-3334/33/2/259 Google Scholar

32. S. Brigadoi et al., “Motion artifacts in functional near-infrared spectroscopy: a comparison of motion correction techniques applied to real cognitive data,” NeuroImage 85, 181–191 (2014).NEIMEF1053-8119 http://dx.doi.org/10.1016/j.neuroimage.2013.04.082 Google Scholar

33. R. J. Cooper et al., “A systematic comparison of motion artifact correction techniques for functional near-infrared spectroscopy,” Front. Neurosci. 6, 147 (2012).1662-453X http://dx.doi.org/10.3389/fnins.2012.00147 Google Scholar

34. S. Lloyd-Fox et al., “Reduced neural sensitivity to social stimuli in infants at risk for autism,” Proc. R. Soc. 280, 20123026 (2013). http://dx.doi.org/10.1098/rspb.2012.3026. Google Scholar

35. J. J. Hox, Multilevel Analysis: Techniques and Applications, Lawrence Erlbaum Associates, Mahwah, New Jersey (2002). Google Scholar

36. J. J. Hox, Multilevel Analysis: Techniques and Applications, 2nd ed., Routledge, New York (2010). Google Scholar

37. H. Sato et al., “Intersubject variability of near-infrared spectroscopy signals during sensorimotor cortex activation,” J. Biomed. Opt. 10(4) 044001 (2005). http://dx.doi.org/10.1117/1.1960907 Google Scholar

38. M. S. Gazzaniga, R. B. Ivry and G. R. Mangun, Cognitive Neuroscience, W. W. Norton & Company, New York (2002). Google Scholar

39. S. W. S. MacDonald, S.-C. Li and L. Bäckman, “Neural underpinnings of within-person variability in cognitive functioning,” Psychol. Aging 24, 792–808 (2009). http://dx.doi.org/10.1037/a0017798 Google Scholar

40. S. W. S. MacDonald, R. S. Stawski, “Intraindividual variability—an indicator of vulnerability or resilience in adult development and aging?” in Handbook of Intraindividual Variability Across the Life Span, , M. Diehl, K. Hooker and M. J. Sliwinski, Eds., pp. 231–57, Taylor & Francis, New York (2015). Google Scholar

41. M. J. Sliwinski, J. Mogle, “Time-based and process-based approaches to analysis of longitudinal data,” in Handbook of Cognitive Aging: Interdisciplinary Perspectives, , S. M. Hofer and D. F. Alwin, Eds., pp. 477–491, Sage, Thousand Oaks, California (2008). Google Scholar

42. M. A. Yücel et al., “Mayer waves reduce the accuracy of estimated hemodynamic response functions in functional near-infrared spectroscopy,” Bio. Opt. Express 7(8), 3078–3088 (2016).BOEICL2156-7085 http://dx.doi.org/10.1364/BOE.7.003078 Google Scholar

43. G. Deco, V. K. Jirsa and A. R. McIntosh, “Emerging concepts for the dynamical organization of resting-state activity in the brain,” Nat. Rev. Neurosci. 12, 43–56 (2011). http://dx.doi.org/10.1038/nrn2961 Google Scholar


Drew W.R. Halliday is a PhD student in clinical neuropsychology at the University of Victoria. His research focuses on the identification of neurophysiological changes that precede the behavioral manifestation of pathology in neurodevelopmental and neurodegenerative conditions. Most recently, he has employed neuroimaging techniques and investigated dispersion in neuropsychological test performance, in order to complement early identification efforts, and to inform targeted intervention strategies with older adults at risk for adverse age-related health outcomes (e.g., MCI, falls).

Bryce P. Mulligan is a neuropsychologist and rehabilitation psychologist at the Ottawa Hospital. He is interested in integrating new techniques in neuropsychological diagnosis, brain physiology, and repeated computerized cognitive testing in an effort to better characterize multitimescale neurocognitive adaptability and improve the early detection of pathological cognitive aging. He has also been involved with meditation and embodiment research, including a randomzied control trial of mindfulness training for older adults with subjective cognitive decline.

Douglas D. Garrett completed his BA degree at the University of Victoria and his PhD at the University of Toronto in 2011, prior to heading the Lifespan Neural Dynamics Group (LNDG) within the Max Planck UCL Centre for Computational Psychiatry and Ageing Research. He was also a fellow with the MPS-UCL Initiative for Computational Psychiatry and Ageing Research from 2011 to 2013, and a senior research scientist/project leader for the intraperson neural dynamics project at the Center for Lifespan Psychology, Max Planck Institute for Human Development in Berlin, Germany.

Stefan Schmidt: Biography is not available.

Sandra R. Hundza received her BSc degree in physiotherapy from the University of Alberta in 1990 and her PhD in neuroscience/kinesiology from the University of Victoria in 2008. She has an extensive clinical background in physiotherapy (neurology, orthopedics and geriatrics). As an associate professor (kinesiology) and adjunct professor (medicine) at the University of Victoria, her research focuses on fall risk as it relates to the interplay between neuromechanical control of gait/balance and cognitive decline in clinical populations including aging and those with neural disorders.

Mauricio A. Garcia-Barrera, PhD, is an associate professor in the Department of Psychology, University of Victoria, Canada. He has been working on developing biologically plausible models to better understand and assess executive control mechanisms. His current research involves three lines of query, examining executive function plasticity associated with sports participation, objective versus subjective effects of sports concussions on executive control, and the neural systems involved in error detection during speech monitoring.

Robert S. Stawski received his PhD in experimental psychology from Syracuse University in 2006. He is an associate professor at the School of Social and Behavioral Health Sciences, College of Public Health and Human Sciences, Oregon State University. His research interest focuses on stress, health, and cognitive aging, as well as research designs and analytic techniques for articulating time-varying and developmental processes. He is a fellow of the Gerontological Society of America

Stuart W. S. MacDonald has helped shape the modern field of cognitive aging by employing advanced measures and methods to characterize patterns and predictors of age-related cognitive decline. Combining expertise in gerontology, cognition, longitudinal methodology, and neuroscience has facilitated his identification of risk factors (e.g., biological, psychological, genetic) that foreshadow cognitive impairment associated with age, dementia, and mortality. His findings have influenced cognitive aging theory and methods, public policy, and personal convictions to age successfully.

© The Authors. Published by SPIE under a Creative Commons Attribution 3.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.
Drew W. R. Halliday, Bryce P. Mulligan, Douglas D. Garrett, Stefan Schmidt, Sandra R. Hundza, Mauricio A. Garcia-Barrera, Robert S. Stawski, Stuart W. S. MacDonald, "Mean and variability in functional brain activations differentially predict executive function in older adults: an investigation employing functional near-infrared spectroscopy," Neurophotonics 5(1), 011013 (26 September 2017). http://dx.doi.org/10.1117/1.NPh.5.1.011013 Submission: Received 7 April 2017; Accepted 29 August 2017
Submission: Received 7 April 2017; Accepted 29 August 2017

Back to Top