1.1.

Overview of the Pipeline

Figure 1 shows the steps to be followed using this pipeline. These steps are detailed below. Briefly, the first step handles preparing and cleaning the digitization data, i.e., estimates of the three-dimensional (3D) positions of the fNIRS optodes on the head. Next, we discuss how to prepare the MRI images. Here, we focus on preparing individual MRI images (see left, top in Fig. 1). For studies using an atlas that already contains a segmented head (see right, top in Fig. 1), the user can skip to setting up and running Monte Carlo simulations to create the light model. Where users have a template MRI which has not yet been segmented, they may wish to use the segmentation steps described in this pipeline, but bypassing the cleaning steps, similar to a user in possession of individual MRIs without background noise in a good orientation. Nonetheless, it is recommended that users read the steps on preparing individual anatomical files to determine whether their files meet the formatting requirements.

Fig. 1

A flow diagram demonstrating the use of this pipeline. Steps on the left indicate steps to be taken with individual MRI images, steps on the right indicate steps to be taken using an existing template image.

Once the MRI data are prepared, the next step is to create the light model. There are multiple software packages that handle this step; we describe one implementation using Homer2¹⁶ ( https://homer-fnirs.org/). It is noteworthy that the pipeline is also fully compatible with Homer3 as well. Then, fNIRS data must be pre-processed. Again, there are multiple options here, and we describe a solution using Homer2. The next step is image reconstruction. Here, we use tools from the NeuroDOT package¹⁷ ( https://github.com/WUSTL-ORL/NeuroDOT_Beta). We then run an individual-level GLM on the data, again, using NeuroDOT tools. Finally, we describe how to transfer the data into a common group space for a group-level analysis. We demonstrate the final data produced by the pipeline using recent data from a study looking at the development of visual working memory (VWM) in children. All code and full documentation for the pipeline steps can be found on the GitHub sites described in each section below.

2.1.

Collecting Digitization Data

To locate the optodes on the head, we collect digitizations with a Polhemus Patriot FCC Class B digitizing device ( https://polhemus.com/scanning-digitizing/digitizing-products/). The device consists of an electromagnetic receiver, placed on top of the head, and a wand used to mark the placement of the tip at each optode. The device default is set to capture the locations in inches from the receiver, but we changed this to centimeters to be compatible with Homer2^14,16 (which we use in subsequent pipeline steps). Typically, digitizations are assumed to be in the format of one point per row, where the first set of points are reference landmarks on the head, followed by the sources, and then the detectors. Figure 2 shows a typical digitization captured with a small amount of error, pulling points out of alignment (in this case due to reaching across the head).

Fig. 2

(a) The full cap layout. (b) A standard digitization with some error on the part of the experimenter, where the area shown is the highlighted region of (a). Detectors 6 and 7 (blue) have been digitized out of alignment with sources D, E, and F (red). c. The same digitization following three-point correction. It is noteworthy that the layout in (c) is the expected layout, matching the highlighted region in (a).

2.2.

Creating Digitization Templates

One of the challenges of using a digitizing pen is that errors in positioning can arise. For instance, in Fig. 2, several points are slightly mis-localized, resulting in a distorted geometry. This can be caused by a slight slippage in the digitizing pen when the participant moves, inducing a slight mis-localization of the pen tip. Such minor errors are difficult to identify on the fly, particularly when the head being digitized is on the body of a squirmy child. To combat errors in mis-localization of points, we developed a package for the R language ( https://www.r-project.org/), digitizeR ( www.github.com/samhforbes/digitizeR), which takes digitized input from multiple participants for a given head size, creates a digitization template, and identifies and fixes mis-localization errors in individual caps.

In particular, once collected, the digitizations for each participant can be read into digitizeR by cap size. For ease of visualizing caps as they change, we also recommend RStudio. The digitizeR package will read in the Polhemus output files (saved in ASCII text format) and create a list for each cap size. Next, digitizeR finds the best caps to align the other caps to for each cap size. This is done using the Kabsch algorithm,¹⁸ a method of finding the optimal rotation between two matrices to reduce the root mean squared deviation. DigitizeR uses this algorithm to align each cap to all the other caps in that size. Since every digitization runs the risk of participant movement throwing off the values for any given point (or the entire cap), each cap is rotated to match every other cap. The assumption here is that all the caps that were digitized well will all align with each other, and the caps that were digitized poorly will fail to align. So, for example, for three caps A, B, and C, if caps A and B are good caps, and C is a poor digitization, A will align with B within a given 3D criterion distance and B will align with A, whereas C will align with neither A nor B. The next step is to use this information about which caps align to create a template. In our hypothetical situation, two caps, A and B, align with each other within the criterion distance, so the template will be defined by calculating the 3D mean of each corresponding point in A and B. Using this process, a template cap is created for each cap size.

Figure 3 shows a template created from a group of caps. In this instance, the template digitization was created from a group of caps with a head circumference of 48 cm. Caps were aligned using the above method where the 3D distance threshold was chosen as 10 cm. This took the reference cap with the most other matches at that cap size where no point was allowed to exceed a 3D distance of 10 cm from the reference; in this instance, the aligned caps had a mean distance of 3.91 cm from the reference cap. The mean of each of the points across all of the matching caps was then calculated to create the template.

Fig. 3

A 3D image of a template created from a group of caps, from digitizeR. Points in blue are head landmarks, points in red are digitized points from the fNIRS cap. Landmarks have been labelled. (a) The cap from the left-hand side; (b) from the front; and (c) from the top.

Once created, templates can be saved in digitizeR and used for all participants in that cap size if desired. This would yield an average acceptable digitization that is unique to each cap size, but captures the data collected across the sample. Alternatively, as we describe in the next section, digitizeR allows using the template to clean and save modified versions of the individual digitizations, allowing the user to preserve as much of the individual-level variation in optode location and cap placement as possible.

2.3.

Correcting Individual Caps

Once templates have been constructed, each participant cap can then be aligned to the template for that cap size, again using the Kabsch algorithm. There are two sources of variation in the digitized locations in 3D space: true variation in cap placement and optode location due to factors such as wrinkles in the cap or the cap being off-center; and errors such as digitization of the incorrect point or movement causing the point to be unrealistically far away from the head. The goal of the correction at the individual cap level is to permit as much of the former variance as possible, which is due to the true placement of the optode on the head, while eliminating the large sources of the latter type of variance.

The digitizeR package allows three different correction methods: a three-step alignment correction method, an iterative method, and a head-wise method. The three-step method allows the user to specify three threshold 3D distances, at which the points will be sequentially tested against the template for that cap size, and those outside those distances will be removed, and then the cap re-aligned without those points. As an example, a three-step alignment on a digitization with distances set at 12, 8, and 6 would first exclude points beyond a 3D distance of 12 (cm), re-align without those points and do a second pass with a threshold of 8 cm, re-align without those points, and then do a third pass with a threshold of 6 cm, removing the points beyond this threshold. After the final distance check, the missing points are replaced with points from the relevant template. The difficulty with the three-step method is selecting the three distances, and this can be dependent on the data that is being used. The data that were cleaned using the template in Fig. 3 were cleaned using threshold distances of 12, 10, and 7 cm, with the end goal of making sure no more than 1/3 of the original data points had been replaced using this method.

The iterative approach takes the same method as the three-step approach, but rather than the user-specified distances, it starts with a 3D distance of 20 cm and replaces points outside of that distance, and then reduces the difference by 1 cm, repeating the procedure of aligning to the template and removing the points outside of the 3D distance each time until the points replaced reach the user-specified proportion of replaced points, or until the points are within half a centimeter of the template. The head-wise method follows a similar logic, but deals with entire caps rather than individual points, so that rather than individual points being replaced, the mean 3D distance from all the points in the cap is calculated against the template, and caps exceeding the threshold distance are replaced with a complete template cap, otherwise, the cap remains as it was.

In this paper, we used the three-step method which uses the R function threestep_alignment(), although the best approach may be dependent on the types of errors being seen in the data. In the case of our dataset, we wanted to minimize the number of points we were replacing while ensuring our caps were as correct as possible. While some trial and error can be necessary to determine the best distances for the caps at hand, the three-point method is a good option for this purpose. The points that were removed as part of this procedure are then filled in with the equivalent point from the template once the cap is aligned to make a complete cap that is aligned in the same space as the template cap.

We tested the three-step method by taking two templates, loading them in as caps and moving the first point [the landmark nasion (NZ)] much further away from the other points. We then used the original templates to repair them. The test of this method is that the three-step alignment correction should correct that point on each of the two templates, but no other points, as all other points derive from the template itself. When this is tested with the three-step method, only these two bad points are corrected, having been replaced with the original template points, and all other points remain the same.

2.4.

Preparing for Use with Homer2/ Homer3

The final versions of the digitized caps after alignment to the template can be saved and exported from R using the built-in functions from the digitizeR package. The save_caps option automatically saves each of the caps in its own folder as a digpts.txt file, saved in the format that is necessary for Homer2/Homer 3, which includes reference, source, or detector labels for each point. In addition, the caps are converted by default to distances in millimeters (from centimeters), making them automatically ready to use for import into Homer2 or Homer3. Full documentation for any of the digitizeR functions can be found at www.github.com/samhforbes/digitizeR.

6.1.

Overview of Image Reconstruction in NeuroDOT

First, the sensitivity volumes for each source-detector pair that were calculated in AtlasViewer were converted into NeuroDOT format. In NeuroDOT, the sensitivity volumes are arranged in a single two-dimensional matrix $A$ with measurements in the first dimension and voxels in the second dimension. A structure variable info contains sub-structures that house meta-data describing the sensitivity volumes (info.tissue), the optode locations (info.optodes), the measurement list (info.pairs), the stimulus paradigm (info.paradigm), and other parameters from data acquisition (info.system). The header that contains the spatial meta-data from the NIFTI *.nii files from AtlasViewer is contained within info.tissue.infoT1. The order of the measurements in $A$ corresponds to that described in the measurement list located in info.pairs. The voxels are ordered as if the entire sensitivity volume is unwrapped into a single vector. To maximize computational efficiency, only voxels with a sensitivity above a given threshold for any of the measurements are included. The indexing of these voxels relative to the entire volume is described in the variable info.tissue.dim. Good_Vox. This step typically saves around a factor of 5 to $10\times$ in disk space and lowers computation time due to lower demand on RAM. For this study, the threshold was set to 0.01 of the maximum sensitivity within the combined volumes of all measurement sensitivity profiles.

Second, the preprocessed optical density data from Homer2 in the *.nirs file (or *.snirf file in the case of Homer3) was converted directly into NeuroDOT format. The NeuroDOT measurement list represented in the table info.pairs is populated from the SD variables in Homer2 (or snirfdata.data.measurementList in Homer3). The table info.pairs contains the source, detector, wavelength, source-detector separation in 2D and 3D, and other information for each measurement. The measurements to be used for reconstruction are listed in the info.MEAS.GI variable. The stimulus paradigm timing information is populated into info.paradigm for use in the post-reconstruction analyses (described below).

To minimize contributions from extra-cerebral tissue in the variance of the measurement data, we performed global signal regression separately for each wavelength. The mean signal across all low-noise measurement channels was regressed from all channels in a manner akin to superficial signal regression as previously described using the NeuroDOT function regcorr.^26–28 Though short measurement channels were not utilized in these data, the mean across all measurement channels provides an estimate of systemic physiology present in all of the data. Removing this global mean increases the contrast to noise and spatial specificity of the variance remaining in the measurement data. The conversion of measurement-wise data to voxel-wise data through reconstruction significantly increases the overall size of the optical data due to the presence of many times more voxels than measurements. To help alleviate the large computational demands, the data were temporally down-sampled to 25 Hz before reconstruction using the NeuroDOT function resample_tts.

Once the optical, sensitivity, and meta-data have been converted into NeuroDOT format, the sensitivity matrix can be inverted separately for each wavelength using the Moore–Penrose generalized inverse with a Tikhonov regularization parameter of ${\lambda }_1=0.01$ and spatially variant regularization parameter of ${\lambda }_2=0.1$ using the NeuroDOT function Tikhonov_invert_Amat.²⁹ The Tikhonov regularization parameter tunes the balance between spatial smoothness and high-spatial frequency noise in the image. Spatially variant regularization helps control for the tendency for diffuse optical tomography (DOT) methods to reconstruct too superficially and is set to optimize spatial correspondence with fMRI based on other studies.³⁰ This process minimizes the objective function: $\min \{{\Vert y_{meas}-Ax\Vert }_2^2+{\lambda }_1{\Vert Lx\Vert }_2^2\}$. We include ${\lambda }_2$ as a spatially variant regularization term, such that: $diag(L)=\sqrt[2]{[diag(A^{T}A)+{\lambda }_2]}$.^26,30,31 Given these regularization steps, the direct inverse that provides the solution, $x=A_{{\lambda }_1{\lambda }_2}^{\#}{y}_{meas}$ is given by the Moore–Penrose generalized inverse: $A_{{\lambda }_1{\lambda }_1}^{\#}=L^{-1}{({\tilde{A}}^T\tilde{A}+{\lambda }_{1}I)}^{-1}{\tilde{A}}^T{y}_{meas}$ with $\tilde{A}=\tilde{A}{L}^{-1}$. Regularization parameters of ${\lambda }_1=0.01$ and ${\lambda }_2=0.1$ work well for a broad range of geometries and array designs. These steps are all carried out within Tikhonov_invert_Amat.m. We also recommend using at least a modest amount of direct spatial smoothing using a Gaussian smoothing kernel to control for speckly noise in the images. Recommended smoothing kernel sizes range from 3 to 10 mm at their full width at half max. The optical data are reconstructed into the voxelated space using the NeuroDOT function reconstruct_img.

Relative changes in oxygenated (${\mathrm{HbO}}_2$), deoxygenated (HbR), and total hemoglobin (HbT) concentrations are obtained from the reconstructed absorption coefficients by spectral decomposition of the extinction coefficients of ${\mathrm{HbO}}_2$ and HbR at the two wavelengths.^32,33 Conversion from relative changes in absorption coefficients to ${\mathrm{HbO}}_2$ and HbR is performed with the NeuroDOT function spectroscopy_img.

6.2.

Sample Image Reconstructed Data

Figure 5 below shows samples of image reconstructed fNIRS data from one participant plotted against the processed channel-based fNIRS data. To create these images, we first placed the centroid of a 2-cm diameter sphere at the maximum value of each channel’s sensitivity volume from the photon migration simulations for each chromophore. Note that we constrained the sensitivity volumes to the cortex to ensure that the maximum value was located in the cortex and not in the skull or surface tissues. Next, we extracted the mean reconstructed fNIRS time series from each sphere, averaging over voxels within the sphere. The red curve in each panel shows the mean image-reconstructed time series, while the black line shows the channel-based data from the associated channel. We show data from six channels in total [see brain images in Figs. 5(e) and 5(f) for mapping of channels to spheres]. We selected these channels because they are located in parietal and ventral occipital regions of interest for our sample data set. As can be seen in the figure, the image-reconstructed data match the channel-based data well for both ${\mathrm{HbO}}_2$ [Figs. 5(a) and 5(b)] and HbR [Figs. 5(c) and 5(d)]. The signals are within the same range, with the image-reconstructed data typically showing a bit more signal intensity.

Fig. 5

Image reconstructed fNIRS data (red lines) plotted against channel-based fNIRS data (black lines). Left panels show HbO time series (in $\mu \mathrm{M}$) for (a) channels 27–29 and (b) channels 32–34; middle panels show HbR time series for (c) channels 27–29 and (d) channels 32–34. (e) The locations of the spheres used to extract data for channels 27–29. It is noteworthy that different colors are used to show the channels across views in (e) as the axial slice does not match the coronal and sagittal views. (f) The locations of the spheres associated with channels 32–34 (axial, coronal, and sagittal views are all aligned).

6.3.

Verifying the Image Reconstruction Approach

To verify the image reconstruction approach, we correlated the average image-reconstructed time series (e.g., red lines in Fig. 5) with the channel-based time series (e.g., black lines in Fig. 5) for all data run for four sample participants. In total, we conducted 504 correlations: seven runs (over four participants) × 36 channels × 2 chromophores (${\mathrm{HbO}}_2/\mathrm{HbR}$). Data from 18 channels were discarded because the channels were pruned in the initial fNIRS processing. The correlation results are shown in Fig. 6. Most of the correlations were quite high and above our minimum acceptable threshold of 0.25. In particular, 473 of 496 correlations were $>0.25$; the mean $r$ value for this subset was 0.8. Thus, only 2.6% of the sample was below the criterion. In these cases, the channel-based fNIRS data were of low quality with poor SNR.

Fig. 6

Verifying the fNIRS image reconstruction data. Each dot shows the correlation between the image-reconstructed fNIRS data averaged over a 2-cm sphere in cortex [see examples in Figs. 5(e) and 5(f)] and the associated channel-based data. Red dots show correlations for HbO; blue dots show correlations for HbR. Data are shown for four participants (see columns). Three participants (06IND012B, 06IND118G, and 06IND239G) had multiple runs, so there are four correlation values per channel (HbO/HbR for run1 and run2). Dashed line shows our target criterion value of 0.25. The vast majority of comparisons were above the criterion value.

9.1.

Participants

Data from seventy-four 4.5-year-old children (34 females, $M$ age = 53.5 months, and $SD=1.3$) who participated in a longitudinal study in Scotland were used to verify the analyses pipeline. The study was conducted in the participants’ homes. Parents gave written informed consent and children gave verbal assent prior to the testing session. The research was approved by the General University Ethics Panel (GUEP 375) at the University of Stirling.

9.2.

Exerimental Task

VWM performance was measured as children took part in a color change detection task.³⁷ The experimenter explained the task to the children using $3\times 3\text{inch}$ flash cards. Once the experimenter was satisfied that the children understood the rules, they moved to the experiment. The experimental task was run in E-prime V.3 software on an HP laptop with a 14-inch screen. In each trial, a fixation appeared on either side of the screen to orient the child to the screen. Then, a memory array of colored squares was presented for 2 s, followed by a blank screen for 1 s (delay period). Then, a test array of colored squares was presented and stayed on the screen until a response was made. After the presentation of the test array, the child was asked whether the colored squares were same or different between the memory and test arrays. The experimenter registered the child’s response using the keyboard. Each trial was followed by an inter-trial interval of 1 s (50% of the trials), 3 s (25% of the trials), or 5 s (25% of the trials). VWM load was manipulated from 1 to 3 items (load 1, load 2, and load 3). There were eight trials for the same condition and eight trials for the different condition for each load. Thus, there were 12 possible conditions for trials in this task, three loads (1, 2, and 3) × 2 trial types (same and different) × accuracy (correct and incorrect).

9.3.

fNIRS Data Acquisition

An NIRX system with eight sources and eight detectors was used to collect neuroimaging data from children as they engaged in the VWM task. Data were collected at 7.81 Hz and the wavelengths of the sources was 760 and 850 nm. Probe geometry with fourteen channels overlaying the bilateral frontal and parietal areas [see Fig. 7(a)] was constructed based on regions of interest from previous fMRI VWM literature.³⁸ The head circumference of the child was measured and a suitable cap was chosen. Four cap sizes were used across all children and source-detector distances were scaled according to cap size (50 cm: 2.5 cm, 52 cm: 2.6 cm, 54 cm: 2.7 cm and 56 cm: 2.8 cm). The distances from the NZ to inion and between the peri-auricular points were measured to ensure that the cap was centrally placed on the child’s head.

Fig. 7

(a) Probe geometry for the task. (b) Significant clusters from the image-reconstructed fNIRS data in the middle frontal gyrus and left angular gyrus (see red clusters in brain image). Time series plots show weighted average ${\mathrm{HbO}}_2$ (solid lines) and HbR (dashed lines) for load 1 (light gray), load 2 (medium gray), and load 3 (black) conditions. (c) Block-averaged activation from the channel-based fNIRS data showing channels overlapping with the red clusters in (b) (see highlighted yellow region in brain image). Plotting conventions as in (b).

9.4.

fNIRS Data Processing

The probe geometry was digitized to obtain the coordinates of sources and detectors for each of the four capsizes. A 4.5-year-old MRI atlas was obtained from the Neurodevelopmental MRI database^39–43 and segmented into four tissue types (scalp, CSF, gray matter, and white matter). Monte Carlo simulations were run with 1 million photons to generate sensitivity profiles for channels for each capsize (step 2, see Fig. 1). The head volume and sensitivity profiles were converted into NIFTI format. Next, channel-based fNIRS raw data were pre-processed, corrected for motion artifacts, filtered, and converted to optical density units. Data were also corrected for superficial scalp artifacts using global signal regression (step 3). These channel-based fNIRS optical density data were integrated with the volumetric sensitivity profiles using image reconstruction as described previously (step 4).

A GLM with 6 regressors (for 12 conditions) was run on each voxel by convolving a modified gamma function from SPM (delay of response = 4; delay of undershoot = 15, dispersion of response = 1; dispersion of undershoot = 1; ratio of response to undershoot = 6; onset = 0; length of kernel = 16)⁴⁴ with a boxcar window of duration 4.5 s (step 5). Note that this gamma function was chosen after visually comparing the duration of response onset, time to peak, and response offset of the generated HRFs and conventional channel-based block averages. Further, this boxcar window was chosen to account for a sample array of 2 s, delay of 1 s and the rest of the time dedicated towards the test array period. Beta coefficient maps were obtained for all 12 conditions in each voxel for each chromophore and participant. Next, these individual-level beta maps were registered to the MNI space (step 6). For Group analyses, the beta maps were entered into a linear mixed effects model (using 3dLME function in AFNI) with within-subjects factors of load, trial type, and chromophore (HbO and HbR). Note that only correct trials were used for our analyses. For the current case-study, we focused on the interaction between load and chromophore to examine whether HbO activation in significant clusters show the canonical trend of increasing activation with increasing VWM load. Only those voxels that contained data from 60% of the subjects were included for further analyses. To control the family-wise error in our data, we used the 3dFWHMx function in AFNI to estimate the empirical ACF in our fNIRS data and fit a mixed ACF model to this function. We then used the mixed-ACF parameters (0.9939 8.3354 5.2485) in 3dClustSim with a voxelwise $p=0.01$, $\text{alpha}=0.05$, and 10,000 iterations. We selected two-sided thresholding with the NN1 option (first-nearest neighbor clustering). The cluster size criterion was 169 voxels.

9.5.

Results

Two clusters showed a significant interaction between load and chromophore: left middle frontal gyrus and left angular gyrus. Figure 7(b) shows the two clusters in red. Block averaged change in HbO/HbR activation extracted from these clusters is also shown. In the left angular gyrus, HbO activation increased with increasing load. In the left middle frontal gyrus, activation was highest for load 2. Activation in both regions show the canonical difference between HbO and HbR activation, with the largest difference visible at the higher loads. Masked sensitivity profiles from channels 5, 6, 7, and 8 in the left frontal cortex and 12, 13, and 14 in the left parietal cortex are shown in yellow. Block-averaged activation from these channels is shown in Fig. 7(c). We see similar trends in activation with increasing load between the block averages from the significant clusters from the image reconstruction approach and block averages from the conventional channel-based approach. Image reconstruction provides more precise localization of the source of activation in the brain. The change in activation with increasing load is more pronounced when localized to the left angular gyrus compared to the changes in activation in the three channels. It is important to note that the goal of employing image reconstruction techniques is not to empirically compare against and arrive at the same result observed in channel space. Instead, it affords the opportunity to obtain spatial localization of activation from multiple channels to specific clusters and compare these activation maps across participant groups and ages while accounting for key sources of spatial variance in the data.