EDP Sciences logo
Submit your research paper
Search
Menu
All issues Volume 22 / No 1 (2026) J. Eur. Opt. Society-Rapid Publ., 22 1 (2026) 48 Full HTML
Open Access
Issue
J. Eur. Opt. Society-Rapid Publ.
Volume 22, Number 1, 2026
Article Number 48
Number of page(s) 11
DOI https://doi.org/10.1051/jeos/2026044
Published online 03 June 2026

© The Author(s), published by EDP Sciences, 2026

Licence Creative CommonsThis is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Near–infrared (NIR) optical tissue imaging is widely used to probe subsurface structures and optical contrast in biological media [1–3], owing to the relatively low absorption and moderate scattering of tissue in this spectral window [4–6]. Accurate interpretation of measured signals depends critically on reliable forward models of photon transport through layered tissue structures [7–10]. Photon propagation in the human head is governed by a complex interplay between tissue optical properties, anatomical geometry, and refractive index mismatches at tissue boundaries [11, 12]. While the influence of bulk tissue composition and macroscopic anatomical variability has been extensively studied [13], the role of fine–scale surface microgeometry in shaping photon transport remains incompletely understood.

The outermost boundary of the tissue domain, corresponding to the scalp surface, is characterized by pronounced microstructural features, most notably skin pores [14, 15]. Such structures introduce localized geometric perturbations and refractive index discontinuities at medium interfaces [9, 16], potentially modifying photon coupling conditions and the angular distribution of injected light. In practice, surface microstructural features such as skin pores are almost universally ignored in the design and modeling of NIR optical imaging systems, where tissue boundaries are treated as smooth planar or gently curved interfaces. This design choice is commonly motivated by the expectation that scattering within superficial tissue layers rapidly randomizes photon trajectories [17–20], effectively suppressing sensitivity to small–scale surface irregularities. From a ray–based propagation perspective, pore–scale structures are similarly assumed to act as weak perturbations relative to the optical source footprint and emission divergence, and are therefore not expected to produce sustained alterations in photon trajectories or coupling efficiency. Accordingly, smooth boundary assumptions are generally considered adequate for predicting optical coupling and power transmission at the tissue surface. However, the validity of these assumptions has not been systematically examined across different photon transport regimes or for physiologically realistic pore geometries.

The present study examines the influence of skin pore microgeometry on NIR photon transport using anatomically realistic, MRI–segmented cadaveric head models and complementary simulation frameworks. Photon transport in biological tissue exhibits a transition from an initial quasi-ballistic regime near sources and interfaces to a diffusive regime as multiple scattering events accumulate, a behavior well established in radiative transfer theory and Monte Carlo (MC) modeling of turbid media [21, 22]. Accordingly, photon propagation across the source–scalp interface and subsequent transport within tissue were modeled using a cascaded simulation framework. Ballistic photon propagation simulations with explicit representation of physiologically relevant pore geometries were first used to resolve surface-scale optical interactions and determine the resulting photon states. These photon states, characterized by their position, direction, and weight, were then used to initialize MC-based diffusive photon transport simulations within anatomically realistic head models. By linking surface-scale geometric optics with subsequent scattering-dominated transport in deeper tissues, this approach enables systematic assessment of how pore-scale surface features influence NIR photon transport and clarifies the conditions under which such features may be neglected or require explicit representation in optical tissue imaging models.

2. Material and methods

The study was conducted using a cascaded simulation workflow, consisting of subsections: Magnetic Resonance Imaging–Derived Cadaveric Head Anatomy (Sect. 2.1), Ballistic Photon Propagation Simulations (Sect. 2.2), Diffusive Photon Transport Simulations (Sect. 2.3), Data Processing and Quantitative Analysis (Sect. 2.4), as detailed in the following subsections.

2.1 Magnetic resonance imaging–derived cadaveric head anatomy

Anatomically realistic cadaveric head geometries used in this study were obtained from a previously published MRI dataset originally acquired for the study DrSVision: A Machine Learning Tool for Cortical Region–Specific fNIRS Calibration Based on Cadaveric Head MRI [23]. The dataset comprises eight adult human cadaveric heads (five male and three female; age range: 67–97 years), selected to exclude cranial trauma, prior surgery, or visible pathology. Cadaveric material was obtained through a licensed anatomical supplier with informed consent provided during donors’ lifetimes, and all procedures were conducted in accordance with ethical and legal guidelines, as detailed in the original study.

High–resolution T1 and T2–weighted 3 Tesla MRI scans enabled segmentation of scalp, skull, cerebrospinal fluid (CSF), and brain tissues relevant to optical modeling. In the present work, the previously segmented geometries and tissue layer thicknesses were reused directly for diffusive photon transport simulations. No new MRI acquisitions or anatomical segmentation procedures were performed.

2.2 Ballistic photon propagation simulations

Ballistic photon propagation at the skin surface was modeled using the COMSOL Multiphysics (v6.2, COMSOL AB, Stockholm, Sweden) [24] platform, employing the geometric optics interface in 2-dimensional (2D) space with explicit representation of physiologically relevant pore geometries with surface roughness to resolve surface-scale optical interactions and determine the resulting photon states. These simulations were performed using representative surface geometries, as the underlying anatomical layers are not involved at this stage. The geometric optics interface operates in the high–frequency limit where the vacuum wavelength is negligible compared to the characteristic dimensions of the geometry. In this regime, Maxwell’s equations are reduced to the Eikonal equation [25], |∇S|2=n2(r), Mathematical equation: $$ {\left|\nabla S\right|}^2={n}^2(r), $$(1)

where S represents the optical phase and n(r) is the local refractive index. This formulation treats photon propagation as a set of discrete rays whose trajectories are determined by integrating the Hamiltonian form of the ray equations, dqdt=∂ω∂k, Mathematical equation: $$ \frac{dq}{dt}=\frac{\partial \omega }{\partial k}, $$(2) dkdt=−∂ω∂q, Mathematical equation: $$ \frac{dk}{dt}=-\frac{\partial \omega }{\partial q}, $$(3)

where q, k, and t denote the position vector, wave vector, and time of flight, respectively. The term ω signifies the angular frequency of the radiation, which acts as the Hamiltonian of the system and remains constant along the ray path. To account for energy transport and the formation of caustics within the absorbing media, the local intensity I was computed by tracking the evolution of the principal radii of curvature (r1, r2) of the wavefront. For a medium characterized by a complex refractive index, the intensity I along the ray path is governed by, I=I0|r1,0r2,0r1r2|exp(−∫‍μads), Mathematical equation: $$ I = I_0 \left| \frac{r_{1,0} r_{2,0}}{r_1 r_2} \right| \exp \left( -\int \mu_a ds \right), $$(4)

where I0, r1,0, and r2,0 are the initial intensity and principal radii of curvature, respectively. The exponential term accounts for attenuation via the Beer–Lambert law, where the absorption coefficient μa is integrated along the incremental path length ds. This curvature–based approach ensures a mesh–independent evaluation of the energy flux, even in regions of high convergence where traditional ray–counting methods fail.

A scalp tissue layer was represented using a rectangular geometry with a thickness of 5.64 mm, corresponding to that of the mean cadaveric head (Supplementary Materials, Table S1, Section A: Magnetic Resonance Imaging–Derived Cadaveric Head Anatomy). Skin pore geometries were explicitly represented using truncated semi–elliptical structures [26, 27] embedded within the top surface of the scalp layer, defined by a fixed diameter–to–depth ratio of 10:1. This semi-ellipsoidal approximation was selected as the optimal first-order geometric model to isolate the influence of surface curvature and aperture-driven diffraction on incident light, as it accurately captures the high aspect ratios and characteristic morphology observed in standard dermatological histology [28]. Furthermore, this geometry accounts for the critical refractive index mismatch at the air–tissue interface, which dictates the initial phase-space distribution of ballistic photons before they enter the diffusive regime [5]. Multiple pore diameters were evaluated across separate simulation runs, with pore diameters of 500, 250, 125, 62.5, 31.25, and 15.625 μm selected to span a physiologically realistic and literature–reported range [29, 30] and to assess scale–dependent ray interactions. To approximate realistic surface irregularities, pore boundary roughness was introduced by perturbing the idealized geometry using Gaussian noise. This was implemented in COMSOL Multiphysics via a stochastic function with normal distribution (mean 0, standard deviation 0.1 μm) and a fixed random seed, which was applied to the parametric representation of the truncated semi–elliptical boundaries. The computational domain was discretized using a free triangular mesh, with minimum and maximum element sizes of 2 μm and 1 mm, respectively. A maximum element growth rate of 1.1, a curvature factor of 0.2, and a narrow–region resolution of 1 were employed to ensure accurate ray–interface interactions while maintaining computational efficiency. Within the geometric optics interface, the intensity computation setting was configured to compute intensity and power in graded media, to correctly account for continuous refractive index variations in heterogeneous tissue [25]. The same wavelength–dependent optical properties used in the MC simulations (Table 1) were also assigned to the scalp tissue layers in the COMSOL Multiphysics model to ensure methodological consistency since the geometric optics interface requires specification of a complex refractive index, the real part n was directly assigned from Table 1, while the imaginary part κ was computed from the absorption coefficient μa according to, κ=μaλ4π, Mathematical equation: $$ \kappa =\frac{\mu_a\lambda }{4\pi }, $$(5)

Table 1

Optical properties of each tissue layer used in the MCX simulations, including absorption coefficient, scattering coefficient, anisotropy factor, and refractive index, evaluated at a wavelength of 735 nm [31].

where λ denotes the operating wavelength 735 nm. A release boundary condition was assigned to the right vertical edge of a 5 × 5 mm square source domain to emulate an unpolarized optical source with an irradiance of 1 W/m, emitting rays in the positive x–direction with the central emission axis aligned horizontally. A 70° half-angle was specified to replicate the divergence characteristics of commonly used NIR optical tissue imaging systems [32]. The source was oriented perpendicular to the scalp layer, vertically (radially) centered, and placed in direct contact with its surface. Ray directions were defined using the conical emission option, with the number of rays per release and the number of rays in wave vector space set to 50 and 100, respectively, resulting in an effective emission of 5000 primary rays per simulation. Unlike stochastic MC methods, geometric optics ray tracing is deterministic; therefore, simulations employing 103–104 rays provide robust solutions, consistent with COMSOL Multiphysics Ray Optics Module guidelines [25]. A wall boundary conditions configured with the disappear option were assigned to the outer perimeters of the computational domain to ensure elimination of rays exiting the model, thereby reducing the number of rays tracked during propagation and improving computational efficiency without influencing internal ray–tissue interactions. Simulations were repeated across 6 distinct skin–pore diameters (500, 250, 125, 62.5, 31.25, and 15.625 μm). A time–dependent ray tracing study was conducted over a 0–500 ps interval with a time step of 10 ps. To evaluate positional sensitivity, a parametric sweep was applied to the vertical location of the skin pore geometry, spanning −1 to 0 mm in increments of 0.1 mm, where 0 mm corresponds to the vertical centerline of the optical source. Simulations for positive vertical positions were not performed, as the computational domain is vertically symmetric about the source centerline; thus, such configurations would yield mirrored photon trajectories and redundant results. For qualitative assessment, 2D ray trajectory maps were generated using a rainbow color scale to visualize ray propagation paths and interface interactions. Simulations were performed on a high performance workstation, with each run requiring approximately 4.36 s. In total, 66 simulation runs were executed, corresponding to 11 vertical pore positions, and 6 pore diameters. Upon completion of the ballistic photon propagation simulations, photon positions, directions, and weights were recorded at the point of exit from the skin pore and exported as .txt files. To transition from the 2D ballistic photon propagation simulations to the 3-dimensional (3D) diffusive photon transport simulations, the 2D photon positions and directions were revolved around the x-axis, effectively distributing photons uniformly around the axis, x3D=y2D·cosθ, Mathematical equation: $$ {x}_{3\mathrm{D}}={y}_{2\mathrm{D}}\cdotp \cos \theta, $$(6) y3D=y2D·sinθ, Mathematical equation: $$ {y}_{3\mathrm{D}}={y}_{2\mathrm{D}}\cdotp \sin \theta, $$(7) z3D=x2D, Mathematical equation: $$ {z}_{3\mathrm{D}}={x}_{2\mathrm{D}}, $$(8)

where θ is the rotation angle around the 2D x-axis and is sampled uniformly from [0, 2π) to distribute photons evenly around the axis, avoiding duplication at θ = 0 and θ = 2π. This procedure generates fully 3D photon states, which were subsequently used as input for MC-based diffusive photon transport simulations.

2.3 Diffusive photon transport simulations

Diffusive photon transport in the eight individual cadaveric head models, as well as in the corresponding mean head model, was simulated using Monte Carlo eXtreme (MCX, v2025.10, Northeastern University, Boston, MA, USA), an open-source MC-based framework for modeling photon transport in 3D turbid media and scattering effects. These simulations were initialized immediately using the photon states obtained from the ballistic photon propagation simulations after photons exited the skin pores. MCX probabilistically tracks photon packets as they undergo scattering and absorption events within tissue [33]. As photons propagate through the domain, the photon packet weight is attenuated via absorption according to the Beer–Lambert law, while scattering events redirect the trajectories. The photon fluence, ϕ, within each voxel is defined as the sum of the weights of all photon packets traversing that voxel, representing the cumulative photon energy. Accordingly, the fluence distribution ϕ (x, y, z) at voxel coordinates (x, y, z) is computed as, ϕ(x,y,z)=∑i=1Nx,y,z‍wi, Mathematical equation: $$ \phi(x,y,z) = \sum_{i=1}^{N_{x,y,z}} w_i, $$(9)

where Nx,y,z denotes the number of photon packets traversing voxel (x, y, z) and wi is the weight of the i-th photon within that voxel. This discrete accumulation of photon weights forms the basis for quantifying optical energy distribution within the tissue domain. The computational domain consisted of a 30 × 70 × 64 voxel grid, with an isotropic voxel size of 1 mm3, providing a balance between spatial resolution and computational efficiency. The voxelated geometry represented the scalp, skull, CSF, and brain tissue layers. Layer thicknesses for non-brain tissues were computed as the mean across the eight cadaveric heads segmented from MRI data [23], with the remaining volume assigned to brain tissue to generate a representative mean-head geometry. Diffusive photon transport simulations were performed both on this mean cadaveric head model and on all individual cadaveric head geometries. For the mean cadaveric head model, simulations were conducted for all six skin pore diameters and 11 vertical pore positions, consistent with the configurations used in the ballistic photon propagation simulations. For the individual cadaveric heads, simulations were also conducted on the same six pore diameters but limited at two vertical skin-pore positions, 1 and 0.5 mm offset from the optical source center, to assess anatomical sensitivity. Each tissue layer was constructed using voxelated box structures, and the corresponding thickness values are provided in the Supplementary Materials (Table S1, Section A: Magnetic Resonance Imaging–Derived Cadaveric Head Anatomy). The optical properties of each tissue layer, including the absorption coefficient (μa), scattering coefficient (μs), anisotropy factor (g), and refractive index (n), were assigned based on values reported in the literature [31] (Table 1).

To emulate the optical emission profile of NIR optical tissue imaging systems, the photon emitter was modeled as a disc source with a diameter of 5 mm and a 70° half-angle, matching the source properties used in the ballistic photon propagation simulations and replicating the divergence characteristics of commonly used NIR optical tissue imaging systems [32]. The detector was modeled as a spherical active area with a diameter of 0.69 mm to represent the effective photosensitive area of typical functional near-infrared spectroscopy (fNIRS) photodetectors [34, 35]. The source and detector were positioned perpendicularly above the scalp layer, centered along the vertical axis, with a fixed 30 mm horizontal separation such that their midpoint aligned with the plane center. Each simulation launched 109 photon packets to minimize statistical fluctuations, producing high-fidelity fluence maps [36, 37]. Refractive index mismatches between tissue layers were incorporated via Fresnel boundary conditions. Simulations were conducted at a wavelength of 735 nm, with initial photon states assigned using .txt files exported from ballistic photon propagation simulations, and a power scaling factor set to unity. 88 simulation runs were performed on the same workstation used for the ballistic photon propagation simulations, with each simulation requiring approximately 480 s. Output data, including volumetric fluence maps and detected photon counts, were stored in .jnii and .jdat formats.

2.4 Data processing and quantitative analysis

Ballistic photon propagation simulation results were post–processed using custom scripts implemented in MATLAB (R2025a, The MathWorks, Inc., Natick, MA, USA) [38]. Photon positions, directions, and weights immediately after exiting the skin pores were extracted, and vertical distance from the optical source center versus photon exit angle scatter plots were generated, with each point color-coded according to the corresponding photon weight using a custom magenta–cyan blue colormap.

Diffusive photon transport simulation results stored in .jnii format were decoded, decompressed, and reshaped into 3D arrays with isotropic voxel dimensions of 1 × 1 × 1 mm3, using custom scripts executed in Python (3.12, Python Software Foundation, Wilmington, DE, USA) [39]. The resulting fluence volumes were visualized from top, side, front, and isometric viewpoints using a custom yellow–green–blue–dark blue colormap with logarithmic intensity scaling to accommodate the wide dynamic range of photon density.

Subsequently, source–detector based DS values were computed in % form by summing fluence within a circular region of interest (ROI) of 10 mm radius, centered at the midpoint between the source and detector locations, at successive 1 mm depth increments along the tissue depth axis perpendicular to the surface. A 10 mm radius ROI was selected to balance spatial specificity with adequate sampling volume, consistent with common practice in optical imaging literature. This scale has been shown to capture relevant optical or hemodynamic responses while minimizing contamination from surrounding regions [40, 41]. DS(z)=∑(x,y)∈R(z)‍ϕ(x,y,z)∑x,y,z‍ϕ(x,y,z)×100%, Mathematical equation: $$ DS(z) = \frac{\sum_{(xy) \in R(z)} \phi(xyz)}{\sum_{x,y,z} \phi(xyz)} \times 100\%, $$(10)

where R(z) denotes the set of voxels within the ROI at depth z. DS versus depth curves were plotted for the mean cadaveric head at each optical source–scalp layer separation. Additionally, detected photon data stored in .jdat format were decoded and used to plot histograms of photon count versus path length for all diffusive photon transport simulations. This post–processing workflow enabled high–resolution spatial and statistical assessment of DS and photon propagation characteristics.

3. Results

Ray–tracing results obtained from the ballistic photon propagation simulations are illustrated in Figure 1 for the mean cadaveric head with a 250, μm diameter skin pore at t = 10 ps. As the vertical position of the skin pore is varied, changes in ray–pore interactions and reflection patterns at the optical source–scalp interface are modeled, with Figures 1b–1c showing increasingly pronounced reflection features.

Thumbnail: Figure 1 Refer to the following caption and surrounding text. Figure 1

Example COMSOL ray–tracing plots at t = 10 ps for the mean cadaveric head with a 250, μm diameter skin pore, zoomed to the source–scalp interface. Skin pore at vertical (a) −1.5, (b) −0.8, and (c) −0.7 mm positions. The color bar represents normalized photon weights (0–1).

The influence of pore geometry on exit trajectories was assessed across 66 combinations of pore diameters and positions. Figures 2a–2c present the most pronounced photon weight reduction and dispersion behavior within the ballistic transport regime, specifically for the 500 μm diameter pore at vertical positions of (a) 1, (b) 0.5, and (c) 0 mm off-axis relative to the optical source center. As the skin pore approaches the optical source axis, the ejected photon angle shifts and photon weight is reduced. At a 1 mm off-axis position (Fig. 2a), there is no pore–photon interaction; therefore, the photon weights remain at their maximum, clustered, and initial photon angles are preserved. However, as the pore alignment shifts toward 0 mm (Figs. 2b–2c), both weight reduction and dispersion behavior are observed for the photon groups positioned at specific radial distance intervals.

Thumbnail: Figure 2 Refer to the following caption and surrounding text. Figure 2

Ballistic photon propagation simulation ejected photon distributions as a function of radial distance from the optical source center and angle (degrees off-axis relative to the optical source center) for a 500 μm diameter pore. Plots correspond to vertical pore positions (a) 1, (b) 0.5, and (c) 0 mm off-axis relative to the optical source center. The color bars represent normalized photon weights.

Spatial fluence color maps for the mean cadaveric head were generated for the same 66 combinations of pore diameters and positions. Additionally, spatial fluence distributions across the eight individual cadaveric heads were simulated for a 500 μm diameter pore at vertical positions of 1 and 0.5 mm off-axis to provide an anatomical sensitivity analysis. Figure 3 presents isometric views of the spatial fluence for the mean cadaveric head at a 30 mm source–detector separation. The 500 μm pore diameter and vertical positions 1, 0.5, and 0 mm (Figs. 3a–3c, respectively) were selected to highlight the maximum observable differences in fluence. Initial photon weights of Diffusive Photon Transport Simulations incorporate the final photon weights from the ballistic photon propagation simulations, thereby these spatial fluence color maps represent the localized photon energy distribution. The transition from off-axis to on-axis pore placement (Figs. 3a–3c) demonstrates a visually observable redistribution of energy density and reduction of photon weights in the vicinity of the source–scalp interface.

Thumbnail: Figure 3 Refer to the following caption and surrounding text. Figure 3

Isometric views of spatial fluence color maps for the mean cadaveric head at a 30 mm source–detector separation: plots correspond to 500 diameter pore at vertical pore positions (a) 1, (b) 0.5, and (c) 0 mm off-axis relative to the optical source center. Black circles and squares indicate the optical source and detector positions, respectively. The color bars represent photon fluence on a logarithmic scale to enhance visibility across a wide dynamic range.

DS profiles as a function of depth were generated for all six pore diameters across eleven vertical positions for the mean cadaveric head, alongside comparisons across the eight individual cadaveric models. Figure 4 presents a representative DS profile for the 500 μm diameter pore at a 30 mm source–detector separation. Analysis of the complete dataset (detailed in Supplementary Materials, Section C: Diffusive Photon Transport Simulations, Figures S149–S156) revealed that the eleven vertical pore positions exhibit virtually identical trends across the entire depth range. This high degree of uniformity indicates that while pore alignment modulates the total energy weight (Fig. 3), the relative spatial distribution of sensitivity is independent of localized surface microtopography. This trend is further corroborated by the comparisons across the eight individual cadaveric heads, which despite measurable variance in decay rates driven by anatomical differences, maintain a consistent baseline sensitivity profile. Together, these observations suggest that the fundamental DS is governed by the macro-anatomical layers of the scalp and skull, remaining robust against both localized geometric perturbations and subject-specific anatomical variations.

Thumbnail: Figure 4 Refer to the following caption and surrounding text. Figure 4

Representative DS profile for a 30 mm source–detector separation with a 500 μm diameter pore (positioned at 0.5 mm off-axis relative to the optical source center). The DS curves obtained across all individual cadaveric models and pore positions were found to be virtually identical, following the same spatial distribution. For visual clarity, this single representative profile is shown to illustrate the system’s sensitivity to underlying tissue volumes, which remains stable despite surface-level intensity fluctuations.

Detailed detected photon data for the mean cadaveric head model at a 30 mm source–detector separation are provided in the Supplementary Materials, Section C: Diffusive Photon Transport Simulations. Specifically, Table S2 summarizes the detected photon counts across all investigated combinations of pore diameters and vertical pore positions. Additionally, the complete set of 66 photon path length histograms is presented in Figures S157–S222. These distributions demonstrate that the temporal characteristics of the detected signal remain consistent with expected diffusive behavior for the specified source–detector separation, regardless of the surface pore geometry.

The statistical distribution of detected photon counts and mean photon weights across the eight cadaveric heads is summarized in Table 2. As shown in Figure 5, shifting the pore position from 1.0 mm to 0.5 mm off-axis resulted in a decrease in the median detected photon count. For both alignment positions, head #3 was identified as a consistent outlier, falling below the lower whisker of the cohort distribution. A parallel trend was observed for the mean photon weights (Fig. 6), where the median weight decreased as the pore was moved closer to the optical source center. Unlike the photon count data, the distribution of mean photon weights exhibited a larger interquartile range across the cohort, and no individual subjects were identified as outliers. These results indicate that while the vertical positioning of the skin pore modulates signal intensity, the underlying anatomical differences across the eight cadaveric heads introduce a more significant variance in these energy metrics than the localized pore alignment itself. This characterizes the variance in signal detection across the MRI-segmented models relative to localized pore-induced energy modulation, establishing that anatomical variability remains the primary determinant of total detected energy.

Table 2

Summary of detected photon statistics across eight cadaveric head models for two vertical pore alignments. The data highlights the inter-subject variability in total photon counts and mean photon weights relative to the systematic shifts induced by sub-millimeter changes in pore proximity.

Thumbnail: Figure 5 Refer to the following caption and surrounding text. Figure 5

Box plots for detected photon counts across eight cadaveric heads at two distinct vertical pore positions: 1, and 0.5 mm off-axis relative to the optical source center.
Thumbnail: Figure 6 Refer to the following caption and surrounding text. Figure 6

Box plots for detected mean photon weight across eight cadaveric heads at two distinct vertical pore positions: 1, and 0.5 mm off-axis relative to the optical source center.

4. Discussion

This study systematically evaluates the influence of skin pore microgeometry on NIR photon transport within anatomically realistic, MRI–segmented cadaveric head models. By employing a cascaded simulation approach, the present work effectively bridges the gap between pore-scale surface realism and macro-scale photon propagation. This sequential integration where ballistic photon propagation simulations through truncated semi–elliptical, rough-walled pores define the initial photon states for subsequent diffusive photon transport simulations provides a robust framework for assessing multi-scale optical interactions. Thereby, study clarifies the conditions under which microscopic surface features can be safely neglected in optical modeling, while simultaneously identifying specific regimes where anatomical and pore-scale variability contribute measurably to the detected signal.

Ballistic photon propagation simulations reveal that skin pore geometry modulates transmitted optical energy in a pore diameter and position-dependent manner, with the effect becoming more pronounced as pore diameter increases and as the vertical distance between the pore and the optical source center decreases. Within this regime, the modeled pore geometries act as localized perturbations at the optical source–scalp interface, redistributing the angular spectrum of incident rays. Examination of the angular trajectories (Figs. 2 and S1–S66) indicates that photons interacting with the rough pore structure tend to acquire larger propagation angles relative to the surface normal. Consequently, rather than remaining confined within a narrow forward-directed cone toward deeper tissues, these photons are redirected outward, increasing beam divergence and reducing the concentration of ballistic energy delivered to deeper layers. As pore diameter increases, this divergence effect influences a larger fraction of the incident beam profile, leading to progressively stronger spatial spreading of transmitted optical energy.

Building upon the initial phase of the study, the final photon states from the ballistic simulations characterized by both pore-induced angular redirection and photon weight reduction were used to initialize the subsequent diffusive photon transport simulations. The results demonstrate a clear bifurcation in how these perturbations propagate: while the high-scattering nature of the scalp and skull induces a rapid loss of directional memory within the first few millimeters [19, 20], the initial attenuation of photon weights persists as a dominant factor (Figs. 3 and S67–S148). Consequently, localized angular perturbations do not propagate into deeper tissue layers [42, 43] and produce no considerable changes in DS profile (Figs. 4 and S149–S156) or path length statistics (Figs. S157–S222). While the stochastic scattering regime effectively randomizes the angular information imparted by localized surface features [17, 18], it cannot recover the optical energy lost during the ballistic phase (Fig. 6). Although this pore-induced weight reduction is an order of magnitude smaller than the inter-subject anatomical variance, it represents a permanent loss of incident energy that scales systematically with pore geometry.

The identification of head #3 as a consistent outlier in detected photon counts (Figs. 5 and 6) provides a mechanistic example of this sensitivity; its exceptionally thick CSF layer (16.1 mm) likely facilitates increased lateral light piping, reducing the fraction of photons returning to the scalp [44]. This reinforces that while pore-scale features introduce systematic perturbations, large-scale anatomical variations remain the primary determinants of the absolute photon budget. From a system design perspective, these findings indicate that while explicit incorporation of skin pores is unlikely to influence the stochastic characteristics of diffusive transport, the resulting energy loss introduces a localized, systematic bias. This effect is most relevant for high-precision topographic mapping [45] or longitudinal fNIRS monitoring, where pore-induced attenuation could be misinterpreted as regional cortical variation or altered sensitivity if not accounted for during probe placement [46]. While this study focuses on reflective-mode optical systems, these findings suggest that surface-induced intensity fluctuations could also be relevant in transmissive-mode configurations, where photons traverse different tissue paths. In such geometries, the initial loss of photon energy at the source–scalp interface may represent a significant factor in the overall link budget. However, as the current study does not explicitly model transmissive geometries, the relative magnitude of these effects remains a hypothesis to be explored in future work.

Certain limitations of the present study should be acknowledged. First, the ballistic photon propagation simulations were performed in 2D using idealized elliptical pore geometries, with surface roughness introduced by applying Gaussian noise to the pore boundaries. While this approach captures key geometric perturbations, it cannot fully represent the 3D complexity of real skin microtopography. Furthermore, while a typical NIR optical source diameter (e.g., 5 mm) simultaneously illuminates a multi-pore array, the present work utilizes a single-pore configuration to isolate the fundamental mechanistic influence of pore geometry. While the aggregate effect of hundreds of pores across a wide source footprint would likely result in spatial averaging of the energy loss, the single-pore model established here serves as a necessary unit analysis to characterize the maximum potential perturbation at the source–scalp interface. Second, during the diffusive photon transport simulations, tissue optical properties were assumed to be homogeneous within each segmented anatomical layer and implemented within a graded-media framework. Consequently, dynamic physiological factors such as variations in skin hydration, local optical coupling conditions, and temporal changes in tissue optical properties were not considered. The inclusion of complex pilosebaceous units such as hair follicles or sebum was excluded from the present study to maintain a focused analysis on the primary source–scalp interface. Unlike surface skin pores, which create discrete geometric perturbations at the immediate point of entry where photons remain ballistic, hair follicles extend deep into the dermal and muscular layers. At these depths, the high scattering coefficient of the surrounding turbid media induces a rapid transition to the diffusive regime, effectively homogenizing the influence of such microstructures [42, 43, 47]. Consequently, the proposed cascaded framework prioritizes the skin pore as the dominant structural variable affecting initial photon angular distribution before the onset of bulk scattering. Despite these simplifications, the use of MRI-segmented cadaveric head models preserves anatomical realism at the macroscopic scale, while the cascaded implementation of ballistic photon propagation and diffusive photon transport regimes provides a mechanistically consistent framework for interpreting how pore-scale surface features influence subsequent photon transport. Finally, direct experimental validation of pore-scale optical coupling effects would require the fabrication of skin-mimicking phantoms with tightly controlled sub-millimeter surface microgeometry and well-defined optical properties. Developing such phantoms represents a substantial experimental challenge and therefore remains beyond the scope of the present computational study.

Future studies should focus on further refining the model by extending the ballistic photon propagation simulations to fully 3D geometries with multiple pore configurations and more realistic surface roughness. This would allow a more faithful representation of the stochastic nature of skin microtopography, provided that pore locations are spatially distributed to intercept the emitted photon bundle within the optical source half-angle of 70°, rather than representing arbitrary multiplicity without optical relevance. Including deep structures like hair follicles, would enable assessment of their complex geometry and contribution to photon transport. Integrating spatially varying optical properties and dynamic surface coupling conditions could further enable assessment of physiological variability under realistic measurement scenarios. Additionally, direct experimental validation remains necessary to confirm pore-scale optical coupling effects.

5. Conclusion

This study investigated the influence of skin pore microgeometry NIR photon transport within anatomically realistic, MRI-segmented cadaveric head models using a cascaded simulation framework that links ballistic photon propagation at the source–scalp interface with subsequent diffusive photon transport in deeper tissues. By explicitly resolving pore-scale surface structures during the ballistic phase and propagating the resulting photon states into diffusive transport simulations, the present work provides a multiscale framework for evaluating how microscopic surface features influence macroscopic photon propagation.

Skin pore microgeometry introduces localized perturbations during the ballistic photon propagation phase, increasing photon angular divergence and reducing forward-directed energy delivery toward deeper tissues. However, due to the strongly scattering nature of scalp and skull tissues, these directional perturbations are rapidly randomized within the first few millimeters of propagation and therefore do not measurably alter DS profiles or photon path statistics in the diffusive regime. In contrast, the pore-induced reduction in photon weights represents a persistent loss of incident optical energy that directly modulates the total photon budget available for transport through deeper tissue layers. At the system level, this attenuation affects both reflective and transmissive-mode optical configurations, although its impact may become more pronounced in transmissive geometries where photons traverse longer tissue paths. Collectively, these findings indicate that pore-scale surface microgeometry primarily influences optical coupling efficiency rather than the stochastic characteristics of photon transport, providing practical guidance for the modeling and design of biomedical optical measurement systems.

Funding

This work was supported by Acibadem Mehmet Ali Aydınlar University General Research Fund with Grant Number FBA20232164.

Conflicts of interest

The authors have nothing to disclose.

Data availability statement

Data associated with this article are available from the corresponding author upon reasonable request.

Author contribution statement

Conceptualization, S.I.Y.; Methodology, S.I.Y.; Software, S.I.Y.; Validation, S.I.Y.; Formal Analysis, S.I.Y.; Resources, M.E.A.; Data Curation, S.I.Y.; Writing–Original Draft Preparation, S.I.Y.; Writing–Review & Editing, S.I.Y. and M.E.A.; Visualization, S.I.Y.; Supervision, M.E.A.; Project Administration, M.E.A.

Ethics approval

This study did not involve new human or animal experiments. All cadaveric MRI data used in this work were acquired as part of a previously published study [23], for which ethical approval and informed consent procedures were obtained and are described in detail in the original publication. The present study exclusively reuses these previously approved and published data.

Supplementary material

Section A: Magnetic Resonance Imaging–Derived Cadaveric Head Anatomy.

Section B: Ballistic Photon Propagation Simulations.

Section C: Diffusive Photon Transport Simulations.

Access Supplementary Material

Acknowledgments

The authors gratefully acknowledge Acibadem Mehmet Ali Aydınlar University for supporting this research. The authors also thank MD Baran Bozkurt (Department of Neurosurgery and Neuroanatomy Laboratory, School of Medicine, Acibadem Mehmet Ali Aydınlar University, Istanbul, Türkiye) for providing the facilities used for the MRI of cadaveric heads acquired in a previous study, which were utilized in the present work.

References

  1. Jöbsis FF, Noninvasive, infrared monitoring of cerebral and myocardial oxygen sufficiency and circulatory parameters, Science 198(4323), 1264–1267 (1977). [CrossRef] [PubMed] [Google Scholar]
  2. Richards-Kortum R, Sevick-Muraca E, Quantitative optical spectroscopy for tissue diagnosis, Ann. Rev.Phys. Chem. 47, 555–606 (1996). [Google Scholar]
  3. Ntziachristos V, Ripoll J, Wang LV, Weissleder R, Looking and listening to light: the evolution of whole-body photonic imaging, Nat. Biotechnol. 23, 313–320 (2005). [Google Scholar]
  4. Anderson RR, Parrish JA, Selective photothermolysis: precise microsurgery by selective absorption of pulsed radiation, Science 220(4596), 524–527 (1983). [Google Scholar]
  5. Jacques SL, Optical properties of biological tissues: a review, Phys. Med. Biol. 58(11), R37–R65 (2013). [CrossRef] [Google Scholar]
  6. Bashkatov AN, Genina EA, Kochubey VI, Tuchin VV, Optical properties of human skin, subcutaneous and mucous tissues in the wavelength range from 400 to 2000 nm, J. Phys. D: Appl. Phys. 38(15), 2543–2555 (2005). [Google Scholar]
  7. Arridge SR, Optical tomography in medical imaging, Inverse Probl. 15(2), R41–R93 (1999). [NASA ADS] [CrossRef] [Google Scholar]
  8. Haskell RC, Svaasand LO, Tsay T, Feng T, McAdams MS, Tromberg BJ, Boundary conditions for the diffusion equation in radiative transfer, J. Opt. Soc. America A 11(10), 2727–2741 (1994). [Google Scholar]
  9. Wang L, Jacques SL, Zheng L, Mcml–monte carlo modeling of light transport in multi-layered tissues, Comp. Meth. Progr. Biomed. 47(2), 131–146 (1995). [Google Scholar]
  10. Kienle A, Patterson MS, Dögnitz N, Bays R, Wagnières G, van den Bergh H, Noninvasive determination of the optical properties of two-layered turbid media, Appl. Opt. 37(4), 779–791 (1998). [Google Scholar]
  11. Dehghani H, Srinivasan S, Pogue BW, Gibson A, Numerical modelling and image reconstruction in diffuse optical tomography, Phil. Trans. R. Soc. A 367(1900), 3073–3093 (2009). [Google Scholar]
  12. Zhang Y, Chen B, Li D, Propagation of polarized light in the biological tissue: a numerical study by polarized geometric monte carlo method, Appl. Opt. 55(10), 2681–2691 (2016). [Google Scholar]
  13. Durduran T, Choe R, Baker WB, Yodh AG, Diffuse optics for tissue monitoring and tomography, Rep. Prog. Phys. 73(7), 076701 (2010). [Google Scholar]
  14. Griffiths C, Jonathan B, Bleiker T, Hussain W, Simpson R, editors, Rook’s Textbook of Dermatology (Wiley-Blackwell, Oxford, UK, 10th edn, 2024). [Google Scholar]
  15. Páez MCZ, Gusmão LVV, Cartell AS, Wortsman X, de Almeida CÁ, Canella Moraes do Carmo CCM, Duque-Estrada B, Anatomy of scalp: clinical, imaging, and histological aspects, in Atlas of Diagnostic Imaging in Dermatology, (Springer Nature 2025), pp. 25–35. [Google Scholar]
  16. Zhou Y, Chan KKH, Lai T, Tang S, Characterizing refractive index and thickness of biological tissues using combined multiphoton microscopy and optical coherence tomography, Biomed. Opt. Exp. 4(11), 38–50 (2012). [Google Scholar]
  17. Chandrasekhar S, Radiative transfer. Classic monograph establishing the foundations of radiative transfer theory (Dover Publications, New York, 1950). [Google Scholar]
  18. Mishchenko MI, Travis LD, Lacis AA, Multiple scattering of light by particles: radiative transfer and coherent backscattering. in Modern systematic treatment of radiative transfer and similarity relations (Cambridge University Press, Cambridge, 2006). [Google Scholar]
  19. Fukui Y, Ajichi Y, Okada E, Monte carlo prediction of near-infrared light propagation in realistic adult and neonatal head models, Appl. Opt. 42(16), 2881–2887 (2003). [Google Scholar]
  20. D. A. Boas, C. Pitris, N. Ramanujam, editors, Handbook of Biomedical Optics (CRC Press, Boca Raton, FL, 2011). [Google Scholar]
  21. Li W, Liu C, Song P. Unified gas-kinetic particle method for frequency-dependent radiation transport (2023). [Google Scholar]
  22. Di Rocco HO, Carbone NA, Iriarte DI, Pomarico JA, Sandoval HFR, Modeling transition diffusive–nondiffusive transport in a turbid media and application to time-resolved reflectance, J. Quant. Spectrosc. Radiat. Transf. 120, 16–22 (2013). [Google Scholar]
  23. Yöner SI, Aksoy ME, Südor HC, İzzetoğlu K, Bozkurt B, Dinçer A, Drsvision: a machine learning tool for cortical region-specific fnirs calibration based on cadaveric head mri, Sensors, 25(20), 6340 (2025). [Google Scholar]
  24. COMSOL AB. COMSOL Multiphysics (COMSOL AB, Stockholm, Sweden, 2024). Version 6.x. [Google Scholar]
  25. COMSOL AB. Ray Optics Module User’s Guide (COMSOL AB, Stockholm, Sweden, 2024). Version 6.2. [Google Scholar]
  26. Born M, Wolf E, Principles of Optics, 7th edn. (Cambridge University Press, Cambridge, 1999). [Google Scholar]
  27. Jeulin D. Analysis and modeling of 3d microstructures, in Mathematical Morphology (John Wiley & Sons, Hoboken, NJ, 2013), pp. 421–444. [Google Scholar]
  28. Shaiek A, Flament F, Francois G, Lefebvre-Descamps V, Barla C, Vicic M, Giron F, Bazin R, A new tool to quantify the geometrical characteristics of facial skin pores. changes with age and a making-up procedure in caucasian women, Skin Res. Technol. 23, 249–257 (2016). [Google Scholar]
  29. Makki S, Barbenel JC, Agache P, A quantitative method for the assessment of the microtopography of human skin, Acta Dermato-Venereologica 59, 285–291 (1979). [Google Scholar]
  30. Omotezako T, Rodrigues MR, Neo Eleanor Y Wang GW, Chong B, Matsubara A, Subvisible microscale texture is present on the crista cutis of the skin and interacts with incident light to create a soft and radiant “shitsukan” appearance, Int. J. Cos. Sci. 47(4), 691–702 (2025). [Google Scholar]
  31. Haeussinger FB, Heinzel S, Hahn T, Schecklmann M, Ehlis AC, Fallgatter AJ. Simulation of near-infrared light absorption considering individual head and prefrontal cortex anatomy: implications for optical neuroimaging, PLoS One 6(10), e26377 (2011). [Google Scholar]
  32. Mahbub P, Leis J, Macka M, Chemometric approach to the calibration of light emitting diode based optical gas sensors using high-resolution transmission molecular absorption data, Analyt. Chem. 90(10), 5973–5976 (2018). [Google Scholar]
  33. Fang Q, Boas DA, Monte carlo simulation of photon migration in 3d turbid media accelerated by graphics processing units, Opt. Exp. 17(22), 20178–20190 (2009). [Google Scholar]
  34. Texas Instruments. OPT101 monolithic photodiode and single-supply transimpedance amplifier data sheet. Technical datasheet (2015). Rev. B. Available from Texas Instruments website. [Google Scholar]
  35. Vishay Intertechnology. BPW34 silicon pin photodiode datasheet (Technical datasheet, 2011). Rev. 2.1. Angle of half sensitivity φ = ±65°; Silicon PIN photodiode. Available from Vishay website. [Google Scholar]
  36. Moradi M, Chen Y, Monte carlo simulation of diffuse optical spectroscopy for 3d modeling of dental tissues, Sensors 23, 5118 (2023). [Google Scholar]
  37. Bürmen M, Pernuš F, Naglic P. Mcdataset: A public reference dataset of monte carlo simulated quantities for multilayered and voxelated tissues computed by massively parallel pyxopto python package, J. Biomed. Opt. 27, 083012 (2022). [Google Scholar]
  38. The MathWorks, Inc., MATLAB (The MathWorks, Inc., Natick, MA, USA, 2024). Version R20xx. [Google Scholar]
  39. Python Software Foundation, Python (Python Software Foundation, Wilmington, DE, USA, 2024). Version 3.12. [Google Scholar]
  40. Das SR, Avants BB, Pluta J, Wang H, Suh JW, Weiner MW, Mueller SG, Yushkevich PA. Measuring longitudinal change in the hippocampal formation from in vivo high-resolution T2-weighted MRI, NeuroImage 60(2), 1266–1279 (2012). [Google Scholar]
  41. Tong Y, Chen Q, Nichols TE, Rasetti R, Callicott JH, Berman KF, Weinberger DR, Mattay VS, Seeking optimal region-of-interest (roi) single-value summary measures for FMRI studies in imaging genetics, PLoS One 11(9), e0151391 (2016). [Google Scholar]
  42. Wyman DR, Patterson MS, Wilson BC, Similarity relations for anisotropic scattering in monte carlo simulations of deeply penetrating neutral particles, J. Computat. Phys. 81(1), 137–150 (1989). [Google Scholar]
  43. Graaff R, Aarnoudse JG, Zijp JR, Sloot PMA, de Mul FFM, Greve J, Koelink MH, Reduced light-scattering properties for mixtures of spherical particles: a simple approximation derived from mie calculations, Appl. Opt. 31(10), 1370–1376 (1992). [Google Scholar]
  44. Okada E, Delpy DT. Near-infrared light propagation in an adult head model. ii. effect of superficial tissue thickness on the sensitivity of the near-infrared spectroscopy signal, Appl. Opt. 42(16), 2915–2921 (2003). [Google Scholar]
  45. Borner K, Blood PD, Silverstein JC, Ruffalo M, Satija R, Teichmann SA, Pryhuber GJ, Misra RS, Purkerson JM, Fan J, Hickey JW, Molla G, Xu C, Zhang Y, Weber GM, Jain Y, Qaurooni D, Kong Y, HRA Team, Buckle A, Herr BW II, Human biomolecular atlas program (hubmap): 3d human reference atlas construction and usage, Nat. Meth. 22, 845–860 (2025). [Google Scholar]
  46. Tak S, Uga M, Flandin G, Dan I, Penny WD, Sensor space group analysis for fnirs data, J. Neurosci. Met. 264, 103–112 (2016). [Google Scholar]
  47. Yodh A, Chance B, Spectroscopy and imaging with diffusing light, Phys. Today 48(3), 34–40 (1995). [Google Scholar]

All Tables

Table 1

Optical properties of each tissue layer used in the MCX simulations, including absorption coefficient, scattering coefficient, anisotropy factor, and refractive index, evaluated at a wavelength of 735 nm [31].

In the text
Table 2

Summary of detected photon statistics across eight cadaveric head models for two vertical pore alignments. The data highlights the inter-subject variability in total photon counts and mean photon weights relative to the systematic shifts induced by sub-millimeter changes in pore proximity.

In the text

All Figures

Thumbnail: Figure 1 Refer to the following caption and surrounding text. Figure 1

Example COMSOL ray–tracing plots at t = 10 ps for the mean cadaveric head with a 250, μm diameter skin pore, zoomed to the source–scalp interface. Skin pore at vertical (a) −1.5, (b) −0.8, and (c) −0.7 mm positions. The color bar represents normalized photon weights (0–1).
In the text
Thumbnail: Figure 2 Refer to the following caption and surrounding text. Figure 2

Ballistic photon propagation simulation ejected photon distributions as a function of radial distance from the optical source center and angle (degrees off-axis relative to the optical source center) for a 500 μm diameter pore. Plots correspond to vertical pore positions (a) 1, (b) 0.5, and (c) 0 mm off-axis relative to the optical source center. The color bars represent normalized photon weights.
In the text
Thumbnail: Figure 3 Refer to the following caption and surrounding text. Figure 3

Isometric views of spatial fluence color maps for the mean cadaveric head at a 30 mm source–detector separation: plots correspond to 500 diameter pore at vertical pore positions (a) 1, (b) 0.5, and (c) 0 mm off-axis relative to the optical source center. Black circles and squares indicate the optical source and detector positions, respectively. The color bars represent photon fluence on a logarithmic scale to enhance visibility across a wide dynamic range.
In the text
Thumbnail: Figure 4 Refer to the following caption and surrounding text. Figure 4

Representative DS profile for a 30 mm source–detector separation with a 500 μm diameter pore (positioned at 0.5 mm off-axis relative to the optical source center). The DS curves obtained across all individual cadaveric models and pore positions were found to be virtually identical, following the same spatial distribution. For visual clarity, this single representative profile is shown to illustrate the system’s sensitivity to underlying tissue volumes, which remains stable despite surface-level intensity fluctuations.
In the text
Thumbnail: Figure 5 Refer to the following caption and surrounding text. Figure 5

Box plots for detected photon counts across eight cadaveric heads at two distinct vertical pore positions: 1, and 0.5 mm off-axis relative to the optical source center.
In the text
Thumbnail: Figure 6 Refer to the following caption and surrounding text. Figure 6

Box plots for detected mean photon weight across eight cadaveric heads at two distinct vertical pore positions: 1, and 0.5 mm off-axis relative to the optical source center.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.