1. Introduction
Functional near-infrared spectroscopy is a method to investigate brain activity non-invasively (Berivanlou et al., 2014; Ferrari & Quaresima, 2012; Hemmati et al., 2012; Rahimpour et al., 2020; Scholkmann et al., 2014). This low-cost imaging method has many applications in the field of neuroscience, including infant cerebral hemodynamic monitoring; classification of chronic diseases; stress level, and mental workload assessment, and intelligence quotient (IQ) estimation (Dadgostar et al., 2018; Firooz & Setarehdan, 2019; Hakimi & Setarehdan, 2018; Jahani et al., 2015; Mirbagheri et al., 2020b, 2020a; Rahimpour et al., 2017; Rahimpour et al., 2018). In this imaging method, a light source and detector are installed on the head; consequently, the re-emission of light from the human skin contains optical information from the human body (Scholkmann et al., 2014). Three types of spectroscopy methods from the brain tissue exist, such as the continuous wave method, time domain, and frequency domain technique (Ferrari & Quaresima, 2012). The first one is simpler and more portable than the other two and can wirelessly send hemodynamic information (Chiarelli et al., 2017; Ferrari & Quaresima, 2012; von Lühmann et al., 2017). Studies and results in this work are focused on continuous wave-functional near-infrared spectroscopy (CW-fNIRS). 
According to Figure 1 brain, activities are associated with the generation of action potentials; consequently, the amount of oxygen supply to that part increases. It is called the hemodynamic response. A single or multiple channels fNIRS can detect these hemodynamic changes.

The single-channel fNIRS helps measure hemodynamic changes based on detected optical modulation. An optical channel is created by setting a light source with a specified distance from the optical receiver. Hemodynamic changes modulate the changes in light received by the detector. It is essential to study the photon propagation profile to obtain the brain tissue’s sensitivity to the cerebral hemodynamic function. The question is how long is the depth of penetration in one single channel fNIRS? Another specific issue that arises is how much the reflectance of the detector is sensitive to depth. The previous study indicates that depth sensitivity decreases exponentially, depending on the source-detector separation (Mirbagheri et al., 2020b; Strangman et al., 2013).
Analytical methods have already been developed to study light emission within simple geometry, such as slab medium (Arridge et al., 1993; Arridge et al., 1992; Carraresi et al., 2001; Cui & Ostrander, 1992; Haselgrove et al., 1992; Patterson et al., 1989; Schweiger et al., 1993). Numerical methods also are used in complex brain models to study light diffusion in tissues (Mansouri et al., 2010; Strangman et al., 2013). The analytical methods take less computation time than statistical approaches, especially when it comes to investigating the effect of several fNIRS channels on depth sensitivity. It is also possible to evaluate the performance of high-density (HD) source and detector design topology on several hemodynamic reconstructions using an analytical approach (Borjkhani & Setarehdan, 2020). In this paper, using a perturbative diffusion equation, the light emission profile for a channel is studied. The photon beam propagation path into the tissue is investigated for both the perturbation theory on diffusion equation (DE) (Carraresi et al., 2001) and the finite element method (FEM). The detector's sensitivity to hemodynamic changes can be determined with the knowledge of the pathway of photons propagation. Both methods improve our understanding of the emission of photons into the human brain. The benefits and drawbacks of each approach will be discussed in this article.
The next section will describe the theory and mathematics governing the model. All the equations are solved in a medium like the human brain for both analytical and FEM methods. In the third section, the photon’s contrast in different X, Y-plane, and at different depths are simulated. The final part will discuss the quantitative comparison between these two methods and further address their benefits and limitations.

2. Material and Methods
Analytical solution
To study the transmission and reflectance of the light between pairs of source and detector concerning the configuration of the source and detectors on the surface of the medium and also to find out the three-dimensional distribution of photon inside it, an appropriate model of photon transport need to be used (Sassaroli et al., 2006). The models developed for this task are based on the radiative transfer theory. The derivatives of the radiative transfer equation (RTE) are stochastic or deterministic. No analytical answer exists to solve this equation. Therefore, simpler models of this equation are extracted. With some assumptions and simplification, DE is derived from RTE (Martelli et al., 2010). DE is a practical mode, and an analytical solution exists for it. The DE solution should be applied to inhomogeneous media similar to brain tissue properties. The photons’ intensity that undergoes many scattering events and is detected by the detector is called reflectance. Reflectance can be obtained by solving the DE Slab geometry. The optical properties of the Slab are considered: μa=0.017 m-1 and μs=1.2 mm-1 and and S=40 mm refractive index nr=1.4. 
In an environment modeled with optical properties identical to the human brain (Shown in Figure 2), an optical source in the location (x0, y0, z0) and a detector in position (x3, y3, z3) are placed at a given distance from the source. By placing an inclusion inside the mediums in (x2, y2, z2), the reflectance concerning the inclusion is calculated. The perturbed reflectance due to inclusion inside the medium can be obtained by (Carraresi et al., 2001):


Where R0(ρ) is the reflectance inside homogeneous media, δRa(ρ) and δRD(ρ) are the absorption and scattering effects of the inhomogeneous sphere, respectively.

Where z3m=-2ms-4mze-zs , z4m=-2ms-4mze+2ze+zs and “s” is the thickness of the Slab and . μa is the distance between source and detector.  is the absorption coefficient, and D is the diffusion coefficient. In steady-state conditions, the perturbation for reflectance due to the inhomogeneity absorption for the channel is:

Where ρ , ze=2AD, D is the diffusion coefficient, and A is dependent on the refractive index:


And also, h and w functions are given (Carraresi et al., 2001):


where  is the effective attenuation coefficient.
The perturbed reflectance in the channel is obtained by sweeping the inhomogeneity in three-dimensional space. The perturbed reflectance due to scattering is considered to be constant for simplicity. The spatial probability distribution profile of photons penetrating tissue at a source spot, scattering into the tissue, and released at an appropriate detector position, represents the spatial sensitivity. 

Finite element method (FEM)
One of the numerical approaches for the solution of DE in the slab medium is FEM. The optical properties of the medium, in this case, is considered the same as section 2.1. The dimension of the Slab is . The source and detector are placed in (40 mm, 60 mm, 0) and (80 mm,60 mm, 0), respectively. The distance between them is 40 mm same as in the previous section. The inclusion with only a 5%-6% change in optical properties of the background medium is inserted under the source and detector’s surface and inside the medium. Figure 3 shows the geometry and mesh for the FEM solution.

The FEM is applied to the steady-state diffusion equation:

Where  is the fluence rate, which is a scalar intensity (in units of W/m2), represents the power of light radiating radially per area at position  and  indicates the power per area from source element at position  (Wheelock et al., 2019). The geometry in Figure. 3 A is discretized into voxels in Figure 3 B. Then Equation 8 is solved for each voxel to calculate the solution of FEM method. The size of the geometry elements controls the size of voxels to reduce the computation time of FEM. Since the size of the source, detector, and inclusion are small, then the voxels beside them are small and become wider at a distance from them. The depth sensitivity for the FEM method is calculated using the defined sensitivity equation (Equation 9):

Where , which is the integral of total power reflected in the detector area in existing of the perturbation in position  with μ absorption coefficient. And  is the integral of the total power reflected in the detector area inhomogeneous medium without any inclusion. Vi refers to the volume of the perturbation. It is better to note that the shape and size of the source and detector are considered a circle with a 1.5 mm radius, and the inclusion is a cubic element with a volume of 2×2×2 mm.
The next section described the simulation results of the spatial sensitivity for both analytical and FEM approaches.

3. Results
For an analytical case, the inclusion has been moved under the source and detector surface to measure the relative perturbation (δRa) ⁄ R0 (Contrast) versus Y-plane. The result of this simulation has been shown in Figure 3 for Y=2 mm, Y=3 mm, Y=4 mm, and Y=5 mm. The result of this simulation indicates that sensitivity is slowly decreasing, and on the other hand, the sensitivity below the source and detector is the highest. Concerning Figure 4 A, the dominant pathway of photons is a banana shape. According to Figure 4 B, Figure 4 C, and Figure 4 D, as the inclusion moves away from the Y=0, the sensitivity is decreasing.
The same procedure is repeated in Figure 5 for Z-plane to achieve an image for spatial sensitivity in depth. According to this simulation, the sensitivity in Figure 5 A is higher (in Z=10 mm) compared to Figure 5 B, Figure 5 C, and Figure 5 D. 

The depth sensitivity of both Analytical and FEM is calculated for the same Slab geometry with the same optical properties. The result of the comparison of sensitivity between these methods is shown in Figure 6.

Regarding this figure, the depth sensitivity reduced gradually in-depth, and the shape of sensitivity for both approaches is almost identical. The depth sensitivity depends on the geometry and size of the source and detector and perturbation size. The volume of perturbation and its associated absorption coefficient directly impact the depth sensitivity. The radius of the source and detector for FEM simulation is 1.5 mm. In the FEM approach, the sensitivity is calculated according to Equation 9. The depth sensitivity of analytical and FEM method for source and detector distance of 40 mm is also compared in Figure 7 (in this comparison, the amount of perturbation for both approaches is δμa×Vi=1mm2).

Finally, Table 1 summarizes the results of the comparison between the two methods. Both of them take advantage of the simple diffusion approximation of RTE. However, a considerable difference is observed in computation times. MATLAB and COMSOL multiphysics are the simulation environments for analytical and FEM simulations, respectively. The software runs on a laptop computer with an Intel Core-i5 processor and 4GB RAM.

4. Discussion
In this paper, perturbation theory is applied for the analytical solution of the DE in slab media based on the method in (Carraresi et al., 2001). The solution of FEM is applied to the diffusion equation and the results of depth sensitivity are compared to both methods. In the FEM solution, the volume and shape of the source, detector, and perturbation can be changed. Even it is possible to calculate the solution for complex geometries (Wheelock et al., 2019).
Using both simulation approaches, one can obtain the trajectory of photon propagation received by the optical detector. Although the FEM can be applied to sophisticated geometries like the human brain and several options to study the effect of source and detector size on simulation, the computation time of FEM is four orders of magnitude higher than analytical simulation. One of the significant contributions of this investigation is to compare the computation time of each method. 
It is assumed that the human brain is only one layer, while realistic results can be obtained by taking several layers (since, in reality, the head model has several layers). The perturbation theory’s accuracy in the article (Carraresi et al., 2001) is compared with the Monte Carlo (MC) analysis results. There is a good agreement between the solution of perturbative DE and Monte Carlo (MC) simulation.
We can examine the depth of penetration and the photon propagation’s shape by regulating the distance between the source and the light detector. This study would help fNIRS researchers design a better measurement setup in the instrumentation part and reconstruct the hemodynamic response concerning spatial sensitivity.
The quantitative comparison concludes that the elapsed simulation time in FEM is higher than analytical. When the number of simulated channels exceeds more than 100 channels, the simulation time changes from a few minutes to a few hours. It is recommended to use the analytical method when the geometry of the head model is simple. For complex geometries, the FEM will be a suitable option despite the computation time challenge. The calculated depth sensitivity for both methods is almost identical. It is expected because, in both procedures, the head model’s optical properties and geometry are similar, and diffusion approximation has been used for both techniques.
The introduced strategies can be applied to study light-tissue interaction and simulation of synthetic fNIRS channels (Bonomini et al., 2015; Borjkhani & Setarehdan, 2020; Torricelli et al., 2005) based on a linear regression model, to statistically estimate the hemodynamic activations in fNIRS data sets. The main concern guiding the algorithm development was the minimization of assumptions and approximations made on the data set for the application of statistical tests. Further, we propose a K-means method to cluster fNIRS data (i.e. channels. The analytical and FEM results will be compared with laboratory results in future work by implementing a near-infrared spectroscopy system and liquid phantom.

Ethical Considerations
Compliance with ethical guidelines
There were no ethical considerations to be considered in this research.

Funding
This research did not receive any grant from funding agencies in the public, commercial, or non-profit sectors.

Authors' contributions
The both authors equally contributed to preparing this article.

Conflict of interest
The authors declared no conflict of interest.

Acknowledgments
The authors would like to thank Alessandro Torricelli, professor of the Department of Physic at the University of Politecnico di Milano for his valuable suggestions, many helpful hints, and scientific supports on developing this manuscript.


References
Arridge, S. R., Schweiger, M., Hiraoka, M., & Delpy, D. T. (1993). A finite element approach for modeling photon transport in tissue. Medical Physics, 20(2 Pt 1), 299–309. [DOI:10.1118/1.597069] [PMID]
Arridge, S. R., Cope, M., & Delpy, D. T. (1992). The theoretical basis for the determination of optical pathlengths in tissue: Temporal and frequency analysis. Physics in Medicine and Biology, 37(7), 1531–1560. [DOI:10.1088/0031-9155/37/7/005] [PMID]
Hemmati Berivanlou, N., Setarehdan, S. K., & Ahmadi Noubari, H. (2014). Evoked hemodynamic response estimation using ensemble empirical mode decomposition based adaptive algorithm applied to dual channel functional near infrared spectroscopy (fNIRS). Journal of Neuroscience Methods, 224, 13–25. [DOI:10.1016/j.jneumeth.2013.12.007] [PMID]
Bonomini, V., Zucchelli, L., Re, R., Ieva, F., Spinelli, L., & Contini, D., et al. (2015). Linear regression models and k-means clustering for statistical analysis of fNIRS data. Biomedical Optics Express, 6(2), 615–630. [DOI:10.1364/BOE.6.000615] [PMID] [PMCID]
Borjkhani, H., & Setarehdan, S. K. (2020). Performance assessment of high-density diffuse optical topography regarding source-detector array topology. Plos One, 15(3), e0230206.[DOI:10.1371/journal.pone.0230206] [PMID] [PMCID]
Carraresi, S., Shatir, T. S., Martelli, F., & Zaccanti, G. (2001). Accuracy of a perturbation model to predict the effect of scattering and absorbing inhomogeneities on photon migration. Applied Optics, 40(25), 4622–4632. [DOI:10.1364/AO.40.004622] [PMID]
Chiarelli, A. M., Libertino, S., Zappasodi, F., Mazzillo, M., Pompeo, F. D., & Merla, A., et al. (2017). Characterization of a fiber-less, multichannel optical probe for continuous wave functional near-infrared spectroscopy based on silicon photomultipliers detectors: In-vivo assessment of primary sensorimotor response. Neurophotonics, 4(3), 035002. [DOI:10.1117/1.NPh.4.3.035002] [PMID] [PMCID]
Cui, W. J., & Ostrander, L. E. (1992). The relationship of surface reflectance measurements to optical properties of layered biological media. IEEE Transactions on Bio-Medical Engineering, 39(2), 194–201. [DOI:10.1109/10.121651] [PMID]
Dadgostar, M., Setarehdan, S. K., Shahzadi, S., & Akin, A. (2018). Classification of schizophrenia using SVM via fNIRS. Biomedical Engineering, 30(02), 1850008. [DOI:10.4015/S1016237218500084]
Ferrari, M., & Quaresima, V. (2012). A brief review on the history of human functional near-infrared spectroscopy (fNIRS) development and fields of application. NeuroImage, 63(2), 921–935. [DOI:10.1016/j.neuroimage.2012.03.049] [PMID]
Firooz, S., & Setarehdan, S. K. (2019). IQ estimation by means of EEG-fNIRS recordings during a logical-mathematical intelligence test. Computers in Biology and Medicine, 110, 218–226.[DOI:10.1016/j.compbiomed.2019.05.017] [PMID]
Hakimi, N., & Setarehdan, S. K. (2018). Stress assessment by means of heart rate derived from functional near-infrared spectroscopy. Journal of Biomedical Optics, 23(11), 1–12.[DOI:10.1117/1.JBO.23.11.115001] [PMID]
Haselgrove, J. C., Schotland, J. C., & Leigh, J. S. (1992). Long-time behavior of photon diffusion in an absorbing medium: Application to time-resolved spectroscopy. Applied Optics, 31(15), 2678–2683. [DOI:10.1364/AO.31.002678] [PMID]
Hemmati, N., Setarehdan, S. K., & Noubari, H. A. (2012). Multi-channel near-infrared spectroscopy (NIRS) system for noninvasive monitoring of brain activity. Paper presented at: Proceedings of 2012 IEEE-EMBS International Conference on Biomedical and Health Informatics, Hong Kong, China, 7 January 2012. [DOI:10.1109/BHI.2012.6211548]
Jahani, S., Berivanlou, N. H., Rahimpour, A., & Setarehdan, S. K. (2015). Attention level quantification during a modified stroop color word experiment: An fNIRS based study. Paper presented at: 2015 22nd Iranian Conference on Biomedical Engineering (ICBME), Tehran, Iran, 27 November 2015. [DOI:10.1109/ICBME.2015.7404124]
Mansouri, C., L'huillier, J. P., Kashou, N. H., & Humeau, A. (2010). Depth sensitivity analysis of functional near-infrared spectroscopy measurement using three-dimensional Monte Carlo modelling-based magnetic resonance imaging. Lasers in Medical Science, 25(3), 431–438. [DOI:10.1007/s10103-010-0754-4] [PMID]
Martelli, F., Del Bianco, S., Ismaelli, A., & Zaccanti, G. (2010). Light propagation through biological tissue. Washington: SPIE press. [DOI:10.1117/3.824746]
Mirbagheri, M., Hakimi, N., Ebrahimzadeh, E., & Setarehdan, S. K. (2020). Quality analysis of heart rate derived from functional near-infrared spectroscopy in stress assessment. Informatics in Medicine Unlocked, 18, 100286. [DOI:10.1016/j.imu.2019.100286]
Mirbagheri, M., Hakimi, N., Ebrahimzadeh, E., & Setarehdan, S. K. (2020). Simulation and in vivo investigation of light-emitting diode, near infrared Gaussian beam profiles. Journal of Near Infrared Spectroscopy, 28(1), 37-50. [DOI:10.1177/0967033519884209]
Patterson, M. S., Chance, B., & Wilson, B. C. (1989). Time resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties. Applied Optics, 28(12), 2331–2336. [DOI:10.1364/AO.28.002331] [PMID]
Rahimpour, A., Dadashi, A., Soltanian-Zadeh, H., & Setarehdan, S. K. (2017). Classification of fNIRS based brain hemodynamic response to mental arithmetic tasks. Paper presented at: 2017 3rd International Conference on Pattern Recognition and Image Analysis (IPRIA), Shahrekord, Iran, 20 April 2017. [DOI:10.1109/PRIA.2017.7983029]
Rahimpour, A., Noubari, H. A., & Kazemian, M. (2018). A case-study of NIRS application for infant cerebral hemodynamic monitoring: A report of data analysis for feature extraction and infant classification into healthy and unhealthy. Informatics in Medicine Unlocked, 11, 44-50. [DOI:10.1016/j.imu.2018.04.001]
Rahimpour, A., Pollonini, L., Comstock, D., Balasubramaniam, R., & Bortfeld, H. (2020). Tracking differential activation of primary and supplementary motor cortex across timing tasks: An fNIRS validation study. Journal of Neuroscience Methods, 341, 108790. [DOI:10.1016/j.jneumeth.2020.108790] [PMID] [PMCID]
Sassaroli, A., Martelli, F., & Fantini, S. (2006). Perturbation theory for the diffusion equation by use of the moments of the generalized temporal point-spread function. II. Continuous-wave results. Journal of the Optical Society of America, 23(9), 2119–2131. [DOI:10.1364/JOSAA.23.002119] [PMID]
Scholkmann, F., Kleiser, S., Metz, A. J., Zimmermann, R., Mata Pavia, J., & Wolf, U., et al. (2014). A review on continuous wave functional near-infrared spectroscopy and imaging instrumentation and methodology. NeuroImage, 85 (Pt 1), 6–27. [DOI:10.1016/j.neuroimage.2013.05.004] [PMID]
Schweiger, M., Arridge, S. R., & Delpy, D. T. (1993). Application of the finite-element method for the forward and inverse models in optical tomography. Journal of Mathematical Imaging and Vision, 3(3), 263-283. [DOI:10.1007/BF01248356]
Strangman, G. E., Li, Z., & Zhang, Q. (2013). Depth sensitivity and source-detector separations for near infrared spectroscopy based on the Colin27 brain template. Plos One, 8(8), e66319.[DOI:10.1371/journal.pone.0066319] [PMID] [PMCID]
Torricelli, A., Pifferi, A., Spinelli, L., Cubeddu, R., Martelli, F., & Del Bianco, S., et al. (2005). Time-resolved reflectance at null source-detector separation: Improving contrast and resolution in diffuse optical imaging. Physical Review Letters, 95(7), 078101. [DOI:10.1103/PhysRevLett.95.078101] [PMID]
von Luhmann, A., Wabnitz, H., Sander, T., & Muller, K. R. (2017). M3BA: A mobile, modular, multimodal biosignal acquisition architecture for miniaturized EEG-NIRS-based hybrid BCI and monitoring. IEEE Transactions on Bio-Medical Engineering, 64(6), 1199–1210. [DOI:10.1109/TBME.2016.2594127] [PMID]
Wheelock, M. D., Culver, J. P., & Eggebrecht, A. T. (2019). High-density diffuse optical tomography for imaging human brain function. The Review of Scientific Instruments, 90(5), 051101.[DOI:10.1063/1.5086809] [PMID] [PMCID]