Mapping human social brain specialisation beyond the neuron using multimodal imaging in human infants

The specialised regional functionality of the mature human cortex partly emerges through experience-dependent specialisation during early development. Our existing understanding of this process is based on evidence from unitary imaging modalities and has thus focused on isolated changes in spatial or temporal precision of neural or haemodynamic activation alone, giving an incomplete picture of the process. We speculate that neural specialisation of function will be underpinned by better coordinated haemodynamic and metabolic changes in a broader orchestrated physiological response. Thus, we present a harmonised framework in which specialisation is indexed by the emergence of coupling between neuronal activity and vascular supply of oxygen and energy. Here, we combine simultaneous measures of coordinated neural activity (EEG), metabolic rate and oxygenated blood supply (broadband near-infrared spectroscopy) to measure emerging specialisation in the infant brain. In 4-to-7-month-old infants, we show that social processing is accompanied by spatially and temporally specific increases in coupled activation in the temporal-parietal junction, a core hub region of the adult social brain. During non-social processing coupled activation decreased in the same region, indicating specificity to social processing. Coupling was strongest with high frequency brain activity (beta and gamma), consistent with the greater energetic requirements and more localised action of high frequency brain activity. We conclude that functional specialisation of the brain is a coordinated activity across neural, haemodynamic, and metabolic changes, and our ability to measure these simultaneously opens new vistas in understanding how the brain is shaped by its environment.


Introduction
Broadband near-infrared spectroscopy (or bNIRS) is a new technique that can be used to 98 quantify the relationship between the neuronal, hemodynamic, and metabolic activity in the 99 infants' brain as it allows the simultaneous and non-invasive acquisition of haemodynamic 100 and metabolic activity concurrently with EEG during functional activation. This technology 101 uses a broad range of optical wavelengths which allows the measurement of the oxidation 102 state of mitochondrial respiratory chain enzyme cytochrome-c-oxidase (CCO), thereby 103 providing a direct measure of cellular energy metabolism (18). CCO is located in the inner 104 mitochondrial membrane and serves as the terminal electron acceptor in the electron transport 105 chain (ETC). It therefore accounts for 95% of cellular oxygen metabolism. In this way, 106 bNIRS allows non-invasive measurement of cellular energy metabolism alongside 107 haemodynamics/oxygenation in awake infants. We recently showed the feasibility of using 108 bNIRS in 4-to-7-month-old typically developing infants (19) and demonstrated the presence 109 of unique task-relevant, regionally specific functional networks where high levels of 110 haemodynamic and metabolic coupling were observed. Here, we integrate this methodology 111 with EEG to identify markers of early brain specialisation with coordinated energetic 112 coupling and neural activity. We develop a novel analysis pipeline to identify localised 113 coupling responses that are modulated by naturalistic social content. We predicted that 114 coupling would be most pronounced in the high-frequency beta and gamma band (20-25) 115 (26), and we hypothesised that we would identify core localised social brain regions with 116 coordinated increases in coupled neural activity, metabolic changes and neurovascular 117 response in the infant brain. 118 119 Results 120 121 Naturalistic social stimuli elicit expected increases in broadband EEG activity: 5-month-122 old infants n=42) viewed naturalistic social and non-social stimuli (Fig 1a) while we 123 concurrently measured EEG and broadband NIRS. Fourier-transform of continuously 124 recorded EEG data from 32 channels (n=35) in one-second segments across the time course 125 of stimulus presentation confirmed robust broadband increases in neural activity in response 126 to social versus non-social stimuli (Fig 1b, replicating (11)). coupled increases in metabolic function and oxygenated blood flow (19). This revealed 141 distinct locations sensitive to social (Fig 2b) and non-social (Fig 2d) processing; the 142 topography of these locations is strikingly similar to the topography of differentiated 143 broadband EEG activity (Fig 2a, Fig 3a). A general linear model (GLM) 157 approach was then used to identify FDR-corrected associations between EEG channels and 158 bNIRS channels that showed significant coupling between metabolic response and 159 oxygenated haemoglobin delivery (Fig 2 b, d). We were looking for bNIRS channels showing 160 the expected patterns of positive associations between EEG and oxCCO and HbO 2 and 161 negative associations with HHb. Figure 3 shows that these associations were primarily 162 concentrated in the beta and gamma bands as predicted (Fig 2 in the supplementary material  163 shows the associations for the theta and alpha bands). Coupled activity was localised to a 164 bNIRS channel (channel 14) positioned over the superior temporal sulcus -temporo-parietal 165 junction region. At this channel, a coupled increase for the social condition and a coupled 166 decrease for the non-social condition was observed (Fig 3 b,

182
Using image reconstruction on the bNIRS data, the spatial sensitivity of the bNIRS location 183 of interest (channel 14) is shown in Figure 4. The method for image reconstruction has been 184 described in detail in the methods section. This indicates that coupled activity was most 185 consistent with the spatial extent of changes in metabolic activity (CCO) and was 186 differentially modulated in the social and non-social conditions. 187 We conducted a multimodal imaging analysis of coordinated neural activation, metabolic 196 demand, and oxygenated haemoglobin delivery in the infant brain. Confirming previous 197 work, naturalistic social and non-social stimuli produce broad haemodynamic changes that 198 can be refined through examining locations with coupled haemodynamic and metabolic 199 activity (19). We and others have also observed broadband differences in EEG responses to 200 social and non-social stimuli (11) that were also observed in the present datasets. However, 201 examining coupling between these two phenomena uncovered a precise pattern in which a 202 specific location at the temporal-parietal junction that differentially responds to both social 203 and non-social stimuli was also coupled with beta and gamma band activity across 204 chromophores in the expected pattern. We contend that this approach allows precision 205 identification of neural specialisation through the coordination of neural, haemodynamic, and 206 metabolic activity. Widespread use of this technique will accelerate our understanding of 207 both the typically and atypically developing brain. 208 209 Our work is consistent with previous studies in identifying increased gamma band activity 210 over temporal and parieto-occipital brain regions during face processing (29-38) (39-42). 211 High-frequency neural firing is associated with localised processing (43) whilst lower-212 frequency activity is associated with larger-scale network changes and transfer of information 213 across systems (44). The increase in lower-frequency activity during social attention also 214 observed here and in other work (11, 45) may support larger-scale connectivity and 215 communication of information through cross-frequency coupling (45). Our work further 216 indicates that measures of metabolic load are a critical nexus in understanding localisation of 217 brain function. Localised high-frequency activity exerts strong metabolic demand (46, 47) 218 and subsequent increases in oxygenated haemoglobin (24,48,49). These increases in 219 metabolic rate are supported by increased activity in the mitochondrial electron transport 220 chain, resulting in the changes in cytochrome-c-oxidase we detected with broadband NIRS. 221 Nitric oxide (which competes with oxygen to bind to cytochrome-c-oxidase) and carbon 222 dioxide (produced as a by-product in the ETC) are key signalling molecule in controlling 223 neurovascular coupling and thus subsequent oxygen delivery (50, 51). Finally, reactive 224 oxygen species produced by the ETC are a key signal in inducing synaptic plasticity (52). 225 Thus, our work is consistent with a model in which social attention induces localised high 226 frequency brain activity in the temporal parietal junction, which increases local metabolic 227 rates, triggering synaptic plasticity and subsequent oxygen delivery to a broader region. 228 229 Our work specifically pinpoints the importance of the temporal-parietal junction in early 230 social brain function. Previous studies measuring haemodynamic activity have identified 231 early sensitivity of this region to social stimuli from at least 4 months (53), alongside a 232 broader network of other regions. Here, we pinpoint this specific location as having coupled 233 neuronal, metabolic, and haemodynamic activity that is modulated in opposite directions by 234 complex social and non-social content. In the adult brain, the temporal-parietal junction has 235 received considerable attention and there are several competing models of its function. It has 236 been linked to mentalising (54, 55) and reorienting attention to behaviourally relevant stimuli 237 (56); it can be viewed as a nexus area where the convergence of attention, language, memory 238 and social processing supports a social context for behaviour ((57) or as a region that is active 239 when awareness of a prediction permits attentional control (58). Intriguingly, recent 240 formulations within the predictive coding framework link the right temporal-parietal junction 241 to a domain-general role in prediction, perhaps representing the precision of priors (59). 242 Predictability has been linked to energy-efficiency, with some computational models showing 243 that energy limitations are the only requirement for driving the emergence of predictive 244 coding (60). Increases in beta/gamma have also been linked to unexpected reward processing 245 (61). Taken together, our results may indicate the early presence of priors for social 246 interaction that are being actively updated (in contrast to the dynamic toys, which may 247 already be more predictable). 248 249 The methods we developed have extensive application in both neurotypical and atypical brain 250 function. Assessing coupling over developmental time will indicate the mechanisms 251 underpinning neural specialisation and constrain theoretical frameworks seeking to explain 252 specialisation in the adult brain. The mechanisms of neurovascular coupling remain unclear 253 in the adult brain (50), and are developing in infancy (17), and novel multimodal and non-254 invasive approaches to their identification could yield significant progress. Computational 255 models could test the role of constraints in energy supply on developing localisation of 256 function. Further, the region identified here also shows atypical haemodynamic 257 responsiveness in infants with later symptoms of autism (62); since mitochondrial 258 dysfunction has become an increasing focus in autism (63) the possibility that atypical 259 coupling may impact specialisation in autism is an important hypothesis to test. Further, our 260 methods have applicability in determining the impacts of early brain injury. Recent work (64) 261 measured both cerebral oxygenation and energy metabolism in neonates with brain injury 262 (hypoxic-ischaemic encephalopathy) and demonstrated that the relationship between 263 metabolism and oxygenation was able to predict injury severity. This therefore provided a 264 clinical, non-invasive biomarker of neonatal brain injury. Indicating applicability across the 265 lifespan (65) simultaneous measurements of cerebral oxygenation, metabolism and neural 266 activity in epilepsy revealed unique metabolic profiles for healthy brain regions in 267 comparison to those with the regions of the epileptic focus. This work demonstrates the 268 strength of combining measurements from multiple modalities to investigate brain states, 269 particularly in clinical populations. 270 271 Our work has several limitations. We used naturalistic stimuli to maximise ecological 272 validity; however, this reduces our ability to probe the function of the temporal-parietal 273 junction across specific stimulus dimensions and this is an important target for future work. 274 Limitations of current technology meant we recorded from the right hemisphere only and 275 thus cannot determine the specificity of our findings to left temporal-parietal junction; 276 engineering advances are required to produce whole-head bNIRS devices. 277 278 Conclusion: Energy metabolism and neural activity are known to be tightly coupled in order 279 to meet the high energetic demands of the brain, both during a task (66, 67) and at rest (68). It 280 has been hypothesised that the level of correspondence between energy metabolism and 281 neuronal activity may be an indicator for brain specialisation (28, 66, 69 Declaration of Interests 320 321 The authors declare that the research was conducted in the absence of any commercial or 322 financial relationships that could be construed as a potential conflict of interest. 323 324 Data availability statement 325 326 The data contains human subject data from minors and guardians provided informed consent 327 to having data shared only with researchers involved in the project, in anonymised form. consisted of a variety of full-colour video clips of actors performing nursery rhymes such as 362 "pat-a-cake" and "wheels on the bus". The non-social videos consisted of dynamic video 363 clips of moving mechanical toys. The visual and auditory components of both social and non-364 social videos was matched. These videos have been used extensively in prior infant studies in 365 both EEG studies (11) and NIRS studies (70, 71). Both social and non-social experimental 366 conditions were presented alternatingly for a varying duration between 8-12 s. The baseline 367 condition consisted of static transport images, for example cars and helicopters, which were 368 presented for a pseudorandom duration of 1 -3 s each for a total of 8 s. Following the 369 presentation of the baseline condition, a fixation cross in the shape of a ball or a flower 370 appeared in the centre of the screen to draw the infant's attention back to the screen in case 371 the infant had become bored during the baseline period. The following experimental 372 condition was then presented once the infant's attention was on the fixation cross. Error! 373 Reference source not found.a depicts the order of stimulus presentation. All infants sat in 374 their parent's lap at an approximate distance of 65 cm from a 35-in screen which was used to 375 display the experimental stimuli. The study began with a minimum 10 s rest period to draw 376 the infant's attention towards the screen during which the infant was presented with various 377 shapes in the four corners of the screen. Following this, the baseline and experimental stimuli 378 were presented alternatingly until the infant became bored or fussy. 379 380 Data acquisition and array placement: bNIRS and EEG data was acquired simultaneously 381 and the bNIRS optodes and EEG electrodes were positioned on the head using custom-built, 382 3-D printed arrays which were embedded within a soft neoprene cap (Neuroelectrics, Spain).  This formed a total of four light sources at the participant-end and each pair of light sources 402 were controlled by a time multiplexing mechanism whereby one pair of light sources was on 403 every 1.4 s. The system also consisted of fourteen detector fibres at the participant-end which 404 were connected to two spectrometers, seven for each spectrometer (in-house developed lens 405 spectrographs and PIXIS512f CCD cameras (Princeton Instruments). The configuration of 406 four light sources and fourteen detectors formed a total of nineteen measurement channels. 407 These were positioned over the occipital cortex and the right hemisphere as shown in Figure  408 5a. The source-detector separation was 2.5 cm. For each infant, intensity counts (or photon counts) from each of the fourteen detectors were 422 used to assess the signal-to-noise (SNR) ratio at each channel and the channels with intensity 423 counts lower than 2000 or higher than 40,000 were excluded (72). If an infant had more than 424 60% of channels excluded, they were excluded from the study. At the group level, five 425 channels over the occipital cortex were excluded due to poor SNR in majority of infants 426 (Channel 3 excluded in 64% of infants, Channel 6 excluded in 83% of infants, Channel 7 427 excluded in 64% of infants, Channel 8 excluded in 79% of infants) and one channel over the 428 right hemisphere was excluded in 100% of infants due to a damaged optical fibre. 429 430 EEG: EEG was used to measure neural activity simultaneously to haemodynamic and 431 metabolic activity using the Enobio EEG system (Neuroelectrics, Spain) which is a wireless 432 gel-based system. The system consisted of 32 electrodes, the locations of which are shown in 433 Figure 5b. The sampling rate of the system was 500 Hz. The experimental protocol in 434 Psychtoolbox sent event markers to both bNIRS and EEG systems using serial port 435 communication which was then used to synchronise the bNIRS and EEG. 436 437 All data were analysed using the EEGlab Toolbox (Schwartz Centre for Computation  438 Neuroscience, UC San Diego, USA) and in-house scripts in Matlab (Mathworks, USA). The 439 raw EEG signal was band-pass filtered between 0.1 -100 Hz and a notch filter (48 -52 Hz) 440 was applied to remove artifacts due to line noise. Following this, blocks of the data were 441 created such that they consisted of the baseline period prior to the stimulus presentation and 442 the entire following stimulus period. These blocks were then segmented into 1 s segments 443 such that for both the baseline and the stimulus, each 8 -12 s presentation of the baseline 444 condition or the stimulus condition yielded 8 -12 x 1 s segments. These 1 s segments 445 consisted of 200 ms of the previous 1 s segment and 800 ms of the current segment and the 446 200 ms was used for baseline correction of each 1 s segment. Segments where the infants 447 were not visually attending to the stimulus were removed. Artifacts were detected using 448 automatic artifact-detection in EEGlab and through manual identification. EEG segments 449 were rejected if the signal amplitude exceeded 200 μ V, or if electro-ocular, movement, or 450 muscular artifacts occurred. Channels with noisy data were interpolated by an algorithm 451 incorporated within EEGlab. Data were then re-referenced to the average reference. 452 453 Within each block (consisting of the baseline period and the stimulus period), each artifact-454 free 1 s segment was subjected to a power analysis to calculate the average root mean square 455 (RMS) power for both low and high frequency bands -theta (3 -6 Hz), alpha (8 -12 Hz), 456 beta (13 -30 Hz) and gamma (20 -60 Hz), within each 1 s segment. This then yielded the 457 average RMS power across the block (baseline period + following stimulus period). Baseline 458 correction was performed by subtracting the average of the 2 s of the baseline period from the 459 entire block. RMS power was chosen as the metric to correlate bNIRS and EEG data as 460 previous studies have demonstrated that task-related BOLD changes are best explained by 461 RMS (75, 76). The blocks were then averaged across trials to obtain an averaged RMS 462 response per participant. A portion of the averaged RMS power was then entered into a GLM 463 analysis described below -this consisted of two seconds of the baseline period and 8 seconds 464 of the stimulus period. Figure 6a provides a visual depiction of how the RMS power was 465 derived from the pre-processed EEG data. For each participant, the RMS power was also 466 averaged across the stimulus period for statistical analysis of the EEG data. For each 467 frequency band, statistical t-tests were performed on this averaged RMS power comparing 468 the social condition versus the baseline (RMS power was averaged during the baseline 469 period), the non-social condition versus the baseline and social versus non-social. The false 470 discovery rate (FDR) procedure using the Benjamin Hochberg method (77) was performed 471 to correct for multiple comparisons. 472 473 Data Analysis: Figure 6b outlines the data analysis pipelines for both bNIRS and EEG data, 474 as well as the procedure for the combined bNIRS-EEG analysis.

490
Combined NIRS-EEG analysis: A GLM (78) approach was employed to investigate the 491 relationship between the bNIRS hemodynamic and metabolic data with the EEG neural data. 492 The canonical GLM typically uses a model of the expected haemodynamic response, i.e. the 493 hemodynamic response function (HRF), to predict the hemodynamic signal. However, given 494 the differences in the haemodynamic response in adults and infants, the standard adult HRF 495 model cannot be assumed for infant data. For example, infants display a delay in their 496 haemodynamic responses (79)(80)(81). In addition, the analogous of the HRF is not established 497 for the metabolic response (i.e. the metabolic response function or MRF). Therefore, the first 498 step of this analysis involved reconstructing the HRF for HbO 2 and HHb and the MRF for 499 oxCCO before combing bNIRS and EEG data. 500 501 The reconstruction of the infant HRF and MRF started with block-averaging the HbO 2 , HHb, 502 and oxCCO signals for social and non-social conditions within each baby. Based on our 503 previous study (19), we selected only the channels that displayed statistically significant 504 responses to the contrast task versus baseline. The single subjects block-averaged responses 505 were averaged across the social and non-social conditions and then across the significant 506 channels. The resulting block-averaged responses were then averaged across the group to 507 obtain a "grand average" HbO 2 , HHb and oxCCO response. 508 509 The grand average was then used in an iterative approach to estimate the HRF and MRF that 510 best fit the HbO 2 , HHb and oxCCO responses. This involved fitting the grand averaged 511 signals with different HRF/MRF models starting from the canonical HRF made of two 512 gamma functions and varying the following parameters: 1) delay of response, 2) delay of the 513 undershoot and 3) ratio of response to undershoot to identify the combination of parameters 514 that best reconstructed the infant HRF/MRF for the social/non-social stimuli. The parameters 515 were varied in increments of 1 s such that the delay of the response was varied from 5 s to 15 516 s from the stimulus onset, the delay of the undershoot was varied from 5 to 20 s and the ratio 517 of the response to the undershoot was varied from 2 to 6 s. All possible combinations of 518 parameters were tested. The grand average responses were fitted with each HRF/MRF in 519 GLM approach, and β -values were obtained for each combination of the HRF/MRF 520 parameters. The β -values were entered into a statistical test and the parameter combinations 521 that yielded the highest, statistically significant β -values (i.e. the model best fitting the data) 522 were selected to reconstruct the infant HRF/MRF. This is approach is similar to those used 523 previously to reconstruct the infant HRF (81) and identified the best fit to be with a 2-s delay 524 of response for HbO 2 and HHb and a 3-s delay of response for oxCCO in comparison to the 525 adult HRF (i.e. 6 s). Moreover, the delay of the undershoot was 9-s earlier for all 526 chromophores and the ratio of the response to the undershoot was 2 for HbO 2 and HHb and 3 527 for oxCCO, in comparison to 6 for the adult HRF. The new reconstructed HRF was then used 528 for the GLM approach to correlate bNIRS and EEG data. The process for estimating the HRF 529 and MRF has been depicted in Figure 7.

534
To constrain the analysis, we chose to investigate coupling of haemodynamic and metabolic 535 with neural activity at specific channels. For this, we used the results from an analysis we 536 described previously that combined bNIRS haemodynamic and metabolic signals (19, 27). 537 The results from this identified task-relevant cortical regions that displayed high levels of 538 haemodynamic and metabolic coupling. The bNIRS channels that displayed significant 539 haemodynamic and metabolic coupling for social and non-social conditions were used here. 540 All EEG channels were used as EEG is not as spatially specific as bNIRS. For each infant, 541 for each chromophore, for each channel and each EEG frequency band, the new infant 542 HRF/MRF that was reconstructed in the previous step was convolved with the events to 543 obtain the "predicted" bNIRS signal. The "predicted" bNIRS signal was then convolved with 544 the EEG RMS power block (consisting only of the data from the stimulus period) at each 545 frequency band to obtain the neural regressor for the bNIRS data, considering both social and 546 non-social conditions together. The design matrix thus included the neural regressor 547 reflecting the increased in EEG activity to the social and non-social stimuli and used to fit the 548 bNIRS data. This was performed for HbO2, HHb, and oxCCO individually for all the 549 channels.
β -values were estimated for each channel and t-tests against 0 were conducted to 550 test whether there was a statistically significant association between bNIRS signals and EEG 551 frequency bands. The false discovery rate (FDR) procedure using the Benjamin Hochberg 552 method (77) was performed to correct for multiple comparisons. The FDR-corrected 553 significant t-values were plotted. This method has been used in numerous studies previously 554 in correlating fMRI BOLD -EEG (20). For each frequency band, FDR-corrected, significant, 555 β -values were also averaged (1) for bNIRS and EEG channels over the right hemisphere and 556 (2) between bNIRS channels in the right hemisphere and EEG channels in the left hemisphere 557 to obtain an estimate of the frequency band where bNIRS and EEG activity associated most 558 strongly within hemispheres and across hemispheres. Only bNIRS channels that displayed 559 significant (prior to FDR correction) haemodynamic and metabolic coupling were used for 560 this analysis (as indicated in Figure and Error! Reference source not found.). For the 561 social condition, these were channels 4, 12, 13 and 14 for HbO 2 , channels 11, 12, 14 and 18 562 for HHb and channels 4, 11, 12, 13, 14 and 18 for oxCCO. For the non-social condition, these 563 were channels 5, 12 and 14 for HbO 2 , channels 5, 12, 14 and 16 for HHb and channels 5, 12, 564 14 and 16 for oxCCO. 565 566 For the bNIRS analysis, data from 25 infants was included while for the EEG analysis, data 567 from 35 infants were included. For the joint bNIRS-EEG analysis, only infants that had both 568 valid bNIRS and EEG data for social and non-social conditions were included and therefore 569 17 infants were included in this analysis. 570 571 Image reconstruction: Image reconstruction was performed on the bNIRS data, at the 572 individual subject level and then averaged across infants to produce a grand average that is 573 shown in Figure 4. For this analysis, three additional long-distance channels were created 574 over the right hemisphere (in addition to the 19 bNIRS channels) that had a source/detector 575 separation of 4.3cm. 576 577 For this analysis, the block averaged attenuation changes at 13 discrete wavelengths (from 578 780 -900 nm at 10 nm intervals) for each infant were selected from the bNIRS data. This 579 was done to reduce the computational burden of the reconstruction (82). A four-layer infant 580 head-model (consisting of the grey matter (GM), white matter (WM), cerebrospinal fluid 581 (CSF) and extra cerebral tissue) was built using averaged MRI data from a cohort of 12-582 month-old infants presented in Shi et al. (83). The Betsurf segmentation procedure (84) was 583 then used to define an outer scalp boundary from the average head MRI template. The 584 voxelised four-layer model was converted to a high-resolution tetrahedral mesh (∼7.8 × 585 10 5 nodes and ∼ 4.7 × 10 6 elements) using the iso2mesh software (Fang & Boas, 2009). The 586 same software was used to create the GM surface mesh (∼5.8 × 10 4 nodes and ∼ 1.2 × 10 5 587 faces), used to visualise the reconstructed images. 588 589 The reconstruction of images of HbO 2 , HHb and Δ oxCCO are described elsewhere (85), 590 using a multispectral approach (86). Wavelength-specific Jacobians were computed with the 591 Toast++ software (87) on the tetrahedral head mesh and projected onto a 50 × 60 × 50 voxel 592 regular grid for reconstruction, using an intermediate finer grid of 100 × 120 × 100 voxels to 593 optimize the mapping between mesh and voxel space. Optical properties were assigned to 594 each tissue type and for each wavelength by fitting all published values for these tissue types 595 (88-90). Diffuse boundary sources and detectors were simulated as a Gaussian profile with a 596 2-mm standard deviation, and Neumann boundary conditions were applied. The inverse 597 problem was solved employing the LSQR method to solve the matrix equations resulting 598 from the minimization and using first-order Tikhonov regularization, with the parameter 599 covariance matrix containing the diagonal square matrices with the background concentration 600 values of the three chromophores (23.7 for HbO 2 , 16 for HHb and 6 for Δ oxCCO) (91, 92) 601 and the noise covariance matrix set as the identity matrix. The maximum number of iterations 602 allowed to the LSQR method was set to 50, and with a tolerance of 10 -5 . The regularization 603 hyperparameter λ was set to 10 -2 . 604 605 The reconstructed images, defined on the same regular grid of the Jacobian, were remapped 606 to the tetrahedral head mesh and then projected to the GM surface mesh, by assigning a value 607 to each node on the GM boundary surface that was equal to the mean value of all the 608 tetrahedral mesh node values within a 3-mm radius. The concentration changes for HbO 2 and 609 HHb were normalised to the maximum concentration change of HbO 2 while Δ oxCCO was 610 normalised to its own maximum change in concentration. 611 612