A compact time-domain diffuse optical tomography system for cortical neuroimaging

Neurocognitive and psychiatric disorders are highly prevalent and treatment expenditures are projected to increase in the coming decades (Mattke et al., 2023; World Health Organization, 2022). Many emerging treatments for these disorders leverage novel mechanisms of action from scientific research for improved clinical results. However, clinical decision making around diagnosis, treatment selection, and treatment monitoring for these complex disorders remains a challenge. Quantitative brain data may be able to help clinicians make informed decisions for personalized care. As an example, evidence is building that non-invasive neuroimaging methods such as functional magnetic resonance imaging (fMRI) (Arns et al., 2023; Cole et al., 2019; Hack et al., 2023; Poirot et al., 2024) or electroencephalography (EEG) (W. Wu et al., 2020) may be instrumental for tailored depression treatments. In their current implementation, non-invasive neuroimaging assays are typically expensive and inaccessible, and may be contraindicated for certain patients. Methods to quantify the brain in a simple, quick, informative way at a point-of-care clinic are urgently needed.

Wearable functional near-infrared spectroscopy (fNIRS) systems have the potential to fill this gap (Li et al., 2023). In the past 15 years, fNIRS systems have become more comfortable and capable, with higher spatial resolution, multimodal measurement capabilities, and standardized data collection and analysis procedures (Ayaz et al., 2022). With these improvements, wearable fNIRS systems enable measurements in situations inaccessible or inconvenient for other forms of neuroimaging, including during free ambulation (Bishnoi et al., 2021), social interactions (Zhao et al., 2024), and at point-of-care clinics.

While these improvements have allowed for an increased range of measurement situations, the underlying technology of most wearable fNIRS systems has remained the same: continuous wave (CW) light. CW systems have the benefit of being relatively inexpensive and straightforward, although systems designed for high-density diffuse optical tomography (HD-DOT) have a large number of optical sources and detectors which can increase cost and system complexity, as well as reduce sampling frequency (Eggebrecht et al., 2014). Alternatively, time-domain (TD-)fNIRS can go beyond traditional HD-DOT systems and interrogate the underlying tissue in greater detail due to its higher depth sensitivity (Ortega-Martinez et al., 2022). TD-fNIRS systems use short pulses of light and detectors capable of measuring single photons to capture a distribution of times of flight (DTOF) of photons. This fine-grained measurement capability allows TD-fNIRS systems to measure absolute optical properties of the underlying tissue including the absorption (μ_a) and reduced scattering (μ_s′ ) coefficients, instead of only measuring relative changes in light intensity as in CW systems. The increased information from TD-fNIRS systems also allows for advanced signal processing methods to more heavily weight photons that arrive late in the DTOF in order to emphasize data from deeper in the head (i.e., from the brain). Despite these benefits, TD-fNIRS systems have limited commercial availability, and the available systems have slow sampling frequencies or limited coverage over the head.

Our team has been working to overcome the technical challenges involved in developing a portable, whole-head coverage, high-density, fast, and scalable TD-fNIRS system. We previously published on the characterization and validation of our prototype system, Kernel Flow1 (Ban et al., 2022). While the prototype system was used in several scientific studies (Castillo et al., 2023; Dubois et al., 2023; Kamar et al., 2024; Ortega-Martinez et al., 2022; Oveisi et al., 2024), it had some limitations: there were gaps in spatial coverage over the head; pre-recorded instrument response functions (IRFs) were not sufficiently stable in time to ensure accuracy of absolute metrics and the system had high-power requirements. As such, we had confined our prior analyses to those performed in channel space and did not fully explore Flow1’s whole-head DOT capabilities.

In this paper, we present the first commercially available whole-head coverage time-domain diffuse optical tomography (TD-DOT) system, Flow2 (Kernel, Culver city, CA, USA). We first show that Flow2 compares favorably with limited channel count research-grade devices on key figures of merit agreed upon by the field, while extending the field-of-view to cover the whole head. We then demonstrate how it can be used to measure absolute concentrations of oxygenated and deoxygenated hemoglobin in the brain (during a breath-hold challenge), and how it can be used to reconstruct focal brain activity (during a sensory and a motor task). Improvements in reconstruction accuracy for TD versus CW data have been theoretically discussed (Arridge & Schweiger, 1997; Selb et al., 2014), and empirically demonstrated with limited coverage systems (Selb et al., 2006). Here we empirically show the benefit of whole-head TD over CW systems. In light of these characterization results and the scalability-by-design of the Flow2 device, we argue that it is well positioned to enable the translation of research findings from both fNIRS and fMRI literature to clinical settings.

The system consists of 40 modules that are arranged in a headset design (Fig. 1A, B). Modules are organized into seven headset superstructure plates that cover the frontal, parietal, temporal, and occipital cortices. Each Kernel Flow2 optical module consists of 3 dual-wavelength sources and 6 detectors located on a 13.5 mm radius from module center. An additional detector located in the center of the module continuously measures the IRF. The 3 source emission points are offset 120 degrees, with a detector located 37 degrees on either side of each source point (Fig. 1C).

Intra-module channel distances are 8.5 mm (6 source–detector pairs), 17.9 mm (6 source–detector pairs), and 26.5 mm (6 source–detector pairs), for a total of 18 dual-wavelength channels within a module. Cross-module channels also provide data, for a total of 2,565 possible channels with a source–detector separation of ≤50 mm over the whole head (3,583 ≤ 60 mm). The actual number of usable channels depends on the light attenuation for a particular participant. Module locations are fixed within a plate, and plates are held together with an adjustable tensioning system. The overall size of the headset fits a range of adult head sizes, with a minimum size of 52-cm circumference and a 32.5-cm bitragion coronal arc. The system weighs 2.4 kg. By comparison, a single detector channel wearable system was reported to weigh 2.5 kg (Lacerenza et al., 2020), although that system includes the weight of the battery while the Kernel Flow2 system is powered over universal serial bus type-C (USB-C) using the USB power delivery (USB-PD) standard.

Each module in the system consists of three major subassemblies: laser assembly, detector assembly, and the optical assembly. All three of these subassemblies are shown together in an exploded view in Figure 1D. The details of each subassembly and the overall system architecture are presented in the following sections.

The Flow2 system has a hierarchical architecture with electronics and wiring harnesses integrated into the headset. The system is cabled over a single USB-C interface that both supplies power (USB-PD) and enables bidirectional communication (USB 2.0) between the data collection computer and the Flow2 system.

The USB-C cable connects to the Flow2 system through the hub sub-assembly that includes an application processor (AP), and a global reference clock. In addition, the hub sub-assembly integrates four electroencephalography (EEG) analog-to-digital converter (ADC) channels that are designed for connecting to active dry electrodes. The hub also handles the primary power negotiation for the USB-PD standard and distribution of power to the rest of the system. Connected to the hub are follower boards that serve as data aggregation points for clusters of modules. Each of these follower boards include a low-power field-programmable gate array (FPGA), additional power conditioning circuitry, a 9-axis IMU, and local USB 2.0 communication interface between the follower and hub.

In total, the Flow2 system supports connection of up to 40 time-domain optical modules and includes 4 active dry EEG channels (approximate locations from 10–10 montage: AF4, AF3, FCz, CPz). The system may be operated with fewer than 40 optical modules, allowing for the removal of unnecessary modules to reduce headset weight or cost. In the current study, we opted to utilize the complete headset configuration to depict activation and deactivation regions across the entire brain during different tasks as described in the following sections.

Each module source location (three per module) has two pulsed edge emitting laser diodes. The paired laser diodes are placed diametrically opposite each other and emit light into a combined micro-prism and lens optical element. The lasers are driven by custom-designed pulse shaping circuitry that efficiently generates laser pulse widths that are <= ~225 ps (690 nm has a pulse width of less than 205 ps, and the 905 nm laser has a pulse width of less than 250 ps). The lasers operate in gain-switched mode, which enables the production of optical pulses that are shorter than the electrical pulse that drives them. In addition to the laser driver circuitry and laser diodes, this sub-assembly contains the power conditioning circuits for the module, as well as the connectivity of the optical module to the follower board which serves as a data aggregation point for clusters of modules.

The detector sub-assembly comprises seven detector application-specific integrated circuits (ASICs) custom designed by Kernel, specifically optimized for conducting time-of-flight measurements in diffuse optical tomography. Notably, the Kernel Flow2 ASICs feature integrated time-to-digital (TDC) circuitry with the photodiodes, ensuring a seamless fusion of key functionalities.

A novel custom dual band-pass coating has been developed and applied to the detector ASIC package. This specialized coating serves as a selective filter, tailored to discriminate against undesired wavelengths within the spectrum, thereby optimizing the detection of photons within the desired spectrum wavelengths of the Flow2 lasers. The implementation of this coating helps to address the challenges associated with ambient light and extraneous signals, ensuring a maximal signal-to-noise ratio and heightened accuracy in TD-fNIRS measurements.

Photon arrival times, derived from on-chip TDCs, are systematically aggregated into histograms and communicated via a serial peripheral interface (SPI) bus to a dedicated FPGA. The synchronization of ASICs across all modules to a 20-MHz global reference clock facilitates the recording of temporally aligned signals between any source and any detector in the system. Furthermore, each detector ASIC incorporates dedicated high-voltage bias circuitry, optimizing the bias for each detector within the system for stable operation at any temperature. Additionally, new circuit architectures were used to reduce power consumption while improving overall TDC performance compared with the Flow1 design.

Within each module, one dedicated detector at the center of the module assumes the role of an on-board IRF detector. This detector captures light, transmitted via a waveguide from the source optical path, yielding a per-pulse waveform that temporally corresponds to each pulse from each respective wavelength. The programmable integration time for constructing histograms on each detector spans from 1 to 800 ms, with the current configuration set at 3.5 ms. Each histogram collected contains signals from only one wavelength. This means our histogram sampling rate is 285.7 Hz, and considering both wavelengths, the system is able to complete spectroscopic measurements at a rate of 142.9 Hz. To avoid optical crosstalk, all lasers are not enabled at the same time. This temporal multiplexing enables lasers in a 38-state pattern, completing a full cycle of data collection for all modules and wavelengths every 76 histograms, corresponding to a system sampling frequency of 3.76 Hz. In the full headset configuration, one source operates at 7.52 Hz, which is a frequency fast enough to accurately capture pulse rate, a complementary measure to the optical properties. This temporal multiplexing provides the additional benefit of bringing the average power per source down to a level that classifies as a class 1 laser device according to the United States Food and Drug Administration Federal Laser Product Performance Standard Code of Federal Regulations Title 21 Section 1040.10 (US FDA FLPPS 21CFR1040.10).

The module optics (Fig. 1D, E) have been carefully designed for several purposes. The optics conduct laser light from the laser diode sources into the scalp, couple light returning from the scalp to the detectors, conform to the curvature of the head, reduce interference from hair by parting it, isolate detected signal to a single detector, and maintain both source output power and detector measurement intensity independent of compression.

There are a total of 9 spring-loaded light pipes contained within the optical module, with 3 source light pipes having an optical exit area of 1.3 mm x 1.3 mm, and the 6 detector light pipes each with a 3 mm diameter. Each of these nine light pipes is optically isolated from one another in separate cavities to prevent optical cross-talk and signal contamination. Both the sources and detectors have a two-lens optical system that maintains optical intensity throughout the compression range of the light pipe springs.

Per each source, there are two laser diodes that are placed offset from each source optic axis and directed into a custom-built optical assembly. The light from the diodes immediately enters the source optic and is redirected by means of an integrated micro-prism with reflective silver coating. The light continues diverging within the first polymethyl methacrylate (PMMA) optical component until it hits the first aspheric surface.

A secondary source optic is a solid body optical component, also made out of PMMA, and is spring loaded with 5 mm of travel. The secondary source optic collects the collimated light from the first aspheric optic, and homogenizes the light prior to exiting the final surface of the light pipe and passing into the user’s scalp. This secondary source optic is also shrouded with an opaque covering to prevent light from leaking from a source directly into a detector light pipe outside of the module housing.

The receiving optics are also designed to use springs to comfortably conform to the user’s head. The input aperture of the detector light pipes is 3 mm, with a rounded edge on the external interface for additional comfort. We have created another two-lens imaging system, consisting of one plano surface, and three aspheric Fresnel surfaces, in order to keep the received optical intensity constant at the detector, regardless of spring compression.

A critical performance metric of a time-domain optical measurement system is the IRF. The IRF is a measure of the uncertainty associated with each timestamp recorded and is reflected in the histogram of accumulated events, when used in time-correlated single-photon counting (TCSPC) applications, as a smearing or broadening of the desired signal being measured.

Each component in the system plays a pivotal role in shaping the system’s IRF, with the laser and detectors emerging as the predominant influencers. The IRF is primarily shaped by the design specifics of the detector process, the laser driver circuitry, and the laser itself. Moreover, the IRF characteristics of these active components exhibit dependencies on the temperature and operating voltages of the electronics and optoelectronics, factors subject to temporal variations and fluctuations during a measurement.

A crucial element in the updated Flow2 module design is the incorporation of a dedicated reference IRF detector within each module (Fig. 1F). This detector is the exact same design as detectors employed in measuring signals from the scalp but is strategically isolated to capture light directly emitted from the lasers without traveling through tissue. This isolation provides a reliable estimate for both the IRF contribution from the detector (owing to its identical design and similar process as the signal-measuring detectors) and the IRF contribution originating from the laser. This IRF measurement is recorded continuously, at the same rate as that of the other detectors. A similar approach has been implemented successfully in a frequency-domain imager (Ban et al., 2016).

The performance of a novel system must first be evaluated, with respect to systems of its class, using established protocols and well-defined figures of merit (Yücel et al., 2021). Three characterization protocols have been internationally accepted for TD diffuse optical systems: (1) the Basic Instrumental Performance (BIP) protocol, (2) the Optical Methods for Medical Diagnosis and Monitoring of Diseases (MEDPHOT) protocol, and (3) the Noninvasive Imaging of Brain Function and Disease by Pulsed Near Infrared Light (nEUROPt) protocol (Pifferi et al., 2005; Wabnitz, Jelzow, et al., 2014; Wabnitz, Taubert, et al., 2014).

The BIP protocol is designed to assess basic hardware performance of time-domain instruments such as detector responsivity, differential non-linearity (DNL), afterpulsing, the system IRF, and system stability (Wabnitz, Taubert, et al., 2014). The protocol is applicable to instruments based on pulsed laser sources with repetition rates of the order of several tens of MHz, fast single-photon detectors, and time-correlated single-photon counting (TCSPC).

The detector responsivity metric in the BIP protocol measures the relative sensitivity of light detection for the system, calculated as the ratio of measured photons exiting from a calibrated phantom to the input illumination (Wabnitz, Taubert, et al., 2014). The calibrated phantom is a diffuse medium, to ensure that the measurement accounts for collection efficiency of the optics as well as the detection efficiency of the detector. The experimental setup consists of a calibrated phantom illuminated on one side with a collimated laser beam emitted from an externally calibrated fiber source and coupled into a lens collimation unit to produce an approximately pencil-sized laser beam and measured on the other side with a Flow2 module detector as shown in Figure 2A. A power meter was placed between the collimated beam and the input face of the phantom to verify the correct power reading prior to each measurement and data collection session. We used an input power of 0.2 mW to yield target counts similar to those reported in the BIP protocol (we reached 0.68 x 10⁶ photons and 1.12 x 10⁶ photons at 690 and 905 nm, respectively, with our system’s default integration time of 3.5 ms).

The DNL measurement in the BIP protocol characterizes variability in the time bin width of the TDCs. Variability in the time bin width results in a non-uniform number of photons in each bin when the TDC is illuminated with a continuous light source. The DNL was measured by applying a uniform illumination source to the detector from a battery-powered CW light source for 100 s, and calculating the relative differences in the number of collected photons per bin. The deviation from the ideal of an equal number of photons in every bin is calculated as the peak-to-peak difference between the maximum and minimum bins, normalized by the mean photon counts over bins:

The IRF is carefully characterized in the BIP protocol to describe the overall time resolution of the system. It is typically measured by coupling the attenuated output of a source directly into a detector. The IRF was measured both at a typical detector and at the dedicated IRF detector. To measure the IRF at a typical detector and because our system does not have the same flexibility of fiber-based systems, a custom fixture was used to capture light in reflectance mode (Fig. 2B). The source beam was adjusted to avoid saturating the detectors with direct illumination, and the light is then reflected off a matte surface to redirect it into the collection optics. To measure the IRF at the dedicated IRF detector, the built-in waveguide that directly couples light from the source optical train to the IRF detection ASIC was used.

The IRF results from a convolution of the laser pulse shape and the temporal response of the detector and associated electronics. Per the protocol (Wabnitz, Taubert, et al., 2014), the IRF was measured by averaging 20 histograms of 1s acquisition time. Because our detector maximum integration time is only 800 ms, we collected these 20 histograms by summing individual 3.5 ms histograms to produce a single 1001 ms histogram. These 20 summed histograms were then averaged to calculate the IRF. The BIP protocol specifies a count rate of 1e6/s, that is, each of the 20 1-s histograms should contain 1e6 photons. As the count rate of our system is much higher than 1e6/s, we also present the results when we match the protocol in terms of photon counts (i.e., much shorter integration time).

From the IRF measurement, we also calculate the afterpulse ratio (RAP), a signal intensity noise source associated with the detector, as defined in the BIP protocol:

where N_mean,bkg and N_mean,dark are the average counts of the background measurement in the tail of the IRF and dark count measurements, respectively, T_laser is the full laser period (1 / repetition rate) and $Δt$ is the time bin width. The afterpulsing ratio was calculated for both 690 nm and 905 nm. The afterpulse ratio is a measure of the increase in the noise floor due to this intensity-dependent noise.

Finally, we measured the stability of the IRF over a 2-h period starting from a cold start. This measurement characterizes the time scale of thermal equilibrium for a single Kernel Flow2 module. We analyze the total intensity, the first moment, and the shape of the IRF as described in the BIP protocol continuously over the 2-h recording period.

The MEDPHOT protocol serves as an evaluative framework for various photon migration instruments, assessing their capacity to accurately retrieve known optical properties within a physiologically relevant range using homogeneous phantoms (Pifferi et al., 2005). Being able to retrieve accurate optical properties at each wavelength used by the system is a prerequisite to being able to retrieve accurate estimates of the concentrations of oxygenated and deoxygenated hemoglobin from the combination of the wavelengths (through the modified Beer–Lambert Law). Importantly, the range of properties represented by these phantoms encompasses properties of human head tissue (Doulgerakis et al., 2019) (skin, skull, cerebrospinal fluid, gray matter, and white matter) (Fig. 4A; Supplementary Table S1). In this study, we employed the Kernel Flow2 modules to conduct measurements on a specific subset of the MEDPHOT kit, comprising 12 solid phantoms (BioPixS, Ireland). These cylindrical phantoms, measuring 50 mm in height and 100 mm in diameter, consist of solid compositions containing titanium dioxide (TiO₂) and absorbing toner at varying concentrations. The phantoms are identified by letters (A, B, C, and D) and numbers (1, 3, and 5), where the letters denote nominal scatter values, and the numbers represent absorption values, as illustrated in Figure 2C.

Utilizing a single Flow2 module, we probed the phantoms from the top within a reflectance geometry setup. The absorption and scattering parameters were estimated by minimizing the disparity between the measured histogram (DTOF) and a predicted histogram, the latter being the outcome of convolving the measured impulse response function (IRF) with an analytical temporal point spread function (TPSF). The TPSF equation was derived from an analytical semi-infinite diffusion model (Patterson et al., 1989) employing an approximation to the Robin boundary conditions (Schweiger et al., 1995). The search for optical properties was carried out using the Levenberg–Marquardt algorithm, focusing on fitting within the range spanning from 80% of the peak on the rising edge to 0.1% of the peak on the falling edge, with a refractive index set to 1.41. We focused on the long within-module channels (26.5 mm source–detector distance).

To ensure the stability of optical property characterization, we subjected the B3 phantom to measurements for a duration of 2 h starting from a cold start, with an independent fit of optical properties from each 1 s of data.

The nEUROPt protocol (Wabnitz, Jelzow, et al., 2014) is employed in the evaluation of devices at the system level through the utilization of optical phantoms designed to replicate brain tissue. We used previously described methods to perform this evaluation (Ban et al., 2022). A liquid phantom was prepared consisting of a mixture of water, India ink, and intralipid emulsion (Intralipid^® 20%), titrated to have optical properties of μ_a = 0.01 and μ_s′ = 1.0 mm^-1 (titrated separately for each wavelength). The liquid phantom tank is made of black anodized aluminum with 0.1-mm thickness mylar windows for the source and detectors, as shown in Figure 2D. Experiments were conducted with a set of black polyvinyl chloride (PVC) cylinders with dimensions (diameter equal to height) of 3.2, 5, and 6.8 mm, corresponding to volumes V_incl/mm³ = 25, 100, and 250. The occlusions were suspended in the titrated solution by thin metal wires (0.4 mm, painted white) and were placed in the path of source–detector pairs formed within the Flow2 module. The target was moved incrementally from a depth of 8 to 36 mm away from the source–detector plane in 2 mm steps. We took measurements at each depth. The original protocol calls for 100 1-s accumulation histograms, at a count rate of 1e6 photons/s. As with the BIP protocol, we matched the original protocol in two ways: with the 1-s integration time (and a higher count rate), and with the target 1e6 photons per histogram (and a lower integration time).

As is customary in TD-fNIRS analysis, we present results after summarizing the full resolution DTOFs in 2 different ways (Lange & Tachtsidis, 2019; Ortega-Martinez et al., 2022; Wabnitz et al., 2020): (1) using the first 3 moments, as described in the previous section; (2) using coarse time gates: here, we first deconvolve the measured IRF from the measured DTOFs, using an FFT deconvolution approach with a regularization factor that dampens higher frequency content (see Wahab & O’Haver, 2023); then we coarsen the data to 500 ps gates (0–500 ps, 500–1000 ps, etc.). For each set of features (moments and gates) and each wavelength, we computed two different metrics: contrast and contrast-to-noise ratio (CNR), as a function of occlusion depth.

Contrast C for a given measurand M is defined as

where M₀ is the reference or baseline measurement and M_i is a measurement made after some i^th change in absorption (in our case an absorbing target) is introduced.

For photon counts, for which absolute contrast values are not as interpretable, we report a relative contrast measurement Cr defined as

Finally, for all measurands, the contrast-to-noise ratio (CNR) is defined as

where σ(M_o) is the standard deviation of the reference measurement across time samples.

Informed consent was obtained from all participants before beginning the study. The study was approved after ethical review by the Advarra IRB (#Pro00074416), and was conducted according to the Declaration of Helsinki. All experimental protocols were approved by the IRB. Participants received monetary compensation for their time, effort, and travel expenses.

Two participants (male, age = 28, Fitzpatrick skin type II, thin hair density, dark blonde, fine width/thickness, curly texture; woman, age = 35, Fitzpatrick skin type III, medium hair density, light brown, medium width/thickness, straight texture) contributed human data. One participant completed two tasks, namely, the finger-tapping and breath-hold tasks, and the other participant completed the auditory task.

This is a commonly used calibration task in the hemodynamic neuroimaging literature, as it produces a large and reproducible response in blood oxygenation and blood flow in the brain (Emir et al., 2008). This task was programmed and presented in Unity. The participant was asked to switch between interleaved periods of holding their breath and periods of paced breathing. During each paced breathing block (30 s), a bright green circle on a black background repeatedly expanded and contracted at a fixed pace (6 s per cycle), and the participant was instructed to use this animation to guide their inhalation and exhalation, respectively (5 breathing cycles were repeated in each block). At the end of paced breathing blocks, the circle changed to yellow, signaling to the participant that this would be their final exhalation and the breath-hold period would occur next. During each breath-hold block (20 s), the fully contracted yellow circle remained on the screen, above which the words “Hold your breath!” were displayed. Below the yellow circle a countdown to zero indicating the time left in the block was displayed. When the countdown timer hit zero, the circle turned back to bright green and paced breathing immediately commenced. This was repeated such that the participant completed a total of six breath-hold blocks and six paced breathing blocks.

This task was programmed and presented in Unity. The task had a block design with two block types: story blocks (n = 8) during which the participant listened to short clips from TED talks and noise blocks (n = 7) during which the participant listened to brown noise. After an initial 10-s rest period, the story and noise blocks (each lasting for 20 s) were presented (via earbuds) in a preset pseudo-randomized order. The participant was asked to keep their eyes open and look at a white fixation cross that was presented on a black background throughout the task.

This task was programmed and presented in Unity. In this task, the participant was asked to sit in a chair with their arms on the armrests such that their palms faced upward, while audio and visual stimuli guided them through randomized periods of left- and right-hand finger tapping (n = 10 blocks per side). Specifically, at the start of each block, the participant was cued in two ways: (1) audibly—brown noise was played through earbuds to either the right or left ear indicating the hand that should be used during the task and (2) visually—a white image of a hand was displayed on a black screen with either an “L” or an “R” inscribed on it, again indicating the left or right hand should be used. Both cues persisted throughout the block (17.3 s). Within each block, the participant was asked to repeatedly tap the thumb of the cued hand to a certain finger on the same hand. A red dot overlaid on a finger of the visual stimulus indicated which finger to tap. Throughout a block, the red dot moved sequentially through each of the four fingers, and each shift to a new finger indicated a new trial (n = 13 trials per block; trial duration = 0.75 s and inter-trial interval = 0.50 s). A brief resting period (20 s) with a white fixation cross on a black screen followed each block.

The data preprocessing procedures have been extensively detailed in our previous studies (Castillo et al., 2023). Initially, we applied a channel selection method based on histogram shape criteria (Ban et al., 2022). Subsequently, histograms derived from the chosen channels were utilized to calculate the moments of the DTOFs, specifically focusing on the sum, mean, and variance moments. The alterations in preprocessed DTOF moments were then translated into changes in absorption coefficients for each wavelength, employing the sensitivities of the various moments to absorption coefficient changes, as outlined in Ortega-Martinez et al. (2022). To determine these sensitivities, a 2-layer medium with a superficial layer of 12 mm thickness was employed. Utilizing a finite element modeling (FEM) forward model from NIRFAST (Dehghani et al., 2008; Doulgerakis-Kontoudis et al., 2017), the Jacobians (sensitivity maps) for each moment were integrated within each layer to assess sensitivities. The changes in absorption coefficients at each wavelength were further converted into alterations in oxygenated hemoglobin and deoxygenated hemoglobin concentrations (HbO and HbR, respectively), employing the extinction coefficients for the respective wavelengths and the modified Beer–Lambert law (mBLL (Huppert et al., 2009)). The HbO/HbR concentrations underwent additional preprocessing through a motion correction algorithm known as Temporal Derivative Distribution Repair (TDDR (Fishburn et al., 2019)). To address spiking artifacts arising from baseline shifts during TDDR, they were identified and rectified using cubic spline interpolation (Scholkmann et al., 2010). Lastly, data detrending was performed using a moving average with a 100-s kernel, and short channel regression was employed to eliminate superficial physiological signals from brain activity (Gagnon et al., 2011), utilizing short within-module channels with a source–detector separation (SDS) of 8.5 mm.

To retrieve absolute optical properties, we used the procedure described in the MEDPHOT methods section (Section 2.2.2). We used a refractive index of 1.4. The absorption coefficient estimates were then converted to HbO and HbR concentrations. A single value for HbO and HbR was obtained by computing the median value across well-coupled long, within-module channels (SDS = 26.5 mm) of two prefrontal modules. The use of a single-layer homogeneous model is of course a simplification, as head tissue is not homogeneous.

A GLM approach was employed in order to elucidate patterns of significant brain activity during the different block conditions of the auditory and finger-tapping tasks. For each task, the activity of each channel (the hemodynamic time course, $y)$ was fitted with a linear model:

where the matrix, $X$⁠, was composed of (i) task relevant regressors (the time course for each block condition represented as a square wave convolved with a canonical hemodynamic response function) and (ii) task irrelevant regressors (namely, drift and low-frequency cosine terms). A least-squares method was used to solve the multiple regression problem, producing the fitted model coefficients, β-weights. These weights quantify the effect that each regressor has on the hemodynamic signal. To evaluate whether the activity of a given channel is modulated across task conditions, β-weights associated with the block types of interest are subtracted (commonly termed as contrasts) and a t-test is employed to statistically compare this difference to zero. The contrasts of interest were story versus noise for the auditory task and left versus right for the finger-tapping task. The resulting test statistics (from each channel) were plotted as a heatmap over the head to visualize patterns of brain activity and regions of interest.

The associated p-values from the GLM analysis were used to identify representative channels that showed significant activation to a given task condition for time course visualization and epoched analyses. These channels were chosen from the regions of interest within each task (i.e., motor and auditory areas for the finger-tapping and passive auditory task, respectively), and to avoid superficial signals, only channels with SDS > 12 mm were considered. The time series for these channels underwent further processing; namely, detrending (with a 100-s kernel) and low-pass filtering (0.01–0.1 Hz finite impulse response—FIR filter). For each task, the time course of HbO and HbR for one representative channel was overlaid on the task design time course (Fig. 7B, E).

Moreover, we performed standard epoching analyses by windowing the channel time courses over blocks and grouping them by block type. For both tasks, the windows started 5 s before a block-start event and extended 30 s beyond (window = [-5 s, 30 s]; where block start = 0). A period preceding the start of each block (-5 to -1 s) was used to baseline the response of each epoch. The average time course (mean) and variability of response (standard error) for each block type was computed over the condition windows.

A finite element model (FEM) of the adult head was developed based on the unbiased non-linear averages of the MNI 152 database (Fonov et al., 2009). The atlas was segmented to five tissue types of skin, skull, CSF, gray and white matter, and discretized into linear tetrahedral elements using NIRFASTSlicer, giving rise to 413,403 nodes and 2,465,366 elements. Optical properties of each tissue layer at each wavelength (690 nm and 905 nm) were assigned based on published values of the adult head (Doulgerakis et al., 2019) (see Supplementary Table S1). The coordinates for each of the 40 modules containing the optical sources and detectors were determined and identified on the surface of the FEM and the time-resolved light propagation model was solved using the diffusion approximation to the light transport equation throughout the domain (Dehghani et al., 2008). The Jacobians (sensitivity functions that map a change in measured data due to a change in optical properties) for the time-resolved data (TPSF; we used a time resolution of 10 ps) for each optical parameter (μ_a and μ_s′) were calculated using the adjoint theorem (Arridge & Schweiger, 1995) at each wavelength and then interpolated to a uniform voxel grid (also known as reconstruction basis) spanning the entire model, with a resolution of 4 × 4 × 4 mm. The use of lower resolution reconstruction basis is crucial for DOT as the problem is highly under-determined: that is the number of measurements is much lower than the number of unknowns. While a high-resolution FEM mesh is needed for the calculation of the time-resolved light propagation to ensure numerical accuracy, a much lower voxel resolution is needed to improve the stability of the inverse problem.

The time-resolved Jacobian for each optical property was then mapped to each data type (intensity, mean time of flight, and variance; see Supplementary Material text for derivation). The Jacobian for each data type was then normalized to the [-1,1] range, and the same normalization was applied to the data, to give all data types an equal weight in the inversion process.

The following equation was solved:

where

-  $J$ is the Jacobian matrix of shape (channels x data types, voxels x 2 [absorption, scattering])

- $X$ is the optical properties vector (voxels x 2, 1)

- $Y$ is the data vector of shape (channels x data types, 1).

A Moore–Penrose pseudoinverse with Tikhonov regularization was used to calculate an approximation of the inverse of the Jacobian to perform a single step linear recovery of the optical properties (Doulgerakis et al., 2019) using the same functional data (DTOFs processed to moments) as outlined earlier. Briefly, the system becomes

Equivalently, and more practically when $J$ is wide (more voxels than measurements)

Here we did not optimize the Tikhonov regularization parameter $λ$⁠, setting it to $0.1||JJT||F,$ (where $||.||F$ designates the Frobenius norm) which in our experience works reasonably well.

Note that we downsampled the data to 1 Hz before performing reconstruction. For each time sample, the previous time sample serves as baseline, meaning that we recover the sample-to-sample changes in absorption assuming that changes are small to allow linearization of the problem. This is a typical approach as the true baseline optical properties are challenging to know experimentally and published literature values are often used (Eggebrecht et al., 2014). The recovered changes in the μ_a within each voxel were mapped to changes to oxygenated/deoxygenated hemoglobin for further processing.

It has been shown that different moments of the DTOFs exhibit different depth sensitivities (Fig. 5A, also see Wabnitz et al., 2020). To assess how the quality of reconstruction varies when using only the sum moment (i.e., intensity—similar to what CW systems measure) versus when additional moments (mean and variance) enabled by TD capabilities are included, we obtained reconstructed HbO/HbR cortical maps in these two scenarios. The time course of reconstructed activity was further analyzed using the same GLM model described above. Lastly, in addition to GLM analyses, we performed an epoched analysis. Here, we considered different ROIs given the task: voxels within 10 mm of the left motor area with the maximum GLM contrast for the finger-tapping task and voxels within 10 mm of the left auditory region with maximum GLM contrast for the passive auditory task. The time course of these ROIs was then epoched and aggregated within each block type for further visualization (Fig. 8).

The responsivity of the Flow2 detectors was calculated to be 9.97 x 10^-9 and 2.28 x 10^-9 m²sr, respectively, at the 2 wavelengths. These are within the range of reported values in research-grade systems (Lanka et al., 2022). Compared with Flow1 at 690 nm, sensitivity was improved by 38% (owing to a combination of improved detector design and optics). The DNL was found to be 0.148. This DNL is relatively high compared with other TD devices (Lanka et al., 2022). However, there is no established threshold of acceptable DNL for TD systems; in fact, DNL can be measured with CW light and compensated for in signal processing if it affects the shape of the time-of-flight histograms. We have not found the DNL of Flow2 to be detrimental to any of the functional characterization results (see next sections, MEDPHOT and nEUROPt), and there is thus no need to perform any signal processing correction.

The afterpulsing ratio (Methods; Eq. 2), a known source of signal-dependent noise, was 0.0024 at 690 nm and 0.0019 at 905 nm. While this metric is not as widely reported in the literature despite being part of the benchmarking protocol, these values are similar to those reported for another system (Re et al., 2013).

The IRF of a TD system is critical to its ability to accurately characterize the media through which light travels. The data recorded by the system for a given source–detector pair (which we refer to as the DTOF) consists of the convolution of the DTOF one expects from the radiative transfer equation (which we refer to as the temporal point spread function, or TPSF) with the IRF of the system (Diop & St. Lawrence, 2013). Therefore, a wide IRF tends to blur information. For our system, the full width at half maximum (FWHM) of the IRF (as measured on the IRF detector, for a representative module and source) was 240 ps at 690 nm and 270 ps at 905 nm, at the beginning of the session (Fig. 3A, B). At the end of the 2-h session, it was 285 ps at 690 nm and 305 ps at 905 nm. This drift with recording time is due to an increase in temperature over the course of the session, which affects several aspects of the electronics. The IRF FWHM is below or similar to the mean as compared with other systems (Lanka et al., 2022). While the IRF is typically characterized by its FWHM, we note that the ability of a TD system to accurately characterize media with higher absorption is determined by the slope of the tail of the IRF. We, therefore, also report the width at other fractions of the maximum, as well as the time constant of the exponential decay (computed over three orders of magnitude). At the beginning of the session: FW10%max = 710 ps (690 nm), 710 ps (905 nm); FW1%max = 1630 ps (690 nm), 1655 ps (905 nm); time decay constant = 292 ps (690 nm), 265 ps (905 nm). After 2 h, FW10%max = 710 ps (690 nm), 725 ps (905 nm); FW1%max = 1740 ps (690 nm), 1715 ps (905 nm); time decay constant = 297 ps (690 nm), 267 ps (905 nm).

A common way to summarize information from time-of-flight histograms is to compute the first three moments of the histogram corresponding to the total counts (sum), mean time-of-flight (first moment), and variance of the times of flight (second central moment) (Liebert et al., 2003; Ortega-Martinez et al., 2022; Wabnitz et al., 2020). Moments have a convenient property: the moments of the DTOF can be obtained from calculating the moments of the TPSF and of the IRF straightforwardly (Liebert et al., 2003). Accordingly, with Flow2, system drift in the DTOF moments can be corrected for, using the internal IRF detector. We demonstrate here that the first three moments of DTOFs (total counts, mean, variance) were consistently stable over the course of the 2-h session after such correction, as shown for a representative channel in Figure 3C.

We found the retrieved optical properties to be very close to the calibrated values for both wavelengths and all 12 phantoms (Fig. 4A, B), with the following average relative errors: 8.4% and 4.3% for μ_a, and 4.2% and 3.8% for μ_s′ at 690 and 905 nm, respectively. These relative errors are similar to what has been reported for other instruments (Lanka et al., 2022). We note that accuracy is poorer at higher absorption levels, which we primarily attribute to a limitation imposed by the slope of the IRF (hence the importance of characterizing it fully). The retrieved and calibrated optical properties showed excellent linearity, with a correlation coefficient over 0.98 for both μ_a and μ_s′ (Fig. 4C, D), and a median relative error with respect to the line of best fit of 6.3% (690 nm) and 3.0% (905 nm) for μ_a, and 2.3% (690 nm) and 3.1% (905 nm) for μ_s′. This performance is in the range of those reported by similar instruments (Lanka et al., 2022).

Despite a noticeable increase in the system temperature (+25 deg C) over a 2-h recording period from a cold start for phantom B3 at 905 nm (Fig. 4E), the retrieved μ_a and μ_s′ were stable (slopes computed over full time range are less than 0.01%/min for μ_a and μ_s′; maximum deviation is 3% in the first 20 min), further confirming that the continuous measurement at the dedicated IRF detector is an effective way to fully correct for system drift. The initial transients are likely attributable to temperature gradients between components as the system heats up from a cold state.

For the moments, we found that both contrast and CNR metrics exhibited a depth-selective profile with higher moments being more sensitive to deeper occlusions (Fig. 5A), as expected from the literature (Wabnitz et al., 2020). Similarly, when using the gating approach (Fig. 5B), later time gates show better contrast and CNR with depth (Selb et al., 2005). It is important to note that with both methods, there is high contrast and CNR at depths >20 mm (which is well into the brain) (M. M. Wu et al., 2022), especially for higher moments and later time gates—demonstrating the advantage of TD-fNIRS systems over CW-fNIRS systems with respect to depth sensitivity. These results are among the best performance of TD systems according to a recent benchmarking study (Lanka et al., 2022).

We observed systematic changes in the absolute HbO and HbR that were locked to the task events (Fig. 6A), and these increases and decreases were consistent with the literature (Emir et al., 2008; MacIntosh et al., 2003). Specifically, HbO revealed a significant increase after the start of the breath-hold period and underwent a decrease when the next breathing period began. Furthermore, both HbO and HbR were modulated by the breathing pattern such that inhaling and exhaling resulted in changes in the signals (inhale: HbO increases, HbR decreases; exhale: opposite). Lastly, we extracted the time-varying heart rate (using the dedicated prefrontal module; Methods), which revealed the presence of strong respiratory sinus arrhythmia (RSA; Fig. 6B), that is, modulation of the heart rate as a function of the phase of the breathing cycle (Tonhajzerova et al., 2016).

A GLM analysis revealed significantly more activity in bilateral auditory cortex during the story blocks than during the noise blocks (Fig. 7A)—a pattern of brain activation consistent with the literature (Belin et al., 2000; Luke et al., 2021).

To confirm the statistical output of the GLM visually, we examined the channel-level oscillations of HbO/HbR concentration over the full time course of the task. When overlaid with the block design of the task, channels that were within the regions of interest (and identified as significant by the GLM analysis) demonstrated the expected hemodynamic activity, effectively tracking the experimental blocks throughout time (representative channel from left auditory cortex shown in Fig. 7B). As expected, the two chromophores fluctuated in opposing directions.

Second, for each channel, we averaged the hemodynamic signal over epochs corresponding to the same block type. Representative epoched responses for single channels located in the right/left auditory cortex are displayed in Figure 7C. These channels exhibited increases in HbO, and corresponding decreases in HbR during the story condition of the task. Notably, they also exhibited negligible or relatively smaller responses to the noise condition and complementary HbO/HbR activity.

During the finger-tapping task, a GLM analysis revealed regions of significant activation in the right motor cortex and regions of significant deactivation in the left motor cortex when comparing left-hand tapping with right-hand tapping (Fig. 7D), as expected from the known organization of motor cortex (Brinkman & Kuypers, 1973).

The time course for a representative channel from the right motor cortex exhibited increases in HbO throughout time that were matched to left blocks and accompanied by decreases in HbR (Fig. 7E). These characteristic hemodynamic responses were even more evident after epoching. The epoch-averaged time courses for representative channels exhibited condition-specific hemodynamic responses. Specifically, a channel in the left motor cortex showed a rise in HbO and a corresponding decrease in HbR during right-tapping blocks (dashed line) (Fig. 7F, left), while a channel in the right motor cortex showed increases and decreases in HbO and HbR, respectively, during left-tapping blocks (Fig. 7F, right). Taken together, these analyses underscore the ability to use Flow2 to record meaningful brain activity during a passive sensory task as well as an active motor task.

GLM activation maps using the sum moment for reconstructing activity in image space exhibited expected patterns for the auditory task: increased activity during story blocks as compared with noise blocks in the auditory region (Fig. 8A). In the finger-tapping task, too, the activation maps revealed elevated contralateral motor cortex activity (Fig. 8B).

When performing reconstruction using all three TD moments, we indeed found deeper sensitivities with higher moments (Supplementary Fig. S2). Consequently, both tasks demonstrated even larger, more significant, statistical contrasts between activations during the conditions in expected regions (Fig. 8C, D). This difference was even more stark when considering the time course of activation in spherical ROIs around the peak of activation (Fig. 8E–H). For example, while the sum moment-based HbO time course in the left motor region during right-tapping blocks exhibited small elevated activity compared with that during left-tapping blocks (Fig. 8F), this contrast in activity was notably larger when using the reconstructed HbO from all moments (Fig. 8H).

In the current study, we introduce and validate the Kernel Flow2 as a wearable TD-DOT system capable of measuring and reconstructing brain activity over the whole head. We have shown the system has performance comparable with limited channel count research-grade devices according to standardized benchmarking protocols, while also extending the field-of-view over the whole head. The system hardware was developed using scalable methods used for consumer electronics such as custom ASICs, plastic molded optics, low-power application processors, and standard communication protocols. This commercially available system improves on the Flow1 prototype (Ban et al., 2022) on many key performance metrics (summarized in Supplementary Table S2). This combination of system performance, cost, and form factor can enable many applications for both neuroscience research and precision neuromedicine.

It is critical to validate the performance of newly developed measurement systems using standardized protocols, before they can be applied for research or clinical purposes. For TD-fNIRS systems specifically, recent multi-laboratory efforts have resulted in comprehensive tests to characterize performance: the BIP, MEDPHOT, and nEUROPt protocols. Each of these focuses on performance from a different angle (Lanka et al., 2022; Pifferi et al., 2005; Wabnitz, Jelzow, et al., 2014; Wabnitz, Taubert, et al., 2014). Briefly, these tests revealed the following for Flow2: (1) the system’s IRF is adequately narrow and stable, which is required for various computations (e.g., computing absolute HbO/HbR concentrations); (2) the system can indeed be used to accurately retrieve optical properties from calibrated optical phantoms with optical properties in the physiological range at both wavelengths, which is required to retrieve accurate concentrations of HbO and HbR; and (3) the measurements from the system have better depth sensitivity than the measurements from CW systems (intensity only). These metrics were similar or better than those reported for limited channel count systems (Lanka et al., 2022). Together, these validations against known ground truth suggest that Flow2 recordings should be of high enough quality to study hemodynamic activity on a human head. We thus proceeded to investigate spatial and temporal patterns of brain activity extracted from neural recordings in vivo on a human subject.

First, we used an easily implemented mild hypercapnic challenge, which consists in a short breath hold (at the end of an exhale). This challenge is known to trigger an increase in arterial PCO₂, resulting in a homeostatic response including cerebral vasodilation and an increase in cerebral blood flow; importantly, this mild challenge does not appear to result in extracranial perfusion effects (MacIntosh et al., 2003). We were able to use curve fitting with the TD data to track the absolute concentrations of oxygenated and deoxygenated hemoglobin, something that traditional CW systems cannot do without making additional assumptions. Because the homeostatic response is limited to cerebral tissue in this hypercapnic challenge, we can demonstrate that our approach of using a homogeneous model is sensitive to absolute oxygenation changes in the brain. Of course, two-layer and even three-layer analytical models could be used instead and may improve CNR. Retrieving absolute properties is of particular interest in many clinical applications of TD-fNIRS such as in stroke assessment (Obrig & Steinbrink, 2011), neurocritical care (Abdalmalak et al., 2017; Thomas et al., 2023), multiple sclerosis (Yang & Dunn, 2015), schizophrenia (Hoshi et al., 2006), and dementia (Viola et al., 2013). Thus far, because of the high cost and bulkiness of these systems, their clinical use cases have been limited to research settings (Lange & Tachtsidis, 2019). Therefore, using a portable system such as Flow2 may present an opportunity to further translate these lines of clinical research into practice.

We then had participants wearing the Flow2 headset engage in two simple tasks: a passive auditory (sensory) task and a finger-tapping (motor) task. Note that both of these tasks are classic neuroimaging paradigms with consensus from the scientific community regarding the foci of functional brain activity (Belin et al., 2000; Luke et al., 2021; Witt et al., 2008). Our analyses showcased localized regions of activation in agreement with the literature: bilateral auditory cortex activity for the passive auditory task and contralateral motor cortex activity for the finger-tapping task. Further, those channels that exhibited the most statistically significant responses to task conditions revealed strong task-locked modulations. It is worth noting that we previously demonstrated the stability of such task-evoked metrics in a larger-scale study with the Flow2 device (Dubois et al., 2024). Taken together, these findings uphold the two pillars of device performance for research and clinical use: validity and reliability.

Individual channels of the Flow2 system yield high-quality hemodynamic signals from the human head. As the Flow2 system is composed of a very large number of channels (over 3,500 with source–detector distance <60 mm) in a wearable helmet-like form factor, it further enables reconstruction of brain activity (DOT) over the whole head. The density of channels is on par with state-of-the-art CW HD-DOT systems (Vidal-Rosas et al., 2023). Furthermore, TD-DOT improved the statistical significance of cortical activations as compared with CW-DOT; we can speculate that this is because of its better depth resolution as demonstrated in our characterization results (e.g., nEUROPt protocol results). We note that the comparison performed here is as targeted and unconfounded as can be, as it was performed using the same input data (TD DTOFs) summarized in different ways. Establishing that in vivo reconstructed activations have better anatomical specificity in the TD-DOT case than in the CW-DOT case would require an independent measurement, for example, using fMRI, to establish a ground truth to which reconstructions could be compared. We focused on using moments of the DTOFs for all analyses, however, there may be analyses based on the arrival times of photons (“gates”) which would further increase sensitivity to deeper brain activity. In addition, to fully establish the increased depth sensitivity of the TD measurements, two-layer phantom measurements could be used to validate the improved depth resolution using a controlled experimental setup.

While the Flow2 system has overcome many challenges of the previous prototype Flow1 system and low-channel-count TD systems, there are still remaining limitations. The full headset remains relatively heavy, as all electronics are built-in. The weight may limit system wear time, thus the headset may not accommodate applications requiring long recordings (although in some cases this can be mitigated with thoughtful study design). Additionally, as with other optical neuroimaging devices, hair can limit signal quality, as it absorbs light which reduces signal at the detectors. We have found that at least 1,000 usable channels are available over the head for most participants (Dubois et al., 2024), there is still room for improvement in mechanical design to ensure good optical coupling by working through all hair types and making contact with the scalp at all locations. Moreover, as with all diffuse optical systems, deep brain structures remain inaccessible. The TD nature of our device does, however, allow for deeper sensitivity in cortical regions as compared with CW systems. Finally, to perform a reconstruction of brain activity, a structural head model must be used which describes the tissue composition (scalp, skull, cerebrospinal fluid, gray matter, white matter), and the precise positioning of all optodes (sources and detectors) must be known. For the reconstructions performed in this work, we used an atlas head model and average optode positions (see Methods). This is the most scalable approach, which does not require individual head measurements or individual structural MRIs, and solely relies on good practices for headset positioning. It is also an approach that makes several assumptions and, therefore, has room for error. There are intermediate ways to improve atlas-based tomography by using individual measurements, such as individually registering optode locations (Mazzonetto et al., 2022) or picking the best atlas from a library (Fonov et al., 2009). There are ongoing efforts to make such methods readily available for Flow2. Additionally, more advanced tomographic reconstruction approaches that fully utilize the DTOFs will be explored and implemented in future iterations (Gao et al., 2002).

Although the application of optical neuroimaging systems in clinical settings has been discussed repeatedly (Bonilauri et al., 2020; Lange & Tachtsidis, 2019; Li et al., 2023; Wheelock et al., 2019), we argue that no system to date has successfully integrated the combination of (1) an easy-to-use form factor, (2) the most advanced fNIRS modality (TD), and (3) high enough density to enable the reconstruction of whole-brain cortical hemodynamic signals. The confluence of these factors in Flow2 facilitates better, larger-scale, and more timely data collection efforts. The combination distinctly positions the device as a premier recording option for both the scientific and clinical communities. Future studies with TD-DOT can translate insights gleaned from neuroimaging research and bring them into practice, therefore, enabling novel advances in precision medicine through scalable neuroimaging-based biomarker identification (Hampel et al., 2023; Nebe et al., 2023; Rossi et al., 2023).

The human data in this study is available at OpenNeuro (https://openneuro.org/datasets/ds005927).

Conceptualization: R.M.F., K.L.P. System design and engineering: Y.C., D.D., R.M.F., V.G., I.O., J.S. Software: G.L., M.J.P., V.S. Formal analysis: J.D., H.D., E.M.K., Z.M.A., I.O. Investigation: V.G., I.O., N.M. Visualization: J.D., E.M.K., Z.M.A., I.O. Writing—original draft: J.D., H.D., E.M.K., Z.M.A., I.O., K.L.P. Writing—review and editing: all authors. Author order was determined alphabetically. Dr. Katherine Perdue agrees to be accountable for all aspects of the work, ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved.

This study was sponsored by Kernel.

The authors of the study, except for H.D., are all employees of Kernel, which sponsored this work. H.D. provided scientific consulting for Kernel. Authors do not receive financial incentives for publications. Some authors are named on granted and pending patents related to this work; however, the patents are held by Kernel and inventors do not receive royalties.

The authors gratefully acknowledge Moriah Taylor for her administrative support for human subjects data collection, and Eli Palmer for software support with task programming.

Supplementary material for this article is available with the online version here: https://doi.org/10.1162/imag_a_00475.