_{small}, RBC

_{large}) 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 RBC

_{small}and RBC

_{larger}, respectively. The pciMC based on the geometric optics can be used to accurately predict photon migration in the optically diffusing, turbid medium.

## 1.

## Introduction

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 Meglinski^{9} 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 brain^{16} 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 pulses^{21, 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$
${\mu}_{s}$
and
[TeX:]
$\mu _a$
${\mu}_{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 Meglinski^{9} 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.

## 2.

##
**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 Takatani^{28} 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.

## 2.1.

### 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

## Eq. 1

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} S = {\rm RN}(0,\!1)\frac{1}{{\mu _c }}\,, \end{equation}\end{document} $$S=\mathrm{RN}(0,1)\frac{1}{{\mu}_{c}},$$## Eq. 2

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \mu _c = \frac{{H(1 - H)}}{{{\rm MCV}}}\sigma _c \,, \end{equation}\end{document} $${\mu}_{c}=\frac{H(1-H)}{\mathrm{MCV}}{\sigma}_{c}\phantom{\rule{0.16em}{0ex}},$$*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$ ${\sigma}_{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}}$ ${m}_{\mathrm{plasma}}$ ). Note that [TeX:] $m_{{\rm plasma}}$ ${m}_{\mathrm{plasma}}$ is approximately equal to that of water ( [TeX:] $m_{{\rm water}}\break = 1.33$ ${m}_{\mathrm{water}}=1.33$ 5), so the photon velocity is derived dividing the velocity in vacuum by [TeX:] $m_{{\rm plasma}}$ ${m}_{\mathrm{plasma}}$ .

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}

## Eq. 3

[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} $$\begin{array}{ccc}\hfill (x,y,z)& =& \left(r\right(\theta )\phantom{\rule{0.16em}{0ex}}\mathrm{sin}\theta \phantom{\rule{0.16em}{0ex}}\mathrm{cos}\phi ,\phantom{\rule{0.16em}{0ex}}r(\theta )\phantom{\rule{0.16em}{0ex}}\mathrm{sin}\theta \phantom{\rule{0.16em}{0ex}}\mathrm{sin}\phantom{\rule{0.16em}{0ex}}\phi ,\phantom{\rule{0.16em}{0ex}}r\phantom{\rule{0.16em}{0ex}}(\theta \left)\phantom{\rule{0.28em}{0ex}}\mathrm{cos}\theta \right),\hfill \\ & & 0\phantom{\rule{0.28em}{0ex}}\leqq \phantom{\rule{0.28em}{0ex}}\theta \phantom{\rule{0.28em}{0ex}}\leqq \phantom{\rule{0.28em}{0ex}}\pi ,\phantom{\rule{1em}{0ex}}0\phantom{\rule{0.28em}{0ex}}\leqq \phantom{\rule{0.28em}{0ex}}\phi \phantom{\rule{0.28em}{0ex}}\leqq \phantom{\rule{0.28em}{0ex}}2\pi \phantom{\rule{0.16em}{0ex}},\hfill \end{array}$$## Eq. 4

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} r(\theta) = 3\sin ^4 \,\theta + 0.75\,. \end{equation}\end{document} $$r\left(\theta \right)=3{\mathrm{sin}}^{4}\phantom{\rule{0.16em}{0ex}}\theta +0.75\phantom{\rule{0.16em}{0ex}}.$$*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 μm

^{3}and its geometric cross section [TeX:] $\sigma _c$ ${\sigma}_{c}$ of 29.5 μm

^{2}using the following approximation:

## Eq. 5

[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} $${\sigma}_{c}=\frac{2{\int}_{0}^{\pi}r\left(\theta \right)\phantom{\rule{0.16em}{0ex}}\mathrm{d}\theta +\pi {\left[r(\pi /2)\right]}^{2}}{2}\phantom{\rule{0.16em}{0ex}}.$$## Eq. 6

[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} $$\begin{array}{ccc}\hfill {P}_{n}& =& \left[r\left({\theta}_{n}\right)\mathrm{sin}{\theta}_{n},0,r\left({\theta}_{n}\right)\mathrm{cos}{\theta}_{n}\right]\hfill \\ \hfill n& =& 1,2,\dots {n}_{\_\mathrm{max}},\phantom{\rule{0.28em}{0ex}}0\phantom{\rule{0.28em}{0ex}}\leqq \phantom{\rule{0.28em}{0ex}}{\theta}_{n}\phantom{\rule{0.28em}{0ex}}\leqq \phantom{\rule{0.28em}{0ex}}\pi /2,\hfill \\ & & {\displaystyle \phantom{\rule{12.0pt}{0ex}}{n}_{\_\mathrm{max}}=\frac{\pi}{2\Delta {\theta}_{n}}.}\hfill \end{array}$$*n_*max becomes 50. The incident photon patterns are required only for [TeX:] $0\; \leqq\; \theta_n\; \leqq\; \pi/2$ $0\phantom{\rule{0.28em}{0ex}}\leqq \phantom{\rule{0.28em}{0ex}}{\theta}_{n}\phantom{\rule{0.28em}{0ex}}\leqq \phantom{\rule{0.28em}{0ex}}\pi /2$ because of an axial symmetry of the model.

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

## Eq. 7

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} m_{\_\!\max} = \frac{{\pi ^2 }}{{\Delta \psi _{m1} \Delta \xi _{m2} }}. \end{equation}\end{document} $${m}_{\_\mathrm{max}}=\frac{{\pi}^{2}}{\Delta {\psi}_{m1}\Delta {\xi}_{m2}}.$$*m*

_{_max}of 100. Therefore, a total number of incident photon pattern becomes

*n*

_{_max}×

*m*

_{_max}= 5000.

## 2.2.

### 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$
${\chi}_{m}$
, the program opens one of the 5000 look-up tables (*n*
_{_max} × *m*
_{_max} = 5000)
[TeX:]
$Q_k (P_n, {\bm \chi} _m)$
${Q}_{k}({P}_{n},{\chi}_{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 parameters^{28} (
[TeX:]
$M_k$
${M}_{k}$
,
[TeX:]
${\hbox{\bi N}}_k$
${\text{N}}_{k}$
,
[TeX:]
$D_k$
${D}_{k}$
). Here, *M*
_{k} are the coordinates of the scattered points,
[TeX:]
${\bm N}_k$
${\mathbf{N}}_{k}$
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$
${M}_{k}$
,
[TeX:]
${\bm N}_k$
${\mathbf{N}}_{k}$
,
[TeX:]
$D_k$
${D}_{k}$
) based on the
[TeX:]
$Q_k (P_n, {\bm \chi} _m)$
${Q}_{k}({P}_{n},{\chi}_{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$
${Q}_{1}\phantom{\rule{0.16em}{0ex}}=\phantom{\rule{0.16em}{0ex}}{R}_{1}$
. If *k* is 2 or greater,
[TeX:]
$Q_k$
${Q}_{k}$
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)$
${Q}_{k}(k\phantom{\rule{0.28em}{0ex}}\geqq \phantom{\rule{0.28em}{0ex}}2)$
as follows:

## Eq. 8

[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} $${Q}_{k}({P}_{n},{\chi}_{m})=(1-{R}_{1})\phantom{\rule{0.16em}{0ex}}(1-{R}_{k})\prod _{i=2}^{k-1}{R}_{i},\phantom{\rule{1em}{0ex}}(k\phantom{\rule{0.16em}{0ex}}\geqq \phantom{\rule{0.16em}{0ex}}2).$$The probability of reflection at the *k*'th interaction point
[TeX:]
$R_k$
${R}_{k}$
depends on the incident angle, refraction angle and the relative refractive index between the RBC (
[TeX:]
$m_{{\rm RBC}}$
${m}_{\mathrm{RBC}}$
) and
[TeX:]
$m_{{\rm plasma}}$
${m}_{\mathrm{plasma}}$
, where
[TeX:]
$m_{{\rm RBC}}$
${m}_{\mathrm{RBC}}$
is expressed by:^{12, 28, 29}

## Eq. 9

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} m_{{\rm RBC}} = m_{{\rm water}} + \beta {\rm MCHC}. \end{equation}\end{document} $${m}_{\mathrm{RBC}}={m}_{\mathrm{water}}+\beta \mathrm{MCHC}.$$^{7, 30}

## Eq. 10

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation}n_i \sin \alpha _i = n_t \sin \alpha _t, \end{equation}\end{document} $${n}_{i}\mathrm{sin}{\alpha}_{i}={n}_{t}\mathrm{sin}{\alpha}_{t},$$## Eq. 11

[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\left({\alpha}_{i}\right)=\frac{1}{2}\left[\frac{{\mathrm{sin}}^{2}\phantom{\rule{0.28em}{0ex}}({\alpha}_{i}-{\alpha}_{t})}{{\mathrm{sin}}^{2}\phantom{\rule{0.28em}{0ex}}({\alpha}_{i}+{\alpha}_{t})}+\frac{{\mathrm{tan}}^{2}\phantom{\rule{0.28em}{0ex}}({\alpha}_{i}-{\alpha}_{t})}{{\mathrm{tan}}^{2}\phantom{\rule{0.28em}{0ex}}({\alpha}_{i}+{\alpha}_{t})}\right]$$For each photon-RBC interaction, the precalculation was pursued until the quantity
[TeX:]
$Q_k$
${Q}_{k}$
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

## Eq. 12

[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} $${Q}_{\mathrm{trap}}({P}_{n},{\chi}_{m})=1-\sum _{k=1}^{k\_\mathrm{max}}{Q}_{k}({P}_{n},{\chi}_{m}).$$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.

## 2.3.

### 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:

## Eq. 13

[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} $$\frac{{P}_{s}}{{P}_{i}}=\mathrm{exp}(-2.303\times e\left(\lambda \right)\times (10\times \mathrm{MCHC})/m\mathrm{Hb}\times {D}_{k}),$$*e*(λ) is the molecular extinction coefficient for the wavelength λ in cm

^{−1}/mol/l, and

*m*Hb 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$
${D}_{k}$
and the intracellular photon velocity based on
[TeX:]
$m_{{\rm RBC}}$
${m}_{\mathrm{RBC}}$
as follows:

## Eq. 14

[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} $$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle {T}_{\mathrm{RBC}}=\frac{{D}_{k}}{c},}\\ {\displaystyle c=\frac{{c}_{0}}{{m}_{\mathrm{RBC}}},}\end{array}\end{array}$$## 2.4.

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

After the pciMC selects the scattered data [TeX:] $M_k$ ${M}_{k}$ , [TeX:] ${\bm N}\!_k$ $\mathbf{N}{k}_{}$ , and [TeX:] $D_k$ ${D}_{k}$ based on the [TeX:] $Q_k$ ${Q}_{k}$ probability distribution, the photon moves from point [TeX:] $P_n$ ${P}_{n}$ to point [TeX:] $M_k$ ${M}_{k}$ by traveling the intracellular distance of [TeX:] $D_k$ ${D}_{k}$ and the propagation vector is switched from [TeX:] ${\bm \chi} _m$ ${\chi}_{m}$ to [TeX:] ${\bm N}\!_k$ $\mathbf{N}{k}_{}$ . 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$ ${\chi}_{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.

## Eq. 15

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \Phi = 2\pi *{\rm RN}(0,\!1) \end{equation}\end{document} $$\Phi =2\pi *\mathrm{RN}(0,1)$$*q*is defined as follows.

## Eq. 16

[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} $$\begin{array}{ccc}\hfill q=a+ix+jy+kz& =& (a;x,y,z)=(a;\mathbf{v}),\hfill \\ \hfill {i}^{2}& =& {j}^{2}={k}^{2}=ijk=-1.\hfill \end{array}$$**is the point or vector in the 3-D coordinate. Although the real number**

*v**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)$ ${P}_{n}=({a}_{1};{\mathbf{v}}_{1})$ and [TeX:] ${\bm \chi} _m = (a_2; {\bm v}_2)$ ${\chi}_{m}=({a}_{2};{\mathbf{v}}_{2})$ , respectively. Next, the incident point and the vector in the simulation coordinate are also defined by [TeX:] $(a_3; {\bm v}_3)$ $({a}_{3};{\mathbf{v}}_{3})$ and [TeX:] $(a_4; {\bm v}_4)$ $({a}_{4};{\mathbf{v}}_{4})$ , 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$ ${\mathbf{v}}_{2}$ and [TeX:] ${\bm v}_4$ ${\mathbf{v}}_{4}$ . Then the description of quatemion [TeX:] $q_1$ ${q}_{1}$ on the rotation from [TeX:] ${\bm v}_2$ ${\mathbf{v}}_{2}$ of [TeX:] ${\bm v}_4$ ${\mathbf{v}}_{4}$ is calculated as

## Eq. 17

[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} $${q}_{1}=\left[(\mathrm{cos}\frac{\alpha}{2};\mathrm{sin}\frac{\alpha}{2}*({\mathbf{v}}_{2}\times {\mathbf{v}}_{4})\right].$$## Eq. 18

[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} $${q}_{2}=\left(\mathrm{cos}\frac{\Phi}{2};\phantom{\rule{0.28em}{0ex}}\mathrm{sin}\frac{\Phi}{2}*{\mathbf{v}}_{2}\right),$$*s*was calculated by the point translating about [TeX:] ${\bm v}_3$ ${\mathbf{v}}_{3}$ . These calculations are as

## Eq. 19

[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={q}_{1}{q}_{2}{s}_{\left|\right|}{\overline{q}}_{2}{\overline{q}}_{1}+{\mathbf{v}}_{3},$$## 3.

## 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 (RBC_{small}) and larger MCV and lower MCHC (RBC_{large}) 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.

## 3.1.

### 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 centrifuged^{32} 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 (RBC_{large}) and the bottom 10% with higher specific gravity as smaller RBCs (RBC_{small}). 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 RBC_{large} and RBC_{small} 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$
${\sigma}_{c}$
and
[TeX:]
$\mu _c$
${\mu}_{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.

pciMC | MCV (fL) | MCHC (g/dL) | H | σc (μm2) | μc (cm−1) |
---|---|---|---|---|---|

RBC_{Large} 1 | 55.1 | 29.3 | 0.17 | 21.915 | 556.442 |

RBC_{Large} 2 | 50.6 | 28 | 0.21 | 20.705 | 654.712 |

RBC_{Large} 3 | 49.2 | 29.2 | 0.2 | 20.322 | 660.864 |

RBC_{Small} 1 | 64.4 | 28.4 | 0.22 | 24.317 | 650.052 |

RBC_{Small} 2 | 58.5 | 27.2 | 0.2 | 22.808 | 623.804 |

RBC_{Small} 3 | 57.3 | 28.1 | 0.21 | 22.495 | 628.129 |

## Table 2

Input parameters and analysis conditions of RWMCS.

Mean MCV | Mean MCV | σs | σs − | μs | σa | μa | |||
---|---|---|---|---|---|---|---|---|---|

RWMCS | (fL) | (g/dL) | Mean H | (μm2) | (μm2) | (cm−1) | g | (μm2) | (cm−1) |

RBC_{Large} (1,2,3) | 51.6 | 28.8 | 0.193 | 35.357 | 0.005 | 1053.856 | 0.99487 | 0.02 | 0.746 |

RBC_{Large} (1,2,3) | 60.1 | 27.9 | 0.21 | 39.903 | 0.006 | 1065.299 | 0.99438 | 0.023 | 0.762 |

## 3.2.

### TSROS

The optical measurements of the RBC_{large} and RBC_{small} 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\_\phantom{\rule{0.16em}{0ex}}r\left(t\right)$
, and the center of the sample signal *s*(*t*) as
[TeX:]
$G\_\,s(t)$
$G\_\phantom{\rule{0.16em}{0ex}}s\left(t\right)$
, the mean time of flight through the sample is given by

## Eq. 21

[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} $$\begin{array}{ccc}\hfill \mathrm{mean}\phantom{\rule{0.16em}{0ex}}\mathrm{time}\phantom{\rule{0.16em}{0ex}}\mathrm{of}\phantom{\rule{0.16em}{0ex}}\mathrm{flight}:\phantom{\rule{0.16em}{0ex}}\Delta t& =& G\_\phantom{\rule{0.16em}{0ex}}s\left(t\right)-G\_\phantom{\rule{0.16em}{0ex}}r\left(t\right)\hfill \\ & =& {T}_{\mathrm{media}}\left(d\right)-{T}_{0}\left(d\right)\hfill \\ & =& T\left(d\right),\hfill \end{array}$$*d*, and [TeX:] $T_0 (d)$ ${T}_{0}\left(d\right)$ is the propagation time through air for the cuvette sample thickness

*d*.

## 3.3.

### Biconcave versus Spherical RBC model and pciMC versus RWMCS

## 3.3.1.

#### 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.

## 3.3.2.

#### 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 Turcu^{10, 11} together with the scattering and the absorption cross sections
[TeX:]
$\sigma _s$
${\sigma}_{s}$
and
[TeX:]
$\sigma _a$
${\sigma}_{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}

## Eq. 22

[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} $${n}_{\mathrm{imag}}=\frac{\mathrm{ln}10}{\pi m\mathrm{Hb}}\lambda \phantom{\rule{0.16em}{0ex}}{\varepsilon}_{\mu}M\xb7\mathrm{MCHC}.$$^{34}the experimentally observed back scattering cross section [TeX:] $\sigma _s ^ -$ ${\sigma}_{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.

## Eq. 23

[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} $$\mathrm{empirical}\phantom{\rule{0.16em}{0ex}}\mathrm{RBC}\phantom{\rule{0.16em}{0ex}}\mathrm{anisotropy}\phantom{\rule{1em}{0ex}}g=\frac{{\sigma}_{s}^{+}-20{\sigma}_{s}^{-}}{{\sigma}_{s}^{+}+20{\sigma}_{s}^{-}}.$$*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:

## Eq. 24

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \tau = \frac{{H\sigma _s }}{{{\rm MCV}}}d, \end{equation}\end{document} $$\tau =\frac{H{\sigma}_{s}}{\mathrm{MCV}}d,$$## Eq. 25

[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} $${g}_{1}\left(\tau \right)={g}^{\tau}\phantom{\rule{1em}{0ex}}\mathrm{if}\tau \phantom{\rule{0.16em}{0ex}}\leqq \phantom{\rule{0.16em}{0ex}}1.8,$$## Eq. 26

[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} $${g}_{2}\left(\tau \right)={g}^{G\left(\tau \right)}\phantom{\rule{0.28em}{0ex}}\mathrm{with}\phantom{\rule{0.28em}{0ex}}G\left(\tau \right)=\frac{(\tau -1){e}^{\tau}+1}{{e}^{\tau}-\tau -1}\phantom{\rule{0.28em}{0ex}}\phantom{\rule{0.28em}{0ex}}\mathrm{if}\phantom{\rule{0.28em}{0ex}}\tau >1.8.$$*d*is the penetrating depth of photons into blood. Although RWMCS is an effective model, Chicea and Turcu reported

^{10, 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:

## Eq. 27

[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} $${g}_{3}\left(\tau \right)={g}^{G(\tau \_\mathrm{max})}\phantom{\rule{2.em}{0ex}}\phantom{\rule{1em}{0ex}}\mathrm{if}\phantom{\rule{0.16em}{0ex}}\tau >\tau \_\phantom{\rule{0.16em}{0ex}}\mathrm{max}.$$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:

## Eq. 28

[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} $$\mu =\frac{1}{2g}\left[1+{g}^{2}+{\left(\frac{1-{g}^{2}}{1+2g\phantom{\rule{0.16em}{0ex}}\mathrm{RN}(0,1)-\mathrm{g}}\right)}^{2}\right].$$*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 10^{6} photons with 10^{4} nodes of the RBC model, showing both OD and the propagation time converge with 10^{4} photons injected into the medium. In Fig. 3, the convergence was tested for the required number of nodes with 10^{4} 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 10^{4} nodes. The mesh condition was described earlier in the Sec. 2.

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

## 4.

## Results

## 4.1.

### 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).

## 4.2.

### Scattering Distribution in pciMC

Figures 4 and 4 show example patterns of photon-RBC interactions inside a RBC having a volume of 86 μm^{3}. 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$
${Q}_{k}$
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$
${Q}_{k}$
reached the maximum interaction number of *k*_max = 100 satisfying the condition of
[TeX:]
$Q_{k\_\max }$
${Q}_{k\_\mathrm{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.

We confirmed that the internal reflections increased as the relative refractive index
[TeX:]
$m_{{\rm RBC}}$
${m}_{\mathrm{RBC}}$
/
[TeX:]
$m_{{\rm plasma}}$
${m}_{\mathrm{plasma}}$
increased, consequently increasing the intracellular optical path lengths.^{28} Concerning the propagation in the suspension of RBC_{small} with an MCHC of 28.8 g/dL and RBC_{large} with 27.9 g/dL for a normalized MCV of 86 μm^{3}, the intracellular path length was slightly different, 2.93 μm in RBC_{small} and 2.92 μm in RBC_{large}, respectively. Figure 5 shows the
[TeX:]
$Q_k$
${Q}_{k}$
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

## Eq. 29

[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} $$\text{pciMC}\phantom{\rule{1em}{0ex}}\text{anisotropy}\phantom{\rule{1em}{0ex}}\mathrm{g}={\int}_{0}^{\pi}Q\left(\theta \right)\mathrm{cos}\theta \mathrm{d}\theta .$$*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.

## 4.3.

### 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 RBC_{large} and in Figs. 7 and 7 for RBC_{small}. 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.

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 RBC_{large} and RBC_{small} suspensions. The ratio of the intracellular migration time to the total time is shown in Fig. 8. For both RBC_{large} and RBC_{small} suspensions, the ratio of the intracellular to the total migration time decreased as the blood layer thickness increased.

## 5.

## Discussion

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 10^{4} 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.

## 5.1.

### 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 μm^{3} 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.

## 5.2.

### 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.

## 5.3.

### 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 RBC_{small} and 11.2±8.5% in RBC_{large} 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 RBC_{small} and 21.7±19.1% in RBC_{large}, respectively. There were statistically significant differences between pciMC and RWMCS with *p* < 0.047 in RBC_{small} for blood layer thickness from 2000 to 4000 μm, and *p* < 0.055 in RBC_{large} 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 RBC_{small} and 2.4±2.4% in RBC_{large}. Those for pciMC simulation were 14.5±3.3 and 9.1±3.4% in RBC_{small} and RBC_{large}, 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 (RBC_{small}) and 0.424±0.049 OD/mm (RBC_{large}), respectively, as compared to 0.431±0.036 OD/mm (RMC_{small}) and 0.4281±0.058 OD/mm (RBC_{large}) obtained in the experiments. The RWMCS also resulted in good agreement with 0.46 OD/mm (RBC_{small}) and 0.42 OD/mm (RBC_{large}). There were no statistically significant differences between pciMC and experiments with the mean deviations between them of 1.4 and 0.9% in RBC_{small} and RBC_{larger}, 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 confirmed^{10, 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.

## 5.4.

### 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. Roggan^{12} 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 (RBC_{small}) and 0.59 OD/mm (RBC_{large}), respectively, than the experiments. Hence the utility of pciMC with sphere model would be lower than that of biconcave model.

## 5.5.

### 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.

## 5.6.

### 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.

## 6.

## Conclusion

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.

## 7.

## Appendix: Computation of Intracellular Photon Pathways

[TeX:] $M_1$ ${M}_{1}$ = [TeX:] $P_n$ ${P}_{n}$ ,

[TeX:] $M_k (k\,\geqq\,2)$ ${M}_{k}(k\phantom{\rule{0.16em}{0ex}}\geqq \phantom{\rule{0.16em}{0ex}}2)$ ,

was decided by

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

2 COND2: [TeX:] $ {\bm U}_{k - 1} \cdot {\bm B} > 0$ ${\mathbf{U}}_{k-1}\xb7\mathbf{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}_{(k-1)x}\phantom{\rule{0.16em}{0ex}}\geqq \phantom{\rule{0.16em}{0ex}}0\cap {M}_{kx}\phantom{\rule{0.16em}{0ex}}\geqq \phantom{\rule{0.16em}{0ex}}0)\cup [{M}_{(k-1)x}\phantom{\rule{0.16em}{0ex}}\leqq \phantom{\rule{0.16em}{0ex}}0\cap {M}_{kx}\phantom{\rule{0.16em}{0ex}}\leqq \phantom{\rule{0.16em}{0ex}}0]$

b On the line equation connecting [TeX:] $M_{(k - 1)x}$ ${M}_{(k-1)x}$ and [TeX:] $M_{(k - 1)z}$ ${M}_{(k-1)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} $$\begin{array}{ccc}\hfill \left|z\right|& =& \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|\phantom{\rule{1em}{0ex}}\hfill \\ & & \times \phantom{\rule{0.16em}{0ex}}\phantom{\rule{0.16em}{0ex}}r\left(0\right).\hfill \end{array}$$4 COND3b: In COND3a, [TeX:] $M_{(k - 1)x}$ ${M}_{(k-1)x}$ and [TeX:] $M_{kx}$ ${M}_{kx}$ are replaced by [TeX:] $M_{(k - 1)y}$ ${M}_{(k-1)y}$ and [TeX:] $M_{ky}$ ${M}_{ky}$ , respectively.

5 COND4: [TeX:] $ {\bm U}_{k - 1} \cdot {{ {\bm B}}}/{{| {\bm B}|}}$ ${\mathbf{U}}_{k-1}\xb7\mathbf{B}/\left|\mathbf{B}\right|$ 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} ]$
${M}_{k-1}=[{M}_{(k-1)x},{M}_{(k-1)y},{M}_{(k-1)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} ]$
${\mathbf{U}}_{k-1}=[{U}_{(k-1)x},{U}_{(k-1)y},{U}_{(k-1)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$
$\mathbf{B}$
is the vector joining from point
[TeX:]
$M_{k - 1}$
${M}_{k-1}$
to the current scanning point
[TeX:]
$M_k$
${M}_{k}$
. 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}$
$\mathbf{B}$
is at parallel with
[TeX:]
${\bm U}_{k - 1}$
${\mathbf{U}}_{k-1}$
.

At the same time of deriving point [TeX:] $M_k$ ${M}_{k}$ , the scattered vector [TeX:] ${\bf N}_k$ ${\mathbf{N}}_{k}$ is also calculated; [TeX:] ${\bf N}_1$ ${\mathbf{N}}_{1}$ and [TeX:] ${\bf N}_k (k\,\geqq\,2)$ ${\mathbf{N}}_{k}(k\phantom{\rule{0.16em}{0ex}}\geqq \phantom{\rule{0.16em}{0ex}}2)$ are based on the reflection law at the incident point and Snell's law at the scattered point, respectively.

## Acknowledgment

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).

## References

**,” Los Alamos Sci., (Special Issue) 125 –130 (1987). Google Scholar**

*The beginning of the Monte Carlo methods***,” Methods and Applications in Neutronics, Photonics and Statistical Physics, 33 –47 Springer-Verlag, New York (1983). Google Scholar**

*MCNP—a general Monte Carlo code for neutron and photon transport***,” Med. Phys., 10 (6), 824 –830 (1983). https://doi.org/10.1118/1.595361 Google Scholar**

*A Monte Carlo model for the absorption and flux distributions of light in tissue***,” Proc. SPIE, IS (5), 102 –111 (1989). Google Scholar**

*A Monte Carlo model of light propagation in tissue***,” Lasers Surg. Med., 9 (2), 148 –154 (1989). https://doi.org/10.1002/lsm.1900090210 Google Scholar**

*Light distributions in artery tissue: Monte Carlo simulations for finite-diameter laser beams***,” Opt. Spectrosc., 72 934 –939 (1992). Google Scholar**

*Light transport in multilayerd scattering media. Monte Carlo modeling***,” Comput. Meth. Progr. Biomed., 47 131 –146 (1995). https://doi.org/10.1016/0169-2607(95)01640-F Google Scholar**

*MCML—Monte Carlo modeling of light transport in multi-layered tissues***,” Phy. Med. Biol., 47 4271 –4285 (2002). https://doi.org/10.1088/0031-9155/47/23/312 Google Scholar**

*Influence of refractive index matching on the photon diffuse reflectance***,” Opt. Express, 13 (23), 9181 –9195 (2005). https://doi.org/10.1364/OPEX.13.009181 Google Scholar**

*New model for light propagation in highly inhomogeneous polydisperse turbid media with applications in spray diagnostics***,” J. Optoelectron. Adv. Mater., 8 (4), 1516 –1519 (2006). Google Scholar**

*Testing a new multiple light scattering phase function using RWMCS***,” Optik, 118 232 –236 (2007). Google Scholar**

*RWMCS—an alternative random walk Monte Carlo code to simulate light scattering in biological suspensions***,” J. Biomed. Opt., 4 (1), 36 –46 (1999). https://doi.org/10.1117/1.429919 Google Scholar**

*Optical properties of circulating human blood in the wavelength range 400-2500 nm***,” Proc. SPIE., 2678 314 –324 (1996). Google Scholar**

*The optical properties of blood in the near infrared spectral range***,” J. Biomed. Opt., 4 (1), 47 –53 (1999). https://doi.org/10.1117/1.429920 Google Scholar**

*Influence of the scattering phase function approximation on the optical properties of blood determined from the integrating sphere measurements***,” Proc. SPIE, 2982 324 –330 (1997). Google Scholar**

*Different phase function approximations to determine optical properties of blood: a comparison***,” Appl. Opt., 36 (1), 21 –31 (1997). https://doi.org/10.1364/AO.36.000021 Google Scholar**

*Theoretical and experimental investigation of near-infrared light propagation in a model of the adult head***,” Phys. Med. Biol., 43 2465 –2478 (1998). https://doi.org/10.1088/0031-9155/43/9/003 Google Scholar**

*Near-infrared optical properties of ex vivo human skin and subcutaneous tissues measured using Monte Carlo inversion technique***,” Physiol. Meas., 23 741 –753 (2002). https://doi.org/10.1088/0967-3334/23/4/312 Google Scholar**

*Quantitative assessment of skin layers absorption and skin reflectance spectra simulation in the visble and near-inrared spectral regions***,” IEEE Trans. Biomed. Eng., 50 (8), 1026 –1033 (2003). https://doi.org/10.1109/TBME.2003.814532 Google Scholar**

*Optical transmission of blood: effect of erythrocyte aggregation***,” Proc. SPIE, 4162 120 –129 (2000). Google Scholar**

*RBC aggregation effects on light scattering from blood***,” Appl. Opt., 18 (20), 3484 –3488 (1979). https://doi.org/10.1364/AO.18.003484 Google Scholar**

*Backscattering of a picosecond pulse from densely distributed scatterers***,” Appl. Opt., 28 (12), 2331 –2336 (1989). https://doi.org/10.1364/AO.28.002331 Google Scholar**

*Time resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties***,” IEEE. J. Quantum Electron., 26 (12), 2186 –2199 (1990). https://doi.org/10.1109/3.64355 Google Scholar**

*Optical reflectance and transmittance of tissues: principles and applications***,” Proc. IEEE, 80 (6), 918 –930 (1992). https://doi.org/10.1109/5.149454 Google Scholar**

*Time-dependent optical spectroscopy and imaging for biomedical applications***,” J. Biomed. Opt., 3 (3), 364 –372 (1998). https://doi.org/10.1117/1.429883 Google Scholar**

*Scattering of light by a red blood cell***,” Appl. Opt., 37 (31), 7410 –7418 (1998). https://doi.org/10.1364/AO.37.007410 Google Scholar**

*Single scattering by red blood cells***,” Appl. Opt., 27 (19), 4027 –4033 (1988). https://doi.org/10.1364/AO.27.004027 Google Scholar**

*Comparison of Mie theory and the light scattering of red blood cells***,” Proc. SPIE, 7573 (757316), 1 –10 (2010). Google Scholar**

*Photon cell interactive Monte Carlo (**pciMC*) model to describe both intracellular and extracellular optical pathways of biconcave red blood cells: phase function and albedo**,” PhD Dissertation, Department of Biomedical Engineering, Case Western Reserve University, Cleveland, OH(1979). Google Scholar**

*On the theory and development of a noninvasive tissue reflectance oximeter***,” J. Quant. Spectrosc. Radiat. Transf., 106 285 –296 (2007). https://doi.org/10.1016/j.jqsrt.2007.01.041 Google Scholar**

*Light scattering by arbitrarily oriented optically soft spheroidal particles: calculation in geometric optics approximation***,” J. Bio. Chem., 48 173 –183 (1943). Google Scholar**

*The absorption spectra of hemoglobin and its derivatives in the visible and near infrared regions***,” Hemorheol. Microcirc., 18 67 –74 (1998). Google Scholar**

*Red blood cell aging and risk of cardiovascular diseases***,” Appl. Opt., 24 (9), 1355 –1365 (1985). https://doi.org/10.1364/AO.24.001355 Google Scholar**

*Flow-cytometric light scattering measurement of red blood cell volume and hemoglobin concentration***,” PhD Dissertation, Dept. of Electrical Engineering, Seattle(1975). Google Scholar**

*Optical diffuse and transmittance from an anisotropically scattering finite blood medium**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

**,” Artif. Organs, 25 (5), 336 –340 (2001). https://doi.org/10.1046/j.1525-1594.2001.025005336.x Google Scholar**

*Computational fluid dynamics as a development tool for rotary blood pumps*