2.1.

Model Structure

2.1.1.

Structure of the fingertip

Details of the anatomical structure of the fingertip can be found in several references.^15–17 The sagittal cross section of the fingertip and skin model has been shown in the top of Fig. 1(a). The bone of the distal phalanx is surrounded beneath by the soft deformable tissue of the pulp and above by the nail unit. The nail bed is richly vascularized with blood flowing from the digital arteries into a network of arterioles, through capillary loops just under the surface, and back out through the venules to the digital veins. Figure 1(b) shows a diagram of the principal arteries and veins. The digital arteries measure is $0.85\pm 0.1\mathrm{mm}$ in diameter.¹⁸ There are approximately three or more arteries that take off from the distal arch measuring $0.58\pm 0.1\mathrm{mm}$.¹⁸ The venous anatomy is significantly variable and consists of a dorsal and palmar system.^19,20

Fig. 1

(a) Structural anatomy of the fingertip and (b) vascular anatomy of the fingertip.

Figure 2 presents the radiographic measurements of the distal phalanx using GE Imagecast PACS, and Table 1 shows the overall mean for measurement.²¹ These data were averaged and further used as the referencing data in the simulation of phantom.

Fig. 2

Radiographic measurements of the distal phalanx using GE Imagecast PACS.

Table 1

Mean dimension by finger.

Finger	Measurement of distal phalanx	Dimension (mm)	Tip—soft tissues of the tip of the distal phalanx (mm)
Index	Width	4.8	$3.84\pm 0.59$
	Height	3.4

2.1.2.

Optical properties

The optical parameters of human skin layers are shown in Fig. 3 and Table 2.²² Evidently not all parameters for skin layers are known yet; some of them were extrapolated. Especially, for anisotropy factor $g$, there is the lack of experimental data.

Fig. 3

Skin model: 1, epidermis; 2, dermis; 3, dermis with plexus superficialis; 4, dermis; 5, dermis with plexus profundus; and 6, subcutaneous tissue.

Table 2

Optical properties of the different layers of the tissue model.

N	Skin layer	λ (nm)	μa (cm−1)	μs′ (cm−1)	g	n	d (μm)
1	Epidermis	830	$0.51\pm 0.14$	$23.37\pm 3.15$	0.86	1.417	100
2	Dermis	830	$0.40\pm 0.11$	$14.02\pm 1.89$	0.86	1.385	200
3	Dermis with plexus superficialis	830	—	—	0.86	—	200
4	Dermis	830	$0.40\pm 0.11$	$14.02\pm 1.89$	0.86	1.385	900
5	Dermis with plexus profundus	830	—	—	0.86	—	600
6	Subcutaneous tissue	830	$1.1\pm 0.1$	$11.62\pm 4.63$	0.86	1.44	500

2.1.3.

Geometry of model

In this paper, a hemispherical volume of tissue is used to represent the tip of a finger on which an oximeter probe is placed. The geometry of the problem is illustrated in Fig. 4. The light source is modeled as a point source emitting monochromatic light at a wavelength in the near-infrared (NIR) region of the spectrum.

Fig. 4

Cross section of hemispherical geometry on which the model is defined.

For mathematical convenience, the point source is represented as a dipole composed of a photon source and a photon sink. The detector receives light at a point opposite to the source (transmission case) or at a point on the same side as the source (reflection case). The finger is treated as a heterogeneous mixture of blood and other tissues containing no blood, and its optical properties are assumed to be determined mostly by dermal and subcutaneous tissues. The pigmented epidermis, which is very thin ($<100\mathrm{pm}$) and scatters light mostly in the forward direction,²³ is neglected in this analysis because it only serves to attenuate the source light before and after it undergoes absorption and scattering in the blood-perfused tissue.

3.1.

Detection Light Intensity

Generally, the detected light intensity of a reflection pulse oximeter is weaker than that of a transmission pulse oximeter with the same source–detector separation. To determine the detectable light intensity of a pulse oximeter, we calculated the detected light for reflection and transmission modes using the analysis of light propagation in tissue.

3.1.1.

Diffusion approximation implementation in numerical methods

Partial differential equation representation of diffusion approximation

The theory of diffusion approximation is based on the radiative transport equation and the conservation of energy. This equation is in general not amenable to an analytic solution for arbitrary light input. Only in specific and restricted geometries are closed-form solutions known.^24,25

The radiative transport equation in a media that scatters and absorbs photons is

Eq. (1)

$$\frac{1}{c}\times \frac{\partial L(r,s,t)}{\partial t}+s·\nabla L(r,s,t)+{\mu }_{t}L(r,s,t)={\mu }_{\mathrm{s}}{\int }_{4\pi }L(r,s^{\prime },t)p(s^{\prime },s)\mathrm{d}{\mathrm{\Omega }}^{\prime }+S(r,s,t),$$

where $L(r,s,t)$ is the average power flux density at point $r$ in the direction of unit vector $s$, at time $t$, in the unit of ($\mathrm{W}.{\mathrm{m}}^{-2}.{\mathrm{sr}}^{-1}$). $C$ is the speed of light in the medium. ${\mu }_{\mathrm{t}}={\mu }_{\mathrm{s}}+{\mu }_{\mathrm{a}}$ is the optical transport coefficient, where ${\mu }_{\mathrm{s}}$ is the scattering coefficient and ${\mu }_{\mathrm{a}}$ is the absorption coefficient. $S(r,s,t)$ is the spatial and angular distribution of the source in the unit of ($\mathrm{W}.{\mathrm{m}}^{-3}.{\mathrm{sr}}^{-1}$). The normalized phase function $p(s^{\prime },s)$ represents the probability of scattering in a direction $s^{\prime }$ from direction $s$.

Various tissue optical properties, such as the absorption coefficient, scattering coefficient, anisotropy index, and refractive index, can be used to calculate the power distributions and trajectories of photons. The scattering behavior is usually described by the Henyey–Greenstein phase function.²⁶ The photons pass that through and interact with scatterers in the biological tissue are subject to a high degree of forward scattering. The phase function describes the power distribution of scattering from the forward direction (0 deg) to the backward direction (180 deg), and its equation is given by

Eq. (2)

$$P_{HG}(\theta )=\frac{\omega (1-g^2)}{{[1+g^2-2g\cos (\theta )]}^{\frac{3}{2}}},$$

where $P_{HG}$ denotes the scattered intensity at a specific angle and $w$ and $\theta$ denote a normalization parameter and the direction of scattering, respectively. Parameter $g$, named the scattering anisotropy parameter, denotes the mean cosine of the scattering angle and is given by

Eq. (3)

$$g=⟨\cos \theta ⟩={\int }_0^{\pi }{P}_{HG}(\theta )\cos \theta 2\pi \sin \theta \mathrm{d}\theta .$$

We use Eq. (2) to obtain the normalized scattering phase function. Figure 5 shows the phase-function diagrams at 656 and 943 nm for different tissues: epidermis, dermis and pulp, bone, and blood. Since scattering anisotropy parameter $g$ is close to 1, the diagrams exhibit a sharp peak at 0 deg, which indicates that light tends to scatter in the forward direction during interactions between photons and biological tissues.

Fig. 5

Normalized angular scattering intensity distributions for epidermis, dermis and pulp, bone, and blood at (a) 656 and (b) 943 nm.

On the other hand, treating the tissue as a heterogeneous mixture, the total absorption coefficient can be written as a sum of the absorption coefficients of the individual tissue components

Eq. (4)

$${\mu }_{\mathrm{a}}=V_{\mathrm{a}}{\mu }_{\mathrm{a}}^{\text{art}}+V_{\mathrm{v}}{\mu }_{\mathrm{a}}^{\mathrm{ven}}+[1-(V_{\mathrm{a}}+V_{\mathrm{v}})]{\mu }_{\mathrm{a}}^{\mathrm{tis}},$$

where $V_{\mathrm{a}}$ and $V_{\mathrm{v}}$ are the volume fractions of arterial and venous blood, respectively, and ${\mu }_{\mathrm{a}}^{\mathrm{tis}}$ is the absorption coefficient of the bloodless tissue. The absorption coefficients of the arterial and venous blood in the above equation are ${\mu }_{\mathrm{a}}^{\text{art}}$ and ${\mu }_{\mathrm{a}}^{\mathrm{ven}}$, respectively.

Like the absorption coefficient, the total scattering coefficient of the idealized heterogeneous tissue can also be written as a sum of the coefficients of the individual tissue components

Eq. (5)

$${\mu }_{\mathrm{s}}=(V_{\mathrm{a}}+V_{\mathrm{v}}){\mu }_{\mathrm{s}}^{\mathrm{b}}+[1-(V_{\mathrm{a}}+V_{\mathrm{v}})]{\mu }_{\mathrm{s}}^{\mathrm{tis}},$$

where ${\mu }_{\mathrm{s}}^{\mathrm{b}}$ and ${\mu }_{\mathrm{s}}^{\mathrm{tis}}$ are the scattering coefficients of whole blood and of the component of the tissue containing no blood, respectively.

The photon fluence is given by

Eq. (6)

$$\mathrm{\Phi }(r,t)={\int }_{4\pi }L(r,s,t)\mathrm{d}\mathrm{\Omega }.$$

The photon flux or current density is given by

Eq. (7)

$$J(r,t)={\int }_{4\pi }sL(r,s,t)\mathrm{d}\mathrm{\Omega }.$$

Both the fluence and the flux are in units of $\mathrm{W}.{\mathrm{m}}^{-2}$.

In vivo experimental light propagates in a highly scattering tissue in which the scattering coefficient is considerably greater than the absorption coefficient and the diffusion approximation is sufficiently accurate. In the $P_1$ approximation, the radiance is expanded in spherical harmonics truncated to the first order, and the radiance equation can be reduced to the time-dependent diffusion equation ^27,28

Eq. (8)

$$\frac{1}{c}\frac{\partial \mathrm{\Phi }(r,t)}{\partial t}+{\mu }_{\mathrm{a}}\mathrm{\Phi }(r,t)-D{\nabla }^2\mathrm{\Phi }(r,t)=S(r,t),$$

where $\mathrm{\Phi }(r,t)$ is the photon fluence rate and $D={\{3[{\mu }_{\mathrm{a}}+(1-g){\mu }_{\mathrm{s}}]\}}^{-1}$ is the diffusion coefficient.

The steady-state diffusion equation for the average diffuse intensity $\mathrm{\Phi }(r)$ in a heterogeneous scattering medium of absorption coefficient ${\mu }_{\mathrm{a}}$ is

Eq. (9)

$$D{\nabla }^2\mathrm{\Phi }(r)-{\mu }_{\mathrm{a}}\mathrm{\Phi }(r)+S(r)=0.$$

Analytical solution of PDE representation of diffusion approximation:

With dimensional parameters defined in Fig. 4, the diffusion approximation equation, by considering the approximation of ${\mu }_{\mathrm{s}}^{\prime }\gg {\mu }_{\mathrm{a}}$ in a highly scattering tissue, is as follows:

Eq. (10)

$$\frac{1}{r^2}\frac{\partial }{\partial r}\left[r^2\frac{\partial \mathrm{\Phi }(r)}{\partial r}\right]-{\alpha }^2\mathrm{\Phi }(r)=-\frac{1}{D}S(r),$$

Eq. (11)

$$D=\frac{1}{3({\mu }_{\mathrm{a}}+{\mu }_{\mathrm{s}}^{\prime })}\approx \frac{1}{3{\mu }_{\mathrm{s}}^{\prime }},$$

where ${\mu }_{\mathrm{s}}^{\prime }$ is the reduced scattering coefficient and $\alpha =\sqrt{\frac{{\mu }_{\mathrm{a}}}{D}}\approx \sqrt{3{\mu }_{\mathrm{a}}{\mu }_{\mathrm{s}}^{\prime }}$ is the effective attenuation coefficient. An explanation of these coefficients in the context of diffusion theory as applied to biological tissues can be found in a review by Cheong et al.²⁹ For a point source, $S(r)$ can be written in terms of the volume delta function $\delta (p)$, as below:

Eq. (12)

$$S(r)=\frac{3P_0}{4\pi }\delta (r),$$

where $P_0$ is the power emitted by the source. We will first solve Eq. (10) for $\mathrm{\Phi }(r)$ within a bounded spherical medium and then derive the solution for the hemispherical geometry as depicted in Fig. 4 from this result using the method of images.³⁰

The general solution to Eq. (10) is

Eq. (13)

$$\mathrm{\Phi }(r)=\frac{A}{r}{e}^{-\alpha r}+\frac{B}{r}{e}^{\alpha r},$$

where $A$ and $B$ are the constants determined by two boundary conditions. The first boundary condition is as follows:

Eq. (14)

$$\underset{r\to 0}{\lim }[\text{flux density}]=\underset{r\to 0}{\lim }[{\int }_{S}J·\mathrm{d}A]=\underset{r\to 0}{\lim }[-4\pi D{r}^2\frac{\partial \mathrm{\Phi }(r)}{\partial r}]=P_0,$$

where $J=-D\nabla \mathrm{\Phi }$ denotes the flux of fluence rate.

The second boundary condition is that all photon flux is absorbed beyond the boundary $r=d$

Eq. (15)

$$\mathrm{\Phi }(r)=0\text{at}r=d.$$

Subsequently, $\mathrm{\Phi }(r)$ is, for bounded spherical geometry

Eq. (16)

$${\mathrm{\Phi }}_{\mathrm{sph}}(r)=\frac{P_0}{4\pi D}[\frac{1}{r(1-e^{-2\alpha d})}{e}^{-\alpha r}+\frac{1}{r(1-e^{2\alpha d})}{e}^{\alpha r}].$$

Also, the AS of the diffusion approximation equation for a point beam incident on the surface of a semi-infinite medium is given by summing the solutions for a pair of positive and negative point sources as follows:

Eq. (17)

$${{\mathrm{\Phi }}_{\mathrm{sph}}(\rho ,z)|}_{z=0}=\frac{P_0}{4\pi D}[\frac{1}{r_1}{e}^{-\alpha r_1}-\frac{1}{r_2}{e}^{-\alpha r_2}],$$

where $z$ is the depth from the surface and $\rho$ is the distance between the point of interest and the $z$-axis, as shown in Fig. 4. $r_1$ and $r_2$ are given as follows:

Eq. (18)

$$r_1={[{(z-z_0)}^2+{\rho }^2]}^{0.5},$$

Eq. (19)

$$r_2={[{(z+z_0+2z_e)}^2+{\rho }^2]}^{0.5}.$$

In these equations, $z_e$ is the distance between the real and extrapolated boundaries for the refractive index mismatched case. $z_0=1/{\mu }_{\mathrm{s}}^{\prime }$ is a small value representing the depth below the surface from which the first-scattered photons emanate. We assume here that the incident photons are converted to scattered photons within one mean scattering length, so $z_e$ is simply given as

Eq. (20)

$$z_e=2DA,$$

where $A$ is related to the internal reflection $r_j$ and can be derived from Fresnel reflection coefficients.³¹ When the boundary has matched refractive indices, $A=1$; otherwise, $A$ is estimated as³²

Eq. (21)

$$A=\frac{(1+r_j)}{(1-r_j)},$$

where

Eq. (22)

$$r_j=-1.44n_{\mathrm{rel}}^{-2}+0.71n_{\mathrm{rel}}^{-1}+0.668+0.0636n_{\mathrm{rel}},$$

where $n_{\mathrm{rel}}$ is the relative index of tissue to air. Now, by applying the method to images, the desired solution for hemispherical geometry can be obtained from Eq. (16)

Eq. (23)

$${\mathrm{\Phi }}_{\mathrm{hemi}}(r)=\frac{-2z_{0}z}{r}\frac{\partial {\mathrm{\Phi }}_{\mathrm{sph}}(r)}{\partial r}=\frac{z_{0}z{P}_0}{2\pi D{r}^2}\left[\frac{1}{(1-e^{-2\alpha d})}{e}^{-\alpha r}\right(\frac{1}{r}+\alpha )+\frac{1}{(1-e^{2\alpha d})}{e}^{\alpha r}(\frac{1}{r}-\alpha \left)\right].$$

The magnitude of the intensity captured by the reflected light detector (labeled “#1” in Fig. 4) is given by

Eq. (24)

$$I_r=-D{\frac{\partial {\mathrm{\Phi }}_{\mathrm{hemi}}(r)}{\partial z}|}_{\genfrac{}{}{0ex}{}{z=0}{r=p}}=-\frac{z_0{P}_0}{2\pi {\rho }^2}\left[\frac{1}{(1-e^{-2\alpha d})}{e}^{-\alpha \rho }\right(\frac{1}{\rho }+\alpha )+\frac{1}{(1-e^{2\alpha d})}{e}^{\alpha \rho }(\frac{1}{\rho }-\alpha \left)\right].$$

The magnitude of the intensity captured by the transmitted light detector (labeled “#2” in Fig. 4) is given by

Eq. (25)

$$I_t=D{\frac{\partial {\mathrm{\Phi }}_{\mathrm{hemi}}(r)}{\partial z}|}_{\genfrac{}{}{0ex}{}{r=d}{\rho =0}}=-\frac{z_0{P}_0}{2\pi d}\left[\frac{1}{(1-e^{-2\alpha d})}{e}^{-\alpha d}\right(\frac{2}{d^2}+\frac{2\alpha }{d}+{\alpha }^2)+\frac{1}{(1-e^{2\alpha d})}{e}^{\alpha d}(\frac{2}{d^2}-\frac{2\alpha }{d}+{\alpha }^2\left)\right].$$

From the above equations, a model of the response from a pulse oximeter can be developed.

Numerical method to solve the diffusion approximation equation

Finite element method:

In this section, the FEM model has been used to perform the forward problem of light propagation in scattering tissue.

Numerical experiments:

Due to the accuracy and efficiency of FEM, we used this method to calculate solutions of the diffusion approximation equation along with the Dirichlet boundary conditions (DBC) and Robin boundary conditions (RBC) presented in Eqs. (26) and (27), respectively,³³

Eq. (26)

$$\mathrm{\Phi }(r)+2\frac{1}{3({\mu }_{\mathrm{a}}+{\mu }_{\mathrm{s}}^{\prime })}\frac{(1+r_j)}{(1-r_j)}\frac{\partial \mathrm{\Phi }(r)}{\partial n}=0,$$

Eq. (27)

$$\mathrm{\Phi }(r)=0,$$

where $n$ is the normal vector.

Assumptions used in the FEM formulation is that $\mathrm{\Phi }(r,t)$ is continuous and approximately obeys the diffusion equation locally within the domain $\mathrm{\Omega }$. To find a piecewise linear approximation and continuous function ${\mathrm{\Phi }}^h(r,t)$ of $\mathrm{\Phi }(r,t)$, $\mathrm{\Omega }$ is divided into $P$ nonoverlapping elements ${\tau }_i(i=1,2,\dots ,P)$, such that $\mathrm{\Omega }=Y_{i=1}^p{\tau }_i$, joined at $D$ vertex nodes. Then, the solution obtained for $\mathrm{\Phi }(r,t)$ is approximated by a piecewise polynomial and continuous function as

Eq. (28)

$${\mathrm{\Phi }}^h(r,t)=\sum {\mathrm{\Phi }}_i(t){\mathrm{\Psi }}_i(r),$$

where the set $\{{\psi }_i(r)\}$ is the basis function set for expansion. Considering these assumptions, the diffusion Eq. (8) is declared as^12,18

Eq. (29)

$$[K(D)+C({\mu }_{\mathrm{a}}c)]\mathrm{\Phi }+M\frac{\partial \mathrm{\Phi }}{\partial t}=S,$$

where

Eq. (30)

$$K_{ij}={\int }_{\mathrm{\Omega }}D(r)\nabla {\mathrm{\Psi }}_j(r).\nabla {\mathrm{\Psi }}_i(r)\mathrm{d}\mathrm{\Omega },$$

Eq. (31)

$$S_j(t)={\int }_{\mathrm{\Omega }}{\mathrm{\Psi }}_j(r){\mathrm{q}}_0(r,t)\mathrm{d}\mathrm{\Omega },$$

Eq. (32)

$$C_{ij}={\int }_{\mathrm{\Omega }}{\mu }_{\mathrm{a}}(r){\mathrm{c}\mathrm{\Psi }}_j(r){\mathrm{\Psi }}_i(r)\mathrm{d}\mathrm{\Omega },$$

Eq. (33)

$$M_{ij}={\int }_{\mathrm{\Omega }}{\mathrm{\Psi }}_j(r){\mathrm{\Psi }}_i(r)\mathrm{d}\mathrm{\Omega },$$

Eq. (34)

$$\mathrm{\Phi }={[{\mathrm{\Phi }}_1(t),{\mathrm{\Phi }}_2(t),\dots ,{\mathrm{\Phi }}_D(t)]}^T,$$

where $\{{\psi }_i(r)\}$ is the basis function and $\{{\psi }_j(r)\}$ is the test function.

By considering the time-independent case, the equation can be written as

Eq. (35)

$$[\nabla .D\nabla -{\mu }_{\mathrm{a}}(r)]\mathrm{\Phi }(r)=-s_0(r).$$

Following the above derivation, the equation was obtained as

Eq. (36)

$$[K(D)+C({\mu }_{\mathrm{a}})]\mathrm{\Phi }=S,$$

where $\mathrm{\Phi }$, $C$, $D$, $K$, and $S$ have the same meaning as mentioned earlier.

The above explanation of the FEM discretization preserves the frequency-domain differential equation (DE). Also, associating the Robin boundary condition, which can be expressed as

Eq. (37)

$$D\nabla \mathrm{\Phi }.\hat{n}=\beta \mathrm{\Phi },$$

where $\beta$ is considered as 0.177 for a tissue–air interface and $\hat{n}$ is the unit vector normal to the boundary, the DE discrete formulation in the frequency domain can be explained as

Eq. (38)

$$[K(D)+C({\mu }_{\mathrm{a}})+\beta A+i\omega c^{-1}M]\mathrm{\Phi }(\omega )=S(\omega ),$$

where $S$ and $\mathrm{\Phi }$ exhibit the source values and nodal solution, respectively, $M$, $K$, and $C$ are system matrices as described above, and $A_{ij}={\int }_{\mathrm{\Omega }}{\mathrm{\Psi }}_j(r){\mathrm{\Psi }}_i(r)\mathrm{d}\mathrm{\Omega }$.

Optical simulation model of fingertip:

The finger base mainly consists of epidermis, phalanx, blood, and dermis and pulp. The fingertip model used in this work has been illustrated in Fig. 6. In this model, the height, thickness, and outer epidermis’ thickness of employed geometry are 15, 10, and 0.1 mm, respectively. Also, the characteristics of distal phalanx have been listed in Table 1. As defined in Sec. 2.1.1, three pairs of arteries with 0.85-mm diameter are considered top and bottom of the distal phalanx. Additionally, there are several veins with 0.85-mm diameter on both sides of the distal phalanx. Table 3 shows the different tissues optical properties at wavelengths 656 and 943 nm.^9–35 In the simulated model, the wavelength of 943 nm and a fixed intensity have been employed for the light point source.

Fig. 6

Cross sections and models of the fingertip.

Table 3

Optical properties used in the tissue simulations: μa (absorption coefficient), μs (scattering coefficient), n (refractive index), and g.

	λ (nm)	μa (mm−1)	μs′ (mm−1)	n	g
Dermis and pulp	656	0.05	8.6	1.38	0.91
	943	0.03	7	1.38	0.91
Epidermis	656	0.05	18.3	1.43	0.81
	943	0.03	12.3	1.42	0.89
Bone	656	0.05	33.1	1.4	0.92
	943	0.02	26.1	1.4	0.94
Arterial blood	656	0.2	77.5	1.363	0.98
	943	0.64	65	1.358	0.99
Venous blood	656	0.72	77.5	1.362	0.98
	943	0.54	65	1.357	0.99

Mesh generation:

In this study, the developed model involved hemispherical geometry and heterogeneous cases. To improve the accuracy of the geometry approximation, six node triangular elements were used instead of three nodes. Also, more accurate answers were obtained by using a quadratic shape function in comparison of three-node linear triangular elements. Figure 7 illustrates the six node triangular elements in the model.

Fig. 7

Six node triangular element.

3.1.3.

Experimental study

The purpose of experiments is to validate and evaluate the performance of the mathematical model by making measurements of steady-state intensity over the tissue geometry. Experiments on 19 healthy volunteers were conducted. All of the volunteers who participated in the study were males with ages ranging between 18 to 43 years with the mean and standard deviation of 27.36 and 7.55, respectively (see Table 4). This study was approved by the ethics committee of the Shahid Beheshti University of Medical Sciences.

Table 4

Biometric detail of the volunteers that was investigated in the study.

Volunteer	Sex	Age	Weight (kg)	BP (mmHg)	Heart rate (bpm)
1	M	18	52	138.84	74
2	M	22	64	132.74	74
3	M	29	75	140.75	80
4	M	20	72	130.65	65
5	M	18	62	130.68	83
6	M	19	61	145.82	78
7	M	23	65	125.55	75
8	M	24	80	120.66	67
9	M	37	83	107.52	62
10	M	27	65	137.86	69
11	M	25	77	125.52	85
12	M	36	80	135.78	75
13	M	31	75	121.65	74
14	M	24	69	115.85	83
15	M	22	68	130.91	75
16	M	29	75	134.80	100
17	M	33	79	110.60	65
18	M	40	81	135.82	66
19	M	43	74	132.65	72

In pulse oximeters, the arterial oxygen saturation ${\mathrm{SaO}}_2$ is calculated as a function of the ac/dc ratio of detected light intensity at systole and diastole. Here, the change of the blood volume due to the heart diastole and heart systole leads to the varying signal of light absorption that is called a photoplethysmograph (PPG). The pulse oximeter can calibrate and estimate ${\mathrm{SaO}}_2$ by the parameter of ratio $R$ obtained from the PPG signal. The defining equation is

Eq. (39)

$${\mathrm{Sa}\mathrm{O}}_2=\frac{[{\mathrm{HbO}}_2]}{[{\mathrm{HbO}}_2]+[\mathrm{Hb}]}\times 100\%=f(\mathrm{ac}/\mathrm{dc}\text{ratio}),$$

where

Eq. (40)

$$\mathrm{ac}/\mathrm{dc}\text{ratio}=R=\frac{(I_{{\text{red}}_{\mathrm{ac}}}/I_{{\text{red}}_{\mathrm{dc}}})}{(I_{{\mathrm{IR}}_{\mathrm{ac}}}/I_{{\mathrm{IR}}_{\mathrm{dc}}})}=\frac{(I_{{\text{red}}_{\max }}-I_{{\text{red}}_{\min }})\times I_{{\mathrm{IR}}_{\max }}}{(I_{{\mathrm{IR}}_{\max }}-I_{{\mathrm{IR}}_{\min }})\times I_{{\text{red}}_{\max }}},$$

where $I_{{\text{red}}_{\mathrm{dc}}}$ and $I_{{\text{red}}_{\mathrm{ac}}}$ denote the DC and AC components of the PPG signal for red light, respectively. $I_{{\mathrm{IR}}_{\mathrm{dc}}}$ and $I_{{\mathrm{IR}}_{\mathrm{ac}}}$ denote the DC and AC components of the PPG signal for NIR light, respectively. The function $f$ is determined by the manufacturer from clinical and in vitro studies.

A PPG finger probe that operates in transmittance and reflectance mode has been used. Eight PPG signals were recorded from the reflectance and transmittance PPG system, consisting of four PPG signals (red DC and AC and infrared DC and AC signals) acquired from each mode that represents ${\mathrm{SpO}}_2$ values. Figures 8 and 9 show a sample of the two infrared and red AC PPG signals, respectively, obtained from the transmittance and reflectance PPG systems of one volunteer.

Fig. 8

Red PPG signals obtained simultaneously from reflectance (green) and transmittance (red) modes.

Fig. 9

Infrared AC PPG signals from all two reflectance (green) and transmittance (red) modes.

4.1.

Simulation Results of Finite Element Method

A point source is present just inside the boundary of the tissue. The radial location of the source was positioned inside the physical boundary by a distance, $z_0=1/({\mu }_{\mathrm{s}}^{\prime })$ for the point source excitation used in the computational algorithm.^39–41 It is worth mentioning that, the resolution of the mesh plays an important role for the stability and accuracy of the FEM model. Meshing is done to divide the entire specimen into nodes where the partial differential equation can be solved. The mesh is intensified where the source is present as it is the area of interest. The mesh of our structure is presented in Fig. 10, we use free mesh parameters with a maximum element size scaling factor = 1, a maximum element growth rate = 1.5, a mesh curvature factor = 0.6, a mesh curvature cutoff = 0.001, and a resolution of narrow region = 0.5.

Fig. 10

Mesh structure for extended boundary FEM calculation.

Figure 11 depicts the propagation of light waves around the point source by the contours. Contours are close together in the proximity of the light source and further toward the other side.

Fig. 11

Contour plot of propagation of light waves.

Point source has a fixed intensity. Fluence values are found throughout the specimen. The values of absorption and reduced scattering coefficients inside the specimen are known (see Fig. 12).

Fig. 12

Surface plot of propagation of light waves.

Figure 13 illustrates the simulated flux density by considering the optimal boundaries. Also, a continuous flux can be seen between the two layers and when moving outward from the two media.

Fig. 13

3-D flux flow between layers.

Figure 14 shows the cross section through the photon density distribution along the $x$-, $y$-, and $z$-axes for the FEM model, on a logarithmic scale. In all three cases, the equivalent isotropic source is placed at position $[0,1/({\mu }_{\mathrm{s}}^{\prime }),-0.004]$.

Fig. 14

(a) Photon density along the $x$-axis from $x=-7.5\mathrm{mm}$ to $x=7.5\mathrm{mm}$ for the FEM model on a logarithmic scale. (b) Photon density along the $z$-axis from $z=0\mathrm{mm}$ to $z=-10\mathrm{mm}$ for the FEM model on a logarithmic scale. (c) Photon density along the $y$-axis from $y=0\mathrm{mm}$ to $y=-14\mathrm{mm}$ for the FEM model on a logarithmic scale.

4.2.

Comparison Between Reflectance and Transmittance Modes and Investigation with Experimental Data

A reflectance oximeter has sensing and emitting diodes on a single surface. Therefore, they can be placed in an anatomical location that does not require two surfaces in close proximity. In addition, their advantage is that, regardless of the patient’s age or size (infants to very large adults), it provides more accurate and stable readings.

The light is emitted by the LEDs, passes through blood vessels and tissue, then reflects off the bone, and passes through the tissues once more, where a significant amount of light will reflect off the skin in the reflectance setup. Thus, reflectance probes have a lower signal-to-noise ratio and a higher offset than the transmittance probes. Also, reflectance setups require a significantly greater amount of light because the detected light intensity of a reflection pulse oximeter is weaker than that of a transmission pulse oximeter with the same source–detector separation (see Fig. 15).

Fig. 15

At the same source–detector separation: (a) light intensity of the reflection pulse oximeter obtained from FEM simulation and (b) light intensity of the transmission pulse oximeter obtained from FEM simulation.

The simulation results and measurement (experimental data) are compared. Figure 16 shows the results of the measurements and corresponding simulations. They show an overall good agreement. As can be seen, the ratio between measured and simulated values is almost constant; hence, the characteristics of the light propagation can be accurately reproduced.

Fig. 16

Results of the simulations and measurements of reflection and transmission modes.

4.3.

Comparison of Monte Carlo, the Analytical Method, and Finite Element Method

In this section, diffuse reflectance AS for fingertip geometry given by Eq. (24) is validated using the MC method. The optical parameters used are presented in Table 3 and under optical simulation model of fingertip. Also, the fluence rates generated by MC, the analytical method, and FEM were compared at a radial distance $r=2\mathrm{mm}$. First, all generated profiles were normalized at a depth of 5 mm; then, they were normalized at a distance of 2 mm, which is out of the region of 1 mean free path (mfp) from the isotropic point source. A reduced scattering albedo (RSA) defined by ${\mu }_{\mathrm{s}}(1-g)/{\mu }_{\mathrm{a}}$ was used to measure the ratio of photon scatting to absorption. The relative refractive index $n_{\mathrm{rel}}$ of the boundary was set at 1.37. Figures 17Fig. 18–19 present a comparison between the FEM and AS for RSAs 45.82, 35.05, and 10.56, respectively. As a result, the relative errors between the photon fluence rate obtained from the FEM model and from the AS were 8.9%, 6.4%, and 3.8% for RSAs of 45.82, 35.05, and 10.56, respectively.

Fig. 17

Comparison between the FEM simulation and AS with $\mathrm{RSA}=45.82$ and $n_{\mathrm{rel}}=1.37$.

Fig. 18

Comparison between the FEM simulation and AS with $\mathrm{RSA}=35.05$ and $n_{\mathrm{rel}}=1.37$.

Fig. 19

Comparison between the FEM simulation and AS with $\mathrm{RSA}=10.56$ and $n_{\mathrm{rel}}=1.37$.

Figures 20 and 21 show for matched and mismatched boundary conditions that AS and MC results do not exactly match near the sources. Furthermore, by comparison of both type of boundary conditions, it can be observed that solutions begin to match after some reduced mfp, that is nearer to the source in case of matched.

Fig. 20

Comparison of AS of diffuse reflectance with the MC method for (a) matched boundary and (b) relative error.

Fig. 21

Comparison of AS of diffuse reflectance with the MC method for (a) mismatched boundary and (b) relative error.

Figures 22 and 23 show the fluence rate profiles and relative errors on a logarithmic scale. Except for a small depth region within 1 mfp of the surface, FEM simulated the distribution of internal fluence rate accurately, compared to MC simulation, with relative errors of less than 3%. For the small depth region close to the light source, FEM simulated the distribution of internal fluence with relative errors of less than 9%, compared to MC simulation. These findings also agree with that of Schweiger et al.⁴²

Fig. 22

(a) Fluence rate along the axis $r=2\mathrm{mm}$ in the tissue simulated by MC and FEM and (b) relative errors with respect to MC simulation.

Fig. 23

(a) Fluence rate obtained using MC, the analytical method, and FEM and (b) their relative errors, compared to MC simulation.

According to the simulation results shown in Fig. 23, all methods underestimated the fluence rate at the boundary by a large error within the region of 1 mfp close to the light source because diffusion approximation was not valid in this region. Furthermore, the analytical method and FEM performed about the same over the entire distances. It can thus be concluded that, overall, the analytical method and FEM have almost identical results compared with MC as the gold standard. The differences are so small between FEM and AS that it is possible to use the diffusion theory for the analysis and prediction of the response of pulse oximeters. Also, it is obvious that the FEM is much more efficient and accurate than the AS method in heterogeneous advanced tissue and has the ability to estimate fluence distribution in highly complex geometries, thus making it a very useful tool in quantifying the optical parameters.

4.4.

Comparison Experimental Results with Finite Element Method and Analytical Solution

Diffuse reflectance defines mainly the light coming out from the medium and diffusely entered into the medium from outside. It depends strongly on the tissue optical properties and is mainly required to determine the optical properties of the medium. Here, a comparative study has been done for diffuse reflectance calculated using FEM and AS for fingertip geometry with respect to different scattering and absorption coefficients. Later, conclusions are presented based on this study. The purpose of experiments is to validate the analytical and FEM models by making measurements of steady-state intensity over the fingertip geometry.

Pulse oximeters use the red and infrared maximum and minimum detected light intensities to calculate ${\mathrm{SaO}}_2$. In this study, our model works in reverse. To obtain experimental data for light intensities, the ${\mathrm{SaO}}_2$ is input and the resulting transmitted and reflected light intensities are calculated to compare with diffuse reflectance values calculated using 3-D FEM and AS. The FEM solution is almost matching with the experiment result. AS is coming as per our expectations, as it has been analyzed that AS is approximately matching with FEM solution (see Fig. 24).

Fig. 24

Comparison of experimental results with FEM and AS.

Also, Fig. 24 shows good agreement between the numerical results and the experimental results. This study helps us to reach a conclusion for the feasibility of the diffusion approximation equation for the investigation of reflectance pulse oximeter in NIR region.