J. Eur. Opt. Society-Rapid Publ. 22, 48( 2026) 473
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 ].
Layers |
l a( mm �1) |
l s( mm �1) |
g |
n |
Air |
0 |
0 |
0 |
1 |
Scalp |
0.016 |
19 |
0.9 |
1.6 |
Skull |
0.018 |
16 |
0.9 |
1.56 |
CSF |
0.004 |
0.3 |
0 |
1.33 |
Brain |
0.09 |
21.5 |
0.9 |
1.4 |
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Þ ¼ XNx; y; z w i; i¼1 ð9Þ
where N x, y, z denotes the number of photon packets traversing voxel( x, y, z) andw i 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 mm 3, 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 skinpore 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( l a), scattering coefficient( l 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 10 9 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 1mm 3, 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