Open Access
Add to BibSonomy
Abstract
Many promising biomedical applications have been proposed for terahertz (THz) spectroscopy and diagnostic imaging techniques. Polarimetric imaging systems are generally useful for enhancing imaging contrasts, yet the interplay between THz polarization changes and the random discrete structures in biological samples is not well understood. In this work, we performed Monte Carlo simulations of the propagation of polarized THz waves in skin and adipose tissues based on the Mie scattering from intrinsic structures, such as hair follicles or sweat glands. We show that the polarimetric contrasts are distinctly affected by concentration, size and dielectric properties of the scatterers, as well as the frequency and polarization of the incident THz waves. We describe the experimental requirements for observing and extracting these polarimetric signals due to the low energy and small angular spread of the back-scattered THz radiation. We analyzed the spatially integrated Mueller matrices of samples in the normal-incidence back-scattering geometry. We show that the frequency-dependent degree of polarization (DOP) can be used to infer the concentrations and dielectric contents of the scattering structures. Our modeling approach can be used to inform the design of the imaging modalities and the interpretation of the spectroscopic data in future terahertz biomedical imaging applications.
© 2024 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
The $30^{th}$ anniversary of open source Monte Carlo (MC) codes for modeling light propagation in tissue was recently celebrated in 2022 [1]. Over this period, MC models describing photon transport in biological tissues have been used to aid the design of diagnostic devices and to explain the observed phenomena. Examples of instrumental designs informed by MC modeling include the source-detector separations [2], angular dependence [3], polarization-sensitive methods [4,5] and coherence-domain techniques [6,7]. Similarly, MC codes are used to predict the response of samples from their bulk optical properties [8] or by use of the inversion techniques [9,10].
During the same period, terahertz (THz) spectroscopy methods have also matured rapidly, which has facilitated the investigation of various biological testing and several in vivo and clinical imaging applications [11–14]. A favored approach in THz biomedical sensing is to extract the dielectric permittivity of the sample from their reflectance or transmittance [15–17]. Further information on sample composition can be gained by applying appropriate effective media models, which relate the THz dielectric spectrum to the fractions and physical parameters of the presumed components. In particular, there has been significant success in the analysis of hydration level of the tissues [18–21], owing to the strong THz absorption by free and bound water. Still, this approach only provides a homogeneous solution, i.e., the loss of THz radiation in a spatial area of the sample is fully described by the Fresnel’s and Beer’s laws [22,23].
The presence of nonuniform dielectric structures can generally produce disturbance of light propagation in an otherwise homogeneous medium [24,25]. However, whether such a disturbance can measurably be confirmed in THz biomedical applications has remained as a point of concern in the scientific community, considering the strong dielectric absorption of the homogeneous medium and reduced scattering at longer wavelengths. Previously, such deviations from bulk optical response of a sample has been experimentally reported in biological tissue with existing THz systems. Feldman et al. have reported that the THz spectra of human skin can be notably affected by the eccrine sweat glands [26], given that active glands have high proton conductivity and are thus analogous to helical antennas [27]. Separately, we have found that the THz response of rat skin with burn injuries is not only accounted by the level of water contents, but also correlates with the density of discrete structures recognized in the sample histology [28–30]. In addition to skin tissue, the idea of combining hydration analysis and correlation to sample non-uniformity has also been proposed for other types of tissue [31,32]. The potential of these investigation comes from two aspects: (1) The presence of nonuniform structures (either intrinsic or extrinsic) opens up new modalities for THz sensing; (2) These nonuniform structures often play an important role in biological processes such as wound healing, diagnostics and therapy [33].
With the enhancement in SNR and expansion of measurement capabilities, new paradigms are emerging in THz biomedical sensing based on alternative contrast mechanisms due to complex heterogeneous, birefringent and stratified sample profiles. For instance, recent THz ellipsometric measurements of in vivo human skin have revealed anisotropic cellular structure of the stratum corneum, birefringence and occlusion effects as well as sensitivity of hydration measurements to surface scattering [34–36]. Although such scattering effects can be mitigated using signal processing methods [37,38], they can also be diagnostically useful and can be measured using highly-accurate and calibrated broadband polarimetry techniques [39,40].
Our objective in this work is to develop predictive models to relate the nonuniform structures in biological samples to THz scattering measurements, and to explain such observations. Similar to the optical regime [41], several approaches can be implemented for THz scattering measurements: (1) angular-resolved magnitude measurements, (2) broadband spectroscopic measurements at a fixed angle, and (3) measurements of coherence or polarization state of the light beam [42]. Angular or frequency dependent THz systems can readily characterize the scattering in granular materials [43–46]. Given the abundance of components and heterogeneity in biological samples, however, it becomes challenging to separate the contribution of scattering from featureless absorption or measurement noise [47]. Meanwhile, THz polarization imaging can potentially provide better contrasts in tissue delineation [48,49]. Still, unlike examples in the optical regime demonstrated by polarization gating [50], selective probing of superficial tissue structures has not been reported in THz frequencies.
In this paper, we preformed MC simulations of the THz radiative transfer in biological tissues as the function of detecting position, frequency and polarization of THz sources. Previously, Ney et al. used similar MC codes to investigate the THz Mie scattering in skin induced by extrinsic nanoparticles [51], whereas here we considered the intrinsic structures that have comparable dimensions to the THz wavelengths. Upon the introduction of scattering particles, we found a diffuse component in the THz reflection from skin, which exhibit different degree of polarization (DOP) as the functions of source frequencies, or the concentration and inside contents of particles. Our approach is motivated by future application of our polarization-sensitive handheld scanner [52,53] for in vivo and longitudinal diagnosis of skin burns [54–56]. Therefore, this paper is focused on similar experimental conditions of instrument and sample type and scattering particles. We hope that the widespread adoption of the simulation algorithms introduced in this work will also be useful for other applications of THz biophotonics [57].
2. Methods
In this section, we provide a detailed description of the numerical models needed for simulating the THz beam propagation in tissue, including the tissue models, the MC algorithms, and data analysis approaches to extract the polarimetric signal contrasts. We have adopted two separate algorithms to model the propagation of THz photons in stratified layers and to analyze the polarization of THz waves after scattering by discrete skin structures. The former issue is addressed by the Monte Carlo Multi-Layered (MCML) codes provided by Wang et al. [58,59], which has gained great success in the optical regime for modeling tissue structures, e.g., in skin and eyeball [60]. The latter is solved by the Meridian Plane method [5,59], which we will refer to as the Polarized-Light MC codes. Further, we made modifications to the original MC codes based on the theory of Mie scattering in absorbing media as described below.
2.1 Modeling of the skin structures
Due to the strong absorption by water, THz radiation has a penetration depth around 0.1–0.5 mm into the skin [12,61], which is less than the typical thickness of the dermis. Therefore, a three-layer [26,51] or two-layer [62,63] approximation is often adopted for simulating the THz spectra of skin. Figure 1 shows the dielectric functions of the epidermis and bulk water calculated using the Double-Debye dispersion model [15]. The THz properties of the dermis tissue are generated by assuming the Maxwell-Wagner mixing law of water and the dry skin components, following the approach described in Ref. [51,64]. For a comparison study, we also considered the adipose tissue as a representative case for less lossy media, since the lipids attenuate THz radiation less strongly than skin. Yanina et al. have reviewed THz measurements of adipose tissue [65], and we have adopted the one reported in Ref. [66], as plotted in Fig. 1.
Fig. 1. The refractive index (a) and absorption coefficient (b) of epidermis, water and adipose tissue in the frequency range between 0.1 and 2 THz.
Another aspect of our skin model is to incorporate discrete particles with size and density resembling observed conditions. If the discrete particle with a radius, $r$, is much smaller than the wavelength ($\lambda$), i.e., the size parameter $x = 2{\pi }r/\lambda \ll 1$, the resultant scattering will be in the Rayleigh regime. For example, between 0.1 and 2 THz, the size parameter $x$ for a 10 $\mu$m diameter cell varies from 0.01 to 0.21. Although large cells may still exhibit Mie scattering at higher THz frequencies, we have focused on skin appendages, which are more likely to induce measurable signal contrasts. The scattering structures proposed by previous THz studies, including sweat glands [62], and hair follicles [67] having an average diameter and length between 100 and 500 $\mu$m [30], are represented by "equivalent spheres" for simplicity [68], which give rise to the size parameter larger than 0.2, in the frequency range between 0.1 and 1 THz [45]. The concentration of particles were calculated from histological sections, and two dielectric constant values were simulated to resemble the skin with both empty (air) and filled (water) sweat glands.
2.2 Monte Carlo simulation algorithms
There are many open-source resources and tutorials available for MC simulation of radiative transport [60,69], including codes and manuals that are readily accessible to the scientific community [59]. Here we focus on providing a brief background necessary for accurate modeling of the skin appendages in the THz frequency range. When a new photon packet is launched in a MC simulation, it travels for a distance of $s$ before experiencing an interaction, such as absorption or scattering. The sampled value of $s$ is related to the optical properties of a turbid medium by,
where $\xi$ is a random number between 0 and 1, and $\mu _{\rm e}$ is the extinction coefficient of the medium. In general, these interactions involve either scattering or absorption. The impact of absorption is described by the albedo, $\alpha$, of the medium: where, and $\mu _{\rm s}$ ($\mu _{\rm a}$) is the scattering (absorption) coefficient of the medium. In each scattering event, the weight of photon (denoting its energy) drops from $W$ to $W\times \alpha$.Subsequently, it is necessary to define the trajectory of the photon after each scattering event. In the MCML code, the photon trajectory is determined based on the Henyey–Greenstein (HG) phase function [58],
which empirically describes the angular pattern of scattering. $P_{\rm HG}$ depends on only the asymetric factor $g$, which is the average cosine of the scattering angle $\theta$.In the Polarized-Light MC codes, this step is instead based on the Mie theory, which derives the relation between incident ($\boldsymbol {\rm E_i}$) and scattered ($\boldsymbol {\rm E_s}$) field amplitudes as,
Equation (5) uses the Jones calculus based on the perpendicular and parallel components to the scattering plane. The scattering matrix only consists of two elements, $S_1(\theta )$ and $S_2(\theta )$, due to the spherical symmetry of particle [70]. The Mueller calculus transform of Eq. (5) is usually adopted so that the operation only involves real values, which is given by,
In addition to the new trajectory, Eq. (6) also determines the polarization change of the photon after each scattering event, which forms the foundation of the Polarized-Light MC codes. On the other hand, the MCML codes make fewer assumptions on the nature of turbid media, and also account for the Fresnel reflections which occur in the stratified structures. We have applied the two methods for different purposes, as described in the logical flowchart of our simulation presented in Fig. 2. While some authors have combined the merits of both methods to develop their own codes [51,71], alternatively we have used the two open-source codes [59] without large modification. Meanwhile, we have incorporated a Mie calculator to generate the required input parameters based on media models, which will be introduced in the following section.
Fig. 2. Flowchart of the input and output data in the steps of numerical simulation for THz radiation transport in skin model.
2.3 Mie scattering in absorbing host media
As introduced earlier, the polarized light MC codes rely on the solutions to Mie scattering as input parameters, in particular the complex functions $S_1$ and $S_2$. Meanwhile, the MCML codes only depend on $\mu _{\rm s}$ ($\mu _{\rm a}$) and $g$, which are often obtained from direct measurements in optical systems. In the THz regime, however, such measurements on the scattering effects in biological tissues are still lacking or inconclusive. Therefore, we need to simulate these properties under Mie scattering assumptions as well.
Conventionally, the solutions to Mie scattering have assumed the host medium to be non-absorbing, i.e., its refractive index $n_{\rm h}$ being a real value. This has been a reasonable approximation in optical regime, where the scattering effects of tissues are usually much stronger than absorption effects. On the contrary, THz waves experience strong absorption due to the water content, and less scattering as the wavelengths become much larger than that of visible light. Also, neglecting the imaginary part of $n_{\rm h}$ would introduce significant error to the boundary conditions at the particle-medium interfaces. Therefore, it is necessary to modify the solutions to Mie scattering to account for the absorbing host media. The derivation of such modifications has been presented in the literature [72–74]. Here, we will briefly introduce key relationships for completeness and to avoid ambiguity in defenition of relevant parameters.
In absorbing host media, absorption occurs in both the host medium and the particle-medium boundary. We have adopted the formulas proposed by Yang et al. [73], where the labels ${\rm p}$ and ${\rm h}$ denote the properties of particle and host medium, respectively. The extinction of a photon packet during the $j$th order of scattering event is expressed by,
Thus, the absorption by host medium is a continuous reduction of $W$ along the optical path ($s_{j-1}$), whereas the absorption of particle is still given by multiplying a ratio of $\alpha$, which is
In the modified Mie theory, the average scattering (absorption) coefficients due to particles are derived as,
Fig. 3. The scattering efficiency $Q_{\rm sca}$ and asymmetric factor $g$ for a spherical particle of 50, 150 , 300 or 1000 $\mu$m diameter, with different particle contents and host media.
3. Results
3.1 Distribution of the THz radiation energy
To study the propagation of THz radiation inside skin, we start with the MCML code, as outlined in Fig. 2, which generates the traces of photon packets based on the bulk optical properties of different tissue layers. Figure 4 shows a point spread function simulation of photon trajectory at 0.6 THz, when the dielectric function of the skin appendages are varied in parts (a) and (b), but under the same assumption of host media, i.e., the thickness ($d$), refractive index ($n$) and absorption coefficient ($\mu _{\rm a,h}$) of the epidermis and dermis layers are the same. We assumed scattering particles of 150 $\mu$m diameter with the concentration of 125 mm$^{-3}$, based on the approach mentioned in section 2.1. using representative skin histology. The dielectric properties of the particles are chosen to represent the two extreme cases of scattering efficiency, i.e., air and water. The former case represents empty sweat glands or hair follicles, while the latter corresponds to structures like filled sweat glands. In both scenarios simulated in Fig. 4, the likelihood of photon absorption (blue lines) are higher than escaped photons upon scattering (red lines). Also, it is evident that the majority of reflected photon packets have only travelled inside the epidermis. Here a thickness of 100 $\mu$m was assumed for epidermis. Since the epidermal thickness varies at different body sites and between individuals [76], there will be more photon packets reaching the dermis if a thinner thickness value is used. Comparing the scattering from air and water particles, the penetration depths appear to be similar, since the attenuation are largely contributed by the absorption in the host media ($\mu _{\rm a,h}$). Still, there are considerable differences in the spatial probability distributions of photons (i.e., the shape of the interaction volume). This is essentially a visual interpretation of the asymmetric factor $g$, whose values are previously calculated in Fig. 3. In other words, when $g$ approaches to 1 a larger degree of scattering in the forward direction is specified. Another noticeable change in Fig. 4(b) is that less photon packets have escaped from skin when the scattering particles are filled with water, which can be attributed to (1) the absorption due to these particles ($\mu _{\rm a,p}$) and (2) smaller mismatch between the refractive index of particles and host media (see Fig. 1). These results suggests that the complex refractive indices have been reasonably incorporated into our MC simulation.
Fig. 4. Simulation of photon propagation, when spherical scatterers consisting of (a) air or (b) water are embedded in skin. Here we used radiation frequency of 0.6 THz, particle diameter of 150 $\mu$m and mean particle distance of 200 $\mu$m (center to center).
Before continuing with the discussion of polarimetric signal contrasts, we need to make some clarifying comments regarding results shown in Fig. 4. First, we have considered a quartz slab above the skin in these simulations, which corresponds to a microscope slide for ex vivo samples or an optical window for in vivo imaging. In the THz regime, the refractive indices of glassy materials (e.g., quartz) are often comparable to or higher than that of skin [77]. If we choose not to apply the quartz layer, the traces of photons change consequently due to total internal reflections of THz waves at the epidermis-air interface. Second, the ratio and spatial distribution of scattered intensities are different from what can be visually estimated from Fig. 4, since the weight information, $W$, is not carried by the ray traces. This information can be discerned, however, from Fig. 5, which shows the probability distributions of the exiting photons as the function of the number of scattering events. The probability distribution functions shown by the histogram (blue bars) suggest that the radiation at 0.6 THz experience different orders of scattering depending on the absorption of the host media. However, the weighted distributions of exiting photons (the energy ratio) reveal that the scattered intensity is mostly due to the single-scattering photons. The red curves show that there is also a minor contribution by the double-scattering photons, whereas the impact of higher-order scattering processes are generally negligible. Note that this conclusion is valid for both weak and high absorption biological tissues, as represented by adipose in Fig. 5(a) and epidermis in Fig. 5(b), respectively. Therefore, the predictions of a single-scattering model will be used to explain other characteristics of the simulation results, namely the polarimetric signatures of the tissue, which will be shown in the following sections.
Fig. 5. The distribution of scattered photons as the function of scattering events in the case of (a) epithelial spheres embedded in adipose, (b) air spheres embedded in epidermis. Inset: assumed optical properties of tissues at 0.6 THz.
3.2 Distribution of the THz polarization states
In addition to the distribution of power, the polarization changes of THz waves occurred during propagation are investigated using the Polarized-Light MC algorithm, as described in Fig. 2 flowchart. This algorithm can generate the spatial maps of the four Stokes parameters ($I, Q, U, V$), while four independent polarization states of the incident light can be selected as the input values. Subsequently, the $4\times 4$ combination of maps can be transformed into the Mueller matrix of the tissue, as shown in Fig. 6 with the normalized values. Our results show that $m_{14}$ and $m_{41}$ are close to zero regardless of the biological tissues, and that certain elements (e.g., $m_{22}$ and $m_{33}$) have symmetric patterns, similar to the results obtained in the optical regime [78,79]. However, it will be difficult experimentally to observe these spatial (angular) patterns with current THz systems, since most of the scattered radiation energy is distributed in a narrow imaging area (angular spread), as shown by the map of $m_{11}$.
Fig. 6. Maps of the normalized Mueller matrix in a $0.2\times 0.2$ cm$^2$ imaging area for 150 $\mu$m epithelial spheres in adipose tissue, which represents the impulse response to a normal incident beam at 0.6 THz.
Instead, we use the pixel sums of these maps to simulate the Mueller matrix captured by an on-axis single-element detector. This matrix M for Fig. 6 is given by,
Here, we note that M is a diagonal matrix and it complies with the additional condition that
which agrees with theoretical predictions resulted from assuming an ensemble of single scatterers, each of them having a reciprocal Mueller-Jones matrix [80,81]. Therefore, this result is consistent with the PDFs shown in Fig. 5. When $m_{22}=-m_{33}$ is also satisfied, the Mueller matrix can be rewritten by substituting $m_{11}$ and $m_{22}$ with symbol $a_1$ and $a_2$, respectively, leading to, where $p$ = $a_2$/$a_1$. This legacy form of the Mueller matrix has been substantiated for the back-scattering from non-chiral particles with a rotational symmetry in a cloud of particles [81,82]. This matrix can be decomposed as a uniform attenuator, i.e., parameter $a_1$, and a nonuniform partial depolarizer, i.e., ${\rm diag}[1, p, -p, 1-2p]$ [83]. Therefore, the polarization changes of the scattered fields ($\boldsymbol {\rm E_s}$) are sufficiently described by the Degree of Polarization (DOP), defined as,According to Eq. (15), the DOP of $\boldsymbol {\rm E_s}$ equals $|p|$ when incident light $\boldsymbol {\rm E_i}$ is linearly polarized (LP), and equals $|1-2p|$ when $\boldsymbol {\rm E_i}$ is circularly polarized (CP). The linear and circular DOP are hence correlated and their order depends on the range of $p$.
Figure 7 presents the DOP of THz waves obtained with the same particle size and density values as previously set in Fig. 4. In general, the linear and circular DOP show similar dependence on the THz frequency, because they are both determined by the the frequency dependence of $p$. In addition, our results show that $p\in (0.5,1)$ is valid under the three simulation conditions, leading to a higher DOP of LP than that of CP ($p>2p-1>0$). This is opposite to the usual outcome for the Mie scattering in tissues in the optical regime [84], where $|m_{33}|<|m_{44}|$. We attribute this obervation in the THz frequencies to the prevalence of single-scattering and the small spread of scattering angle $\theta$, considering that singly scattered radiation stays fully polarized at $\theta =180^\circ$.
Fig. 7. (a) Diffuse refletance $R_d$ and (b) degree of polarization (DOP) for the total THz radiation scattered by epithelial particles in adipose tissue. (c) and (d) are obtained for air voids in skin. (e) and (f) are obtained for water particles in skin. All simulations use the same particle diameter of 150 $\mu$m and mean particle distance of 200 $\mu$m.
3.3 Concentration dependence of the scattering signal contrasts
In Eq. (5), the total intensity of $\boldsymbol {\rm E_s}$ is proportional to that of $\boldsymbol {\rm E_i}$ by a ratio of $a_{1}$, given by Eq. (15), which is also sometimes referred to as the diffuse reflectance. As shown in Fig. 7, the diffuse reflectance, $R_d$, calculated by the MC simulations is in good agreement with the following relation derived by Morel & Prieur [85],
where absorption coefficient $a$ was approximated by $\mu _{\rm a,h}+\mu _{\rm a,p}$, and the backscattering coefficient $b_b$ was calculated from the integral of $S_{11}(\theta )$ (Eq. (6)) in the backward hemisphere. As shown in Fig. 3, the asymmetric factor $g$ approaches to 1 as the frequency reaches 2 THz, indicating less scattering in the backward direction. The scattering anisotropy, together with the increasing absorption ($a$), explain why $R_d$ is diminished in the higher frequency region.It can be inferred from Eq. (17) that $R_d$ is proportional to $b_b$ (when $b_b\ll {a}$) and thus to the particle concentration $\rho$. Figure 8(a) and (b) shows the variation of $R_d$ and circular DOP when diluting the concentration of air voids in the epidermis, with the red curves ($\rho =125$ mm$^{-3}$) representing the same data as those in Fig. 7(c) and (d). The decrease in $R_d$ are similar at different frequencies, and a linear relation between $R_d$ and concentration can be confirmed at 0.6 THz (in Fig. 8(c)). Meantime, the circular DOP approaches to unity as the concentration decreases (in Fig. 8(d)). This figure provides a reference for estimating the SNR requirement of THz polarimetric systems in order to characterize the scattering structures in skin.
Fig. 8. Variations of (a) the diffuse reflectance $R_d$ and (b) the circular DOP in 0.1–2 THz as changing the concentration of air voids in the epidermis. (c) and (d) compare the concentration dependence of $R_d$ and circular DOP at a single frequency (0.6 THz) between the scatters made by air and water.
4. Discussion
Comparing the results of the THz radiation transfer model in skin with adipose tissue, the absorption in our models rapidly increases in the presence of water content in either the scattering particles or host media, leading to low values of the single-scattering albedo $\alpha$. This broadband behavior can be observed over the range of 0.1–2 THz, over which both the scattering and absorption coefficients rise with higher frequencies. Consequently, the low albedo increases the amount of single scattering relative to multiple scattering. This relation is rather straightforward in the MC simulation. The weight of photons drop by a fixed ratio of $\alpha$ during each interaction, therefore the energy ratio of exiting photons having undergone 1, 2, …, n scattering events decay geometrically with the orders [86,87]. Such geometric series can be seen in Fig. 5, with $\alpha =0.11$ in (a) and $\alpha =0.23$ in (b). The discrete structures of interests can be represented with different particle shape, size or concentration, however, it is unlikely to cause dramatic increase in $\alpha$ and hence the simulation results would still favor single scattering.
At $\theta =180^\circ$, the single-scattering by spherical particles stay as fully polarized [70]. With normal incidence, the specular reflections by the layered interfaces also propagate in this direction. Therefore, the THz waves reflected by skin will be measurably depolarized only if sufficient off-axis scattered rays are collected by the detector. In section 3.2 we investigated the pixel sums in an area wider than the point spread function, simulating a $2\pi$ collection solid angle. The contrast in DOP will be more significant if the off-axis components can be selectively resolved by the detector, which is yet a challenging task since the angular spreads of the THz reflection are rather narrow. Lastly, the angular patterns of scattering are also affected by assuming oblique incidence or non-spherical particles, which is worth further investigation.
5. Conclusion
We have performed MC simulations for polarized THz radiation transfer in the absorbing host media doped with spherical particles, representing the potential Mie scattering by the accessory structures in biological tissues. The tissue models exhibit low single-scattering albedo when water contents are introduced in either the particles or host media, leading to the prevalence of single scattering in the scattered THz radiation energy. In the geometry of normal incidence and total backscatter, the skin models behave as a uniform attenuator with a nonuniform partial depolarizer. We have investigated the diffuse reflectance and DOP of scattered THz waves as a function of frequency with different source polarization states, and in different conditions of the concentration and dielectric constant of scattering particles. These results suggest alternative signal contrasts that are distinctly associated with the skin appendages, and also provide estimation of the SNR requirement for observation of these effects using THz polarimetric instruments. Future work may include validation of the simulation results with controlled measurements using phantoms or ex vivo tissue samples. The simulation models and algorithms developed in this paper will enable and expedite the discovery of future applications of THz biophotonics.
Funding
Stony Brook University; National Institute of General Medical Sciences (GM112693); U.S. Army Medical Research Acquisition Activity (W81XWH-21-1-0258).
Acknowledgments
The authors thank Dr. Mahmoud E. Khani for useful discussions.
This work was supported in part by the U.S. Army Medical Research Acquisition Activity (USAMRAA) through the Military Burn Research Program (MBRP) under Award No. W81XWH-21-1-0258, the National Institute of General Medical Sciences (GM112693), and Stony Brook University.
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
References
1. S. L. Jacques, “History of Monte Carlo modeling of light transport in tissues using mcml.c,” J. Biomed. Opt. 27(08), 083002 (2022). [CrossRef]
2. T. J. Farrell, M. S. Patterson, and B. Wilson, “A diffusion theory model of spatially resolved, steady-state diffuse reflectance for the noninvasive determination of tissue optical properties in vivo,” Med. Phys. 19(4), 879–888 (1992). [CrossRef]
3. J. R. Mourant, J. Boyer, A. H. Hielscher, et al., “Influence of the scattering phase function on light transport measurements in turbid media performed with small source–detector separations,” Opt. Lett. 21(7), 546–548 (1996). [CrossRef]
4. S. Bartel and A. H. Hielscher, “Monte Carlo simulations of the diffuse backscattering mueller matrix for highly scattering media,” Appl. Opt. 39(10), 1580–1588 (2000). [CrossRef]
5. J. C. Ramella-Roman, S. A. Prahl, and S. L. Jacques, “Three Monte Carlo programs of polarized light transport into scattering media: part i,” Opt. Express 13(12), 4420–4438 (2005). [CrossRef]
6. D. A. Boas and A. G. Yodh, “Spatially varying dynamical properties of turbid media probed with diffusing temporal light correlation,” J. Opt. Soc. Am. A 14(1), 192–215 (1997). [CrossRef]
7. M. Kirillin, I. Meglinski, V. Kuzmin, et al., “Simulation of optical coherence tomography images by Monte Carlo modeling based on polarization vector approach,” Opt. Express 18(21), 21714–21724 (2010). [CrossRef]
8. I. V. Meglinski and S. J. Matcher, “Quantitative assessment of skin layers absorption and skin reflectance spectra simulation in the visible and near-infrared spectral regions,” Physiol. Meas. 23(4), 741–753 (2002). [CrossRef]
9. C. R. Simpson, M. Kohl, M. Essenpreis, et al., “Near-infrared optical properties of ex vivo human skin and subcutaneous tissues measured using the Monte Carlo inversion technique,” Phys. Med. Biol. 43(9), 2465–2478 (1998). [CrossRef]
10. A. Roggan, M. Friebel, K. Doerschel, et al., “Optical properties of circulating human blood in the wavelength range 400-2500 nm,” J. Biomed. Opt. 4(1), 36–46 (1999). [CrossRef]
11. E. P. J. Parrott, Y. Sun, and E. Pickwell-MacPherson, “Terahertz spectroscopy: Its future role in medical diagnoses,” J. Mol. Struct. 1006(1-3), 66–76 (2011). [CrossRef]
12. X. Chen, H. Lindley-Hatcher, R. I. Stantchev, et al., “Terahertz (THz) biophotonics technology: Instrumentation, techniques, and biomedical applications,” Chem. Phys. Rev. 3(1), 011311 (2022). [CrossRef]
13. A. Leitenstorfer, A. S. Moskalenko, T. Kampfrath, et al., “The 2023 terahertz science and technology roadmap,” J. Phys. D: Appl. Phys. 56(22), 223001 (2023). [CrossRef]
14. A. G. Markelz and D. M. Mittleman, “Perspective on terahertz applications in bioscience and biotechnology,” ACS Photonics 9(4), 1117–1126 (2022). [CrossRef]
15. E. Pickwell, B. E. Cole, A. J. Fitzgerald, et al., “In vivo study of human skin using pulsed terahertz radiation,” Phys. Med. Biol. 49(9), 1595–1607 (2004). [CrossRef]
16. P. C. Ashworth, E. Pickwell-MacPherson, E. Provenzano, et al., “Terahertz pulsed spectroscopy of freshly excised human breast cancer,” Opt. Express 17(15), 12444–12454 (2009). [CrossRef]
17. M. E. Khani, Z. B. Harris, O. B. Osman, et al., “Triage of in vivo burn injuries and prediction of wound healing outcome using neural networks and modeling of the terahertz permittivity based on the double debye dielectric parameters,” Biomed. Opt. Express 14(2), 918–931 (2023). [CrossRef]
18. D. B. Bennett, W. Li, Z. D. Taylor, et al., “Stratified media model for terahertz reflectometry of the skin,” IEEE Sens. J. 11(5), 1253–1262 (2011). [CrossRef]
19. Z. D. Taylor, J. Garritano, S. Sung, et al., “Thz and mm-wave sensing of corneal tissue water content: In vivo sensing and imaging results,” IEEE Trans. Terahertz Sci. Technol. 5(2), 184–196 (2015). [CrossRef]
20. A. A. Gavdush, N. V. Chernomyrdin, G. A. Komandin, et al., “Terahertz dielectric spectroscopy of human brain gliomas and intact tissues ex vivo: double-debye and double-overdamped-oscillator models of dielectric response,” Biomed. Opt. Express 12(1), 69–83 (2021). [CrossRef]
21. A. Chen, Z. B. Harris, A. Virk, et al., “Assessing corneal endothelial damage using terahertz time-domain spectroscopy and support vector machines,” Sensors 22(23), 9071 (2022). [CrossRef]
22. K. I. Zaytsev, A. A. Gavdush, N. V. Chernomyrdin, et al., “Highly accurate in vivo terahertz spectroscopy of healthy skin: Variation of refractive index and absorption coefficient along the human body,” IEEE Trans. Terahertz Sci. Technol. 5(5), 817–827 (2015). [CrossRef]
23. W. Cheong, S. Prahl, and A. Welch, “A review of the optical properties of biological tissues,” IEEE J. Quantum Electron. 26(12), 2166–2185 (1990). [CrossRef]
24. S. Alekseev and M. Ziskin, “Distortion of millimeter-wave absorption in biological media due to presence of thermocouples and other objects,” IEEE Trans. Biomed. Eng. 48(9), 1013–1019 (2001). [CrossRef]
25. F. Théberge, M. Châteauneuf, J. Dubois, et al., “Spectral artifacts from non-uniform samples analyzed by terahertz time-domain spectroscopy,” Opt. Express 17(13), 10841–10848 (2009). [CrossRef]
26. I. Hayut, A. Puzenko, P. Ben Ishai, et al., “The helical structure of sweat ducts: their influence on the electromagnetic reflection spectrum of the skin,” IEEE Trans. Terahertz Sci. Technol. 3(2), 207–215 (2013). [CrossRef]
27. Y. Feldman, A. Puzenko, P. Ben Ishai, et al., “Human skin as arrays of helical antennas in the millimeter and submillimeter wave range,” Phys. Rev. Lett. 100(12), 128102 (2008). [CrossRef]
28. M. H. Arbab, T. C. Dickey, D. P. Winebrenner, et al., “Characterization of burn injuries using terahertz time-domain spectroscopy,” in Advanced Biomedical and Clinical Diagnostic Systems IX, vol. 7890A. Mahadevan-Jansen, T. Vo-Dinh, W. S. Grundfest MD, eds., International Society for Optics and Photonics (SPIE, 2011), p. 78900Q.
29. M. H. Arbab, T. C. Dickey, D. P. Winebrenner, et al., “Terahertz reflectometry of burn wounds in a rat model,” Biomed. Opt. Express 2(8), 2339–2347 (2011). [CrossRef]
30. M. H. Arbab, D. P. Winebrenner, T. C. Dickey, et al., “Terahertz spectroscopy for the assessment of burn injuries in vivo,” J. Biomed. Opt. 18(7), 077004 (2013). [CrossRef]
31. I. Ozheredov, M. Prokopchuk, M. Mischenko, et al., “In vivo THz sensing of the cornea of the eye,” Laser Phys. Lett. 15(5), 055601 (2018). [CrossRef]
32. A. S. Kucheryavenko, N. V. Chernomyrdin, A. A. Gavdush, et al., “Terahertz dielectric spectroscopy and solid immersion microscopy of ex vivo glioma model 101.8: brain tissue heterogeneity,” Biomed. Opt. Express 12(8), 5272–5289 (2021). [CrossRef]
33. S. A. Eming, P. Martin, and M. Tomic-Canic, “Wound repair and regeneration: Mechanisms, signaling, and translation,” Sci. Transl. Med. 6(265), 265sr6 (2014). [CrossRef]
34. X. Chen, Q. Sun, J. Wang, et al., “Exploiting complementary terahertz ellipsometry configurations to probe the hydration and cellular structure of skin in vivo,” Adv. Photonics Res. 2(1), 2000024 (2021). [CrossRef]
35. J. Wang, H. Lindley-Hatcher, X. Chen, et al., “THz sensing of human skin: A review of skin modeling approaches,” Sensors 21(18), C1 (2021). [CrossRef]
36. H. Lindley-Hatcher, R. I. Stantchev, X. Chen, et al., “Real time THz imaging—opportunities and challenges for skin cancer detection,” Appl. Phys. Lett. 118(23), 230501 (2021). [CrossRef]
37. M. H. Arbab, D. P. Winebrenner, E. I. Thorsos, et al., “Application of wavelet transforms in terahertz spectroscopy of rough surface targets,” in Terahertz Technology and Applications III, vol. 7601L. P. SadwickC. M. M. O’Sullivan, eds., International Society for Optics and Photonics (SPIE, 2010), p. 760106.
38. M. H. Arbab, A. Chen, E. I. Thorsos, et al., “Effect of surface scattering on terahertz time domain spectroscopy of chemicals,” in Terahertz Technology and Applications, vol. 6893K. J. LindenL. P. Sadwick, eds., International Society for Optics and Photonics (SPIE, 2008), p. 68930C.
39. K. Xu, E. Bayati, K. Oguchi, et al., “Terahertz time-domain polarimetry (thz-tdp) based on the spinning e-o sampling technique: determination of precision and calibration,” Opt. Express 28(9), 13482–13496 (2020). [CrossRef]
40. K. Xu, M. Liu, and M. H. Arbab, “Broadband terahertz time-domain polarimetry based on air plasma filament emissions and spinning electro-optic sampling in GaP,” Appl. Phys. Lett. 120(18), 181107 (2022). [CrossRef]
41. A. K. Ray, A. Souyri, E. J. Davis, et al., “Precision of light scattering techniques for measuring optical parameters of microspheres,” Appl. Opt. 30(27), 3974–3983 (1991). [CrossRef]
42. J. Pearce, Z. Jian, and D. M. Mittleman, “Spectral shifts as a signature of the onset of diffusion of broadband terahertz pulses,” Opt. Lett. 29(24), 2926–2928 (2004). [CrossRef]
43. L. M. Zurk, B. Orlowski, D. P. Winebrenner, et al., “Terahertz scattering from granular material,” J. Opt. Soc. Am. B 24(9), 2238–2243 (2007). [CrossRef]
44. O. B. Osman and M. H. Arbab, “Mitigating the effects of granular scattering using cepstrum analysis in terahertz time-domain spectral imaging,” PLoS One 14(5), e0216952 (2019). [CrossRef]
45. M. E. Khani, O. B. Osman, and M. H. Arbab, “Diffuse terahertz spectroscopy in turbid media using a wavelet-based bimodality spectral analysis,” Sci. Rep. 11(1), 22804 (2021). [CrossRef]
46. M. E. Khani, Z. B. Harris, M. Liu, et al., “Multiresolution spectrally-encoded terahertz reflection imaging through a highly diffusive cloak,” Opt. Express 30(18), 31550–31566 (2022). [CrossRef]
47. A. Kucheryavenko, I. Dolganova, A. Zhokhov, et al., “Terahertz-wave scattering in tissues: Examining the limits of the applicability of effective-medium theory,” Phys. Rev. Appl. 20(5), 054050 (2023). [CrossRef]
48. C. S. Joseph, R. Patel, V. A. Neel, et al., “Imaging of ex vivo nonmelanoma skin cancers in the optical and terahertz spectral regions optical and terahertz skin cancers imaging,” J. Biophotonics 7(5), 295–303 (2014). [CrossRef]
49. P. Doradla, K. Alavi, C. S. Joseph, et al., “Detection of colon cancer by continuous-wave terahertz polarization imaging technique,” J. Biomed. Opt. 18(9), 090504 (2013). [CrossRef]
50. S. P. Morgan and I. M. Stockford, “Surface-reflection elimination in polarization imaging of superficial tissue,” Opt. Lett. 28(2), 114–116 (2003). [CrossRef]
51. M. Ney and I. Abdulhalim, “Ultrahigh polarimetric image contrast enhancement for skin cancer diagnosis using inn plasmonic nanoparticles in the terahertz range,” J. Biomed. Opt. 20(12), 125007 (2015). [CrossRef]
52. Z. B. Harris and M. H. Arbab, “Terahertz phasr scanner with 2 khz, 100 ps time-domain trace acquisition rate and an extended field-of-view based on a heliostat design,” IEEE Trans. Terahertz Sci. Technol. 12(6), 619–632 (2022). [CrossRef]
53. K. Xu, Z. B. Harris, and M. H. Arbab, “Polarimetric imaging of back-scattered terahertz speckle fields using a portable scanner,” Opt. Express 31(7), 11308–11319 (2023). [CrossRef]
54. O. B. Osman, Z. B. Harris, J. W. Zhou, et al., “In vivo assessment and monitoring of burn wounds using a handheld terahertz hyperspectral scanner,” Adv. Photonics Res. 3(5), 2100095 (2022). [CrossRef]
55. M. E. Khani, Z. B. Harris, O. B. Osman, et al., “Supervised machine learning for automatic classification of in vivo scald and contact burn injuries using the terahertz portable handheld spectral reflection (PHASR) scanner,” Sci. Rep. 12(1), 5096 (2022). [CrossRef]
56. M. E. Khani, O. B. Osman, Z. B. Harris, et al., “Accurate and early prediction of the wound healing outcome of burn injuries using the wavelet Shannon entropy of terahertz time-domain waveforms,” J. Biomed. Opt. 27(11), 116001 (2022). [CrossRef]
57. O. Smolyanskaya, N. Chernomyrdin, A. Konovko, et al., “Terahertz biophotonics as a tool for studies of dielectric and spectral properties of biological tissues and liquids,” Prog. Quantum Electron. 62, 1–77 (2018). [CrossRef]
58. L. Wang, S. L. Jacques, and L. Zheng, “MCML—Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Programs Biomed. 47(2), 131–146 (1995). [CrossRef]
59. S. A. Prahl and S. L. Jacques, “Monte Carlo light scattering programs,” https://omlc.org/software/mc/. Accessed on 8 December 2022.
60. V. Periyasamy and M. Pramanik, “Advances in Monte Carlo simulation for light propagation in tissue,” IEEE Rev. Biomed. Eng. 10, 122–135 (2017). [CrossRef]
61. E. Pickwell-MacPherson, “Practical considerations for in vivo thz imaging,” Terahertz Sci. Technol. 3(4), 163–171 (2010). [CrossRef]
62. S. R. Tripathi, E. Miyata, P. B. Ishai, et al., “Morphology of human sweat ducts observed by optical coherence tomography and their frequency of resonance in the terahertz frequency region,” Sci. Rep. 5(1), 9071 (2015). [CrossRef]
63. J. Wang, Q. Sun, R. I. Stantchev, et al., “In vivo terahertz imaging to evaluate scar treatment strategies: silicone gel sheeting,” Biomed. Opt. Express 10(7), 3584–3590 (2019). [CrossRef]
64. Y. Feldman, A. Puzenko, P. B. Ishai, et al., “The electromagnetic response of human skin in the millimetre and submillimetre wave range,” Phys. Med. Biol. 54(11), 3341–3363 (2009). [CrossRef]
65. I. Y. Yanina, V. V. Nikolaev, O. A. Zakharova, et al., “Measurement and modeling of the optical properties of adipose tissue in the terahertz range: Aspects of disease diagnosis,” Diagnostics 12(10), 2395 (2022). [CrossRef]
66. E. Pickwell and V. P. Wallace, “Biomedical applications of terahertz technology,” J. Phys. D: Appl. Phys. 39(17), R301–R310 (2006). [CrossRef]
67. N. Chopra, K. Yang, Q. H. Abbasi, et al., “Thz time-domain spectroscopy of human skin tissue for in-body nanonetworks,” IEEE Trans. Terahertz Sci. Technol. 6(6), 803–809 (2016). [CrossRef]
68. T. C. Grenfell and S. G. Warren, “Representation of a nonspherical ice particle by a collection of independent spheres for scattering and absorption of radiation,” J. Geophys. Res. 104(D24), 31697–31709 (1999). [CrossRef]
69. S. L. Jacques, “Tutorial on Monte Carlo simulation of photon transport in biological tissues [invited],” Biomed. Opt. Express 14(2), 559–576 (2023). [CrossRef]
70. H. C. van de Hulst, Light Scattering by Small Particles (Dover Publications, 1981).
71. Y. Zhang, B. Chen, and D. Li, “Propagation of polarized light in the biological tissue: a numerical study by polarized geometric Monte Carlo method,” Appl. Opt. 55(10), 2681–2691 (2016). [CrossRef]
72. Q. Fu and W. Sun, “Mie theory for light scattering by a spherical particle in an absorbing medium,” Appl. Opt. 40(9), 1354–1361 (2001). [CrossRef]
73. P. Yang, B.-C. Gao, W. J. Wiscombe, et al., “Inherent and apparent scattering properties of coated or uncoated spheres embedded in an absorbing host medium,” Appl. Opt. 41(15), 2740–2759 (2002). [CrossRef]
74. J. Yin and L. Pilon, “Efficiency factors and radiation characteristics of spherical scatterers in an absorbing medium,” J. Opt. Soc. Am. A 23(11), 2784–2796 (2006). [CrossRef]
75. C. Matzler, “Matlab functions for mie scattering and absorption,” Tech. Rep. No. 2002–08, Institute of Applied Physics, University of Bern (2002).
76. J. Sandby-Møller, T. Poulsen, and H. C. Wulf, “Epidermal thickness at different body sites: relationship to age, gender, pigmentation, blood content, skin type and smoking habits,” Acta Derm. Venereol. 83(6), 410–413 (2003). [CrossRef]
77. M. Naftaly and R. Miles, “Terahertz time-domain spectroscopy: A new tool for the study of glasses in the far infrared,” J. Non-Cryst. Solids 351(40-42), 3341–3346 (2005). [CrossRef]
78. P. Yang, H. Wei, G. W. Kattawar, et al., “Sensitivity of the backscattering mueller matrix to particle shape and thermodynamic phase,” Appl. Opt. 42(21), 4389–4395 (2003). [CrossRef]
79. J. C. Ramella-Roman, S. A. Prahl, and S. L. Jacques, “Three Monte Carlo programs of polarized light transport into scattering media: part ii,” Opt. Express 13(25), 10392–10405 (2005). [CrossRef]
80. C. J. R. Sheppard, A. Bendandi, A. Le Gratiet, et al., “Characterization of the mueller matrix: Purity space and reflectance imaging,” Photonics 9(2), 88 (2022). [CrossRef]
81. O. Arteaga and R. Ossikovski, “Mueller matrix analysis, decompositions, and novel quantitative approaches to processing complex polarimetric data,” in Polarized Light in Biomedical Imaging and Sensing: Clinical and Preclinical Applications, (Springer, 2023), pp. 25–60.
82. M. I. Mishchenko and J. W. Hovenier, “Depolarization of light backscattered by randomly oriented nonspherical particles,” Opt. Lett. 20(12), 1356–1358 (1995). [CrossRef]
83. S.-Y. Lu and R. A. Chipman, “Interpretation of mueller matrices based on polar decomposition,” J. Opt. Soc. Am. A 13(5), 1106–1113 (1996). [CrossRef]
84. V. Sankaran MD, J. T. Walsh Jr., and D. J. Maitland, “Comparative study of polarized light propagation in biologic tissues,” J. Biomed. Opt. 7(3), 300–306 (2002). [CrossRef]
85. A. Morel and L. Prieur, “Analysis of variations in ocean color,” Limnol. Oceanogr. 22(4), 709–722 (1977). [CrossRef]
86. A. Morel and B. Gentili, “Diffuse reflectance of oceanic waters: its dependence on sun angle as influenced by the molecular scattering contribution,” Appl. Opt. 30(30), 4427–4438 (1991). [CrossRef]
87. A. D. Code and B. A. Whitney, “Polarization from Scattering in Blobs,” Astrophys. J. 441, 400 (1995). [CrossRef]
