In this chapter, we shall briefly describe the Monte Carlo algorithm adapted for biomedical optics applications and discuss several implementations of this algorithm for modeling of glucose sensing with OCT, spatial resolved reflectometry, time domain, and frequency domain techniques.
In this chapter, we shall briefly describe the Monte Carlo algorithm adapted for biomedical optics applications and discuss several implementations of this algorithm for modeling of glucose sensing with OCT, spatial resolved reflectometry, time domain, and frequency domain techniques.
Key words: Monte Carlo simulation, optical parameters, biotissues, glucose sensing, time domain technique, frequency domain techique, optical coherent tomography.
Monte Carlo-based simulations of photon transport in living tissues have become the “gold standard” technique in biomedical optics. Indeed, Monte Carlo enables to numerically solve the nonlinear radiative transfer equation that can not be solved analytically but for simplest particular cases which do not describe real biological objects. With Monte Carlo, light transport in biological objects of high complexity can be studied in a wide range of optical parameters and experimental geometries, including different modes of object illumination and signal detection. It allows accounting for optical phenomena that can take place when light interacts with biotis-sues: absorption, scattering, fluorescence, reflection and refraction at interfaces.
As in all other types of modeling, the adequacy of results obtained with Monte Carlo simulation heavily depends upon the choice of basic parameters that govern light propagation in the tissue. In the frames of the discrete particle model of biotissues and representation of the light flux as a train of single photons or photon packets, which are typically used to implement the Monte Carlo algorithm in the biomedical optics field, the process of light propagation in tissue is simulated by a chain of successive random events of scattering, absorption, reflection and, in some cases, emission of fluorescence photons forming random photon trajectories inside the tissue. Some of the photons that do not get absorbed in the tissue may reach a photodetector (photodiode, CCD camera, etc.) after leaving the tissue and contribute to the signal.
The basic parameters that are used in Monte Carlo simulation of photon trajectories and, consequently, the output signal of the detector are: absorption and scattering coefficients of the medium, scattering particle phase function, and the particle scattering anisotropy factor defined as a mean cosine of the scattering angle. These parameters determine all the important characteristics of light propagation in the particulate media: mean free path and transport pathlength of the photon, reduced scattering coefficient of the medium. The geometric structure of the modeled object or medium and the boundary conditions are also accounted for. When the Monte Carlo algorithm is implemented, the random values of the photon free path and two scattering angles are calculated at each scattering event using random number generators thus yielding a three-dimensional trajectory for each photon. The number of photon trajectories that need to be calculated mainly depends upon the statistical precision that has to be achieved in the given numerical experiment which in its turn depends on the number of photons that reach the detector and contribute to the signal and insure the needed signal-to-noise ratio. These numbers heavily depend on the optical properties of the tissue under study at the wavelength that is used. Typically from 105to 1010 input photons are used to insure the statistical validity of the simulated signal.
Monte Carlo simulation is used to solve both the direct and inverse optical problems. In the first case the output signal formed by the detected photons is the final goal of the simulation. In the second case this signal is further used to calculate the optical parameters of the tissue by comparing it with experimentally measured signal and/or theoretical predictions. Usually this is performed as an iterative procedure. Actually the basic parameters that were mentioned above (absorption and scattering coefficients of the medium, and the particle scattering anisotropy factor) can not be measured experimentally. They are usually calculated from the physically measurable values: diffuse backscattering of light from the tissue and diffuse and collimated transmission of light through the tissue. These values can be determined numerically with the help of Monte Carlo simulation of photon trajectories.
The anisotropy factor can also be calculated from the particle scattering phase function, which can be either defined from general theoretical considerations or calculated for certain types of scattering particles (e.g., basing on the precise Mie theory in the case of spherical particles , on approximations of the Mie theory or more complicated numerical solutions of the Maxwell equations in the case of non-spherical shapes ). Another option is to measure the particle scattering phase function experimentally with a goniophotometer. This option has become feasible with the development of optical trapping technique.
The origin of light scattering in biotissues is in the refractive index mismatch between their neighbouring components (cells, nuclei, intercellular and extracellular liquids, etc.). The amount of scattering is determined by the sizes, shapes and concentrations of the components and by the level of the refractive index mismatch. Being dissolved in water, blood plasma or interstitial fluid, glucose changes their refractive indexes, the value of this change being proportional to the concentration of the glucose in solution.
The refractive index mismatch between intra- and extracellular liquids affects the cell’s phase function and scattering cross-section. The effect of glucose dissolved in a solution on the refractive index of the latter was studied quantitatively by Tuchin et al. . It was shown that the glucose-induced change in the refractive index is d n = 2.73 × 10−5 mM−1. Larin et al. report on the value d n = 2.55 × 10−5 mM−1 .
The effect of glucose concentration on the optical parameters µs, µ a , g, and n of the media was discussed in Refs. [5, 6]. The effect of glucose on the suspension of polystyrene particles was considered in Ref. . It was shown that for this solution the changes in the scattering coefficient, µs, are of the order of −0.02% mM−1 glucose concentration. The anisotropy factor, g, also changes with the increase in the glucose concentration by +0.0007% mM−1. It was shown  that an increase in the glucose concentration in the physiological range (3–30 mM) may also decrease the scattering coefficient by 0.22% mM−1 due to cell volume change. The authors considered the rabbit ventricular myocytes as an example of the scattering media, however they generalize the results to other solutions. The mentioned values for the effect of glucose on optical properties were further used by other authors, e.g., in papers [7–9].
The basic idea of the Monte Carlo method is obtaining the required solution via numerous repetition of random tests and further statistical analysis of the obtained results. When applied to solving the problem of light propagation in a scattering medium, the Monte Carlo method is based on numerous calculations of random motion trajectories based on the preset medium characteristics. The Monte Carlo method can be used for the study of the propagation of photons [10, 11], electrons , protons  and other moving particles and quasi-particles. It is widely applied for comparison of experimental or theoretical data with results of computer simulation and further interpretation of the obtained results [14–17].
The application of the Monte Carlo method is based on using the macroscopic properties of a medium that are supposed to be uniform in a certain area. The simulation does not account for the radiation energy redistribution within a single scatterer. Known algorithms allow one to consider light propagation in complex multilayer media with various optical properties , different geometries , finite size of the incident beam [19, 20] and light polarization [21–23].
Theoretical solution of the problem of light propagation in scattering media implies a solution of the radiative transfer equation (RTE) which can be written as follows :
where L(r,s) is radiance at point r to the direction s; p(r,s,s′) is the scattering phase function characterizing also spatial distribution of scatterers; dW′ is the unit solid angle to the direction s′. Eq. (3.1) is valid for the case of absence of light sources inside the medium.
The Monte Carlo method is a convenient tool for simulating the light propagation in biotissues [25–28] because the direct analytical and numerical solution of RTE is strictly limited by the complexity of boundary conditions. The main limitation of the Monte Carlo method is the requirement for significant calculation powers, but it becomes less important due to modern development of computers.
Simulation of the photon transport in single- and multilayer media with various geometries can be performed implementing the Monte Carlo algorithm. Each layer is basically described by the following optical parameters: scattering coefficient µs, absorption coefficient µ a , anisotropy factor g or scattering phase function p(s,s′), where s,s′ are directions of photon propagation before and after scattering, refractive index n, thickness and boundary shape. The refractive index of the medium surrounding the object is also taken into account.
The photon propagation in the object is described in Cartesian coordinates. The current position of the photon is determined by the coordinates (x, y, z) and its current movement direction is determined by the orienting cosines: γ x = e x · r, γ y = ey · r, γ z = e z · r, where r is a unit vector of the photon velocity, e x , e y , e z are unit vectors of the coordinate axes.
Reflection and refraction of the photon on the boundary between the media with different refractive indices are calculated according to the Fresnel law for non-polarised radiation:
where α i and α t are the incident and refraction angles of the photon, ni and nt are the corresponding refractive indices of the media. The refraction angle α t is calculated in accordance with Snell’s law
The calculation algorithm used in the simulations is schematically represented in Fig. 3.1. Let us consider an iteration of this algorithm in more details. The impinging point of a photon and its initial direction are determined by a random generator in accordance to the preset parameters of the incident beam describing the spatial and angular power distributions. The technique of generation of random numbers with a given probability distribution is described in details in . Further, basing on the optical parameters of the upper layer (the only one in the case of l-layer medium) the photon free path length is calculated using the following probability density function:
where the average free path length is defined as
The random free path is determined according to the following formula:
where ξ is a random number uniformly distributed in the interval [0,1] and generated by a random number generator program.
Further, a new photon propagation direction is calculated according to the scattering phase function:
Figure 3.1 Monte Carlo algorithm used in the simulations.
The scatterers are usually supposed to be spherically symmetrical and, consequently, the polar angle φ is supposed to be uniformly distributed in the interval [0,2π], while the azimuthal angle θ is calculated in accordance with the phase function for a single scatterer implementing the algorithm described in . The directional cosines of the velocity vector at each scattering event are modified as follows:
If the incident angle is close to normal (i.e., | γ z | > 0.99999 ) the modification of the directional cosines is calculated according to the following formulae:
After the calculation of the free path length and the directional cosines, the new coordinates are calculated as follows:
where x 0, y 0, z 0 are the initial photon coordinates. After this calculation, the processing of one scattering act is finished and all these steps are repeated to process the next scattering act.
Accounting for the absorption of radiation is performed in the following way. In order to increase the statistics of the calculation it is supposed that each photon is rather a photon packet characterized by initial weight [11, 20], the latter being decreased at each scattering event by the value of
where P0 is the current photon weight. Alternative way of absorption account is a modification of the photon weight in accordance with the Lambert-Beer law at each scattering event:
where l is a pathlength between two consequent scattering events. However the latter variant averages the absorption over space making it separated from the exact scattering events. In the simulations the former variant of absorption account is usually used.
Calculation of the new coordinates and propagation direction is repeated until the photon reaches a boundary of the layer or becomes absorbed. In the case when the photon reaches a boundary the probability of its reflection or refraction is determined by the Fresnel’s coefficients. If the photon is reflected from the boundary the calculation of the propagation parameters is repeated; otherwise these calculations are performed based on the optical properties (µ s , µ a , g or p(s, s′), n) of the other layer into which the photon penetrated. In the case when the photon leaves the considered medium or becomes absorbed, the final parameters of the photon (final position, direction, weight, etc.) are processed and stored and the tracing of a new random photon trajectory begins.
The processing of the output photon data consists in saving the information about the photon trajectory for collecting the necessary statistics and further generalization of the data. For each simulated measurement technique the specific output data set is important. The photons collected are selected according to the preset detection conditions (detector size and position, numerical aperture).
Optical coherence tomography (OCT) is a modern method of noninvasive imaging of the inner structure tissues, tissue phantoms and other strongly scattering media based on the principles of low-coherent interferometry [30, 31].
The basis of any OCT setup is a Michelson interferometer with a low-coherent light source (superluminiscent light emitting diode (SLD) or femtosecond laser) placed in one of its arms (Fig. 3.2). The studied object is situated in the sample arm, a scanning mirror moving with constant speed is placed in the reference arm, a photodetector detecting the interference signal originated from optical mixing of the light coming from the sample and the reference arm is placed in the fourth (detector) arm.
The in-depth scanning (in z-axis direction) is performed by uniform movement of the scanning mirror in the reference arm of the interferometer (so-called A-scan). The amplitude of the signal detected by the photodetector is proportional to local value of the scattering coefficient at the corresponding probing depth while the signal frequency is determined by the velocity of the mirror.
Presence of the local areas with optical properties different from the average ones in the object under study leads to an alteration in the detected signal amplitude. Hence, from the dependence of the signal amplitude on the scanning time during a scan one can obtain the in-depth distribution of scattering and absorption coefficients inside the object. In transversal direction the scanning is performed by consecutive shifting of the probing beam axis at a definite step, for example, using an electromechanical system.
Figure 3.2 Schematics of the OCT setup (BS is the beam splitter).
When simulating the OCT signal, the photon tracing inside the object under study is performed by Monte Carlo technique. The simulation technique was developed in papers [27, 28, 32–37]. For the calculation of the OCT signal only the photons fitting the given detection conditions (detector position, size and numerical aperture) are selected.
For the calculation of the OCT signal the process of interference signal formation resulting from the interaction of waves reflected from the reference mirror and backscattered from the object should be simulated. The beam splitter reflection coefficient is supposed to be 0.5 and the numbers of photons injected into the reference and sample arm of the OCT setup are supposed to be equal.
When simulating the photon propagation in the object its optical pathlength is calculated. After completing the simulations for a large number of photons the distribution of the detected photons over their optical pathlengths is obtained as the output data. This distribution can be considered as an envelope of the OCT signal with very narrow (δ-function) coherence function, because it characterizes the in-depth distribution of object optical properties [34, 38].
The OCT signal itself can be calculated in accordance with the formula for the interference part :
where Ir and Is are the amplitudes of waves scattered from reference and sample arms correspondingly; Δl is their optical pathlength mismatch. Consequently, accounting for the low coherence of the probing radiation, the interference signal can be calculated as follows:
where Nr (t, Δli ) and Ns (t, Δli ) are the numbers of photons detected from the reference and the object arm, respectively, C(Δli,l coh) is the coherence function, and l coh is the coherence length of the probing radiation. This time-dependence can be recalculated into depth-dependence given the scanning mirror velocity is known. In our simulations the coherence function is supposed to be Gaussian, and the OCT signal is calculated according to the following formula:
In Ref.  the Monte Carlo simulation of OCT signal was performed without accounting for the coherence function. In this study the OCT signal was obtained as a square root of the distribution of the detected photons over their optical path-lengths. Accounting for the photons phase shifts when propagating in a medium in the presence of moving scatterers allows one to simulate the optical coherence Doppler tomography (OCDT) signal .
However, due to the fact that the photon packets are considered in the simulations it is necessary to account for the independent interference of each detected photon packet with the reference radiation and the simulated OCT signal should be calculated as a superposition of these partial interference fringe patterns. The number of injected photons needed for the simulation of one A-scan is usually chosen in regard to the required calculation time and statistical accuracy.
Accounting for the speckle formation is an important problem for the simulation of OCT signals. The main sources of speckles in OCT are fluctuations of the position of the studied object and multiple scattering in the medium causing random phase shifts. The first type of the speckles is not taken into account in the simulation because it is assumed that the object is fixed and this type of speckles can be neglected. Moreover, modern OCT devices utilizing the fiber optics allow one to fix the optical fiber directly to the studied object avoiding this kind of speckles.
It is worth mentioning that the speckle effect is present in the simulated OCT signal because the latter is calculated as a sum of partial interference fringe patterns. The summation of the fringe patterns with random phase shifts results in a speckle structure of the signal. In the experiment this speckle structure is caused by a random position of the scatterers inside the scattering medium. To suppress the speckle effect, the spatial or temporal averaging of the OCT signal is usually implemented. In Monte Carlo simulation of the OCT signal, the calculation of the output signal as a sum of envelopes of the partial interference signals is used instead of spatial averaging. The formula for the interference signal calculation in this case will be the following:
Another problem that may appear when simulating the propagation of a broadband radiation is the dependence of the medium optical properties on the wavelength leading to dispersion effects. However, the optical properties of biotissues and their phantoms in the NIR range exhibit weak wavelength dependence and these effects can be neglected in simulations related to these tissues.
Alternative approaches to Monte Carlo simulations of OCT signals were proposed in papers [39, 40], in which the simulations were based on the extended Huygens-Fresnel principle allowing one to account accurately wave effects. In paper  a fixed-particle Monte Carlo method was considered allowing one more accurate calculation of speckle patterns.
The comparison of the experimental and simulated OCT signals from intralipid solution with various concentrations of glucose was performed in paper . The solution under study was positioned in a plain 1-mm thick cuvette. Intralipid solution is the medium widely used for mimicking the scattering properties of dermis and epidermis. The basic idea of the Monte Carlo simulations of the glucose effect on the signals obtained in various optical diagnostics techniques including OCT is the changing of the optical parameters of the media in which the light is propagated according to the relations mentioned in section 2.
The slope of the OCT signal was chosen as criterion for glucose sensing because the changes in glucose level affect the scattering properties of the medium affecting the slope of the signal. Increasing the concentration of glucose reduces the refractive index mismatch between scattering particles and the solution media and hence decreases also the scattering coefficient. As a consequence, the signal attenuation decreases and more optical power can be detected from larger depths, which results in a change of the OCT signal slope in accordance with the Lambert-Beer’s law:
where µ t = µ a + µ s and z is the depth co-ordinate in the media. The comparison of experimental and Monte Carlo results for 5% Intralipid solution are presented in Fig.3.3 with solid and empty circles, respectively.
The results for 2% Intralipid solution have larger deviations from the fitting line compared to the case of 5% solution because the total amount of the backscattered photons, contributing to the level of the signal, decreases with the decrease in concentration, and, hence, we need more accurate intensity detection in the experiment and larger statistics in the simulations.
The results for 5% Intralipid solution exhibit a 4.4% signal slope change in the experiment and about 7% change in the simulations, when the amount of glucose is varied from 0 mg/dl to 1000 mg/dl. The results for 2% Intralipid solution (see Fig. 3.4) exhibit about 13% signal slope change in the simulations, when the amount of glucose changes from 0 to 1000 mg/dl.
The difference in experimental and simulated glucose-induced changes may be caused by the assumption that the glucose effect on µ s equals to −0.22% mM−1, as stated for blood samples in refs [4–6]. However, for Intralipid solutions this effect may be different and needs additional studies.
In conclusion, it is worth saying that both the simulation and experimental results show the sensitivity of the OCT signals to glucose concentration. Although quantitatively they differ due to details not accounted in the model, qualitatively they look alike which proves that Monte Carlo simulations can be effectively used for qualitative estimation of the OCT technique in glucose sensing.
Figure 3.3 Glucose-induced changes in the OCT signal slope for 5% Intralipid TM , obtained from the experiment and the simulation; error bars indicate ±1 standard deviation .
Figure 3.4 Glucose-induced changes in the OCT signal slope for 2% Intralipid TM , obtained from the simulation.
One of the methods that can potentially be implemented for noninvasive measurements of blood glucose level is spatial resolved reflectometry (SRR) [9, 42]. The essence of this technique is in measuring the dependence of the diffuse backscattered light intensity on the source-detector separation. Such measurements are usually performed with the help of an array of detectors or a stepwise moving detector. The SRR technique was also demonstrated to be effective for blood oxygenation detection .
In this section we shall consider the numerical study of SRR potentialities for glucose level detection implementing Monte Carlo technique in the case of a three-layer biotissue phantom and discuss the choice of the measurement parameters. The schematic layout of the simulated experiment is shown in Fig. 3.5.
Figure 3.5 Schematic layout of the simulated experiment.
In our calculations, the embedding depth L 1 of the blood layer varied from 100 to 300 µ m which corresponds to the range of depths of the upper plexus in human skin. The thickness of the blood layer L 2 was fixed and equal to 200 µ m. The thickness of the lower skin layer was chosen so that the total thickness of the sample resulted in 10 mm allowing to consider the lower skin layer as semi-infinite according the values of its optical properties. The optical parameters of the skin-mimicking layers were averaged from the values reported in papers [18, 44–46], while for the blood-mimicking layer the values were chosen basing on the data reported in [14, 47]. The values used for simulations are shown in Table 3.1. The wavelength was chosen as 820 nm which belongs to the so-called diagnostic transparency window (600–1300 nm) and is widely used in the noninvasive diagnostics of tissues. In the table the transport length in both media l tr = 1/(µ a + µ ′ s ) characterizing the chaotization of the photon movement direction is also shown.
µ s , mm−1
µ a , mm−1
For the simulations, two values of the anisotropy factor characterizing skin (0.9) and water solution of Intralipid (0.7) were used. Because this parameter plays a significant role in the process of light propagation in a scattering medium, the analysis of its effect on the obtained results is of importance. The variation of the glucose level is considered in the range from 0 to 500 mg/dl, the extreme low and high values corresponding to strictly pathological cases for a human organism. In simulation, it was supposed that the change in glucose concentration does not affect the optical properties of the superficial skin layer because in human organism the glucose level changes at first in blood and then in tissues neighboring to blood vessels, while in superficial layer the changes are late and non-significant.
The number of launched photons in the simulation was chosen as 109. To increase the calculation statistics in the case of axial-symmetrical medium the dependence of the backscattered intensity on source-detector separation is calculated not along a definite direction, but in all directions which allows obtaining the result accounting for more photon trajectories.
The signals obtained for the medium mimicking the skin with g = 0.9 for embedding depth of the blood layer of 200 µ m and glucose concentration of 0 mg/dl and 500 mg/dl are shown in Fig. 3.6(a). Figure 3.6(b) shows the same signals scaled for the region where the relative signal difference is maximal.
The presented results are normalized by the detector area and the probing beam power with an assumption that the numerical aperture of the detector NA=0.24 which corresponds to the detection angle of 28◦. Thus, the magnitude I in Fig. 3.6 has the dimension of mm−2. Fresnel reflection from the detector surface is about 1% of the signal level; therefore it was not taken into account in the simulation. Fresnel refraction can be also reduced by administering an index-matching liquid onto the surface of the object (optical clearing ). The scattered light power P registered on the detector surface can be calculated from data in Fig. 3.6 according to the formula:
Figure 3.6 a) SRR-signals for the 3-layered medium at the embedding depth of blood layer L 1 = 200 µ m for glucose concentration of 0 mg/dl (solid line) and 500 mg/dl (dashed line); b) the same signals for the region of maximal relative difference (shown in Fig. 3.6a with a rectangle).
where P 0 is the power of the probing radiation, S d is the detector area, R is Fresnel reflection from the surface of medium. It is obvious that the calculation according to this formula is correct only in the case when the linear size of the detector is less than the distance over which the change of the SRR-signals is significant.
An important measurement parameter is the intensity of the probing radiation. In the case of glucose detection the limitations for this intensity are quite strict, because even an increase in the light-induced heating temperature of the probed medium by 1◦C can cause changes in the optical properties of the object leading to an error in the determination of the glucose level . According to Ref. , the maximal admissible intensity of the probing radiation for the considered problem is around 1 W/mm2. The spatial resolution in the SRR-signal calculation corresponding to the linear size of the detector (or optical fiber diameter) is 20 µ m. In this case the detector area is about 400 µ m2. At the source-detector separations ranging from 0 to 8 mm, for the maximal admissible intensity of the probing radiation the detected power varies from 10−4 to 10−8 W and the corresponding intensity varies from 10−1 to 10−5 W/mm2. At the wavelength used, the radiation power of such orders can be detected with modern photo detectors with high accuracy. To raise the power of the detected signal we can increase the detector area, as the selected value of 20 µ m is basically a parameter of simulation rather than the real fiber diameter.
For optimizing the position of the detector relative to the source, an estimation of relative sensitivity S of the SRR-signal to the glucose concentration in dependence on the source-detector separation r can be performed according to the following expression :
where I 0 and I 500 are SRR, signals corresponding to the glucose concentration of 0 and 500 mg/dl. The results obtained for three considered embedding depths of the blood layer are shown in Fig. 3.7. From this figure one can see that for g = 0.9 (Fig. 3.7a) the relative sensitivity S has a local maximum at the source-detector separations below 0.5 mm. The minimum at r above 2 mm is due to the fact that the difference I 0 − I 500 changes its sign; however the magnitude S is defined as nonnegative. Consequently, at the condition I 0 − I 500 = 0 it has a minimum.
Figure 3.7 Relative sensitivity of the SRR-signal to the change of glucose level from 0 to 500 mg/dl in the blood layer and the deeper skin layer for the embedding depths of the blood layer of 0.1, 0.2 and 0.3 mm; the anisotropy factor g = 0.9 (a) and 0.7 (b) .
One can see that as the source-detector separation increases the relative sensitivity S reaches again the value corresponding to the local maximum around r = 7 mm. At this distance from the source, the power of the detected signal decreases by three orders of magnitude in comparison to its value at r = 0.4 mm, which makes reason to perform the measurements at r = 0.4 mm. Moreover, as it will be shown further, the photons forming the backscattering signal at r = 7 mm reach the depths comparable with r in the considered medium. However, the human skin thickness is only 2– 3 mm although in this model we have considered a thicker phantom.
For the case g = 0.7 (Fig. 3.7b) the relative sensitivity S also has a local maximum around 0.5 mm. However, as one can see from this figure, in this case the dependence of the maximum position and maximum value on the blood layer embedding depths is stronger. Moreover, S for small r values is lower than in the case of g = 0.9, which is related to an increase in randomization of the photons’ directions before they reach the layers where the optical parameters are sensitive to glucose (blood layer and below). Although for the case g = 0.9 the photon transport length significantly exceeds the embedding depths of the blood layer, for g = 0.7 it is only 0.3 mm (see Table 3.1), which is comparable with the blood layer embedding depths. Thus, at the maximal value of the considered blood layer embedding depths, the regime when the photon direction is already chaotic before entrancing the glucose sensitive layers takes place. In this case, higher sensitivity to glucose is observed at larger r values, where the backscattering is quite diffuse.
For a more detailed analysis of the photon trajectories forming the SRR-signal, the scattering maps describing the density of photon scattering events in the medium under study can be considered . The scattering maps shown in Fig. 3.8 were obtained for the embedding depth of the blood layer L 1 = 200 µ m and for three different values of r: 0.4, 0.8, and 1.2 mm.
The scattering maps show the trajectories distribution which has higher density and lower width in the points of emission and detection and lower density and higher broadening width between these points. The density of the trajectories is proportional to the local brightness which is higher in the region corresponding to the blood layer because the maps show the density of scattering events of the photons forming the signal and the scattering coefficient of blood is approximately five times higher than that of skin.
For the case g = 0.9 the trajectory distributions are less broadened than for the case g = 0.7, which is caused by higher randomization of the photon directions in the latter case. With an increase in the source-detector separation the overlap of the picture parts with the highest brightness in the regions of the source and the detector is reduced.
Finally, let us analyze the dependence of the SRR-signal on glucose concentration in blood and deep skin layers for the embedding depth of blood layer L 1 = 100 µ m at 2 combinations of parameters: (a) skin anisotropy factor g = 0.9 at r = 0.4 mm; (b) skin anisotropy factor g = 0.7 at r = 0.3 mm. As it was shown before, these combinations are close to optimal. The simulation results obtained in Ref.  are shown in Fig. 3.9a and 3.9b, respectively.
It is seen that in the range of glucose concentrations from 0 to 500 mg/dl there exists a clear correlation between the detected SRR-signal and the glucose concentration in the medium. In the case of a fiber with the cross section area of 400 µ m2 and the maximal input intensity of 1 W/mm2, which were mentioned above, the detected power would change from 1.28 to 1.7 µ W at the change of the glucose concentration from 0 to 500 mg/dl. Such measurement would demand quite high accuracy of the detection. However, as noted before, the mentioned value of detecting fiber cross section area is a parameter of simulation rather than the real characteristic of the measuring system. It is obvious that with an increase in the detector area the detected power will increase too, which reduces the requirements to the sensitivity of the photodetector. For instance, for the ThorLabs detector APD210 the sensitivity is 105 V/W in the considered wavelength range, which is sufficient for the reliable measurements with the detector fiber cross section area about 0.01 mm2. In this case the detected power would vary from 32.1 to 34.7 µ W, which would lead to a change in the detected voltage from 3.21 to 3.47 V. Hence, the sensitivity to glucose changes would be around 0.52 mV/(mg/dl).
Figure 3.8 Scattering maps of photons forming the SRR-signal from a 3-layered skin model for the embedding depth of the blood layer L 1 = 200 µ m and the skin anisotropy factor g = 0.9 (a) and 0.7 (b). Distances along X and Z axis are measured in mm .
Figure 3.9 Dependence of signal intensity at a) g = 0.9 and r = 0.4 mm, b) g = 0.7 and r = 0.3 mm on glucose concentration in blood and deep skin layers at the embedding depth of blood layer L 1= 100 µ m .
The effect of glucose concentration on propagation of ultrashort laser pulses in a model light-scattering suspension was studied in Ref. . It was shown that it is possible to estimate the glucose level in the physiological range of concentrations (100–500 mg/dl) by measuring the peak intensity and total energy of the detected pulses in backscattering geometry. In this section we shall discuss the applicability and advantages of the time domain (time-of-flight) technique in the problem of noninvasive glucose sensing taking into consideration a three-layer model of the skin. Figure 3.10 schematically represents the arrangement of the simulated experiment. The geometry and the optical properties of the skin phantom as well as the parameters of the detecting fibers are similar to those described in the previous section.
Figure 3.10 Schematic layout of the simulated experiment.
The only difference is the type of a light source and the arrangement of the probe. The latter consists of 5 adjacent optical fibers placed in line on the top of a phantom as shown in Fig. 3.10. In this case, we consider the incoming (probing) signal as a δ-pulse in time domain.
The output pulses for all five fibers calculated at L1 = 100 µ m and two glucose concentrations are presented in Fig. 3.11. From this figure one can see that for all considered fibers the time profiles of the scattered pulses corresponding to different glucose concentrations slightly differ. This fact allows one to consider the time-offlight (TOF) technique as potentially applicable for glucose sensing. The essential difference between g = 0.9 and g = 0.7 cases is that in the former case the scattered pulse has one peak while in the latter case it has two peaks. This effect is caused by the fact that the amount of light backscattered from a definite layer is characterized by the value of (1 − g). In this connection, the mismatch between these values for the skin and blood layers is larger in the case g = 0.7 causing larger difference in the amounts of backscattered light and, hence, in the amplitude of the detected pulse. For g = 0.9 this difference can be seen only for the first fiber.
Once the difference in pulse shapes for different glucose concentrations is discovered one should find the proper quantitative characteristics for recovering the glucose level value from the measured pulse shape.
Figure 3.11 TOF signals for two different glucose concentrations: 0 (solid line) and 500 (dashed line) mg/dl for L1 = 100 µ m, anisotropy factor g = 0.7 (a) and 0.9 (b).
TOF measurements allow to use more parameters (e.g., total pulse energy, peak intensity, peak position) for signal characterization as compared to spatial resolved measurements and thus seem to be more informative for glucose sensing, however, more expensive. It is also possible to increase the sensitivity implementing the time gating (see below).
Let us consider the peak value of the measured pulse, the time-integral value (total energy of the pulse) and the partial time-integrals in the intervals from zero to some prefixed time (i.e., the integration time). Dependencies of the peak value and the total energy of the pulse on glucose concentration for the first fiber are presented in Fig. 3.12. One can see that both chosen values exhibit fairly linear dependence on the glucose concentration. However, we should mention that the numerical simulations are performed for an ideal case and when applying the TOF method experimentally one should take into account the measurement error, effect of possible noises and limited dynamic range of the measuring system.
Figure 3.12 TOF signals for two different glucose concentrations: 0 (solid line) and 500 (dashed line) mg/dl for L 1 = 100 µ m, anisotropy factor g = 0.7 (a) and 0.9 (b).
Now let us consider the partial pulse energies calculated under the condition that the photons are collected only within the definite time intervals. In our simulation we chose these intervals to be 5, 10, 20, 30, 40, 50 and 130 ps long. The relative sensitivities of the partial pulse energies to glucose concentration for all five fibers for L1 = 100 µ m and g = 0.9 are presented in Fig. 3.13.
From this figure one can see that the maximal sensitivity is characteristic to the second and the third detectors placed at distances corresponding to the distances of maximal sensitivity obtained with the SRR technique. The dependence of the sensitivity on the integration time differs for different fibers. However for the second and the third fibers the sensitivity decreases with an increase in the integration time. From this fact one can conclude that given the detecting system allows such integration, the shorter integration time intervals (around 5 or 10 ps long) are preferable.
Figure 3.13 Relative sensitivities of the partial pulse energies to glucose concentration change from 0 to 500 mg/dl for L1 = 100 µ m and g = 0.9 for different fibers (the numbers are shown above the curves) versus the integration time interval.
In this section we will consider the applicability and advantages of the frequency domain technique in the problem of noninvasive glucose sensing. It has close connection with the time-resolved (time domain) measurements considered in the previous section. The hardware or software Fourier-analysis of the scattered pulses allows to simultaneously obtain the amplitude and phase responses of the medium for the continuous set of harmonics.
This method is based on the registration of dynamic response of intensity of scattered light to the modulation of the incident laser beam intensity within a wide frequency modulation band (0.1 - 10 GHz). The measuring parameters in the case of frequency domain measurements are the modulation depth of scattered intensity m and the phase shift between modulation phases of incident and detected radiation Df as shown schematically in Fig. 3.14. In comparison with the time-of-flight method, the frequency domain method is simpler for implementation, easier for the interpretation of the results and has greater noise immunity.
In the frequency domain, the signal formation can be interpreted in terms of so-called photon density waves. These waves resemble the heat-waves arising due to the absorption by a medium of the modulated laser radiation. Each photon makes a random walk inside the scattering medium; however all together they form a photon density wave, spreading from the source, with a frequency equal to the modulation frequency of the source. The photon density waves have all the properties typical of the waves of other types, e.g., refraction, interference, diffraction, and dispersion . In highly scattering media with small absorbance far from the borders, light sources and detectors, the light propagation can be considered as a diffusion process described by the time-resolved diffusion equation for the photon density function Φ(r,t).
Figure 3.14 Schematic representation of the reflectance (a) and transmittance (b) modes of the frequency domain measurement technique.
For a point light source located in the point r = 0 with the harmonic modulation of the intensity
the solution of the diffusion equation for the homogeneous infinite medium can be presented in the following form :
where Φ DC (r), Φ AC (r) and ϕ (r) are the constant (direct current) component, amplitude and phase of the photon density wave, respectively, measured at a distance r between the source and detector; ω = 2πν, where ν is the frequency of modulation. The modulation depth of the detected signal is defined by the relation:
The variable component of the obtained solution represents an outgoing spherical wave with the centre in the point r = 0, which oscillates at the modulation frequency ω and has a phase shift relative to the phase value at the point r = 0. The measurements of νΦ(r,w), and ϕ(r,ω) values allow to determine the reduced scattering coefficient µ′ s = µs(1 − g), absorption coefficient µ a and the distribution of these parameters inside the medium under study. Measurements of the reduced scattering and absorption coefficients using the frequency domain technique with the accuracy high enough for glucose sensing have been reported in .
Monte Carlo simulation of scattering signals in frequency domain is based on the fact that the values of the phase shift and modulation depth of detected radiation can be obtained for the continuous modulation frequency band just with the help of Fourier-transform of the medium transfer function h(t) (the medium response to a δ-pulse) :
where m 1 and m 0 is the modulation depth of the detector and source, respectively, ϕ is the phase shift.
Figure 3.15a shows the calculated responses of semi-infinite homogeneous scattering media with scattering coefficients µ s = 1 and 5 mm−1 to a δ-pulse for the case of back reflectance measurements. Absorption coefficient was assumed to be 0, the anisotropy factor g = 0.7, the detector radius was 0.5 mm and the distance between source and detector was 6 mm.
Figure 3.15 Responses of the scattering media with different values of scattering coefficient to a δ-pulse (a); frequency dependences of phase shift between source and detector (b) and modulation depth (c) for the same values of the scattering coefficient obtained with help of Fourier-transform.
Figures 3.15b and 3.15c show the frequency dependencies of the phase shift between source and detector and the modulation depth calculated for the same values of the scattering coefficient obtained with help of Fourier-transform. One can see that a decrease in the scattering coefficient of a homogeneous semi-infinite medium leads to a decrease in the phase shift and an increase in the modulation depth of the detected signal at the certain value of modulation frequency.
Let us consider a semi-infinite layer of a medium with optical properties similar to those of Inralipid solution of 2% concentration  and a 1 mm-thick layer with optical properties close to those of blood [14, 47] embedded at a depth of 1 mm. Parameters of the detector and source-detector separation are the same as mentioned above. Fig. 3.16 shows the calculated phase response of the medium with blood glucose concentrations of 0 and 500 mg/dl.
Figure 3.16 Frequency-resolved dependences of the phase shifts between the source and detector for the blood layer with glucose concentration of 500 mg/dl (dashed line) and without glucose (solid line) embedded into the scattering medium (2% Intralipid solution).
From this figure one can conclude that even the maximal possible physiological change of blood glucose concentration causes but changes in the phase shift. This fact demands that special measures should be taken to raise the sensitivity of the measuring device. Figure 3.17 shows that with an increase in the modulation frequency the difference between the phase shifts presented in Fig. 3.16 gradually increases. In general this result is clear since with an increase in the modulation frequency the sensitivity of the frequency domain method to changes of the scattering coefficient should increase. However, the nonmonotonous character of this dependence needs additional interpretation.
It has been experimentally tested that the uncertainty of phase measurements for the majority of devices is about 0.1 deg and does not depend on the modulation frequency. According to this estimate, the detection of the maximal glucose changes is possible at the modulation frequencies higher than ∼400 MHz. However in real practice it is necessary to detect smaller changes in glucose concentration resulting in the phase shift differences lower than in the above considered case. Thus we can conclude that such measurements should be performed with higher modulation frequencies (in the gigahertz range).
Figure 3.17 Frequency-resolved dependence of difference of phase shifts for blood layer with glucose concentration of 500 mg/dl and without glucose embedded into the scattering medium (2% intralipid solution) (boxes) and its approximation (solid line).
The Monte Carlo method of stochastic numerical simulation of light propagation in strongly scattering media is a powerful technique that can be efficiently implemented for solving various biomedical problems. In particular, it can be used to predict the effect of glucose on the output signal in glucose sensing experiments implementing different measurement techniques. In this chapter, we have shown that the signals of optical coherence tomograph, as well as reflectometers performing measurements in spatial resolved, time resolved, and frequency domain modes can be simulated with the Monte Carlo method. The simulation results that we described demonstrate the sensitivity of the signals obtained in application of the considered techniques to glucose concentration changes in the physiological range. However, even in the case of rather simple biotissue models which do not account for local structural nonhomogeneities and dynamic instabilities typical for live objects, e.g., human skin, this sensitivity, as of today, is not so high that would satisfy practical implementations for in vivo measurements in human patients. Evidently, more sophisticated models of human skin and experimental measurement techniques should be designed to get further progress in this field. Monte Carlo simulations can be used on all stages of their development to evaluate the attained sensitivity to glucose.