NIR light propagation in a digital head model for traumatic brain injury (TBI)

: Near infrared spectroscopy (NIRS) is capable of detecting and monitoring acute changes in cerebral blood volume and oxygenation associated with traumatic brain injury (TBI). Wavelength selection, source-detector separation, optode density, and detector sensitivity are key design parameters that determine the imaging depth, chromophore separability, and, ultimately, clinical usefulness of a NIRS instrument. We present simulation results of NIR light propagation in a digital head model as it relates to the ability to detect intracranial hematomas and monitor the peri-hematomal tissue viability. These results inform NIRS instrument design specific to TBI diagnosis and monitoring. the opposite side of the head as the baseline, and if there is a difference in the optical density on one side of the head versus the other, then there is likely a hematoma. Our simulated optical density measurements for swollen and externally bruised heads show that this may not be a valid baseline. The simulations and results presented in this paper guide the selection of wavelengths, source-detector separations, and processing algorithms to design new optical instruments for TBI detection and monitoring.


Introduction
Between 2007 and 2010, there was an increase from 1.75 million to 2.5 million confirmed cases of traumatic brain injury (TBI) in the United States [1,2]. This increase is due in part to increased awareness of TBI and improved diagnostic equipment, but it also indicates the importance of TBI and the importance of diagnosing TBI correctly. TBI can cause hemorrhage and acute changes in blood flow, oxygenation, and edema that may potentially be detected by one or more clinical imaging modalities: near infrared spectroscopy (NIRS), magnetic resonance imaging (MRI) [3][4][5], magnetic resonance spectroscopy (MRS) [5], diffusion tensor imaging (DTI) [6], computed tomography (CT) [7], single-photon emission CT (SPECT) [8], diffuse correlation spectroscopy (DCS) [9,10], and others [11]. Multiple recent review papers and meeting summaries provide excellent overviews of the successes of different imaging modalities for TBI and provide references for further study. In 2007, Tisdall and Smith presented the prevailing methods of imaging TBI and concluded that a multi-modal approach offers the best potential for improving patient outcome, but multimodal image analysis and display is challenging [12]. Combining information from two or more different imaging modalities complicates instrument design, software analysis, and deciding which actionable information to present to trauma surgeons. In 2012, the Pediatric TBI Neuroimaging Workgroup evaluated over ten emerging imaging modalities and discussed the challenges associated with standardizing data and image processing [13]. NIRS and other optical imaging modalities are favored over modalities that require ionizing radiation for infants. The complexities of pediatric neuroimaging is different than those found in adult neuroimaging due to hair density, and scalp and skull absorption thicknesses. This paper focuses on optical imaging of TBI in adults.
In 2014, at a TBI workshop in Canada, clinicians, radiologists, and other health care professionals met to discuss the roadmap for TBI neuroimaging and concluded that extensive clinical studies are required to understand how patient prognosis is comprehended by imaging [14]. Patient recruitment for extensive controlled studies is difficult because TBI patients are injured away from the hospital, so it would be valuable for these studies to be conducted with a portable instrument. Also in 2014, Davies, et al, with the UK NIH specifically reviewed NIRS monitoring of adult TBI and concluded that NIRS has the potential for adding physiological insight in conjunction with accepted methods of neurological imaging, but current "state-of-the-art" NIRS technology cannot replace CT or invasive intra-cranial pressure monitoring [15]. Optical imaging offers advantages of simpler, non-ionizing instrumentation and can image physiology Functional near infrared spectroscopy (fNIRS) is a functional neuroimaging method that quantifies absolute or relative changes in tissue function, based on HbO 2 , Hb, H 2 O, lipids, or other chromophores. Ferrari and Quaresima in their review of fNIRS estimate that nearly 700 fNIRS units were commercially available for clinical and research use in 2012 [16]. Translating fNIRS from the research community to wide use in the diagnostic and treatment community is a critical next step in the technology development. The three main methods of fNIRS in increasing levels of instrumentation complexity are continuous wave (CW), frequency domain (FD), and time domain (TD). CW fNIRS is relatively simple in instrumentation design with a continuous wave illuminator and a detector to measure the total amount of signal attenuated by the tissue, but it is limited to only detecting relative chromophore concentration and has low spatial resolution. In a review of CW fNIRS technology, Scholkmann, et al, present the transition from single measurements, to 2D imaging, to 3D imaging capability of fNIRS in the past two decades [17]. Fundamental limitations of CW fNIRS are difficult to overcome without high density optode arrangements and innovative signal processing. FD fNIRS requires slightly more complex instrumentation with a frequency modulated source and phase-locked detector, but it is able to quantify absolute concentration of chromophores and can image deeper than CW [18,19]. TD fNIRS requires the most complex instrumentation due to the need for picosecond time gating on each detector, but it also has the most promise for high resolution imaging of chromophore concentration [20,21]. In a review of TD fNIRS technology, Torricelli, et al, present the benefits of TD for imaging depth and resolution [22]. Advanced optical techniques related to TD fNIRS focus on the distribution of time of flight (DTOF) [23] and the correlation between diffuse time signals [9,10]. These advanced techniques have proven useful in non-invasive neurocritical monitoring. As source and detector technology improves in the next decades, TD instruments may have smaller footprints but they will likely remain more expensive than CW instruments.
Monte Carlo simulations have been used extensively to model the transport of light in NIRS applications for the head [19,24-31], but to our knowledge it has not been related to the detection of intracranial hematomas or the monitoring of peri-hematomal progression of TBI. We use Monte Carlo simulations to explore the potential of fNIRS as a pre-hospital screening tool to detect hematomas or as a continuous bedside monitor of TBI progression. NIRS can be implemented as either tethered or untethered [32] and monitoring in ambulances and medical helicopter transports is possible as long as the optodes are isolated from ambient light [33]. NIRS technology is also becoming smaller and more transportable for field use [34,35]. NIRS imaging has previously been developed for detecting intracranial hematomas based on the global optical density differences from a region on the head where there is a hematoma compared to a region where there is no hematoma [35][36][37][38][39][40][41]. One issue with this method is that in blunt trauma brain injury, when a patient is hit on one side of the head, a bruise or swelling appears on the side of the head that is hit and an intracranial hematoma may occur on the opposite side of the head, where the brain collides with the inside of the skull [42]. In this case, and a CT has not yet been performed, a false negative could be recorded by NIRS if the hematoma and the bruise absorb the same amount of the incident light.
For continuous bedside monitoring of TBI progression, it is clinically relevant to image ischemia and water content in the peri-hematomal region so that physicians can predict tissue viability and make informed decisions about treatment options. This is a more challenging proposition than intracranial hematoma detection because it requires quantifying changes in blood oxygenation and water content over time. TD NIRS as a continuous bedside monitor for cerebral blood flow has been explored previously [9,10,23,43], and recent advances in TD fNIRS instrumentation allows for estimation of oxygenation in the brain cortex [21].
In this study, we have run Monte Carlo simulations with multiple wavelengths in a digital phantom with and without bruising, swelling, and hematomas for both CW and TD NIRS. We first show whether NIRS is able to detect a hematoma without a non-hematoma reference. We also explore the capability of CW and TD NIRS with respect to peri-hematomal imaging.

Multi-layered digital phantoms
Monte Carlo simulations of light propagation in the digital phantoms were run to aid in understanding the design space for resolving subdural and epidural hematomas and functional imaging of the tissue perimeter surrounding a hematoma. Two multi-layered head models consisting of (1) scalp, (2) skull, (3) CSF, (4) brain (gray matter and white matter), are constructed as digital phantoms in FRED (Photon Engineering) and are illustrated in Fig. 1. The first was a parallel plate model with normal thicknesses of each layer (1) 3 mm, (2) 7 mm, (3) 2 mm, and (4) 100 mm [44].
The second digital phantom was a full head model constructed of bi-cubic mesh surfaces extracted from an MRI scan of an adult that was readily available with the SPM12 opensource software package [45] with a resolution of 1 mm x 1mm x 3 mm. The voxelated brain in MATLAB was converted into XYZ nodes at the boundaries between air, scalp, skull, CSF, and brain by edge filtering of the MRI intensity in each slice. XYZ nodes were loaded into FRED as an XY matrix with Z heights for each layer boundary in the FRED bi-cubic mesh surface definition. In both digital phantoms, the tissue optical properties for four selected wavelengths were defined according to Table 1. Wavelengths between 600 nm and 1000 nm are typical for NIRS instrumentation due to high scattering below 600 nm and high absorption above 1000 nm [15]. Absorption in the scalp and skull was not higher at 1050 nm than 950 nm, so it is more likely that fNIRS instrumentation operates at wavelengths shorter than 1000 nm because of the roll-off in quantum efficiency of silicon detectors at 1000 nm.

Continuous wave (CW) simulation methods
To study the ability to resolve hematomas located between the scalp and the cortex, we repeated multiple CW simulations of a 0.5 mm collimated source with four wavelengths: 750 nm, 850 nm, 950 nm, and 1050 nm. The wavelengths were chosen to cover spectral attenuation differences between Hb, HbO 2 , and water. An intracranial hematoma, with optical properties defined as blood (SpO 2 > 98%), was placed in the CSF layer of the parallel plate digital phantom to model a hyperacute subdural/epidural hematoma. In the first set of simulations, the hematoma was located 15 mm away from the source and increases in diameter (0 mm, 5 mm, 15 mm, 20 mm, 25 mm, and 35 mm) resulting in a 0 cc, 0.2 cc, 1.4 cc, 2.5 cc, 3.9 cc, and 7.7 cc hematoma volumes. The hematoma volumes were chosen to cover a range of small hematomas associated with mild TBI. The source was placed at the surface midpoint injecting light into the phantom, and consisted of 4.9x10 8 initial rays with Monte Carlo scattering of each ray in the tissue with loss of power at each scatter event until it: (a) escaped the scalp-air boundary, (b) scattered more than 5.0x10 4 consecutive times in one tissue layer, or (c) is reduced in power by 10 −30 . If a ray escapes the scalp-air boundary its irradiance was counted by an absorbing plane detector. The absorbing plane detector was 201x201 0.5 mm pixels centered at the source location. The number of source rays was sufficient so that multiple rays were detected in each 0.5 mm square area up to 50 mm from the source on the scalp-air boundary.

Time domain (TD) simulation methods
To simulate TD methods for hematoma monitoring, we measure the optical path length of all detected photons in the MRI-based digital head model. We present a new method of simulating the TD process by tracking optical path length in the FRED ray-traces. Since the bulk refractive index of all tissue is modeled as the same value (n = 1.37), optical path length is directly proportional to time of flight (1 mm = 4.57 ps). Refractive index of CSF is slightly lower than solid tissues, but the product of cross-sectional area and the refractive index difference is negligible compared to the contribution from solid tissues, so we can approximate a constant speed of light within the digital phantoms.
A 2.5 cc ellipsoidal hematoma, with optical properties defined as blood (SpO2 > 98%), is located ~15 mm in depth below the scalp-air surface, and ~15 mm away from the source location radially. In FRED, all detected photons with greater than 80 mm optical path length are saved and stored in 1 mm wide bins to accumulate a distribution of time of flight (DTOF). A rectangular integration of 5 bins is used to compile pixelated spatial maps of the diffuse reflectance escaping the scalp in a moving 230 ps time window starting 365 ps after the laser pulse, which is assumed to be instantaneous in this model. In instrument design, both the source pulse duration and the detector integration window should be considered and can be measured and deconvolved by signal processing. This paper focuses only on understanding the detector integration window required for maximum contrast.

CW simulations on an intracranial hematoma model
The CW simulations described in Section 2 result in spatial profiles of the irradiance on the array of detectors for each simulation run. A sample of simulation runs with different hematoma sizes outlined in red is shown in Fig. 2. The hematoma is placed 15 mm away from the source, which is located near the center of peak contrast as shown later in Fig. 6. Additionally, a vertical line profile through the source and hematoma is sampled and plotted in Fig. 3. These simulations support previous conclusions that there is an optical density difference between the uninjured brain and brain with a hematoma when the source-detector separation (S-D) is greater than 30 mm [36-38]. A 0.2 cc hematoma is slightly distinguishable as the diffuse reflectance profile begins to separate from the baseline "no hematoma" profile, and a 1.4 cc hematoma is definitely distinguishable from the baseline. The lack of separation between the 3.9 cc and 7.7 cc profiles in Fig. 3 show that as the hematoma increases above a critical size, the contrast in the diffuse reflectance stops increasing. In this paper, contrast is defined as the difference in reflectance between the "no hematoma" case (blue) and the "hematoma" case (green), divided by the "no hematoma" case, i.e. contrast = (bluegreen)/blue. The results are presented to inform instrument design for detector sensitivity and dynamic range requirements based on different S-D and wavelength. Though the plots versus wavelength are qualitatively similar, the quantitative differences are important for multiwavelength instrument design. Fig. 2. Irradiance on a 100 mm x 100 mm detector array for a hematoma located 15 mm below the source (outline shown in red circle). There is a decrease in irradiance at the bottom of the detector array correlated to the increase in hematoma size, and almost no change at the top of the array. Additional analysis from the first round of simulations shows that resolvable contrast decreases when the detector noise floor increases (Fig. 4). Resolvable contrast is considered as the point where the contrast between baseline and a 2.5 cc hematoma is 50%, and the signals below the noise floor are calculated based relative attenuation from a source with intensity set for the maximum permissible exposure (MPE) in accordance with ANSI Z.136.1 standards for skin exposure for each wavelength. For example, the MPE for 750 nm is 2.4 J/cm 2 or 0.24 W/cm 2 for a 10 second duration exposure. Figure 4 is important for instrument design because it shows how contrast increases as S-D increases, but only when the detector noise floor is low enough to support large S-D distances. The noise floors plotted are higher than what is typically achieved by commercial CW NIRS instrumentation. This is intentional to drive toward an instrument design where less sensitive detectors may be implemented. Fig. 4. Resolvable contrast for a 2.5 cc hematoma as a function of detector noise floor (shown in legend). This assumes the hematoma is centered 15 mm from the source. If detector noise floor is higher than 10 nW, then 50% contrast may not be achievable, and there will be no contrast if the detector is spaced 30 mm away from the source.

Effects of bruising and swelling
Baseline measurements of diffuse cranial reflectance may be difficult to attain if bruising or swelling exists, which is the case with most patients suffering from blast or collision. To help understand the impact of bruising and swelling, the second set of simulations were run with a 2.5 cc hematoma in a bruised digital phantom and a swollen digital phantom. In this paper, we use the term "bruise" to refer to only to blood in the scalp and skull and "hematoma" to refer to intracranial bleeding. The bruised digital phantom was constructed by adjusting the tissue properties of the scalp and skull to include 10% blood (µ scalp,new = 0.9µ scalp + 0.1µ blood and µ skull,new = 0.9µ skull + 0.1µ blood , respectively for absorption and scattering coefficients). The swollen digital phantom was constructed by increasing the thickness of the scalp layer to 6 mm -twice its original thickness. Only two wavelengths were chosen (750 nm and 950 nm) because the previous results showed enough similarity between the four wavelengths. Figure  5 compares the diffuse reflectance profiles for bruised, swollen, and normal simulations. The bruised case greatly decreases the diffuse reflectance for all S-D at both wavelengths, but more drastically at 950 nm where the blood absorption coefficient is higher. At S-D greater than 30 mm, the "Bruised" lines seem to plateau which is likely due to higher absorption, which decreases the number of detected rays in the bruised digital phantom. Comparing the "Bruised" and "Swollen" lines to the "Normal+Hematoma" line indicates that using a baseline reading from a bruised or swollen region of a patient's cranium may appear to degrade the signal in the same way as a hematoma. This could result in false negative readings when detecting hematomas using bi-lateral symmetry for a baseline. "Normal" phantom is with tissue properties as described in Table 1, "Bruised" phantom assumes 10% extra blood in scalp and skull tissues, and "Swollen" phantom assumes a scalp with double thickness. The hematoma for the dotted lines is 2.5 cc located 15 mm from the source.

CW simulations for locating hematomas
Locating and determining the size of hematomas may be possible with multiple sourcedetector pair spaced over the patients cranium. To model requirements for optode spacing, a 2.5 cc hematoma located between the cortex and the scalp is moved axially from −40 mm to + 40 mm in relative to the source location in 5 mm increments. The irradiance incident on a 5 mm detector spaced + 30 mm from the source is calculated for each location of the hematoma and the results are presented in Fig. 6. The results show that a 20 mm diameter hematoma broadens to approximately 40 mm in width (FWHM of the green curves) when measured at the surface of the scalp, assuming a grid of 5 mm spaced sources and detectors. Both wavelengths show similar contrast curves even though the total reflectance decreases as wavelength increases. This implies that the wavelengths are equivalent in imaging depth and contrast, therefore any differences in reflectance by wavelength can truly be attributed to functional changes (i.e. ischemia or water content).

TD simulations for hematoma detection
Results of the TD simulations were analyzed as DTOF for each wavelength and for multiple source detector separations out to 40 mm. DTOF is a standard way to view the detected laser pulse as it is broadened by scattering and attenuated by absorption in the biological tissues. Figure 7 shows the DTOF for a detector located 25 mm from the source in the two cases: with and without hematoma present. The "No Hematoma" curve matches analytical results from other publications [22] and the "Hematoma" curve represents the attenuation caused by the 2.5 cc hematoma located ~15 mm deep and ~15 mm away from the source. The DTOF signal when the hematoma is present is nearly the same for all wavelengths, whereas the DTOF signal when no hematoma is present is shifted in time and intensity across wavelengths. The figure shown is for a single S-D separation, but the DTOF peak shifts to the left in time as S-D decreases and shifts to the right as S-D increases.
To evaluate the improvements in contrast achievable by TD processing, the detected reflectance for increasing integration time windows is presented in Fig. 8. All time windows are centered about 1,170 ps which is the DTOF peak for a 30 mm S-D separation. In the left of Fig. 8, at short S-D separations (<10 mm) the magnitude of reflectance increases significantly as the integration time window increases. In the right of Fig. 8, the contrast between "Hematoma" and "No Hematoma" is slightly higher for shorter integration time windows at medium S-D separations (10 mm < S-D < 25 mm). However, at large S-D separations (>25 mm), it appears that there is no significant difference in contrast between short and long time windows. This indicates that at large S-D separations, TD NIRS may not offer more contrast enhancement in the detection of hematomas as compared to CW NIRS. However, the separation in the 230 ps and 1,840 ps curves in the contrast plot of Fig. 8 indicates that TD NIRS may show increased contrast with medium S-D separations as compared to CW NIRS.

Discussion
Monte Carlo simulation results are presented to inform the design of NIRS instrumentation specifically for aiding in the diagnosis and treatment of TBI. Blast, blunt force, and collision induced TBI typically results in external contusions as well as internal contusions, sometimes on opposite sides of the head. For this reason, we simulated the effects of swelling and bruising on the CW NIRS signal for multiple wavelengths at multiple source detector distances. Due to these effects, improved instruments for detecting hematomas caused by TBI should look at more than a single S-D spacing for bilateral optical density differences. To facilitate placement of sources and detectors on patient heads, many companies have developed caps with discrete spacing between sources and detectors [46]. It may be valuable to design caps with more flexible light coupling to the scalp at a continuum of S-D separations to monitor peri-hematomal tissue function in addition to detecting and locating hematomas.
Hematoma sizes modeled in this paper are on the small end of those seen in mild to severe TBI patients. In a parallel study, CT imaging on four TBI patients found subdural hematomas of 1 cc -55 cc in size. Monte Carlo simulations (Fig. 3) indicate that for S-D separations out to 35 mm, there is no contrast advantage between a 3.9 cc (25 mm diameter) hematoma and a 7.7 cc (35 mm diameter) hematoma when they are both located between the source and detector. A 35 mm diameter hematoma is more likely to be detected, however, than a 25 mm hematoma though because of the probability is higher that a significant portion of the hematoma is between the source and detector. That means that detection of small hematomas likely requires more sources and detectors or more measurement location over the entire head to get adequate coverage for detecting mild TBI. From a clinician or emergency medical technician perspective, it is preferred to have a single instrument with multiple source/detector pair to improve the speed of scanning a patient's entire head.
Time domain approaches have advantages in imaging depth for NIRS, and we have developed a method to use commercial non-sequential raytracing software to simulate TD NIRS in 3D head models. Optical design files can be easily imported into the same software, which makes this new method of simulation very useful for instrument design and analysis. Comparison of TD simulations in an MRI-based head model showed that TD with short integration time windowss (230 ps) only has minimal advantages over long integration time windows (1,840 ps) for detecting the presence or absence of hematomas. Contrast is invariant from 1 nanosecond windows and longer. This means that cheaper electronics and larger detectors (with higher SNR) may be implemented in future TBI instrumentation.
There is still a challenge of having baseline measurements taken from a healthy patient to compare with measurements taken after the patient has experienced a potentially traumatic event. Previous approaches have been to use the same cranial location on the opposite side of the head as the baseline, and if there is a difference in the optical density on one side of the head versus the other, then there is likely a hematoma. Our simulated optical density measurements for swollen and externally bruised heads show that this may not be a valid baseline. The simulations and results presented in this paper guide the selection of wavelengths, source-detector separations, and processing algorithms to design new optical instruments for TBI detection and monitoring.