Accuracy of homogeneous models for photon diffusion in estimating neonatal cerebral hemodynamics by TD-NIRS

: We assessed the accuracy of homogenous (semi-infinite, spherical) photon diffusion models in estimating absolute hemodynamic parameters of the neonatal brain in realistic scenarios (ischemia, hyperoxygenation, and hypoventilation) from 1.5 cm interfiber distance TD NIRS measurements. Time-point-spread-functions in 29- and 44-weeks postmenstrual age head meshes were simulated by the Monte Carlo method, convoluted with a real instrument response function, and then fitted with photon diffusion models. The results show good accuracy in retrieving brain oxygen saturation, and severe underestimation of total cerebral hemoglobin, suggesting the need for more complex models of analysis or of larger interfiber distances to precisely monitor all hemodynamic parameters.


Introduction
Brain injuries are an important cause of morbidity and mortality in neonates. Cerebrovascular disturbances, such as impaired cerebral hemodynamics and autoregulation (CA), are involved in the pathogenesis of brain injury [1][2][3][4]. CA is the ability to maintain constant the cerebral blood flow, despite changes in perfusion pressure [5,6]; and it may be impaired in many neonatal conditions such as prematurity [7], hypoxic-ischemic encephalopathy [8], intracranial hemorrhage [9], congenital heart disease [10]. Failure of CA leaves brain tissues unprotected against strong variations of blood flow and could increase the risk of brain injury and mortality, making cerebral perfusion and oxygenation monitoring techniques necessary in Neonatal Intensive Care Units. Nowadays, many of the techniques clinically used [11,12] such as positron emission tomography, contrast enhanced magnetic resonance imaging (MRI), Xenon clearance are bulky, need patient transport and immobility, and some of them are invasive or make use of radioactive tracers. Transcranial Doppler ultrasound is widely used in neonatology to measure cerebral blood flow velocity in target artery, but its results are limited to the macro-vasculature [13], and in particular clinical derangement the macro-vasculature blood flow may not reflect the one in cerebral region [14]. Moreover, it is operator-dependent [15] and it does not allow continuous monitoring cerebral hemodynamic parameters. In this scenario, Near Infrared Spectroscopy (NIRS) appears to be an interesting alternative: it is a non-invasive technique which allows for monitoring brain hemodynamics and CA, at bed side [16,17].
NIRS techniques exploit the low absorption of near infrared light by the main biological tissues' constituents, allowing light to penetrate a few cm inside tissues. Differences in absorption spectra of tissues' constituents (such as oxygenated (O 2 Hb) and deoxygenated (HHb) hemoglobin) permit to evaluate their concentration variations in the probed area.
In continuous wave NIRS [18] (CW-NIRS) light of constant intensity is injected in tissues and by measuring backscattered light amplitudes, at different wavelengths, variations of hemoglobin concentration in the probed tissue are retrieved. It is the most widespread NIRS technique, frequently used in neonatology thanks to its lower level of complexity and its stronger presence in the market. As draw back, it does not allow to directly estimate absolute concentrations of HHb and O 2 Hb, when a single interfiber distance is used, but only their variations.
Frequency domain NIRS (FD-NIRS) [19] exploits variations, in amplitude and phase of modulated light, due to the interaction of light with biological tissues, whereas time domain NIRS (TD-NIRS) [20] profits from variations in shape and amplitude of injected light pulses with short (i.e. picosecond) duration. FD-NIRS and TD-NIRS are more complex techniques as compared to CW-NIRS, but allow to non-invasively estimate absolute concentrations of O 2 Hb and HHb, therefore providing estimates also for total hemoglobin (tHb = O 2 Hb + HHb) and tissue oxygen saturation (S t O 2 = O 2 Hb /tHb). Moreover, both FD-NIRS and TD-NIRS have superior depth sensitivity compared to CW-NIRS [21]. NIRS estimates of the hemodynamic parameters are influenced by probed tissues' structures and geometry. Indeed, many studies have been published on the effects of extra cerebral tissues (ECT) on NIRS measurements in the adult head [22][23][24].
When NIRS techniques are used for monitoring neonatal brain hemodynamics, a precise estimation of vital parameters such as brain tHb and S t O 2 is essential. Thus, a study on the influence of superficial tissues and head structure is fundamental. The differences in size and curvature between adults and neonates head geometry suggest that separate studies need to be performed. Concerning neonatal head structure, previous works have been addressed to determine the influence of superficial layers and the accuracy in estimating optical properties when CW-NIRS [25,26] and FD-NIRS [27] are employed. Regarding the accuracy in estimating cerebral hemodynamic parameters in neonates, the only study that have been published up to now is the one by Barker et al. [28] in 2014. They investigated the accuracy of FD-NIRS in estimating brain S t O 2 and tHb. Few TD-NIRS studies have been conducted to determine the optical and hemodynamic properties of neonatal head [29,30] and some TD-NIRS clinical studies have been performed to monitor neonatal cerebral hemodynamics [31,32]. To our knowledge no systematic study on the TD-NIRS accuracy in estimating neonatal brain hemodynamics has been performed. In our opinion, precise estimation of critical parameters as tHb and S t O 2 is fundamental to adopt TD-NIRS as monitoring technique of neonatal brain health.
Therefore, in this work we assessed the accuracy of homogeneous models for photon diffusion in estimating neonatal cerebral hemodynamics by TD-NIRS. We simulated light propagation inside realistic meshes of 29 and 44 weeks' postmenstrual age (PMA) neonates' head structures, through Monte Carlo (MC) method. Absorption coefficient of brain tissue was varied to model realistic hemodynamic scenarios common in preterm neonates (ischemia, hyperoxygenation and hypoventilation). We mimicked realistic photon distribution of time of flight (DTOF) by convolving the simulated time point spread function (TPSF) with the instrument response function (IRF) of the BabyLux device [31]. Finally, the DTOF were analysed with two homogeneous models for photon diffusion: semi-infinite and spherical.

Mesh-based MC method
TD-NIRS measurements on neonates' heads were simulated with the MC method. In particular, a mesh-based MC tool (MMCLAB [33]) was used to simulate photon migration inside a realistic atlas of three dimensional (3D) head structure [34]. The head models used in our simulations were selected from an online dataset, containing meshes of preterm and term neonates' heads from 29 to 44 weeks PMA [34]. We simulated photon propagation in the two extreme cases, 29 and 44 weeks PMA, to maximize the influence that head structures, dimensions, and shapes would have on the retrieved optical and hemodynamic properties.
The head structures were divided into three domains (see Fig. 1) with different optical properties (absorption coefficient µ a , and reduced scattering coefficient µ' s ): ECT (comprehensive of scalp and skull), cerebrospinal fluid (CSF), and brain (composed by white matter, grey matter, cerebellum and brainstem).

Optical and hemodynamic parameters
The optical parameters used as baseline in our simulations were the one selected by Barker et al. [28] and are summarized in Table 1. For the ECT, the optical parameters used are the one of the scalp, indeed skull thickness is extremely thin [35], suggesting us to approximate its contribution as negligible. The anisotropy parameter g of head tissues was considered homogeneous and set to 0.89 [36]. The refractive index n was assumed constant over wavelength. The absorption coefficient of brain tissue was modified to mimic three hemodynamic scenarios common in preterm neonates: reduction of brain S t O 2 from 70% to 40% (scenario 1: ischemia); increase of brain S t O 2 from 70% to 100% (scenario 2: hyperoxygenation); reduction of brain S t O 2 from 70% to 40% with simultaneous increase of tHb from 53 µM to 93 µM (scenario 3: hypoventilation). The values of O 2 Hb, HHb, tHb, S t O 2 and the correspondent absorption coefficient (derived by Beer's law [28]) simulated in the brain are reported in Table 2. In deriving the hemodynamic parameters a water fraction of 90% was assumed in the brain [37]. No hemodynamic variations were simulated in ECT and CSF, and their optical properties were maintained constant for each one of the considered scenarios. The inter-subject variability of tissue optical and hemodynamic properties was reproduced by repeating all simulations considering other eight baseline conditions: i) two different values of µ' s for ECT (±25%); ii) two different values of µ' s for brain (±25%); iii) two different values of µ a for the ECT (±25%); iv) two different values of µ a for the brain (modeling tHb variations ±25%). A scheme of the performed simulations and their use for estimating the error in the optical and hemodynamic parameters is reported in Table 3.

MC simulations
In each simulation, 10 9 photons were injected in the Afp7 position of 10-20 International system of EEG electrode placement, and the detection point was positioned 15 mm away measured along the scalp surface (Fig. 2). This value was selected to mimic the behavior of BabyLux device [31]: TD-NIRS device designed for monitoring cerebral hemodynamics of preterm neonates.
In our simulations, µ a values were initially set to 0 cm −1 for all tissue types, and the absorption effect was added a posteriori [38] by decreasing the photon weight by exp(-µ a,i L i ), where µ a,i and L i are respectively the absorption coefficient and the pathlength traveled by photons in the i-th tissue type. To model a realistic DTOF, the simulated TPSFs were convolved with the IRF of the BabyLux device [31], after IRF normalization to unit area. The IRF (Fig. 3) used were Table 3. Scheme of the optical and hemodynamic parameters retrieved with the corresponding variations of the baseline, for the three scenarios (Table 2).

Simulation #
Change with respect to Table 1 Used for estimating error in µ a brain , tHb brain , and S t O 2 brain Used for estimating error in µ' s brain 0  characterized by full with at half maximum (FWHM) of 190 ps, and time bin of 25 ps. Shot noise was added, as Poisson distribution [39], while other noise sources (thermal noise and dark current) were considered negligible in the case of hybrid photomultiplier tube, detector type used in BabyLux device [31]. The collected photons were approximately 10 5 and 5·10 4 in the case of 29 and 44 weeks PMA mesh, respectively. To reduce statistical noise, all simulations were repeated 10 times, and the results averaged. All simulations were performed at 690 nm and 830 nm.

DE model for the fitting procedure
Absorption and reduced scattering coefficients were retrieved by fitting the simulated DTOF with the solution of the Diffusion Equation for both a homogeneous semi-infinite model [40] and a homogeneous spherical model [41]. In the latter case, to take into account the effect of the finite volume of the head, we approximated the neonates' head to a sphere, considering the distance between Inion and Nasion positions (Iz and Nz) as the diameter of the sphere: d = 2r = (Iz-Nz). The source detector distance ρ = 15 mm was measured along the sphere surface, and the angular distance θ between injected light and photon collection was computed as the angle underlying the circumference arc of length equal to the inter-fiber distance (θ =ρ/r). The partial current boundary condition (PCBC) was applied to solve the DE for both the semi-infinite and the spherical homogeneous model.
The accuracy of hemodynamic parameters estimated with the semi-infinite and spherical homogeneous model for photon migration was evaluated through the mean absolute percent error (MAPE) computed over the 10 repetitions of simulation 0 in Table 3: where X is the desired hemodynamic parameter. In Sec.3, the mean MAPE for the three studied scenarios (Table 2) is reported, estimated as the average over the 70 simulations done for each of these conditions (the 10 repetitions of all the 7 conditions, reported in Table 2, for each scenario of simulation 0 of Table 3). Variations of optical parameters (simulations from 1 to 6 in Table 3) generate a dispersion (D X ) of estimated MAPE x , which was computed as: Where, X sim 0 is the hemodynamic parameters retrieved for simulation 0 reported in Table 3, and X sim #n is the hemodynamic parameters retrieved for simulation #n reported in Table 3, varying n from 1 to 6.

Optical parameters
The reduced scattering coefficients retrieved by analyzing simulated DTOF with the semi-infinite homogeneous model at 690 nm and 830 nm are reported in Fig. 4, for the three scenarios described in Table 2. The corresponding absorption coefficients can be founded in the Supplement 1 (Sec. S2.1.1). In all cases the estimated reduced scattering coefficient shows an underestimation with respect to the nominal values (Fig. 4), and that is more evident in case of 29 weeks PMA mesh. The effect of inter subject variability (variations of ±25% of µ a ECT , µ a brain and µ' s ECT ) is reported as shadows determined by the maximum and minimum values of the estimated parameter. The parameter which mostly influences the range of variation of µ' s was µ' s ECT . Variations of ±25% of all the other optical properties (µ a ECT , µ a brain ) weakly affect the retrieved µ' s .

Hemodynamic parameters
Hemodynamic parameters tHb and S t O 2 for the three analyzed conditions in Table 2 are presented in Fig. 5. HHb, O 2 Hb, and the corresponding absorption coefficients, can be founded in Supplement 1 (Sec. S2.1.1 and Sec. S2.1.2). The shadows represent the range of variations due to intersubject variability of optical properties (±25% of µ a ECT , µ' s brain , and µ' s ECT ). tHb estimated with the semi-infinite homogeneous model is lower than the nominal values in the brain. Comparing the results obtained for the two meshes, it is evident that tHb retrieved for 29 weeks PMA mesh is underestimated with respect to 44 weeks PMA mesh. Indeed, MAPE tHb is higher and it reaches its maximum value when scenario 3 was simulated (Table 4). Conversely, S t O 2 estimates are close to the values simulated in the brain, and the MAPE StO2 is much lower than MAPE tHb (Table 4). Also, the range of variations due to intersubject variability is reduced compared to tHb, as the reduction of dimensions of blue and red shadows in Fig. 5 suggests.

Spherical homogeneous model analysis
To reduce the error related to surface curvatures, and consider the finite volume of the neonates' heads, the simulated data were analyzed with a spherical homogeneous model for photon migration. The radius of the spheres was computed as half of the distance between Inion and Nasion and resulted to be 41 mm and 59 mm for 29 and 44 weeks meshes, respectively.  Table 2.

Optical parameters
The estimated reduced scattering coefficient, at 690 nm and 830 nm, are shown in Fig. 6. The values reported are the one estimated for simulation 0 in Table 3, and the shadows represent the range of variability due to intersubject variability (± 25% µ' s ECT , ± 25% µ a ECT , ± 25% tHb brain ). The results are like those discussed in 3.1.1. Also, in this case the estimated scattering coefficients is lower than the nominal one in the brain, and in case of 29 weeks PMA mesh the underestimation is stronger. Moreover, as described in paragraph 3.1.1, µ' s ECT is the parameter which mostly influences the estimates. Comparing the results shown in Fig. 6 and Fig. 4, the scattering coefficient retrieved with the spherical homogeneous model are higher than the one estimated with the semi-infinite homogeneous model, therefore closer to the nominal values.  Table 2.

Hemodynamic parameters
Hemodynamic parameters estimated with the spherical model for photon migration for the three scenarios of Table 2 are summarized in Fig. 7. Also in this case the reported values are the one estimated for simulation 0 of Table 3, and the shadows represent the range of variations due to intersubjact variability (± 25% µ' s ECT , ± 25% µ' s brain , ± 25% µ a ECT ). Obtained results are like those reported in 3.1.2 for the semi-infinite homogeneous model. tHb is underestimated, and it appears more evident in case of 29 weeks PMA mesh. The tHb retrieved with spherical model shows lower values than the one obtained with the semi-infinite homogeneous model analysis ( Fig. 5 and Fig. 7), and the corresponding MAPE tHb are higher (Table 5). Conversely, as in 3.1.2 the estimated S t O 2 shows values like the nominal ones in the brain. Slopes of the curves reproducing estimated S t O 2 values (red and blue curves in Fig. 7 b, d, f) are higher than in semi-infinite homogeneous model, and more similar to the nominal S t O 2 values in the brain (grey curves in Fig. 7 b, d, f).  (Table 3)

Discussions
In this work we aimed to assess the accuracy of homogeneous models for TD-NIRS in measuring brain optical properties and hemodynamic parameters in neonates. Mesh-based MC tool was used to simulate photon propagation inside realistic atlases of head structures. Neonates' fast growth, and consequent change of head structure, was considered by analyzing TD-NIRS signals in two meshes relative to 29 and 44 weeks PMA neonates' head. The interfiber distance was set to 15 mm to mimic the behavior of the BabyLux device. The choice of the 15 mm interfiber distance was a tradeoff between average penetration depth through neonate head tissues and bulkiness of the probe. The thickness of superficial tissues in the region (5 mm of radius around injection and detection position) of our simulations were: for ECT 2.7 mm and 3.5 mm, for CSF 2.5 mm and 3.9 mm, in 29 and 44 weeks PMA meshes, respectively. The thin layer of ECT and the small absorption and scattering coefficients of CSF allow detected light to penetrate in the deep tissue of the brain. Indeed, the percentages of pathlength travelled by photon in the three domains (for simulation 0 of Table 3) are: ECT 48% (55%), CSF 17% (11%), brain 35% (34%) for 29 (44) weeks PMA mesh. Simulations results were then analyzed with two analytical models for photon migration: semi-infinite homogeneous model, which is one of the most widespread models of analysis because of its low level of complexity; and spherical homogeneous model, considered to reduce the effects of head curvature and finite volume. The radii of the spheres used to represent neonates' head were chosen, as reported in Sec. 2.3, as a half of the distance between Inion and Nasion position, allowing to adapt this model of analysis also in case of real life measurements. In this section the results obtained, presented in Sec. 3, will be discussed.

Effects of superficial layers
The reduced scattering coefficient shows a systematic underestimation with both the semi-infinite (Fig. 4) and the spherical (Fig. 6) homogeneous model. The error is higher in case of spherical homogeneous model and 29 weeks PMA mesh. It is well known as reported by Comelli et al. [42], that, when a homogeneous model is used, the reduced scattering coefficient is strongly influenced by superficial layers. The presence of CSF, which is a low scattering media, could affect the retrieved µ' s , in both models, and it is probably the cause of the µ' s underestimation. Ancora et al. [24] observed no peak shift on simulated DTOF when CSF thickness was increased, which suggests that an increase of CSF does not influence the retrieved scattering coefficient. However, they simulated an adult head structure, with mean ECT thickness of 12.5 mm. Early photons, which are the ones that mostly affect peak shift, have very low probability to travel inside CSF layer. In case of our work, instead, photon propagation was simulated on newborns' head structures, which are characterized by a lower ECT thickness (about 3 mm of thickness), making the probability of early photons to travel inside CSF layer not negligible, and with the consequent result of affecting the retrieved scattering coefficient. Head dimensions and its structure strongly affect sensitivity of TD-NIRS in retrieving brain optical properties. Indeed, the differences obtained in the estimated scattering coefficients of the two meshes for both model of analysis (see Fig. 4 and Fig. 6) could be related to the different thickness of ECT layer in 29 and 44 weeks PMA head structures. ECT is thicker in 44 weeks PMA mesh, with respect to 29 weeks PMA mesh (3.5 mm and 2.7 mm in case of 44-29 weeks PMA mesh, respectively), leading to a stronger influence of that tissue in the estimation of reduced scattering coefficient. This could be the cause for the higher value of µ' s retrieved in case of 44 week PMA mesh. Indeed, being µ' s ECT higher than µ' s CSF and µ' s brain , an increment of ECT thickness could explain the enhancement in the estimation of reduced scattering coefficient for 44 weeks PMA mesh. Thus, the smaller underestimation of scattering coefficients in the 44 weeks PMA head structures is likely due to a combination of the two effects discussed in the previous rows related to the presence of both ECT and CSF: the presence of the CSF in very shallow layer causes a strong reduction in the estimated µ' s ; on the other side the presence of the ECT causes an increase of the estimated µ' s and compensates the CSF effects. The higher the ECT thickness, the stronger is its influence, and its compensation of the CSF effect. Indeed, the percentage pathlength travelled by detected photons in the CSF is higher in 29 weeks PMA mesh (17%) than in 44 weeks PMA mesh (11%), while the percentage pathlength travelled in the ECT presents an opposite behaviour: it is higher in 44 weeks PMA mesh (55%), and lower (48%) in 29 weeks PMA mesh. Therefore, in 44 weeks PMA mesh the estimation error is lower, even if the superficial layer thickness is higher.
Moreover, as reported in Secs. 3.1.1 and 3.2.1, the parameter which mostly influences the range of variation (red and blue shadows in Fig. 4 and Fig. 6) of estimated µ' s is µ' s ECT . Combining the two discussed effects, that is higher thickness of ECT for 44 weeks PMA mesh, and the strong influence of µ' s ECT variations on retrieved scattering coefficient, we can justify the larger dimension of the shadows in case of 44 (red shadows in Fig. 4 and Fig. 6) with respect to the case of 29 weeks PMA mesh (blue shadows in Fig. 4 and Fig. 6).

Effects of head curvature
Comparing the results obtained with semi-infinite and spherical homogeneous model for photon migration, the estimated values of reduced scattering coefficient are smaller when semi-infinite homogeneous model is applied (see Fig. 4 and Fig. 6). This result is in accordance with previous findings: Sassaroli et al. [41] investigated the accuracy of semi-infinite, cylindrical and spherical homogeneous models of analysis in retrieving optical properties of MC simulations on spherical geometry. They observed an underestimation of scattering coefficient, which increases with source-detector distance, when the semi-infinite homogeneous model was applied. They also observed small errors when wrong radius of curvature was used in the analysis, which empower the solidity of our approximation in measuring sphere radius. This result highlights the influence of head curvature on retrieved optical properties, and the importance of considering an appropriate model to take in to account head curvature.

Effects of superficial layers
As in case of reduced scattering coefficient, also total hemoglobin concentration was strongly underestimated (see Fig. 5 and Fig. 7). In general, optical properties are affected by the presence of superficial layers [42], and variations of absorption coefficients due to ECT and CSF layer are reflected in retrieved hemodynamic parameters. Also in this case, the presence of CSF layer affects hemoglobin concentration, and could be the cause of its underestimation. As reported from Ancora et al. [24], by analyzing the behavior of late-photons, a CSF layer thickness increase causes DTOF curves to exhibit less steep slopes, suggesting a reduction in the absorption coefficients retrieved from DTOF curves. Reduction of the absorption coefficient causes tHb underestimation, which is in accordance with previous findings. For instance, Barker et al. [28] simulated multidistance FD-NIRS measurements on neonate's head structure, considering different source detector distances (10-25 mm, 15-30 mm, 20-35 mm, and 25-40 mm), and obtained lower concentration of tHb compared to that simulated in the brain. Discrepancies with respect to brain tHb were lower when source detector distance was increased, confirming that the possible cause could be the influence of superficial layers.
As in case of reduced scattering coefficient, results for 44 weeks PMA mesh show lower MAPE tHb , due to the thicker ECT layer and a less important effect of the CSF layer.
It is worth noting that S t O 2 is weakly influenced by ECT and CSF layers, probably O 2 Hb and HHb underestimation compensate in some way when saturation is computed, and a better sensitivity to brain saturation is obtained (with MAPE StO2 much lower than in case of tHb, see Table 4 and Table 5).

Effects of head curvature
Comparing the results presented in Fig. 5 and Fig. 7, it appears evident that values retrieved for tHb are lower when spherical homogeneous model is used to analyze simulated DTOF. This result is a consequence of the lower absorption coefficient obtained with spherical homogeneous model, than with semi-infinite homogeneous model, and it agrees with literature. Indeed, Sassaroli et al. [41] observed an overestimation of absorption coefficient when semi-infinite homogeneous model was used to analyze simulations on spherical (and cylindrical) geometries.
Concerning S t O 2 , the main difference among the two models is the slope of the curves representing the estimated S t O 2 (Table 6), which is higher in case of spherical homogeneous model. It appears more evident in case of Ischemia ( Fig. 5(b) and Fig. 7(b)), where the angular coefficient of the linear regression curve fitting the retrieved S t O 2 rises from 0.77, for semi-infinite homogeneous model, to 0.98 for spherical homogeneous model: this fact allows to better retrieve S t O 2 variations. On the contrary, a slight increase of MAPE StO2 can be observed in case of spherical homogeneous model. Table 6. Angular coefficient of the linear regression curve that best fit S t O 2 estimated from simulation 0 (Table 3) with the semi-infinite and spherical homogeneous model, averaged over 10 repeated simulations of each scenario ( Differences between the two models of analysis are stronger when simulations on 29 weeks PMA mesh are considered, probably as consequence of the smaller radius of the sphere (41 mm) used to model head structure, compared to the one used in case of 44 weeks PMA mesh (59 mm). Indeed, the higher is the radius, the better is the semi-infinite approximation as model of analysis.

Effects of the increase of tHb
As presented in Table 2, S t O 2 values simulated in the brain for scenario 1 and 3 were the same (from 70% to 40%), however tHb is constant in case of scenario 1, and increases in case of scenario 3. As a consequence, brain absorption coefficients at the two wavelengths were higher in case of scenario 3 ( Table 2). As suggested by Fig. 5(b) and (f) (and in case of spherical homogeneous model Fig. 7(b) and (f)), an increase of absorption coefficient of the brain affects not only the estimation of tHb, but also the retrieved S t O 2 . The slope of the S t O 2 curves retrieved for scenario 3, in Fig. 5(f) (and Fig. 7(f)), is smaller than that one in Fig. 5(b) (and Fig. 7(b)) for scenario 1, suggesting that an increase of brain tHb (and absorption coefficient) reduces the sensitivity of TD-NIRS in measuring brain saturation.
Moreover, the tHb increase leads to shadows enlargement. In our opinion, it could be related to an increase of the absorption coefficient of the deeper tissue (i.e. brain), which causes a reduction of detected photons that travelled for longer time inside the brain. Thus, the sensitivity to superficial layers (ECT and CSF) increases, together with the dispersion of retrieved tHb values: as a matter of fact, the parameter which mostly influences the dimension of shadow areas is µ a ECT .

Conclusions
In this work we studied the sensitivity of TD-NIRS in assessing neonatal brain hemodynamic parameters, when homogenous models of analysis are employed. We considered differences in head morphology between preterm and term neonates, simulating light propagation in realistic 3D head structures of 29 and 44 weeks PMA newborns, using MC method. Real life measurements were modelled, varying brain absorption coefficients according to three common conditions in neonates: ischemia (scenario 1), hyperoxygenation (scenario 2) and hypoventilation (scenario 3). Simulated DTOF were analyzed with two analytical homogeneous models for photon migration based on Diffusion Equation: semi-infinite and spherical. Effects of superficial layers and head curvature on retrieved optical and hemodynamic parameters were analyzed, highlighting differences between the two meshes and the two models adopted for analysis. The results presented in Sec. 3 indicate an underestimation of brain total hemoglobin concentration, probably related to the influence of superficial layers, which was more evident in the case of 29 weeks PMA mesh and spherical homogeneous model. Indeed, the CSF strongly influenced the retrieved optical and hemodynamic parameters, and due to its low scattering and absorption coefficients (two order of magnitude lower than ECT and brain), it causes a reduction of estimated brain parameters. On the contrary the ECT is characterized by higher optical properties than CSF and brain, and its presence causes an increase of the retrieved brain optical properties. Thus, the effect of ECT partially balances the one of CSF: the higher the thickness of ECT, the stronger is the balance, and the lower is the underestimation. Conversely, the estimated S t O 2 exhibits good accuracy to brain values, and low values of MAPE StO2 . These results were in accordance with a previous work [28], where sensitivity of multi-distance FD-NIRS on neonates hemodynamic shows similar behavior. Head curvature influenced hemodynamic parameters: as consequence, we observed that spherical homogeneous model of analysis allows to better retrieve brain saturation variations.
The presented results empower the ability of TD-NIRS technique in monitoring neonatal brain hemodynamics but highlight the need for more complex methods of analysis which allow to reduce the influence of superficial layers and head curvatures: a possibility could be the use of MC or other numerical fitting procedures that exploits generalized atlases of three dimensional head structures, or the use of layered model of analysis.