1 November 2010 Photon-cell interactive Monte Carlo model based on the geometric optics theory for photon migration in blood by incorporating both extra- and intracellular pathways
Author Affiliations +
J. of Biomedical Optics, 15(6), 065001 (2010). doi:10.1117/1.3516722
A photon-cell interactive Monte Carlo (pciMC) that tracks photon migration in both the extra- and intracellular spaces is developed without using macroscopic scattering phase functions and anisotropy factors, as required for the conventional Monte Carlos (MCs). The interaction of photons at the plasma-cell boundary of randomly oriented 3-D biconcave red blood cells (RBCs) is modeled using the geometric optics. The pciMC incorporates different photon velocities from the extra- to intracellular space, whereas the conventional MC treats RBCs as points in the space with a constant velocity. In comparison to the experiments, the pciMC yielded the mean errors in photon migration time of 9.8±6.8 and 11.2±8.5% for suspensions of small and large RBCs (RBCsmall, RBClarge) averaged over the optically diffusing region from 2000 to 4000 μm, while the conventional random walk Monte Carlo simulation gave statistically higher mean errors of 19.0±5.8 ( p < 0.047) and 21.7±19.1% (p < 0.055), respectively. The gradients of optical density in the diffusing region yielded statistically insignificant differences between the pciMC and experiments with the mean errors between them being 1.4 and 0.9% in RBCsmall and RBClarger, respectively. The pciMC based on the geometric optics can be used to accurately predict photon migration in the optically diffusing, turbid medium.
Sakota and Takatani: Photon-cell interactive Monte Carlo model based on the geometric optics theory for photon migration in blood by incorporating both extra- and intracellular pathways



For analyzing optical propagation through scattering and absorbing turbid media, the Monte Carlo (MC) method has become a powerful technique to track photons as they migrate through a complex medium. The entire course of photon trajectories from the entry into to exit from the medium can be tracked in the computer and can even solve for the parameters of the medium. Apparently, the most famous early use of MC goes back 1930 when Fermi introduced the MC to calculate the properties of newly discovered neutrons.1 Since the 1950s, MC has been actively used for various applications relating to development of the nuclear/hydrogen bomb, atmospheric optics, geophysics. In the 1980s, the progress in computer technology helped model more complex, randomly inhomogeneous systems including biological media.2, 3, 4, 5, 6, 7 Among them, the MCML (MC maximum likelihood) developed by Wang 7 simulated migration process of photons in multilayered tissues that comprised of homogeneous layers with different optical properties, where the internal reflection at the plane boundaries was expressed using geometric optics. Similarly, Churmakov 8 in 2002 verified the importance of refractive index matching at the boundary to model diffuse reflectance by splitting the photon packets into reflection and refraction at the boundary based on the geometric optics. For randomly inhomogeneous polydisperse scattering media, such as sprays, aerosols, smoke, fog, and foams, Berrocal and Meglinski9 in 2005 characterized the medium by different extinction coefficients and phase functions. In a low-hematocrit and thin layered blood medium, random walk Monte Carlo simulation (RWMCS), introduced by Chicea and Turcu,10, 11 accurately represented the multiple light scattering anisotropies by adjusting the empirical scattering phase function.

To derive the parameters of the blood medium such as blood flow, hematocrit, and oxygenation, Roggan 12 utilized the inverse MC technique combined with the double-integrating sphere measurements to study properties of flowing blood including hematocrit, shear rate, osmolarity, hemolysis, and oxygenation that affect optical propagation in the visible and near-IR wavelength regions. Other studies include the effect of the scattering phase function on the measured optical properties of red blood cells (RBCs) by Yaroslavsky, 13, 14, 15 noninvasive diagnostics of organ circulation in brain16 and in subcutaneous tissues.17, 18

We are interested in studying the optical properties of flowing blood in the extracorporeal circulation system to noinvasively extract vital information such as hematocrit and hemoglobin content. The optical properties of flowing blood vary depending on the number, the shape, the volume, the intracellular hemoglobin content, the distribution, the orientation and the aggregation of RBCs in the circuit.19, 20 The use of time-resolved analysis using picosecond or femtosecond ultrashort pulses21, 22, 23, 24 can be combined with MC to quantify the propagation time through the blood medium by solving the inverse problem. The MC codes that have been developed for blood medium by the previous workers, however, treat RBCs as points in the space represented by the macroscopic optical properties such as scattering and absorption constants [TeX:] $\mu _s$ μs and [TeX:] $\mu _a$ μa ; and an anisotropy factor g(− 1 ≦ g ≦ 1); and a scattering phase function derived using the Mie, Rayleigh-Gans, or straight-ray theory.25, 26, 27 Hence, the simulation using the conventional MC code does not represent actual interaction between photons and cells in the space. In reality, some of the photons injected into blood medium propagate inside the cells, and some are reflected at the boundary. The velocity of light differs from intracellular to extracellular space due to the difference in refractive indices. The approach by Berrocal and Meglinski9 is promising to model inhomogeneous distribution of RBCs with varying sizes in space, but it requires various RBC characteristics already mentioned connected to optical parameters including scattering phase functions.

To overcome the shortcomings of the conventional MC code and to obtain quantitative measures of interaction between photons and blood cells, in this paper we propose a photon-cell interactive Monte Carlo (pciMC) code that tracks photon migration in both the extra- and intracellular spaces without requiring macroscopic parameters such as scattering phase function and blood cell anisotropy.28 The interaction of photons at the plasma-cell boundary of a randomly oriented three-dimensional (3-D) biconcave RBC is modeled using the geometric optics. A single photon-cell interaction model of Lugovtsov, 29 who applied the Snell's and Fresnel's laws to model interactions at the boundary of an optically soft spheroidal particle, and the MC code of Churmakov, 8 who compiled the internal reflection profile at the plane boundary, are utilized to model interactions between photons and biconcave RBCs. We developed MC code pciMC without requiring macroscopic scattering phase function and the anisotropy factor g to compute photon migration time and optical density for a given blood medium by accounting for both the intra- and extracellular migration pathways.

This paper examines the validity and utility of the pciMC for quantifying optical density and migration time of photons through a turbid medium such as blood without requiring empirical anisotropy value g and the scattering phase function, as required for conventional MC simulation. The accuracy of pciMC in predicting optical density and migration time is discussed in comparison with the conventional RWMCS.


pciMC Model

Figure 1 shows a flow chart of the pciMC program. The tracking of photons starts by defining a step size, followed with extraction of scattering information stored in the look-up tables. Inside the broken line in Fig. 1 are the look-up tables that contain the precalculated escape probability data of photons from the RBC based on the geometric optics. This look-up table approach was introduced to reduce the calculation cost of pciMC code, as reported by Sakota and Takatani28 in 2010. When the preset number of photons have been tracked, the program returns with the extra- and intracellular propagation time, and the net sum of photon weights reflecting the optical density (OD) through the medium.

Fig. 1

Flow chart of the pciMC algorithm.



RBC Model and Incident Point

We assumed that the incident light is a collimated beam of light having the emission profile of a standard normal Gaussian distribution from a 1-mm-diam optical fiber.

As the first step of simulation, the pciMC program calculates the step size S defined by


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} S = {\rm RN}(0,\!1)\frac{1}{{\mu _c }}\,, \end{equation}\end{document} S=RN(0,1)1μc,


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \mu _c = \frac{{H(1 - H)}}{{{\rm MCV}}}\sigma _c \,, \end{equation}\end{document} μc=H(1H)MCVσc,
where RN(0, 1) is the uniform random number taking the values between 0 and 1, [TeX:] $\mu _c$ μc is the geometric interaction constant in inverse centimeters, and H is the volume fraction of RBCs called hematocrit (0 < H < 1). The mean corpuscular volume (MCV) is the RBC volume in cubic centimeters, and [TeX:] $\sigma _c$ σc is the mean geometric cross section of the RBC in square centimeters. The propagation time for the step distance is calculated assuming the extracellular refractive index to be that of plasma ( [TeX:] $m_{{\rm plasma}}$ mplasma ). Note that [TeX:] $m_{{\rm plasma}}$ mplasma is approximately equal to that of water ( [TeX:] $m_{{\rm water}}\break = 1.33$ mwater=1.33 5), so the photon velocity is derived dividing the velocity in vacuum by [TeX:] $m_{{\rm plasma}}$ mplasma .

Traveling the step size S, the photon interacts with the 3-D biconcave RBCs. The photon-RBC interaction takes place in the interaction coordinate, which is defined independently of the MC simulation coordinate. The 3-D biconcave RBC is located at the center of the interaction coordinate shown in Fig. 2. The surface coordinates of the RBC model are defined in reference to Hammer 26


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} (x,y,z) &=& (r(\theta)\,\sin \theta \,\cos \phi, \,r(\theta)\,\sin \theta \,\sin \,\phi, \,r\,(\theta)\;\cos \theta),\nonumber\\ &&0\;\leqq\; \theta\; \leqq\; \pi, \quad 0\;\leqq\; \phi \;\leqq\; 2\pi \,, \end{eqnarray}\end{document} (x,y,z)=(r(θ)sinθcosφ,r(θ)sinθsinφ,r(θ)cosθ),0θπ,0φ2π,


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} r(\theta) = 3\sin ^4 \,\theta + 0.75\,. \end{equation}\end{document} r(θ)=3sin4θ+0.75.
Equation (3) expresses the 3-D biconcave RBC surface coordinates x, y, and z, respectively. Equation 4 is the biconcave RBC surface function, where θ is the angle formed by the position vector and the z axis, and ϕ is the angle formed by the projection of the position vector and the x axis to the xy plane. The surface data are discretized with the number of total nodes becoming 100×100 by cutting the model with an incremental angle of Δθ = 1.8 deg and Δϕ = 3.6 degree, respectively. Should the unit of Eq. 4 be micrometers the volume of the RBC model would become 86 μm3 and its geometric cross section [TeX:] $\sigma _c$ σc of 29.5 μm2 using the following approximation:


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \displaystyle\sigma _c = \displaystyle\frac{{2\int_0^\pi {r(\theta)\,{\rm d}\theta + \pi \left[ {r({\pi }/{2})} \right]^2 } }}{2}\,. \end{equation}\end{document} σc=20πr(θ)dθ+πr(π/2)22.
Upon interaction of photons with RBCs, the pciMC program randomly selects the incident point [TeX:] $P_n$ Pn and the vector [TeX:] ${\bm \chi }_m$ χm , assuming random orientation of RBCs. Note that [TeX:] $P_n$ Pn and [TeX:] ${\bm \chi }_{m_{} }$ χm are the interaction coordinates stored in the look-up table. Also [TeX:] $P_n$ Pn are defined by


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} P_n &=&[ {r(\theta _n)\sin \theta _n, 0,r(\theta _n)\cos \theta _n }]\nonumber\\ n &=& 1,2,\ldots {n}_{\rm \_max},\; 0\;\leqq\; \theta_n \;\leqq\; \pi/2,\nonumber\\ &&\displaystyle\hspace*{1pc}{ n}_{\_\!\max} = \frac{\pi }{{2\Delta \theta _n }}. \end{eqnarray}\end{document} Pn=[r(θn)sinθn,0,r(θn)cosθn]n=1,2,n_max,0θnπ/2,n_max=π2Δθn.
In Eq. 6, [TeX:] $\Delta \theta _n$ Δθn is the same sampling angle of 1.8 deg used in Eq. 3, so n_max becomes 50. The incident photon patterns are required only for [TeX:] $0\; \leqq\; \theta_n\; \leqq\; \pi/2$ 0θnπ/2 because of an axial symmetry of the model.

Fig. 2

Relationships between MC coordinate and interaction coordinate.


Now, [TeX:] ${\bm \chi} _m$ χm are defined as follows

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} {\bm \chi} _m = \sin {\bm \psi} _{m1} \,\,\cos \,\,\xi _{m2} {\bm \tau} + \sin \psi _{m1} \sin \xi _{m2} {\bm \mu} + \cos \psi _{m1} {\bm \nu}, \nonumber \end{eqnarray}\end{document} χm=sinψm1cosξm2τ+sinψm1sinξm2μ+cosψm1ν,
[TeX:] $m=1,2,\ldots1,1\,{\rm m}_{\rm max},\quad \pi \,\,\leqq\,\, \psi _{m1} < 3 \pi /2,$\break$ 0 \,\, \leqq\,\, \xi_{m2} \,\,\leqq 2 \pi,$ m=1,2,1,1mmax,πψm1<3π/2,0ξm22π,


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} m_{\_\!\max} = \frac{{\pi ^2 }}{{\Delta \psi _{m1} \Delta \xi _{m2} }}. \end{equation}\end{document} m_max=π2Δψm1Δξm2.
Here, [TeX:] ${\bm \nu}$ ν is the normal unit vector perpendicular to the surface at [TeX:] $P_n$ Pn , and [TeX:] $\bm\mu$ μ and [TeX:] $\bm \tau$ τ are the tangential unit vectors, respectively, where they are orthogonal to each other: and [TeX:] $\Delta \psi _{m2}$ Δψm2 and [TeX:] $\Delta\xi _{m2}$ Δξm2 are 9 and 36 deg, respectively, yielding m _max of 100. Therefore, a total number of incident photon pattern becomes n _max ×m _max = 5000.


Photon-RBC Interaction Model Based on Geometric Optics

Upon photons interacting RBCs at the incident point defined by P n and [TeX:] ${\bm \chi} _m$ χm , the program opens one of the 5000 look-up tables (n _max × m _max = 5000) [TeX:] $Q_k (P_n, {\bm \chi} _m)$ Qk(Pn,χm) specific to the incident point that contain the precalculated probability distribution of photon's escaping from the intracellular space of RBCs expressed by Eq. 8, and returns output parameters28 ( [TeX:] $M_k$ Mk , [TeX:] ${\hbox{\bi N}}_k$ Nk , [TeX:] $D_k$ Dk ). Here, M k are the coordinates of the scattered points, [TeX:] ${\bm N}_k$ Nk are the propagation vectors at the scattered points, and D k are the propagation distances from the entry points to the k-th interaction points inside the RBC. The pciMC program selects output parameters ( [TeX:] $M_k$ Mk , [TeX:] ${\bm N}_k$ Nk , [TeX:] $D_k$ Dk ) based on the [TeX:] $Q_k (P_n, {\bm \chi} _m)$ Qk(Pn,χm) probability distribution.

At each interaction point denoted by a k value, escape probabilities of photons from RBCs as precalculated using Eq. 8 are stored in the look-up table. For example, if the table contains only k = 1, reflection occurs at the entry point and the photon does not enter the intra-cellular space giving [TeX:] $Q_1\,=\,R_1$ Q1=R1 . If k is 2 or greater, [TeX:] $Q_k$ Qk indicates the probability of photon escaping from the intracellular space at the k'th interaction point. The internal reflection profile reported by Churmakov 8 is applied to calculate [TeX:] $Q_k (k\; \geqq\; 2)$ Qk(k2) as follows:


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} Q_k (P_n, {\bm \chi_m}) = (1 - R_1)\,(1 - R_k)\prod\limits_{i = 2}^{k - 1} {R_i },\quad (k\,\geqq\, 2). \end{equation}\end{document} Qk(Pn,χm)=(1R1)(1Rk)i=2k1Ri,(k2).

The probability of reflection at the k'th interaction point [TeX:] $R_k$ Rk depends on the incident angle, refraction angle and the relative refractive index between the RBC ( [TeX:] $m_{{\rm RBC}}$ mRBC ) and [TeX:] $m_{{\rm plasma}}$ mplasma , where [TeX:] $m_{{\rm RBC}}$ mRBC is expressed by:12, 28, 29


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} m_{{\rm RBC}} = m_{{\rm water}} + \beta {\rm MCHC}. \end{equation}\end{document} mRBC=mwater+βMCHC.
In Eq. 9, β is the empirically determined quantity equal to 0.001942, and MCHC is the mean corpuscular hemoglobin concentration in grams per deciliter. The refraction angle and reflection probability were calculated using Eq. 10 (Snell's law) and Eq. 11 (Fresnel's probability), respectively, where Lugovtsov applied the same method for analysis of single interaction between a photon and a spheroidal RBC model7, 30


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation}n_i \sin \alpha _i = n_t \sin \alpha _t, \end{equation}\end{document} nisinαi=ntsinαt,


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} R(\alpha _i) = \frac{1}{2}\left[ {\frac{{\sin ^2\; (\alpha _i - \alpha _t)}}{{\sin ^2\; (\alpha _i + \alpha _t)}} + \frac{{\tan ^2\; (\alpha _i - \alpha _t)}}{{\tan ^2\; (\alpha _i + \alpha _t)}}} \right] \end{equation}\end{document} R(αi)=12sin2(αiαt)sin2(αi+αt)+tan2(αiαt)tan2(αi+αt)
In Eq. 10, [TeX:] $\alpha _i$ αi is the incident angle, and [TeX:] $\alpha _t$ αt is the refracted angle. If [TeX:] $\alpha _i$ αi is larger than the critical angle (possible only when the refractive index of incident side [TeX:] $n_i$ ni is lager than that of the refracted side [TeX:] $n_t$ nt ), the internal reflection occurs, and [TeX:] $R(\alpha _i)$ R(αi) is assigned the value of 1.0. Equation 11 is the average of the reflectance for the two orthogonal polarization directions.

For each photon-RBC interaction, the precalculation was pursued until the quantity [TeX:] $Q_k$ Qk became less than 0.001% or the number of interaction reached the maximum limit set by k_max < 100. The result indicated that out of 5000 precalculated Q k tables there were some cases when the interaction number exceeded 100 with Q k values larger than 0.001% due to occurrence of total reflection inside the RBC. Such a case was defined as “trapping” and the probability of photons trapped indefinitely inside the RBC was calculated as


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} Q_{{\rm trap}} (P_n, {\bm \chi} _m) = 1 - \sum\limits_{k = 1}^{k\_\max } {Q_k (P_n, { {\bm \chi}} _m)} . \end{equation}\end{document} Qtrap(Pn,χm)=1k=1k_maxQk(Pn,χm).
Should the trapping event occur, photons do not exit the RBC. The photon's weight completely diminishes due to absorption by the intracellular hemoglobin, and hence the program terminates such a photon. The other situation is “reincidence,” where the photons hit the same RBC again immediately after their exit because of the biconcave structure of RBCs. This is defined as the “reincidence” and ignored.

In deriving scattered point coordinates M k, the following four conditions were assured for optimizing the accuracy and computation time. The condition 1 (COND1) required that the vector at M k point was in the same direction as that at k-1 point. The condition 2 (COND2) ensured that the photon vector at the k-1 point made an acute angle with the vector joining M k-1 and M k points. The third condition (COND3) was that the line joining the two points M k-1 and M k remained inside the RBC. The occurrence of such a situation was prevented by the COND3. Finally, the condition 4 (COND4) stated that the photon vector at M k-1 was expected to be parallel with the vector joining two points M k-1 and M k with the dot-product between them to be close to 1.0. The mathematical expressions of COND1 through COND4 are provided in the the appendix.


Photon Weight Loss Inside RBC and Propagation Time

Using the propagation distance D k inside the RBC, the photon weight loss was estimated by the Beer-Lambert law:


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \frac{{P_s }}{{P_i }} = \exp (- 2.303 \times e(\lambda) \times (10 \times {\rm MCHC})/m{\rm Hb} \times D_k), \end{equation}\end{document} PsPi=exp(2.303×e(λ)×(10×MCHC)/mHb×Dk),
where [TeX:] $P_i$ Pi and [TeX:] $P_s$ Ps are the incident and scattered photon weights, respectively; e(λ) is the molecular extinction coefficient for the wavelength λ in cm−1/mol/l, and mHb is the molecular weight of hemoglobin equal to 64,500 g/mole. We used the molecular extinction data of human hemoglobin obtained by the rigorous investigation of Horecker.31 For the wavelength of 651 nm used in the experiment, e(651 nm) is equal to 362.4 cm−1/mol/ l.

The intracellular migration time T RBC was calculated using [TeX:] $D_k$ Dk and the intracellular photon velocity based on [TeX:] $m_{{\rm RBC}}$ mRBC as follows:


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \begin{array}{c} \displaystyle T_{\rm RBC} = \frac{{D_k }}{c}, \\[7pt] \displaystyle c = \frac{{c_0 }}{{m_{{\rm RBC}} }}, \\ \end{array} \end{eqnarray}\end{document} TRBC=Dkc,c=c0mRBC,
where c is the photon velocity within RBC, and [TeX:] $c_0$ c0 is the photon velocity within vacuum.


Translation from the Photon-RBC Interaction Coordinate to the pciMC Simulation Coordinate

After the pciMC selects the scattered data [TeX:] $M_k$ Mk , [TeX:] ${\bm N}\!_k$ Nk , and [TeX:] $D_k$ Dk based on the [TeX:] $Q_k$ Qk probability distribution, the photon moves from point [TeX:] $P_n$ Pn to point [TeX:] $M_k$ Mk by traveling the intracellular distance of [TeX:] $D_k$ Dk and the propagation vector is switched from [TeX:] ${\bm \chi} _m$ χm to [TeX:] ${\bm N}\!_k$ Nk . Since these scattered data are in the photon-RBC interaction coordinate shown in Fig. 2, the data must be converted to those of the MC simulation coordinate. These two coordinates are interrelated by the incident point and the vector to satisfy the conditions in all the cases so that when the interaction coordinate is rotated around the incident vector [TeX:] ${\bm \chi} _m$ χm , the incident vector in the photon-RBC interaction coordinate fits the simulation coordinate. Therefore, we must decide this rotational angle Φ. In this study, Φ is randomly selected as follows.


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \Phi = 2\pi *{\rm RN}(0,\!1) \end{equation}\end{document} Φ=2π*RN(0,1)
Φ is the azimuthal angle, where in pciMC, the azimuthal scattering is the same as other simulations such as MCML (Ref. 7) and RWMCS (Refs. 10 and 11). Equation 15 means that the RBC is randomly oriented around the incident vector. Although the MCML and RWMCS likewise showed the calculation method about the scattered vector, we employed the method called the “quatemion,” which is frequently used in the computer graphics. Here, simple description of coordinate translation using quatemion is provided. First, quatemion q is defined as follows.


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} q = a + ix + jy + kz &=& (a;x,\!y,\!z) =(a;{\bm v}), \nonumber\\ i^2 &=& j^2 = k^2 = ijk = - 1. \end{eqnarray}\end{document} q=a+ix+jy+kz=(a;x,y,z)=(a;v),i2=j2=k2=ijk=1.
Here, v is the point or vector in the 3-D coordinate. Although the real number a is arbitrary, a = 0 is the most simple description. The incident point and the vector in the interaction coordinates are defined by [TeX:] $P_n = (a_1; {\bm v}_1)$ Pn=(a1;v1) and [TeX:] ${\bm \chi} _m = (a_2; {\bm v}_2)$ χm=(a2;v2) , respectively. Next, the incident point and the vector in the simulation coordinate are also defined by [TeX:] $(a_3; {\bm v}_3)$ (a3;v3) and [TeX:] $(a_4; {\bm v}_4)$ (a4;v4) , respectively. These are already known prior to the photon-RBC interaction event. The scale of the interaction coordinate is conformed to that of the simulation coordinate. We calculate the angle α formed by [TeX:] ${\bm v}_2$ v2 and [TeX:] ${\bm v}_4$ v4 . Then the description of quatemion [TeX:] $q_1$ q1 on the rotation from [TeX:] ${\bm v}_2$ v2 of [TeX:] ${\bm v}_4$ v4 is calculated as


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} q_1 = \left[ {(\cos \frac{\alpha }{2};\sin \frac{\alpha }{2}*({\bm v}_2 \times {\bm v}_4)} \right]. \end{equation}\end{document} q1=(cosα2;sinα2*(v2×v4).
If [TeX:] ${\bm v}_2$ v2 and [TeX:] ${\bm v}_4$ v4 are not linearly independent, the alternative independent vector is used in place of [TeX:] ${\bm v}_2$ v2 for Eq. 17. The description of quatemion [TeX:] $q_2$ q2 for random rotation about the axis of [TeX:] ${\bm v}_2$ v2 is also calculated as


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} q_2 = \left(\cos \frac{\Phi }{2};\;\sin \frac{\Phi }{2}*{\bm v}_2\right), \end{equation}\end{document} q2=cosΦ2;sinΦ2*v2,
where Φ was already shown in Eq. 15. The scattered point in the interaction coordinate was [TeX:] $M_k = (a_5; \;{\bm v}_5)$ Mk=(a5;v5) . To calculate the scattered point in the simulation coordinate, [TeX:] $(a_1; \;{\bm v}_1)$ (a1;v1) was translated to the original point. Then the scattered point became [TeX:] $s_{||} = (a_5 - a_1; \; {\bm v}_5 - {\bm v}_1)$ s||=(a5a1;v5v1) . According to Eq. 18, [TeX:] $s_{||}$ s|| was rotated by [TeX:] $q_2$ q2 , and additionally rotated by Eq. 17. Likewise, [TeX:] $q_1$ q1 was executed. Finally, the scattered point in the simulation coordinate s was calculated by the point translating about [TeX:] ${\bm v}_3$ v3 . These calculations are as


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} s = q_1 q_2 s_{||} \overline q _2 \overline q _1 + {\bm v}_3, \end{equation}\end{document} s=q1q2s||q¯2q¯1+v3,
where [TeX:] $\overline q _1$ q¯1 and [TeX:] $\overline q _2$ q¯2 are the conjugations of [TeX:] $q_1$ q1 and [TeX:] $q_2$ q2 , respectively. Simultaneously, the scattered vector in the simulation coordinate was also calculated using [TeX:] ${ {\bm N}}_k = (a_6; \,{ {\bm v}}_6)$ Nk=(a6;v6) as


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} {\bm s} = q_1 q_2 {\bm N}_k \overline q _2 \overline q _1 . \end{equation}\end{document} s=q1q2Nkq¯2q¯1.


Materials and Methods

The accuracy of the pciMC model for predicting OD and migration time of photons through blood medium were validated against the measurements made using a specially designed time and space resolved optical spectrometer (TSROS) described in Sec. 3.2. As an experimental medium, fresh porcine blood obtained from a local slaughterhouse was utilized, where the blood was centrifuged to separate the RBCs into RBCs with smaller MCV and higher MCHC (RBCsmall) and larger MCV and lower MCHC (RBClarge) according to their specific gravity. Then, each RBC samples were suspended in the phosphate-buffered saline (PBS) solution to provide the experimental medium with different MCV and MCHC. The MCV and MCHC were measured using a Coulter counter whose accuracy was verified against the microcapillary centrifuge method. The measured MCV and MCHC values were used in the pciMC and RWMCS, where their predictions were compared against the experimental OD and migration time of respective RBC suspension measured by a time and space resolved optical system described in the following.


Preparation of Blood Samples

Fresh porcine blood obtained from a local slaughterhouse was anticoagulated using 100 cc of 4% sodium citrate and 10,000 units of heparin per 1 L of blood. The hemoglobin oxygen saturation was adjusted to the level greater than 98% using an oxygenator. To investigate the effects of MCV and intracellular hemoglobin concentration MCHC on the optical characteristics in terms of OD and migration time, the blood sample was first centrifuged for 15 min with 2600 G and the plasma portion in the top layer was discarded. The resultant packed RBCs with the hematocrit of around 80 to 90% were further centrifuged32 for 60 min with 15000 G and 30°C. After the second centrifugation, approximately 10% of the top layers with lower specific gravity were extracted as relatively larger RBCs (RBClarge) and the bottom 10% with higher specific gravity as smaller RBCs (RBCsmall). These cells were washed in the PBS (PBS:300 mosmol/L, pH 7.4) to remove free plasma hemoglobin and to adjust the hematocrit level to approximately 20% by suspending the RBCs in the PBS. The MCV and the MCHC of the RBCs-PBS suspension were then measured using a Coulter counter (CelltacαMEK-6358, Nihon Kohden Inc., Tokyo, Japan).

Table 1 summarized the parameters of the RBClarge and RBCsmall as measured using a Coulter counter. The accuracies of MCV and MCHC measurements were 1.91 and 6.0%, respectively, calibrated against a microcapillary hematocrit method. Geometric cross section and interaction constants [TeX:] $\sigma _c$ σc and [TeX:] $\mu _c$ μc of each suspension were derived by inserting the MCV, MCHC, and H values in Eqs. 2, 5. Table 2 listed the mean values of MCV and MCHC obtained from the data in Table 1, and the optical constants were derived using the Mie theory and used for RWMCS simulation.

Table 1

Input parameters and analysis conditions of pciMC.

pciMCMCV (fL)MCHC (g/dL)Hσc (μm2)μc (cm−1)
RBCLarge 250.6280.2120.705654.712
RBCLarge 349.
RBCSmall 164.428.40.2224.317650.052
RBCSmall 258.527.20.222.808623.804
RBCSmall 357.328.10.2122.495628.129

Table 2

Input parameters and analysis conditions of RWMCS.

Mean MCVMean MCV σs σs − μs σa μa
RWMCS(fL)(g/dL)Mean H(μm2)(μm2)(cm−1)g(μm2)(cm−1)
RBCLarge (1,2,3)51.628.80.19335.3570.0051053.8560.994870.020.746
RBCLarge (1,2,3)



The optical measurements of the RBClarge and RBCsmall suspended in PBS were performed using a TSROS. The TSROS comprises of a picosecond light pulsar PLP-02 C4725, a laser head (74 mW, 53 ps pulse width, 651 nm wavelength), a streak camera C4334, and a delay unit C1097-01, all made by Hamamatsu Photonics Inc. (Hamamatsu, Japan). The ultrashort light pulses from the laser head were first split into a reference and an excitation signal, the latter was carried in one optical glass fiber (beam diameter: PHI = 1 mm) to a collimator lens (diameter = 30 mm, 60-mm focal length, DLB-30-60PM, Sigma Koki Inc., Japan), guided to a space-resolved cuvette and then illuminated to the RBC suspension through a 0.15-mm-thick microscope cover glass. A specially designed optical cuvette enabled fine adjustment of the blood sample thickness d with the resolution of 10 μm. The immerging light from blood sample through a 1.35-mm-thick microscope cover glass was then guided to a streak camera using another optical fiber (outer diameter = 2.4 mm, detection area diameter = 1.0 mm) axially aligned to the incident fiber. Both the photons transmitted through a RBC suspension and those through the reference fiber were brought to a streak camera to analyze their intensity levels using the software HPD-TA.7.1.0 (Hamamatsu Photonics Inc., Japan).

The transmittance through a sample was expressed in terms of the intensity of the transmitted intensity s(t) normalized to the reference signal intensity. The time of flight through the sample was computed using Eq. 21. With the center of the reference signal r(t) as [TeX:] $G\_\,r(t)$ G_r(t) , and the center of the sample signal s(t) as [TeX:] $G\_\,s(t)$ G_s(t) , the mean time of flight through the sample is given by


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \hspace*{-10pt}{\rm mean} \,{\rm time}\, {\rm of}\, {\rm flight}:\,\Delta t &=& G\_\,s(t) - G\_\,r(t)\nonumber\\ &=& T_{\rm media} (d) - T_0 (d) \nonumber\\ &=& T(d), \end{eqnarray}\end{document} meantimeofflight:Δt=G_s(t)G_r(t)=Tmedia(d)T0(d)=T(d),
where [TeX:] $T_{{\rm media}} (d)$ Tmedia(d) is the propagation time through blood layers of thickness d, and [TeX:] $T_0 (d)$ T0(d) is the propagation time through air for the cuvette sample thickness d.


Biconcave versus Spherical RBC model and pciMC versus RWMCS


Biconcave versus spherical RBC model

The accuracy of the pciMC with the biconcave RBC model was first compared against the spherical RBC model. The number of nodes and the analysis conditions were kept the same between the two RBC models.


pciMC versus RWMCS

The computation time, and predictions for OD and propagation time through a blood medium derived by pciMC were compared against those by the RWMCS model. For RWMCS simulation, the model proposed by Chicea and Turcu10, 11 together with the scattering and the absorption cross sections [TeX:] $\sigma _s$ σs and [TeX:] $\sigma _a$ σa derived by Mie theory were used to compute OD and propagation time. As for the refractive index, the real part was defined by Eq. 9, and that of the imaginary part was calculated as.33


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} n_{{\rm imag}} = \frac{{\ln 10}}{{\pi m{\rm Hb}}}\lambda \,\varepsilon _\mu M\cdot{\rm MCHC}{\rm .} \end{equation}\end{document} nimag=ln10πmHbλɛμM·MCHC.
In Eq. 22, λ is the wavelength in micrometers, and [TeX:] $\varepsilon _\mu M$ ɛμM is the micromolar extinction coefficient of hemoglobin at a specific wavelength λ in square centimeters per micromole. According to Reynolds,34 the experimentally observed back scattering cross section [TeX:] $\sigma _s ^ -$ σs of human RBCs at 633 nm is approximately 20 times that of Mie theory, so the anisotropy factor g used in RWMCS was expressed as follows.


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} {\rm empirical}\, {\rm RBC}\, {\rm anisotropy}\quad g = \frac{{\sigma _s ^ + - 20\sigma _s ^ - }}{{\sigma _s ^ + + 20\sigma _s ^ - }}. \end{equation}\end{document} empiricalRBCanisotropyg=σs+20σsσs++20σs.
Here, [TeX:] $\sigma _s ^ +$ σs+ and [TeX:] $\sigma _s ^ -$ σs are the forward and backward scattering cross sections calculated by the Mie theory. Moreover, for appropriate execution of RWMCS, the anisotropy g is usually expressed as a decreasing function of the optical depth [Eqs. 24 through 26] due to multiple scattering effects in the blood medium:


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \tau = \frac{{H\sigma _s }}{{{\rm MCV}}}d, \end{equation}\end{document} τ=HσsMCVd,


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} {\rm }g_1 (\tau) = g^\tau \quad {\rm if \tau }\,\leqq\, {\rm 1}{\rm .8,} \end{equation}\end{document} g1(τ)=gτifτ1.8,


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} g_2 (\tau) = g^{G(\tau)}\; {\rm with }\;G(\tau) = \frac{{(\tau - 1)e^\tau + 1}}{{e^\tau - \tau - 1}}\;\;{\rm if\; \tau > 1}{\rm .8}{\rm .} \end{equation}\end{document} g2(τ)=gG(τ)withG(τ)=(τ1)eτ+1eττ1ifτ> 1.8.
In Eqs. 24 through 26, τ is the optical depth, and d is the penetrating depth of photons into blood. Although RWMCS is an effective model, Chicea and Turcu reported10, 11 its applicability was limited for the region of multiple scattering with the optical depth of τ< 2. The anisotropy value g becomes too low in deeper medium where optical diffusion usually takes place. Assuming that the anisotropy value converges to a constant level in the diffusion region,29 in this study Eq. 27 was introduced to limit the change in g value:


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} g_3 (\tau) = g^{G(\tau \_\max)} \qquad\quad{\rm if\, \tau > \tau \_\,max}. \end{equation}\end{document} g3(τ)=gG(τ_max)ifτ> τ_max.
As for the analysis of RWMCS, the τ_max value was selected to curve fit the simulation results of both OD and propagation time to the experiments.

The azimuthal angle required for RWMCS was calculated using Eq. 15, while the scattering angle θ was derived as θ = cos −1 μ, where μ was given by Henyey-Greenstein phase function:


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \mu = \frac{1}{{2g}}\left[ {1 + g^2 + \left({\frac{{1 - g^2 }}{{1 + 2g\,{\rm RN}(0,1) - {\rm g}}}} \right)^2 } \right]. \end{equation}\end{document} μ=12g1+g2+1g21+2gRN(0,1)g2.
To tracking photons in both pciMC and RWMCS, there was no limitation in the xy space with respect to the z direction of blood layer thickness. In the pciMC, the photon was terminated when the conditions for Eq. 12 were satisfied. The photon weight threshold of 0.001 was used to terminate tracking photons based on the detection limit of the streak camera. In simulating the experimental condition, the blood layer was sandwiched between two thin glass plates with a refractive index of 1.5. The incident point was selected by assuming a standard normal Gaussian distribution within the 1.0-mm diameter of the incident beam. Detection area was also adjusted to be a 1.0-mm diameter.

The calibration of the instrument was conducted for the blood thickness of 100 μm because the nonblood conditions were too difficult for the streak camera. Prior to analysis, we also confirmed the RBC mesh size and the number of photons tracked for convergence in pciMC for the blood thickness of d = 4000 μm, the hematocrit of 20%, the MCV of 55.1 fL, and the MCHC of 29.3 g/dL. Figure 3 shows the errors normalized to that with 106 photons with 104 nodes of the RBC model, showing both OD and the propagation time converge with 104 photons injected into the medium. In Fig. 3, the convergence was tested for the required number of nodes with 104 photons injected into the medium. The photon trajectories inside the RBC suspension were shown for 100, 1600, and 10,000 RBC nodes. The errors for both OD and propagation time converged with 104 nodes. The mesh condition was described earlier in the Sec. 2.

Fig. 3

Convergence requirement for MC simulation : (a) number of photons required and (b) number of nodes required.


The RWMCS was programmed using the MCML source code.7 In all the analysis, the simulation tracked 104 photons. In the RWMCS simulation, the photon velocity remained constant in the medium with the refractive index of 1.335.




Computation Cost of pciMC

The computation time of pciMC compared against the RWMCS (τ_max = 2) using the source code of MCML with an MCV of 55.1 fL, an MCHC of 28.8 g/dL, an H of 0.193, and a blood thickness d of 4000 μm. The pciMC computation time was approximately 1.375 times that of RWMCS (τ_max = 2).


Scattering Distribution in pciMC

Figures 4 and 4 show example patterns of photon-RBC interactions inside a RBC having a volume of 86 μm3. The dark lines in the figures are the intracellular photon paths. Figure 4 is a typical case of the intracellular photon trajectory. The incident point coordinate was (x, y, z) = (3.75 μm, 0 μm,−1.64×10−7 μm) with an incident vector of (0.710,−0.672,−0.219). The parameter [TeX:] $Q_k$ Qk was pursued until the number of interaction k became 4. The photons almost refracted through the RBC membrane for k = 2 because of the very small difference in the refractive index from the intra- to extracellular medium. In contrast, Fig. 4 shows a very rare case when [TeX:] $Q_k$ Qk reached the maximum interaction number of k_max = 100 satisfying the condition of [TeX:] $Q_{k\_\max }$ Qk_max 0.001% or lower. In Fig. 4, the incident point coordinate was (x, y, z) = (3.62766 μm, 0 μm, 0.45828 μm) with the incident vector of (−0.284947, 0.939347, 0.190871). In the case of Fig. 4, the incident angle at each interaction point became large causing repeated reflection. The indefinite intracellular reflection took place so that the photon could not escape from the intracellular space.

Fig. 4

Example of photon-RBC interactions showing photon trajectory inside the RBC model: (a) general case of photon trajectory from the number of interactions k = 1 to 4 and (b) rare case of photon trajectory of k = 100.


We confirmed that the internal reflections increased as the relative refractive index [TeX:] $m_{{\rm RBC}}$ mRBC / [TeX:] $m_{{\rm plasma}}$ mplasma increased, consequently increasing the intracellular optical path lengths.28 Concerning the propagation in the suspension of RBCsmall with an MCHC of 28.8 g/dL and RBClarge with 27.9 g/dL for a normalized MCV of 86 μm3, the intracellular path length was slightly different, 2.93 μm in RBCsmall and 2.92 μm in RBClarge, respectively. Figure 5 shows the [TeX:] $Q_k$ Qk distributions of scattering predicted by Eq. 8 as a function of scattering angle for biconcave and spherical RBC models with an MCHC of 27.9 g/dL in comparison to Henyey-Greenstein phase function. The anisotropy value was calculated as the mean cosine from


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \hbox{pciMC anisotropy } {\rm g} = \int_0^\pi {Q(\theta)\cos \theta {\rm d}\theta }. \end{equation}\end{document} pciMCanisotropyg=0πQ(θ)cosθdθ.
In Fig. 5, g = 0.9819 for the biconcave model, while g = 0.9805 of the spherical model was slightly smaller. The Henyey-Greenstein function in Fig. 5 was plotted for g = 0.9819 of the same value of biconcave model. The curves showed simillar trend with the spherical model.

Fig. 5

Photon scattering distribution as a function of scattering angle from 0 to 15 degree for biconcave and spherical RBC model and those by Henyey-Greenstein phase function.



OD and Propagation Time by pciMC, RWMCS, and Experiments

The OD and propagation time for the sample thickness varying from d = 500 to 4000 μm derived by pciMC for biconcave and spherical RBC models, RWMCS models for τ_max = 2, and experiments using porcine blood samples were summarized in Figs. 6 and 6 for RBClarge and in Figs. 7 and 7 for RBCsmall. The RWMCS predictions were optimized to best fit the experiments by selecting the appropriate τ_max values so as to yield relative measure of the pciMC validities.

Fig. 6

(a) Optical density and (b) propagation time versus blood thickness of RBCsmall suspension.


Fig. 7

(a) Optical density (b) propagation time versus blood thickness of RBClarge suspension.


Concerning the overall trend for OD, both pciMC and RWMCS showed lager errors for the blood layer thickness of 500 μm, with errors decreasing with increases in the blood layer thickness. The pciMC using spherical RBC model revealed a strong tendency of forward scattering compared to the pciMC with biconcave model, revealing faster migration time in the former than the experiments.

Figure 8 shows the intracellular photon migration time and number of RBCs interacted with photons in the RBClarge and RBCsmall suspensions. The ratio of the intracellular migration time to the total time is shown in Fig. 8. For both RBClarge and RBCsmall suspensions, the ratio of the intracellular to the total migration time decreased as the blood layer thickness increased.

Fig. 8

(a) Intracellular propagation time and the number of RBCs interacting with photons and (b) percentage of the intracellular propagation time normalized to the total propagation time predicted by the pciMC for RBC suspensions.




The newly developed photon-RBC interactive MC pciMC code provided for the first time the analysis of intra- and extracellular photon migration in the blood, where the photon-RBC interaction was modeled using geometric optics. The interaction constant referred as scattering constant in RWMCS was defined in terms of mean geometric cross section of the biconcave RBC model and its concentration. At the cell-plasma boundary, the law of geometric optics was applied to determine probability and direction of scattering. The RBC model represented by a 3-D biconcave shape having a finite volume and intra-cellular hemoglobin comprised of 104 nodes and had a total of 5000 interaction patterns with incident photons.

The advantages of the pciMC are that, first, it required neither macroscopic empirical scattering phase functions nor anisotropy values as required for the conventional MC code. Second, the pciMC enabled us to track photon migration through blood medium reflecting microscopic changes caused by alterations in RBC orientation, cell size, and cell shape. Third, the model yielded for the first time the intracellular photon migration from which optical properties of individual cells could be derived. The computation time was optimized by the use of look-up tables to 1.375 times that of the RWMCS (τ_max = 2) code. The discussion here focuses on validities, peculiarities and limitations of the pciMC model in comparison to the conventional empirical RWMCS simulation.


Microscopic Geometric Optics in Blood

We first discuss the validity of geometric optics utilized to express photon-cell interaction at the cell boundary. Concerning applicability of geometric optics, Lugovtsov 30 reported that a single RBC scattering can be expressed using geometric optics approximation, provided the size parameter value defined as 2πr/λ is larger than 35, where r is the radius of particles, and λ is the wavelength of light. In our study, the size parameter value turned out to be 23 for an MCV of 60 μm3 with a wavelength of 651 nm. Shorter wavelengths could be used to increase the size parameter value, but with higher absorption of light in blood. Although the size parameter value of 23 is marginal in comparison to 35 proposed by Lugovtsov, the pciMC with biconcave RBC model showed fairly good agreement with the experiments for blood layer thickness from 2000 to 4000 μm, where optical diffusion may be present. In thinner blood layers from 500 to 1500 μm, coherency of the incident photon direction is maintained and interference takes place between the incident and scattered photons, increasing the errors between the simulation and experiments. As photons go through multiple interactions with RBCs, the coherent components are rapidly lost, resulting in incoherent propagation. In the incoherent field in thicker blood layers, the interference is diminished. The pciMC is better suited to describe photon-cell interaction in an optically diffusing medium.


Errors Sources

Before discussing the validity of pciMC, sources of errors are discussed first. A question was raised concerning accuracy of blood parameters (Tables 1, 2) used for predictions of OD and migration time by pciMC and RWMCS. Accuracies of MCV and MCHC measured by a Coulter counter were verified against a microcapillary hematocrit method, yielding that the MCV was accurate to 1.91%, while the MCHC was accurate to 6.0% due to measurement in two quantities hemoglobin content and blood hematocrit. We believe the effects of blood parameter accuracies are far less than those of the time- and space-resolved measurement system and simulation.

The error rates of the pciMC in comparison to the experiments decreased as the blood layer thickness increased. In addition to the interference effects in the thinner blood layers, as mentioned in Sec. 5.1, the blood-layer-thickness-dependent changes in the errors resort to limitation in a streak camera utilized for measurement of both intensity level and propagation time. In the thin blood layer from 500 to 1000 μm for short propagation time, it was difficult to differentiate the differences in propagation delay due to 1 ps resolution of the streak camera, while in the thick blood layer, intensity resolution of 0.001 affected the accuracy. In the thick layer, photons that had traveled longer distances due to multiple scattering effects were not registered because of low intensity, reducing the contribution of photons having longer propagation time, hence skewing the propagation time distribution toward shorter side. In the simulation for both pciMC and RWMCS, when the photon's weight became less than 0.001 of the incident photons as assumed from the streak camera limitation, they were automatically removed from the analysis. Hence, larger errors seen for the blood layer thickness of 4000 μm for both pciMC and RWMCS simulations are due to (1) not counting the photons with their weight less than 0.001 of the incident photons and (2) detection limit of the streak camera of lower intensity light.

In addition, the mismatch in the detector aperture between the simulation and experiments widened the discrepancies. For example, by narrowing the detector aperture, propagation time became faster with the increase in the OD because smaller numbers of ballistic or sneak photons having faster propagation time with reduced optical intensity thus increased OD were captured by the optical fiber.


Validity of the pciMC Biconcave RBC Model

To evaluate the validity of pciMC, both the OD and migration time should be considered, but the OD is a dimensionless number, ratio between the detected light intensity and the reference signal level, while the migration time is an absolute number, hence the latter should be prioritized to validate the pciMC. The pciMC yielded the mean errors in the migration time of 9.8 ± 6.8% in RBCsmall and 11.2±8.5% in RBClarge averaged over the blood layer thickness from 2000 to 4000 μm, while the RWMCS gave considerably higher mean errors of 19.0±5.8% in RBCsmall and 21.7±19.1% in RBClarge, respectively. There were statistically significant differences between pciMC and RWMCS with p < 0.047 in RBCsmall for blood layer thickness from 2000 to 4000 μm, and p < 0.055 in RBClarge for blood layer thickness from 2500 to 4000 μm, respectively. These large differences between the two models may resort to the fact that the pciMC accounted for different photon velocities from outside to inside the cells, while the RWMCS treated RBCs as points in the space with the photon velocity assumed constant and estimated using the refractive index of the plasma.

As for accuracy of predicting OD, in the blood layer thickness from 2000 to 4000 μm, the RWMCS (τ_max = 2) optimized through proper selection of anisotropy factor g yielded the mean error of 6.2±2.6% in RBCsmall and 2.4±2.4% in RBClarge. Those for pciMC simulation were 14.5±3.3 and 9.1±3.4% in RBCsmall and RBClarge, respectively, higher than the RWMCS. The larger deviations from the experiments in the pciMC were carried over from 500- to 1500-μm thickness where multiple scattering and interference effects created uncertainties in simulation. Regardless of larger mean errors from 2000- to 4000-μm blood layer thickness, the slope of OD changes with respect to blood layer thickness was 0.437±0.052 OD/mm (RBCsmall) and 0.424±0.049 OD/mm (RBClarge), respectively, as compared to 0.431±0.036 OD/mm (RMCsmall) and 0.4281±0.058 OD/mm (RBClarge) obtained in the experiments. The RWMCS also resulted in good agreement with 0.46 OD/mm (RBCsmall) and 0.42 OD/mm (RBClarge). There were no statistically significant differences between pciMC and experiments with the mean deviations between them of 1.4 and 0.9% in RBCsmall and RBClarger, respectively.

In an optically diffusing turbid medium, the light loss as reflected in the OD occurs due to absorption of light. The gradient of OD with respect to sample thickness in an optically diffusing medium indicates absorption properties of the medium.35 Based on this rational, we can say that the fairly constant OD gradients expressed by the pciMC and the experiments for the blood layer thickness from 2000 to 4000 μm reflect absorption properties of the medium. Since the photon absorption in the pciMC model occurs inside the RBCs relating to intracellular propagation time and MCHC of RBCs, and since in the experiments it occurs when photons collide with RBCs, the excellent agreement in OD gradients between the pciMC and the experiments supports the validity of pciMC for predicting propagation time of photons in the turbid, optically diffusing blood medium.

The pciMC is better suited to describe photon-cell interaction in an optically diffusing medium as mentioned in Sec. 5.1. The RWMCS, on the other hand, can describe the multiple scattering anisotropies with the accuracy confirmed10, 11 for the optical depth τ< 2. In the boundary region between the multiple scattering and the optical diffusion region, from 500- to 1000-μm medium thickness, with the optical depth τ of around 66 to 132, the use of RWMCS and pciMC will result in errors, because both coherent and incoherent components of light are present.


pciMC Biconcave versus Spherical RBC Model

Concerning the cell shape, the pciMC with biconcave RBC model yielded better agreement with the experiments than the spherical RBC model. As expected, the spherical RBC model showed higher scattering probability in the forward direction around a 0-deg angle (Fig. 5). Figure 5 also showed a comparison with the Henyey-Greenstein phase function. The spherical model showed similar trend to the Henyey-Greenstein scattering phase function in the 0- to 5-deg range. Roggan12 commented that the Henyey-Greenstein scattering phase function overestimates the forward scattering by RBCs. Roggan's speculation was verified in our study, because the pciMC with sphere RBC model also estimated faster propagation time similar to Henyey-Greenstein phase function. In contrast, pciMC with a biconcave model resulted in lower forward scattering with better accuracy in photon migration.

Regarding OD predictions by pciMC for a spherical model, the OD gradients in the 2000- to 4000-μm blood layers were larger with 0.53 OD/mm (RBCsmall) and 0.59 OD/mm (RBClarge), respectively, than the experiments. Hence the utility of pciMC with sphere model would be lower than that of biconcave model.


Intracellular Photon Migration Time

The pciMC simulation provided the intracellular photon migration based on the geometric optics. This analysis is the estimation that cannot be verified experimentally using the current or other setups. The intracellular photon migration time increased almost exponentially as the blood layer thickness increased [Fig. 8], while the percentage of intracellular migration time to the total propagation time decreased with the blood layer thickness [Fig. 8]. The numbers of photon-RBC interaction increased as the blood layer thickness increased. As the photon-RBC interaction increases, the scattering angle from the main axis increases with the multiply scattered photons going out of the detector site. The reduction in the percentage of intracellular migration time can be speculated such that as the blood layer thickness increases, the detection percentage of photons that have interacted less numbers of RBCs propagating in the forward small angles increase. Since those photons that have interacted with fewer RBCs have a shorter intracellular migration time, the percentage of intracellular migration time decreases [Fig. 8]. Also, the cutoff threshold of 0.001 used in the simulation and the detection capability of the streak camera contributed to the reduced percentage of the intracellular migration time for photons that have traveled a longer distance with a weight of less than 0.001 were eliminated from the analysis.


Usage of pciMC Model

Now that a pciMC model has been developed, we can propose very interesting analysis of microscopic interaction between photons and RBCs. We have been engaged in noninvasive and real-time diagnosis of blood in terms of hematocrit, hemoglobin content, and oxygen saturation during extracorporeal circulation support. To achieve our goal, the microscopic light scattering and absorption model of RBCs under blood flow must be developed. Since the conventional MC simulation expresses the RBC as a point in space using a macroscopic scattering expression such as an empirical anisotropy value and a phase function, it cannot quantify scattering changes as a function of RBC shape, orientation, and distribution in the flow. We believe the pciMC can provide direct understanding of these effects. Also, the computational fluid dynamics (CFD) for the blood flow within the cardiovascular devices such as rotary blood pumps are required for the optimal development of devices.36 By connecting CFD and pciMC, the relationship between the blood stream and the RBC distribution may be understood. Thus, pciMC may contribute to blood rheology and the development of cardiovascular devices. In diagnosing and treating various diseases, the pciMC may be useful to differentiate normal and abnormal blood cells through quantification of the optical properties such as transmission, reflection, and intracellular propagation time in individual RBCs without relying on any empirical parameters. The newly developed pciMC will contribute to early diagnosis and treatment of the diseases through quantification of optical parameters of blood cells.



The pciMC model treating optical interactions at the plasma-cell boundary based on the geometric optics quantified for the first time the OD and migration time of photons by tracking photons both extra- and intracellular spaces. The pciMC did not require empirical quantities such as anisotropy factors and scattering phase functions as required for the conventional MCs. The pciMC based on the geometric optics can be best utilized to quantify photon migration in the optically diffusing medium such as blood whose feature may be useful in diagnosing status of blood cells.


Appendix: Computation of Intracellular Photon Pathways

[TeX:] $M_1$ M1 = [TeX:] $P_n$ Pn ,

[TeX:] $M_k (k\,\geqq\,2)$ Mk(k2) ,

was decided by

  • 1 COND1: [TeX:] ${\bm U}_{k - 1} \cdot {\bm A}$ Uk1·A > 0,

  • 2 COND2: [TeX:] $ {\bm U}_{k - 1} \cdot {\bm B} > 0$ Uk1·B> 0 ,

  • 3 COND3a: a or b was satisfied.

    • a [TeX:] $(M_{(k - 1)x}\, \geqq\, 0 \cap M_{kx}\, \geqq\, 0) \cup [M_{(k - 1)x} \,\leqq\, 0 \cap M_{kx} \,\leqq\, 0]$ (M(k1)x0Mkx0)[M(k1)x0Mkx0]

    • b On the line equation connecting [TeX:] $M_{(k - 1)x}$ M(k1)x and [TeX:] $M_{(k - 1)z}$ M(k1)z , when x = 0, the value |z| < r(0). Hence,

    [TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray*} \hspace*{-29pt}| z | &=& \left| {\frac{{ \!-\! M_{(k - 1)x} [M_{kz} \!-\! M_{(k - 1)z} ] \!+\!\! M_{(k - 1)z} [M_{kx} \!-\! M_{(k - 1)x} ]}}{{M_{kx} \!\!-\! M_{(k - 1)x} }}} \right|\quad \nonumber\\ &&\times \,\,r(0). \end{eqnarray*}\end{document} |z|=M(k1)x[MkzM(k1)z]+M(k1)z[MkxM(k1)x]MkxM(k1)x×r(0).

  • 4 COND3b: In COND3a, [TeX:] $M_{(k - 1)x}$ M(k1)x and [TeX:] $M_{kx}$ Mkx are replaced by [TeX:] $M_{(k - 1)y}$ M(k1)y and [TeX:] $M_{ky}$ Mky , respectively.

  • 5 COND4: [TeX:] $ {\bm U}_{k - 1} \cdot {{ {\bm B}}}/{{| {\bm B}|}}$ Uk1·B/|B| should be close to 1.0

In the preceding equation, [TeX:] $M_{k - 1} = [M_{(k - 1)x}, M_{(k - 1)y},\break M_{(k - 1)z} ]$ Mk1=[M(k1)x,M(k1)y,M(k1)z] is the scattered point on the (k-1)' th interaction, and [TeX:] ${\bm U}_{k - 1} = [U_{(k - 1)x}, U_{(k - 1)y}, U_{(k - 1)z} ]$ Uk1=[U(k1)x,U(k1)y,U(k1)z] is the vector for photons traveling inside RBC on (k-1)' th interaction. Here A is the normal vector at the scanning points, and [TeX:] $\bm B$ B is the vector joining from point [TeX:] $M_{k - 1}$ Mk1 to the current scanning point [TeX:] $M_k$ Mk . COND3 is required to prevent the B vector from cutting across the biconcave section and reentering the other side. COND4 ensures that [TeX:] ${\bm B}$ B is at parallel with [TeX:] ${\bm U}_{k - 1}$ Uk1 .

At the same time of deriving point [TeX:] $M_k$ Mk , the scattered vector [TeX:] ${\bf N}_k$ Nk is also calculated; [TeX:] ${\bf N}_1$ N1 and [TeX:] ${\bf N}_k (k\,\geqq\,2)$ Nk(k2) are based on the reflection law at the incident point and Snell's law at the scattered point, respectively.


This work was supported by Grants-in-Aid under the Research Fellowships from the Japan Society for the Promotion of Science for Young Scientists (No. 19.10735, recipient Daisuke Sakota, April 2007 to the present) and the Basic Research B (No. 18300149, Prof. Setsuo Takatani PI).


1.  N. Metropolis, “The beginning of the Monte Carlo methods,” Los Alamos Sci. (Special Issue), 125–130 (1987). Google Scholar

2.  R. A. Forester and T. N. K. Godfrey, “MCNP—a general Monte Carlo code for neutron and photon transport,” in Methods and Applications in Neutronics, Photonics and Statistical Physics, R. Alcouffe, R. Dautray, A. Forster, G. Ledanois, and B. Mercier, Eds., pp. 33–47, Springer-Verlag, New York (1983). Google Scholar

3.  B. C. Wilson and G. Adam, “A Monte Carlo model for the absorption and flux distributions of light in tissue,” Med. Phys. 10(6), 824–830 (1983). 10.1118/1.595361 Google Scholar

4.  S. A. Prahl, M. Keijzer, S. L. Jacques, and A. J. Welch, “A Monte Carlo model of light propagation in tissue,"Proc. SPIE IS(5), 102–111 (1989). Google Scholar

5.  M. Keijzer, S. L. Jacques, S. A. Prahl, and A. J. Welch, “Light distributions in artery tissue: Monte Carlo simulations for finite-diameter laser beams,” Lasers Surg. Med. 9(2), 148–154 (1989). 10.1002/lsm.1900090210 Google Scholar

6.  I. V. Yaroslavsky and V. V. Tuchin, “Light transport in multilayerd scattering media. Monte Carlo modeling,” Opt. Spectrosc. 72, 934–939 (1992). Google Scholar

7.  L. Wang, S. L. Jacques, and L. Zheng, “MCML—Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Meth. Progr. Biomed. 47, 131–146 (1995). 10.1016/0169-2607(95)01640-F Google Scholar

8.  D. Y. Churmakov, I. V. Meglinski, and D. A. Greenhalgh, “Influence of refractive index matching on the photon diffuse reflectance,” Phy. Med. Biol. 47, 4271–4285 (2002). 10.1088/0031-9155/47/23/312 Google Scholar

9.  E. Berrocal and I. Meglinski, “New model for light propagation in highly inhomogeneous polydisperse turbid media with applications in spray diagnostics,” Opt. Express 13(23), 9181–9195 (2005). 10.1364/OPEX.13.009181 Google Scholar

10.  D. Chicea and I. Turcu, “Testing a new multiple light scattering phase function using RWMCS,” J. Optoelectron. Adv. Mater. 8(4), 1516–1519 (2006). Google Scholar

11.  D. Chicea and I. Turcu, “RWMCS—an alternative random walk Monte Carlo code to simulate light scattering in biological suspensions,"Optik 118, 232–236 (2007). Google Scholar

12.  A. Roggan, M. Friebel, K. Dorschel, A. Hahn, and G. Muller, “Optical properties of circulating human blood in the wavelength range 400-2500 nm,” J. Biomed. Opt. 4(1), 36–46 (1999). 10.1117/1.429919 Google Scholar

13.  A. N. Yaroslavsky, I. V. Yaroslavsky, T. Goldbach, and H.-J. Schwarzmaier, “The optical properties of blood in the near infrared spectral range,” Proc. SPIE. 2678, 314–324 (1996). Google Scholar

14.  A. N. Yaroslavsky, I. V. Yaroslavsky, T. Goldbach, and H.-J. Schwarzmaier, “Influence of the scattering phase function approximation on the optical properties of blood determined from the integrating sphere measurements,” J. Biomed. Opt. 4(1), 47–53 (1999). 10.1117/1.429920 Google Scholar

15.  A. N. Yaroslavsky, I. V. Yaroslavsky, T. Goldbach, and H.-J. Schwarzmaier, “Different phase function approximations to determine optical properties of blood: a comparison,"Proc. SPIE 2982, 324–330 (1997). Google Scholar

16.  E. Okada, M. Firbank, M. Schweiger, S. R. Arridge, M. Cope, and D. T. Delpy, “Theoretical and experimental investigation of near-infrared light propagation in a model of the adult head,” Appl. Opt. 36(1), 21–31 (1997). 10.1364/AO.36.000021 Google Scholar

17.  C. R. Simpson, M. Kohl, M. Essenpreis, and M. Cope, “Near-infrared optical properties of ex vivo human skin and subcutaneous tissues measured using Monte Carlo inversion technique,” Phys. Med. Biol. 43, 2465–2478 (1998). 10.1088/0031-9155/43/9/003 Google Scholar

18.  I. V. Meglinski and S. J. Matcher, “Quantitative assessment of skin layers absorption and skin reflectance spectra simulation in the visble and near-inrared spectral regions," Physiol. Meas. 23, 741–753 (2002). 10.1088/0967-3334/23/4/312 Google Scholar

19.  L. D. Shvartsman and I. Fine, “Optical transmission of blood: effect of erythrocyte aggregation,” IEEE Trans. Biomed. Eng. 50(8), 1026–1033 (2003). 10.1109/TBME.2003.814532 Google Scholar

20.  L. D. Shvartsman and I. Fine, “RBC aggregation effects on light scattering from blood,” Proc. SPIE 4162, 120–129 (2000). Google Scholar

21.  K. Shimizu, A. Ishimaru, L. Reynolds, and A. P. Bruckner, “Backscattering of a picosecond pulse from densely distributed scatterers,” Appl. Opt. 18(20), 3484–3488 (1979). 10.1364/AO.18.003484 Google Scholar

22.  M. S. Patterson, B. Chance, and B. C. Wilson, “Time resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties,” Appl. Opt. 28(12), 2331–2336 (1989). 10.1364/AO.28.002331 Google Scholar

23.  B. C. Wilson and S. L. Jacques, “Optical reflectance and transmittance of tissues: principles and applications,” IEEE. J. Quantum Electron. 26(12), 2186–2199 (1990). 10.1109/3.64355 Google Scholar

24.  B. C. Wilson, E. M. Sevick, M. S. Patterson, and B. Chance, “Time-dependent optical spectroscopy and imaging for biomedical applications,” Proc. IEEE 80(6), 918–930 (1992). 10.1109/5.149454 Google Scholar

25.  A. G. Borovoi, E. I. Naats, and U. G. Oppel, “Scattering of light by a red blood cell," J. Biomed. Opt. 3(3), 364–372 (1998). 10.1117/1.429883 Google Scholar

26.  M. Hammer, D. Schweitzer, B. Michel, E. Thamm, and A. Kolb, “Single scattering by red blood cells,” Appl. Opt. 37(31), 7410–7418 (1998). 10.1364/AO.37.007410 Google Scholar

27.  J. M. Steinke and A. P. Shepherd, “Comparison of Mie theory and the light scattering of red blood cells,” Appl. Opt. 27(19), 4027–4033 (1988). 10.1364/AO.27.004027 Google Scholar

28.  D. Sakota and S. Takatani, “Photon cell interactive Monte Carlo (pciMC) model to describe both intracellular and extracellular optical pathways of biconcave red blood cells: phase function and albedo,” Proc. SPIE 7573, 757316 1–10 (2010). Google Scholar

29.  S. Takatani, “On the theory and development of a noninvasive tissue reflectance oximeter,” PhD Dissertation, Department of Biomedical Engineering, Case Western Reserve University, Cleveland, OH (1979). Google Scholar

30.  A. E. Lugovtsov, A. V. Priezzhev, and S. Yu. Nikitin, “Light scattering by arbitrarily oriented optically soft spheroidal particles: calculation in geometric optics approximation,” J. Quant. Spectrosc. Radiat. Transf. 106, 285–296 (2007). 10.1016/j.jqsrt.2007.01.041 Google Scholar

31.  B. L. Horecker, “The absorption spectra of hemoglobin and its derivatives in the visible and near infrared regions," J. Bio. Chem. 48, 173–183 (1943). Google Scholar

32.  M. V. Kameneva, K. O. Garrett, M. J. Watach, and H. S. Borovetz, “Red blood cell aging and risk of cardiovascular diseases,” Hemorheol. Microcirc. 18, 67–74 (1998). Google Scholar

33.  D. H. Tycko, M. H. Metz, E. A. Epstein, and A. Grinbaum, “Flow-cytometric light scattering measurement of red blood cell volume and hemoglobin concentration,” Appl. Opt. 24(9), 1355–1365 (1985). 10.1364/AO.24.001355 Google Scholar

34.  L. O. Reynolds, “Optical diffuse and transmittance from an anisotropically scattering finite blood medium,” PhD Dissertation, Dept. of Electrical Engineering, University of Washington, Seattle (1975). Google Scholar

35.  A. Ishimaru, Wave Propagation and Scattering in Random Media, IEEE/OUP Series on Electromagnetic Wave Theory, IEEE Press, New York and Oxford University Press, Oxford, Tokyo, Melbourne (1997). Google Scholar

36.  G. W. Burgreen, J. F. Antaki, Z. J. Wu, and A. J. Holmes, “Computational fluid dynamics as a development tool for rotary blood pumps,” Artif. Organs 25(5), 336–340 (2001). 10.1046/j.1525-1594.2001.025005336.x Google Scholar

Daisuke Sakota, Setsuo Takatani, "Photon-cell interactive Monte Carlo model based on the geometric optics theory for photon migration in blood by incorporating both extra- and intracellular pathways," Journal of Biomedical Optics 15(6), 065001 (1 November 2010). https://doi.org/10.1117/1.3516722

Back to Top