
==== Front
bioRxiv
BIORXIV
bioRxiv
2692-8205
Cold Spring Harbor Laboratory

39229084
10.1101/2024.08.20.608859
preprint
1
Article
Revisiting equivalent optical properties for cerebrospinal fluid to improve diffusion-based modeling accuracy in the brain
http://orcid.org/0000-0001-6826-8925
Lewis Aiden Vincent a
http://orcid.org/0000-0003-0805-935X
Fang Qianqian ab*
a Northeastern University, Department of Bioengineering, 360 Huntington Avenue, Boston, USA, 02115
b Northeastern University, Department of EECS, 360 Huntington Avenue, Boston, USA, 02115
* Qianqian Fang, q.fang@neu.edu
21 8 2024
2024.08.20.608859https://creativecommons.org/licenses/by-nd/4.0/ This work is licensed under a Creative Commons Attribution-NoDerivatives 4.0 International License, which allows reusers to copy and distribute the material in any medium or format in unadapted form only, and only so long as attribution is given to the creator. The license allows for commercial use.
nihpp-2024.08.20.608859.pdf
Significance:

The diffusion approximation (DA) is used in functional near-infrared spectroscopy (fNIRS) studies despite its known limitations due to the presence of cerebrospinal fluid (CSF). Nearly all of these studies rely on a set of empirical CSF optical properties, recommended by a previous simulation study, that were not selected for the purpose of minimizing DA modeling errors.

Aim:

We aim to directly quantify the accuracy of DA solutions in brain models by comparing those with the gold-standard solutions produced by the mesh-based Monte Carlo (MMC), based on which we derive updated recommendations.

Approach:

For both a 5-layer head and Colin27 atlas models, we obtain DA solutions by independently sweeping the CSF absorption (μa) and reduced scattering (μs′) coefficients. Using an MMC solution with literature CSF optical properties as reference, we compute the errors for surface fluence, total brain sensitivity and brain energy-deposition, and identify the optimized settings where the such error is minimized.

Results:

Our results suggest that previously recommended CSF properties can cause significant errors (8.7% to 52%) in multiple tested metrics. By simultaneously sweeping μa and μs′, we can identify infinite numbers of solutions that can exactly match DA with MMC solutions for any single tested metric. Furthermore, it is also possible to simultaneously minimize multiple metrics at multiple source/detector separations, leading to our new recommendation of setting μs′ = 0.15 mm−1 while maintaining physiological μa for CSF in DA simulations.

Conclusion:

Our new recommendation of CSF equivalent optical properties can greatly reduce the model mismatches between DA and MMC solutions at multiple metrics without sacrificing computational speed. We also show that it is possible to eliminate such a mismatch for a single or a pair of metrics of interest.

fNIRS
Cerebrospinal fluid
Monte Carlo simulations
Diffusion approximation
Brain imaging
Photo-biomodulation
==== Body
pmc1 Introduction

In the wavelength range between 600 nm and 1,100 nm, near-infrared (NIR) light can penetrate several centimeters of biological tissues with a highly scattering trajectory, as a result of relatively low tissue absorption. Multiple imaging techniques, such as functional near infrared spectroscopy (fNIRS) and diffuse optical tomography (DOT), capitalize upon this behavior to non-invasively monitor hemodynamics in cortical tissue resulting from brain activities. Similarly, in photobiomodulation (PBM) applications, clinicians use this phenomenon to deliver light to deep tissues for therapeutic purposes. Due to the highly scattering and stochastic nature of light-tissue interactions, researchers rely on quantitative modeling techniques to predict light dosages and analyze their measurements. Two widely used numerical techniques for quantitatively modeling light-tissue-interactions are the diffusion approximation (DA)1 and the Monte Carlo (MC) method.2

MC is a stochastic solver to the radiative transfer equation (RTE), a differential-integral equation known to be accurate for modeling light transport in general random media including biological tissues. MC is also relatively easy to implement and can be easily parallelized. However, the primary challenge MC faces is its high computational cost, as it requires to launch large numbers of photon packets to achieve a stable solution with acceptable stochastic noise. Over the past two decades, the widespread use of graphics processing units (GPU) has drastically reduced MC modeling time from several hours3 to tens of seconds.4 MC-based models have also been extended to accommodate increasingly complex tissue shapes, growing from infinite layered media5 to spatially heterogeneous voxel-based,3, 6 mesh-based7 or hybrid shape representations.8, 9 As a result, MC solutions have been increasingly seen in routine data analysis aside from serving its transitional role of providing gold-standard solutions.10

In comparison, DA solves a simplified version of the RTE by ignoring the ballistic behaviors of photons near a collimated source or in void/low-scattering regions. This results in a simpler elliptic partial different equation (PDE) that can be conveniently solved using numerical techniques such as finite-element (FE) or finite-difference (FD) methods. Typical solution time for a DA forward solution using an FE solver, such as NIRFAST11 or Redbird,12 is on the scale of a fraction of a second. This is significantly faster than MC solutions, even with GPU accelerations, and the output is deterministic (i.e. free of stochastic noise). Because of the high computational efficiency, the DA has been actively used in DOT image reconstructions especially when arrays of sources and detectors are used. However, when modeling light transport in the brain, a number of previous studies had demonstrated that the presence of low-scattering tissues such as cerebrospinal fluid (CSF) can produce erroneous solutions.13–18 A number of hybrid methods have been proposed to properly model voids and low-scattering tissues, such as CSF, lung, synovial fluid, cysts etc, however, these methods have received only limited adoption due to increased complexity.

The CSF layer is generally known to have low scattering, however, its literature absorption (μa) and reduced scattering coefficients (μs′) show a wide range of values, ranging between 0.0004 mm−1 and 0.004 mm−1 for μa,19, 20 and between 0.001 mm−1 and 3.0 mm−1 for μs′ 10, 21 due to diverse measurement methods and modeling assumptions. The CSF layer also occupies a complex anatomical space, filling primarily the subarachnoid space bounded by the arachnoid mater at the outer surface and the pia mater at the inner surface,22 as well as the folding space on the cerebrum surface, known as sulci. The CSF in the sulci are predominantly transparent, filling a complex folding geometry. Sulci’s widths and depths are highly dependent on locations, ranging from 0.2 to 2.3 mm for width and 5.7 to 13 mm for depth in young adults, varying further between age groups and genders.22 Multiple approaches for modeling CSF in the brain exist, including treating it entirely as a translucent fluid,20, 21 assigning separate optical properties to the subarachnoid space and sulci,10 and combining it with the cortical tissue using empirically derived bulk properties.23

A widely cited approach for extending DA in modeling light transport in the brain is described by Custo et al.20 In this work, the CSF layer is treated entirely as a diffusive medium with a recommended empirical reduced scattering coefficient of 0.3 mm−1 – determined by the typical inverse line-of-sight distance of the CSF layer.20 Because of the simplicity, this recommmendation has been widely adopted to justify the use of DA in modeling the CSF in the brain tissues.24, 25 Diffusion solvers such as NIRFAST11 and NeuroDOT26 also include this recommended CSF scattering property in many built-in examples for brain related data analyses, and received wide adoption among fNIRS research.27–29 However, most of the works utilizing this approach took the recommended values as the optimal solution without further scrutinizing the limitations on how such a recommendation was derived.

We want to highlight that it is particularly important to understand the conditions upon which the recommended CSF optical properties in Custo et al. were derived. First of all, this recommendation was drawn entirely based on comparing between MC solutions at varying CSF reduced scattering coefficients, instead of directly comparing between MC and DA solutions. Secondly, the chosen value μs′ = 0.3 mm−1 was determined as the upper-bound beyond which the MC solutions start to show large deviations from the respective ground-truth simulations; this “upper-bound” criterion was also quite different from the “optimal value” that best approximates DA with MC in the CSF that most use-cases of this recommendation commonly assumed. Thirdly, the physical quantity studied in the previous work is specifically limited to fluence and partial-path-lengths; the impact to other types of optical measurements, such as brain sensitivity – desired when solving fNIRS brain activity recovery and image reconstructions – and energy deposition – desired in many PBM related analyses – were not discussed. Lastly, the previous work only examined the impact of varying μs′; the impact of simultaneously altering CSF absorption coefficient μa was not considered.

Here, we would like to revisit this widely adopted recommendation by addressing the afore-mentioned limitations. Specifically, we aim to directly compare DA with MC solutions and seek to derive more appropriate CSF equivalent optical properties based on minimizing their differences. In addition to sweeping its reduced scattering coefficient, we also allow the CSF’s absorption coefficient to change, adding a new degree-of-freedom to help reduce the model mismatch. Moreover, we expand the comparison between DA and MC to include total brain sensitivity and energy deposition, extending the new recommendation towards broader optical brain imaging/therapy techniques. Furthermore, we compare our DA and MC solutions on both a simplified 5-layered head model as well as a more complex adult brain atlas – Colin27 30 – at two common wavelengths, seeking further generalization of our findings. Here we use our extensively developed mesh-based MC (MMC), known for its high accuracy among various MC solvers,7 with graphics processing unit (GPU) acceleration4 to provide the reference solution. To further remove the confounding systematic differences due to varying spatial discretization, we apply an identical set of tetrahedral meshes for use in both the FE DA solver and MMC.

In the remainder of this paper, we present our methodologies and results for minimizing DA modeling errors comparing to results obtained from MC. First, we describe the layered head and atlas anatomical models used, as well as the numerical solvers used for DA and MC respectively. Then, we describe the metrics we derive from the DA and MC solutions, such as brain sensitivity, surface fluence, and GM energy deposition; these metrics are used to quantify the errors caused by using DA in fNIRS and PBM applications. Finally, we discuss the results from both the layered and brain-atlas models, demonstrating that use of updated optical properties for CSF in DA models can significantly reduce, or even eliminate, mismatch against the “gold standard” MC models.

2. Methods

2.1 Anatomical models and tetrahedral mesh generation

We perform MC and DA simulations of light transport in two brain anatomical models frequently seen in literature: 1) a 5-layer head model and 2) Colin27 brain atlas.31 The 5-layered head model is created with layer thicknesses based on the average thickness of the atlas layers.31 Note that the scalp and skull are treated with one set of optical properties as in literature.31 A tetrahedral mesh of the Colin27 atlas is generated using the Brain2Mesh toolbox31 and is derived from the Colin27 magnetic resonance imaging (MRI) atlas, with four layers: combined scalp and skull, CSF, gray matter (GM), and white matter (WM). The physiological values for each layer, used in MC simulations serving as the ground-truth, are described in Table 1. We want to note here that all simulations in Custo et al.20 used an assumed CSF μa value of 0.004 mm−1. However, this value is 10× larger than the physiological μa value of 0.0004 mm−1 reported in Strangman et al.,19 which was also cited as the source of the optical properties. For consistency in the rest of our analysis, we use the μa values from the upstream source of Strangman et al.19 in all our reference simulations.

2.2 Forward models

We use our GPU-accelerated, mesh-based Monte Carlo (MMC) photon transport simulator4 – an open-source software that has been widely validated and disseminated among the biophotonics community33–38 – to create the “ground-truth” solutions. Briefly, MMC uses tetrahedral meshes to produce MC simulations calculating fluence at every node7 or element. The use of tetrahedral meshes allows simulations to consider more realistic biological tissues with curved and complex boundaries. All MC solutions are produced with a relatively large number (109) of photon packets to ensure stable results. For DA, we apply another in-house MATLAB toolbox named “Redbird-m” to provide diffusion solutions using a finite-element method (FEM).12 Redbird-m was developed from our various previous works, extending from optical breast imaging39 and structural-prior guided reconstructions.40, 41 To properly approximate a collimated source, such as a coupling fiber used in fNIRS probes, in Redbird-m, we sink both sources and detectors by 1/μs′ from the tissue-air boundary along the light incident direction.42

In all simulations reported below, both MMC and Redbird-m produce solutions over the same tetrahedral mesh generated by our Iso2Mesh mesh generator.31 This allows us to minimize discrepancies due to different discretization strategies. In addition, both solvers have implemented normalization methods to produce solutions that correspond to the Green’s function of the respective mathematical models, therefore, they can be directly compared at all nodal positions.

To simulate a typical fNIRS probe configuration, a pencil beam source is placed on the top surface of the chosen head model; a linear array of disk-shaped detectors with a radius of 1.5 mm, with a geodesic distance to the source ranging between 2.0 cm to 3.5 cm with an increment of 0.5 cm, are placed on one or both sides of the source. An additional near-separation detector was placed at a 0.84 cm geodesic distance from the source.

2.3 Metrics for accuracy assessments

To quantify the mismatch between MC and DA, we compute common metrics relevant to evaluating performance in PBM and fNIRS applications. For PBM, the total energy deposition within the GM layer (Egm) is computed to characterize the dosage of light energy that reaches the brain for therapeutic usage. For fNIRS, multiple metrics are used: 1) spatially resolved fluence,  Φ(r→), and detected fluence values,  Φ(r→s, r→d), sampled at detector locations (r→d) for any given source at rs are extracted to assess the similarity of optical measurements between MC and DA,13 2) the spatially resolved μa Jacobian, J(r→), i.e. sensitivity of μa at each location (r→), is computed to assess the loss of sensitivity caused by using DA,43 3) the total brain sensitivity (Sgm), computed by summing the Jacobian within the GM region, to estimate the overall impact of forward model accuracy to the recovery of brain hemodynamics, and finally, 4) fraction of GM sensitivity over total sensitivity (Egm) is used to quantify the relative impact of the forward models to fNIRS signal recovery. The definitions of each of the above metrics are detailed below.

In each of the used head models, a forward solution is produced using MMC and Redbird-m, respectively, by placing a pencil beam over the source position, with an incident direction along the normal direction of the surface. Both MC and DA simulations produce forward solutions of normalized fluence, Φ, as a Green’s function defined at each spatial location  r→∈Ω where Ω denotes the simulated domain. To calculate the total energy deposition to the brain, Egm, an integral of the forward fluence solution,  Φ(r→, r→s), multiplied by μar→ is performed within the GM region, Ωgm.

(1) Egm=∫r→∈ΩgmΦr→,r→sμar→dr→

To calculate the fluence at various detectors,  Φ(r→d, r→s), a simple linear interpolation is applied to obtain the normalized fluence at the exact coordinate r→d of the detector based on the forward fluence solution  Φ(r→, r→s) within the enclosing tetrahedron.

We apply the adjoint method1, 44 to compute the Jacobian matrix in both DA and MC. This is done by multiplying the forward solution simulated from the source,  Φ(r→, r→s), by an adjoint solution,  Φ(r→, r→d), obtained by simulating a source at the location of the detector r→d. To compute the total brain sensitivity, Sgm, we integrate the spatially resolved Jacobian over the entire GM region as (2) Sgm=−∫r→∈ΩgmΦr→,r→sΦr→,r→dΦr→d,r→sdr→.

Finally, we derive the fraction of Jacobian sensitivity over the total sensitivity, Fgm, by taking a ratio of the total brain sensitivity over the Jacobian integrated through the entire head model Ω as (3) Fgm=Sgm−∫r→∈ΩΦr→,r→sΦr→,r→dΦr→d,r→sdr→.

For each of the above metrics (M), a signed percentage error of DA relative to MC is computed as (4) ϵDA=MDA−MMCMMC×100%,

where ϵDA<0 represents a case where DA underestimates the ground truth, and ϵDA>0 represents overestimation.

3 Results

A tetrahedral mesh of 332,430 nodes and 2,015,332 elements is created using Iso2Mesh for the 5-layered head model; another tetrahedral mesh containing 151,097 nodes and 930,046 elements is produced for the Colin27 atlas. Both mesh models are shared between the MMC and Redbirdm solvers. All MMC simulations are launched using an NVIDIA RTX 2080 SUPER GPU on a Ubuntu 20.04 Linux server. DA solutions computed using Redbird-m are obtained using an AMD Threadripper 3990X processor on the same server.

3.1 Assessing accuracy of diffusion approximations using literature recommendations

Cross-sectional contour plots across all source and detector positions are shown in Fig. 1 comparing DA with the respective MC reference solutions at various conditions. In Figs. 1(a) and 1(d), we artificially set the CSF to a diffusive medium of μa = 0.0026 mm−1 and μs′ = 1 mm−1 with a goal of validating Redbird-m DA solver against the reference MMC solutions. With no surprise, the Redbird-m and MMC solutions are closely aligned across the entire simulation domain in both tested head models. The excellent agreement is also indicated by the relatively uniform and low relative errors (color shade) in most of the brain regions. In Figs. 1(b) and 1(e), we set the CSF to the anticipated physiological values (μa = 0.0026 mm−1, μs′ = 0.001 mm−1).19, 20 The significantly elevated errors in CSF, GM and WM layers as well as at larger source-detector separations once again verify the inability of DA to model low-scattering tissues directly.

In Figs. 1(c) and 1(f), we set the CSF’s optical properties to μa = 0.0026 mm−1 and μs′ = 0.3 mm−1) as recommended by Custo et al.20 It is clear from these results that this approach significantly reduces the overall mismatch between DA and MC. However, notable spatially-resolved errors ranging between 10% to 50% can still be observed within the CSF, GM and WM regions, as well as in large source-detector separations.

3.2 Optimization of optical properties in a 5-layer head model

To identify the optimal CSF equivalent optical properties for approximating MC simulations, in this section, we compute the DA forward solutions by sweeping CSF absorption (μa) and reduced scattering coefficient (μs′) in a large search space that encompasses the typical values seen in biological tissues. The search space for μa ranges between 0 mm−1 and 0.04 mm−1 with a step size of 0.002 mm−1; that for μs′ ranges between 0 mm−1 and 0.4 mm−1 in increments of 0.02 mm−1. The only exception is at μa=μs′ = 0 mm−1, where we had to set μs′ to a small value (0.001 mm−1) to allow the DA solver to produce valid solutions.

The signed errors (ϵDA) for Egm, Φ, Sgm and Fgm as defined in Section 2.3, at a source-detector separation of 35 mm, are computed and plotted in Figs. 2(a) to 2(d), respectively. The color scales of all error contour plots are normalized to be between −100% and 100%, with positive errors (where DA overestimates the metric compared to MC) are shown in red and negative errors (where DA underestimates MC) are shown in blue. A red-colored square indicates the literature recommended value at μa = 0.0026 mm−1 and μs′ = 0.3 mm−1. A star is used to mark the physiological value that we used to run the MC simulation. A dashed line is plotted on each panel indicating the “zero-error contour” line where DA estimates exactly match those from MC (i.e. no error). Based on these plots, we found that the literature recommended CSF optical properties can result in −8.7% error in Egm, −35% error in Φ, −48% error in Sgm and −52% error in Fgm (negative errors suggest that DA underestimates MC solutions).

We want to highlight that, for any of the given metrics, it is possible to perfectly match DA with MC values, thanks to the extra degree-of-freedom when allowing both μa and μs′ to vary simultaneously. In fact, there are an infinite number of such solutions, indicated by the continuous zero-error contour line. In other words, every combination of μa and μs′ values along this line would allow the DA solution to exactly match the expected value computed by the MC model.

3.3 Optimization of optical properties in an adult atlas model

We repeat the above computation over the Colin27 atlas, and the error contour plots between the DA estimates (Egm, Φ, Sgm and Fgm) at all combinations of μa and μs′ values and the respective ground-truth values obtained using MC are plotted in Figs. 2(e)-2(h). Again, we only show the plots using optical properties at 830 nm as an example; those at 690 nm look similar and are not shown.

Comparing to the plots derived from the layered head model, we found a few notable differences. First, while the error plots for Φ and Sgm still present the zero-error contours; Egm and Fgm report underestimated estimates when using DA over all tested μa and μs′ value combinations. In both cases, a minimal error property combination, marked by a green disk in Figs. 2(e) and 2(h) is computed for both of these metrics. Another notable difference is that the contour lines of the error surfaces for the 4 selected metrics show different shapes between the layered and atlas brain models, suggesting that the choice of brain anatomical models does have notable impact to forward solutions. Nonetheless, in both cases, the overall error surfaces from all tests demonstrate a smooth monotonic trend. The desired CSF equivalent optical properties that minimize the DA modeling errors, as either dashed line or green disks, can be readily identified from Fig. 2.

3.4 Simultaneously minimizing errors in two or more optical metrics

The results shown in Fig. 3 demonstrate that it is possible to completely eliminate or, in some cases, minimize modeling errors for a given metric when using DA in brain simulations by choosing a specific set of equivalent CSF optical properties. However, it is generally desirable to recommend a set of CSF equivalent optical properties that can simultaneously minimize two or more optical metrics. To achieve this goal, in Figs. 3(a)-3(b), we first show the overlay of the zero-error contour lines for Sgm and Φ, respectively, for the Colin27 atlas at 830 nm across 4 different separations to illustrate basic rationales when minimizing two metrics simultaneously. From Fig. 3(a), it appears that the zero-error contours for Sgm at various separations intersect each other in a compact region, indicated by a shaded circle. It is clear that the μa and μs′ values at any intersection point of two zero-error contours is able to completely eliminate the DA modeling error for both separations. For example, μa = 0.00173 mm−1 and μs′ = 0.138 mm−1 are the optimal CSF optical properties when one aims to minimize the error caused by DA for brain sensitivity at both 30 mm and 35 mm source-detector separations, with the 830 nm source. The μa and μs′ values near the center of the cluster of the intersection points, roughly located at μa = 0.0026 mm−1 and μs′ = 0.14 mm−1, is suitable to minimize the error for all 4 tested source-detector separations despite that it can not eliminate the error like those values at the exact intersection points. Similar optimization can be made on the error contour plots for fluence measurements Φ shown in Fig. 3(b). In this case, the intersection points are less clustered than those for Sgm. Nonetheless, optimal CSF property settings can be found for every pair of separations.

To generalize our findings, in Figs. 3(a)-3(d), we overlay the zero-error contour lines computed from all tested metrics and source-detector separations. We show such collective zero-error contour plots for both the 5-layer (c, e) and Colin27 atlas (d, f) at either 830 nm (c, d) or 690 nm (e, f). In these plots, red lines represent errors for Φ; green lines show the errors for Sgm, and blue lines show those for Fgm. When a zero-error contour is not found in the search space, a marker of the corresponding color is shown marking the μa and μs′ values that minimize the error of the respective metric.

It is clear that there is no single solution that can simultaneously eliminate DA modeling errors in all metrics and separations. However, we would like to highlight some general observations, from which we attempt to offer an updated recommendation for the equivalent CSF optical properties in DA.

First of all, nearly all zero-error contour lines are located below the previously recommended values (red squares), suggesting that lowering the recommended μs′ values could potentially reduce the overall modeling errors. Secondly, the brain sensitivity (Sgm, shown in green) generally shows a more clustered intersection distribution than other metrics, with the optimal μa value in the vicinity of the physiological μa value. Thirdly, simultaneously minimizing surface fluence (Φ) at multiple separations generally requires a larger μa value, but the μs′ of these intersection points is comparable to those for minimizing errors of other metrics, which is about 1/2 to 1/3 of the literature recommendation. To better guide the interpretation of these findings, in Figs. 3(e)-3(d), we draw a yellow-shaded circle on each plot indicating the rough clustered location of the zero-error contours and their intersections.

Despite the fact that there is not a single μa and μs′ combination that could minimize all metrics, we still feel strongly that recommending a single set of μa and μs′ equivalent CSF properties is still quite helpful, especially considering the wide adoption the similar recommendation from Custo et al.. Consolidating our findings described above, we suggest to lower the μs′ from the previously recommended 0.3 mm−1 to 0.15 mm−1 while maintaining μa to match the respective physiological values.

3.5 Verification of error reduction at optimized CSF property values

To verify that the updated recommendation of CSF equivalent properties for DA can lead to significantly lower errors over the previously recommended values, we recomputed the fluence distributions and the Jacobians for a source-detector separation of 30 mm, similar to those shown in Figs. 4(a) and 4(c), using the updated recommendations and show side-by-side comparisons in Fig. 4 before and after this optimization. The spatial distributions of the DA (gray) and MC (white) solutions are indicated as contour lines (the closer the match, the better); the spatially-resolved percentage errors (the lower the better) of the fluence and Jacobian are also plotted as the color map in these plots.

From fluence and Jacobian distributions in both 5-layer and atlas models, the new recommendation significantly improves the match with the ground-truth MC solutions across the domain, with particularly notable improvement in the CSF, GM and WM layers. The error distributions in all plots also become significantly more uniform across the simulated domain while shifting towards the low-error end. The remaining mismatch between DA and MC largely aggregates in the CSF region, while the peak error is significantly lowered when using the new recommendation. It is worth highlighting that the mismatch within the GM region has been improved dramatically. Because the GM region is particularly important in fNIRS data analysis, our updated recommendation will likely result in enhanced fNIRS analysis accuracy.

4. Discussions

To the best of our knowledge, this work represents the first systematic comparison between DA and MC in the handling of the low-scattering CSF tissues in brain/full-head light transport simulations. We focus on revisiting a set of equivalent CSF optical properties in DA based recommended by a widely cited work by Custo et al., understanding its limitations and seeking to revise it to achieve improved modeling accuracy. Despite the relatively straightforward methodology, we believe that our updated recommendation for modeling CSF in DA is highly significant and could have a broad impact given the widespread use of the previously recommended CSF optical property values.

Based on the results presented in the above section, we want to highlight a number of key findings. First, we demonstrate that in most of the investigated optical metrics, it is possible to exactly match DA with MC when a particular optical metric is of interest. As a matter of fact, in most tested metrics, there are an infinite number of such solutions. We believe that this is a result of the additional degree-of-freedom offered by allowing both μa and μs′ to vary. In comparison, most previous works were focused on optimizing μs′ only. Secondly, a unique combination of CSF μa and μs′ often exists to simultaneously match DA with MC when two optical metrics are considered, indicated as the intersection between any pair of zero-error contours shown in Fig. 3. Thirdly, exactly matching DA with MC solutions at more than 2 optical metrics becomes impossible, however, clusters of intersection points of zero-error curves exist and could be utilized to minimize modeling errors for a specific subset of the desired metrics. Moreover, from Figs. 2 and 3, it is clear that the previously recommended CSF optical properties tend to underestimate nearly all tested metrics when used with DA. This is indicated by the blue-colored regions where the red-square markers are located in Fig. 2 and the fact that most zero-error contours shown in Fig. 3 are below the red-square markers along the y-axis. Finally, using the distributions of the zero-error contours from all metrics, we reduce these complex configurations into a simple updated recommendation: we recommend lowering the CSF equivalent μs′ from 0.3 mm−1 in previous recommendation to 0.15 mm−1, while keeping the absorption coefficient at the physiological value.

Another advance made through this work compared to previous works is the extension from only matching the surface fluence between DA and MC to a number of fNIRS/PBM relevant metrics, including total GM sensitivity (Sgm) and brain energy deposition (Egm). On the one hand, the distinctive distributions of the zero-error contour plots for each of these metrics suggest that choosing different optical metrics to optimize could lead to different optimal CSF property settings. On the other hand, comparing to the previously recommended CSF optical properties, the optimal CSF properties across various metrics all seem to require a lower μs′ value. Similarly, our systematic benchmarks also extend the simulation domain from an atlas model in Custo et al. to also consider a layered brain model, which has also been frequently used in the literature.21, 31 While the error contour plots and zero-error curves show notable visual differences, the observations between both head models, as summarized above, are generally similar.

We recognize that there are a few limitations of this current study. First of all, our primary interest is to characterize the errors between DA with MC solutions under the influence of CSF optical properties, while assuming all other settings are identical – including the geometry of the tested brain model, the assumed optical properties of all other brain layers, even the discretization methods (i.e. meshes). One should understand that exploring variations of these other simulation assumptions, which is beyond the scope of this work, would result in different optimal CSF optical properties. However, one could apply the same methodology as presented here to optimize CSF settings for specific application needs. Secondly, the updated CSF optical properties with μs′ = 0.15 mm−1 and μa at the literature physiological value is a compromise between simplicity and desired accuracy. On the one hand, the widespread use of the the previous recommendation by Custo et al. clearly demonstrates the need for a simple and easy-to-use approximation. On the other hand, as shown in Figs. 2 and 3, minimizing various optical metrics is a complex problem and there is not a simple solution. Thirdly, the study reported here specifically focuses on minimizing the DA-to-MC errors in selected forward and reconstruction related metrics. Despite that our updated recommendation can yield a significant reduction in error, including complete elimination in many cases, whether such metric-wise improvement can also make a significant impact on fNIRS/PBM data analyses depends on the type of data analysis and application. In general, many fNIRS studies currently focus on relative changes between stimuli and baseline conditions. The systematic modeling errors revealed in this work may be partially alleviated due to the use of ratiometric measurements. Regardless, for any application where the Custo et al. CSF DA properties were useful, the new recommendation should be readily applicable and is anticipated to achieve improved results. Finally, a recent work by Hirvi et al.10 explored further refinement of CSF by separately modeling the nearly transparent CSF2 region in the inner brain and the remaining semi-diffusive CSF space. Here we limit our domain to 5-layer or Colin27 model with a single CSF region and leave the modeling of CSF2 for future works.

5 Conclusion

In summary, we have systematically revisited a widely adopted recommendation on CSF equivalent optical properties to enable fast DA modeling in fNIRS data analysis with improved accuracy. By directly comparing DA with reference solutions computed by mesh-based MC, we performed comprehensive characterizations of the errors between DA and MC across various fNIRS/PBM relevant metrics, including surface fluence (Φ), total GM energy deposition (Egm), total GM sensitivity (Sgm) and brain sensitivity fraction (Fgm), at common source-detector separations across two brain models and two wavelengths. We demonstrated that by allowing simultaneous adjustments of μa and μs′, one can exactly match many single metrics between DA and MC, with an infinite number of optimal solutions along the zero-error contours. Unique optimal μa and μs′ value combinations also exist among many pairs of metrics, indicated by the intersection between the two respective zero-error contours. After reviewing the overall distributions of the zero-error contours and their intersections, we found that the previously recommended CSF properties tend to underestimate many of the fNIRS/PBM relevant metrics due to their relatively high μs′ value and consequently higher attenuation.

To offer the community a convenient set of CSF optical properties for DA, we consolidated our findings across all combinations of the tested metrics, separations, brain models and wavelengths, and suggest a revised recommendation: one should set the CSF’s equivalent μs′ to 0.15 mm−1 and maintain μa at its physiological value. Compared with the previously suggested recommendation, we demonstrated that DA simulations using the revised recommendation not only significantly reduce the errors of particular fNIRS/PBM metrics, but also reduce the spatially resolved error throughout our anatomical models, especially in the GM region. In addition, we carefully described our computational protocol for this optimization; one could follow this protocol to rederive optimal CSF properties if additional constraints arise. Given the widespread use of the previously recommended CSF properties, we anticipate that our updated recommendation could generate a broad impact and improve the accuracy of fNIRS data analysis.

Acknowledgments

This research is supported by National Institutes of Health (NIH) grants R01-GM114365, R01-EB026998, and U24-NS124027.

Fig 1 Comparison between diffusion approximation (DA) and Monte Carlo (MC) models in fluence distributions (contour lines) and percentage errors (color maps) in (a-c) a 5-layered head model and (d-f) the Colin27 atlas with a 830 nm source. The cerebrospinal fluid (CSF) optical properties are set to (a, d) a non-physiological diffusive medium of μa = 0.0026 mm−1 and μs′ = 1 mm−1, (b, e) the assumed physiological values of μa = 0.0026 mm−1, μs′ = 0.001 mm−1, and (c, f) the literature recommended equivalent values for DA at μa = 0.0026 mm−1 and μs′ = 0.3 mm−1) while MC utilizes CSF’s physiological values.

Fig 2 Error contour plots for DA computed at a range of CSF values compared to the MC reference solutions computed with CSF’s physiological values (μa = 0.0026 mm−1, μs′ = 0.001 mm−1) in a layered head model (a-d) and atlas model (e-h). We report (a, e) GM energy deposition (Egm), (b, f) detector fluence (Φ) at 35 mm separation, (c, g) total GM sensitivity (Sgm) at 35 mm separation, and (d, h) fraction of GM sensitivity (Fgm) at 35 mm separation.

Fig 3 Aggregated zero-error contours and minimum error points for DA metrics computed at a range of CSF optical properties with MC as the reference solution. We show error minimization for GM sensitivity with an 830 nm source (a), detector fluence with an 830 nm source (b), and all metrics with 830 nm (c) and 690 nm (d) sources in the atlas model. In the layered-head model, we show error minimization of all metrics with an 830 nm source (e) and 690 nm source (f).

Fig 4 Cross sections of fluence and the 30 mm source-detector separation Jacobian in DA and MC models with a 830 nm source using previously recommended values (μa = 0.0026 mm−1, μs′ = 0.3 mm−1) and our recommended values (μa = 0.0026 mm−1, μs′ = 0.15 mm−1) with absolute error (color maps) and log-scale contours of DA (brown) and MC (white). We show a direct comparison for fluence (a, b) and sensitivity (c, d) for the layered head model, as well as fluence (e, f) and sensitivity (g, h) in the atlas model.

Table 1 Assumed absorption (μa) and reduced scattering coefficients (μs′), both in mm−1, based upon literature used to obtain the ground-truth results using Monte Carlo simulations.

	Skull/Scalp19	Cerebrospinal fluid19, 20	Gray matter32	White matter32	
Properties (mm−1)	μa	μs′	μa	μs′	μa	μs′	μa	μs′	
690 nm	0.0159	1.00	0.0004	0.001	0.02	0.88	0.07	6.00	
830 nm	0.0191	0.86	0.0026	0.001	0.03	0.70	0.09	4.29	

Disclosures

No conflicts of interest, financial or otherwise, are declared by the authors.

Code, Data, and Materials Availability

Our open-source DA solver, Redbird-m can be accessed at https://github.com/fangq/redbird-m; our mesh-based MC solver can be accessed at https://github.com/fangq/mmc. The raw data describing relative error for all metrics and combinations of CSF μa and μs′ is provided in https://neurojson.org/db/cotilab/CSF_Neurophotonics_2024.
==== Refs
References

1 Arridge S. R. and Hebden J. C. , “Optical imaging in medicine: II. Modelling and reconstruction,” Physics in Medicine & Biology 42 , 841 (1997).9172263
2 Zhu C. and Liu Q. , “Review of Monte Carlo modeling of light transport in tissues,” Journal of Biomedical Optics 18 , 050902 (2013). Publisher: SPIE.
3 Boas D. A. , Culver J. P. , Stott J. J. , , “Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head,” Opt. Express 10 , 159–170 (2002).19424345
4 Fang Q. and Yan S. , “Graphics processing unit-accelerated mesh-based Monte Carlo photon transport simulations,” Journal of Biomedical Optics 24 , 115002 (2019). Publisher: SPIE.31746154
5 Wang L. , Jacques S. L. , and Zheng L. , “MCML—Monte Carlo modeling of light transport in multi-layered tissues,” Computer Methods and Programs in Biomedicine 47 (2 ), 131–146 (1995).7587160
6 Fang Q. and Boas D. A. , “Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,” Opt. Express 17 , 20178–20190 (2009).19997242
7 Fang Q. , “Mesh-based Monte Carlo method using fast ray-tracing in Plücker coordinates,” Biomedical Optics Express 1 , 165 (2010).21170299
8 Yan S. and Fang Q. , “Hybrid mesh and voxel based Monte Carlo algorithm for accurate and efficient photon transport modeling in complex bio-tissues,” Biomedical Optics Express 11 , 6262–6270 (2020). Publisher: Optica Publishing Group.33282488
9 Yuan Y. , Yan S. , and Fang Q. , “Light transport modeling in highly complex tissues using the implicit mesh-based Monte Carlo algorithm,” Biomedical Optics Express 12 , 147–161 (2021). Publisher: Optica Publishing Group.33520382
10 Hirvi P. , Kuutela T. , Fang Q. , , “Effects of atlas-based anatomy on modelled light transport in the neonatal head,” Physics in Medicine & Biology 68 , 135019 (2023).
11 Dehghani H. , Eames M. E. , Yalavarthy P. K. , , “Near infrared optical tomography using NIRFAST: Algorithm for numerical model and image reconstruction,” Communications in numerical methods in engineering 25 , 711–732 (2008).20182646
12 Fang Q. , Carp S. A. , Selb J. , , “A multi-modality image reconstruction platform for diffuse optical tomography,” in Biomedical Optics (2008), paper BMD24, BMD24, Optica Publishing Group (2008).
13 Okada E. , Firbank M. , Schweiger M. , , “Theoretical and experimental investigation of near-infrared light propagation in a model of the adult head,” Applied Optics 36 , 21–31 (1997). Publisher: Optica Publishing Group.18250644
14 Arridge S. R. , Schweiger M. , Hiraoka M. , , “A finite element approach for modeling photon transport in tissue,” Medical Physics 20 (2 ), 299–309 (1993).8497214
15 Hielscher A. H. , Alcouffe R. E. , and Barbour R. L. , “Comparison of finite-difference transport and diffusion calculations for photon migration in homogeneous and heterogeneous tissues,” Physics in Medicine & Biology 43 , 1285 (1998).9623656
16 Fukui Y. , Ajichi Y. , and Okada E. , “Monte Carlo prediction of near-infrared light propagation in realistic adult and neonatal head models,” Applied Optics 42 , 2881–2887 (2003). Publisher: Optica Publishing Group.12790436
17 Dehghani H. , Arridge S. R. , Schweiger M. , , “Optical tomography in the presence of void regions,” JOSA A 17 , 1659–1670 (2000). Publisher: Optica Publishing Group.10975376
18 Heiskala J. , Hiltunen P. , and Nissila I. ¨, “Significance of background optical properties, time-resolved information and optode arrangement in diffuse optical imaging of term neonates,” Physics in Medicine & Biology 54 , 535 (2009).19124950
19 Strangman G. , Franceschini M. A. , and Boas D. A. , “Factors affecting the accuracy of near-infrared spectroscopy concentration calculations for focal changes in oxygenation parameters,” NeuroImage 18 , 865–879 (2003).12725763
20 Custo A. , Wells W. M. III , Barnett A. H. , , “Effective scattering coefficient of the cerebral spinal fluid in adult head models for diffuse optical imaging,” Applied Optics 45 , 4747 (2006).16799690
21 Okada E. and Delpy D. T. , “Near-infrared light propagation in an adult head model. I. Modeling of low-level scattering in the cerebrospinal fluid layer,” Applied Optics 42 , 2906–2914 (2003). Publisher: Optica Publishing Group.12790439
22 Kochunov P. , Mangin J.-F. , Coyle T. , , “Age-related morphology trends of cortical sulci,” Human Brain Mapping 26 (3 ), 210–220 (2005). eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/hbm.20198. 16161162
23 Farina A. , Torricelli A. , Bargigia I. , , “In-vivo multilaboratory investigation of the optical properties of the human head,” Biomedical Optics Express 6 , 2609–2623 (2015). Publisher: Optica Publishing Group.26203385
24 Dehghani H. , White B. R. , Zeff B. W. , , “Depth sensitivity and image reconstruction analysis of dense imaging arrays for mapping brain function with diffuse optical tomography,” Applied Optics 48 , D137–D143 (2009). Publisher: Optica Publishing Group.19340101
25 Eggebrecht A. T. , White B. R. , Ferradal S. L. , , “A quantitative spatial comparison of high-density diffuse optical tomography and fMRI cortical mapping,” NeuroImage 61 , 1120–1128 (2012).22330315
26 Eggebrecht A. T. and Culver J. P. , “NeuroDOT: An extensible Matlab toolbox for streamlined optical functional mapping,” in Diffuse Optical Spectroscopy and Imaging VII (2019), paper 11074 26, 11074 26, Optica Publishing Group (2019).
27 Zhan Y. , Eggebrecht A. T. , Culver J. P. , , “Image Quality Analysis of High-Density Diffuse Optical Tomography Incorporating a Subject-Specific Head Model,” Frontiers in Neuroenergetics 4 (2012). Publisher: Frontiers.
28 Defenderfer J. , Forbes S. , Wijeakumar S. , , “Frontotemporal activation differs between perception of simulated cochlear implant speech and speech in background noise: An image-based fNIRS study,” NeuroImage 240 , 118385 (2021).34256138
29 Collins-Jones L. H. , Cooper R. J. , Bulgarelli C. , , “Longitudinal infant fNIRS channel-space analyses are robust to variability parameters at the group-level: An image reconstruction investigation,” NeuroImage 237 , 118068 (2021).33915275
30 Holmes C. J. , Hoge R. , Collins L. , , “Enhancement of MR Images Using Registration for Signal Averaging,” Journal of Computer Assisted Tomography 22 , 324 (1998).9530404
31 Tran A. P. , Yan S. , and Fang Q. , “Improving model-based functional near-infrared spectroscopy analysis using mesh-based anatomical and light-transport models,” Neurophotonics 7 , 015008 (2020). Publisher: SPIE.32118085
32 Yaroslavsky A. N. , Schulze P. C. , Yaroslavsky I. V. , , “Optical properties of selected native and coagulated human brain tissues in vitro in the visible and near infrared spectral range,” Physics in Medicine & Biology 47 , 2059 (2002).12118601
33 Wu M. M. , Chan S.-T. , Mazumder D. , , “Improved accuracy of cerebral blood flow quantification in the presence of systemic physiology cross-talk using multi-layer Monte Carlo modeling,” Neurophotonics 8 , 015001 (2021). Publisher: SPIE.33437846
34 Hochuli R. , Powell S. , Arridge S. , , “Quantitative photoacoustic tomography using forward and adjoint Monte Carlo models of radiance,” Journal of Biomedical Optics 21 , 126004 (2016). Publisher: SPIE.27918801
35 Dehaes M. , Grant P. E. , Sliva D. D. , , “Assessment of the frequency-domain multi-distance method to evaluate the brain optical properties: Monte Carlo simulations from neonate to adult,” Biomedical Optics Express 2 , 552–567 (2011). Publisher: Optica Publishing Group.21412461
36 Brigadoi S. , Aljabar P. , Kuklisova-Murgasova M. , , “A 4D neonatal head model for diffuse optical imaging of pre-term to term infants,” NeuroImage 100 , 385–394 (2014).24954280
37 Wu M. M. , Horstmeyer R. W. , and Carp S. A. , “scatterBrains: an open database of human head models and companion optode locations for realistic Monte Carlo photon simulations,” Journal of Biomedical Optics 28 , 100501 (2023). Publisher: SPIE.37811478
38 Wu M. M. , Perdue K. , Chan S.-T. , , “Complete head cerebral sensitivity mapping for diffuse correlation spectroscopy using subject-specific magnetic resonance imaging models,” Biomedical Optics Express 13 , 1131–1151 (2022). Publisher: Optica Publishing Group.35414976
39 Fang Q. , Selb J. , Carp S. A. , , “Combined optical and x-ray tomosynthesis breast imaging,” Radiology 258 (1 ), 89–97 (2011).21062924
40 Fang Q. , Moore R. H. , Kopans D. B. , , “Compositional-prior-guided image reconstruction algorithm for multi-modality imaging,” Biomed. Opt. Express 1 , 223–235 (2010).21258460
41 Deng B. , Brooks D. H. , Boas D. A. , , “Characterization of structural-prior guided optical tomography using realistic breast models derived from dual-energy x-ray mammography,” Biomed. Opt. Express 6 , 2366–2379 (2015).26203367
42 Haskell R. C. , Svaasand L. O. , Tsay T.-T. , , “Boundary conditions for the diffusion equation in radiative transfer,” JOSA A 11 , 2727–2741 (1994). Publisher: Optica Publishing Group.7931757
43 Hernandez-Martin E. and Gonzalez-Mora J. L. , “Diffuse optical tomography in the human brain: A briefly review from the neurophysiology to its applications,” Brain Science Advances 6 , 289–305 (2020). Publisher: SAGE Publications Ltd.
44 Yao R. , Intes X. , and Fang Q. , “Direct approach to compute Jacobians for diffuse optical tomography using perturbation Monte Carlo-based photon “replay”,” Biomedical Optics Express 9 , 4588–4603 (2018). Publisher: Optica Publishing Group.30319888
